Biases in regional carbon budgets from Open Access

Abstract. Recent advances in atmospheric transport model inversions could significantly reduce uncertainties in land carbon uptake through the assimilation of CO 2 concentration measurements at weekly and shorter timescales. The potential of these measurements for reducing biases in estimated land carbon sinks depends on the strength of covariation between surface fluxes and atmospheric transport at these timescales and how well transport models represent this covariation. Daily to seasonal covariation of surface fluxes and atmospheric transport was estimated in observations at the US Southern Great Plains Atmospheric Radiation Measurement Climate Research Facility, and compared to an atmospheric transport model inversion (CarbonTracker). Covariation of transport and surface fluxes was stronger in CarbonTracker than in observations on synoptic (daily to weekly) timescales, with a wet year (2007) having significant covariation compared to a dry year (2006). Differences between observed and CarbonTracker synoptic covariation resulted in a 0.3 ppm CO 2 enhancement in boundary layer concentrations during the growing season, and a corresponding enhancement in carbon uptake by 13% of the seasonal cycle amplitude in 2007, as estimated by an offline simplified transport model. This synoptic rectification of surface flux variability was of similar magnitude to the interannual variability in carbon sinks alone, and indicates that interannual variability in the inversions can be affected by biases in simulated synoptic rectifier effects. The most significant covariation of surface fluxes and transport had periodicities of 10 days and greater, suggesting that surface flux inversions would benefit from improved simulations of the effects of soil moisture on boundary layer heights and surface CO 2 fluxes. Soil moisture remote sensing could be used along with CO 2 concentration measurements to further constrain atmospheric transport model inversions.


Introduction
Land and ocean ecosystems take up more than half of the carbon dioxide emitted by anthropogenic sources, but observational networks lack the coverage needed to identify the regional ecosystem processes responsible for these sinks (Dolman et al., 2010).Although spatially sparse, concentration measurements are temporally dense (Keppel-Aleks et al., 2011) and have been assimilated in transport models at higher temporal resolutions as computing power has increased, from monthly (Baker et al., 2006) to weekly (Peters et al., 2007) and diurnal timescales (Gourdji et al., 2012).The number of high-frequency and continuous surface and airborne measurements has also increased in recent years (Patra et al., 2008;Masarie et al., 2011).Whether this highfrequency data assimilation arrives at more accurate estimates of regional carbon sinks will depend on how strongly surface CO 2 fluxes covary with atmospheric transport and mixing at these frequencies, and how well land surface and transport models simulate this covariation in data assimilation systems (Patra et al., 2008).
Solar radiation absorbed at the surface promotes both boundary layer turbulence and photosynthesis, driving diurnal and seasonal covariation of surface fluxes and boundary layer depth.Vegetation draws CO 2 from deeper boundary layers in summer and respires into shallower boundary layers in winter, enhancing boundary layer concentrations relative to the free troposphere.These concentration gradients Published by Copernicus Publications on behalf of the European Geosciences Union.I. N. Williams et al.: Regional carbon budgets at synoptic timescales persist in annual averages, even for annually balanced but seasonally varying ecosystem fluxes, analogous to an electrical rectifier that converts alternating to direct current (Denning et al., 1995).While deeper boundary layers dilute concentrations exchanged at the surface, transport by the divergent circulation determines how efficiently these concentrations are mixed with those of the overlying free troposphere.These circulations become increasingly important in maintaining vertical concentration gradients at longer timescales (Williams et al., 2011).
The timescale dependence of concentration variability explains why concentration measurements often reflect an approximate balance between surface fluxes and boundary layer mixing at diurnal timescales (Yi et al., 2004), as well as between surface fluxes and vertical transport at seasonal timescales (Bakwin et al., 2004;Helliker et al., 2004;Lai et al., 2006).Yet neither of these approximations applies to concentration variability between diurnal and seasonal timescales, referred to here as synoptic variability.The synoptic covariation of boundary layer heights and surface fluxes can create persistent concentration gradients similar to seasonal rectifier effects.Deeper boundary layers and stronger photosynthetic uptake coincide with enhanced insolation following the passage of synoptic weather systems, weakening the depletion of boundary layer concentrations relative to the free troposphere in summer (Corbin and Denning, 2006).Alternatively, stronger vertical transport and mixing might accompany increased cloud cover and weaker photosynthetic uptake at some sites, which would magnify the depletion of boundary layer concentrations in summer.
Failure to adequately simulate synoptic rectifier effects could bias regional carbon sink estimates.Transport model inversions estimate surface CO 2 fluxes by adjusting modeled fluxes to minimize differences between observed and modeled concentrations.The timescale over which these differences are minimized determines which model processes are constrained by observations.Weekly timescales might be too coarse for optimization schemes to correct for model biases in synoptic rectifier effects.For example, the onset of vegetation water stress could occur over a period of a few days (e.g., Daly and Porporato, 2005), which would reduce carbon uptake and increase sensible relative to latent heating, resulting in deeper boundary layer heights and increased mixing with the free troposphere.
Evidence for biases in transport model inversions includes aircraft vertical profiling campaigns that found systematic biases in modeled vertical profiles of CO 2 , where models that more closely matched the observed profiles inferred weaker northern and stronger tropical land carbon uptake (Stephens et al., 2007).Furthermore, clear-sky biases in satellite observations will require transport models to assimilate concentration measurements with meteorology representative of the time and location of each remotely sensed sample (Corbin and Denning, 2006;Parazoo et al., 2011).Surface flux biases could therefore result from assimilating satellite data in transport models that fail to adequately resolve the synoptic covariation of surface fluxes and atmospheric transport.
Previous studies inferred synoptic rectifier effects indirectly, from large variability in concentrations at synoptic scales during field campaigns (e.g., Lin et al., 2006) and from synoptic storm systems passing over CO 2 observing sites (Hurwitz et al., 2004).Full transport model simulations have demonstrated the effects of covariation in transport and surface fluxes on CO 2 concentrations at seasonal timescales (Chan et al., 2008).Here we quantify the impact of synoptic timescale covariation on surface flux inversions by forcing a simplified transport model with surface and atmospheric data without covariation, or with covariation specified at each timescale between daily and seasonal.We show that synoptic rectifier effects significantly impact surface flux inversions in some years and are stronger in models than in observations over the Southern Great Plains (SGP), indicating the potential for further improvement in surface flux inversions through improved surface and boundary layer simulations, and higher-frequency data assimilation.

Stochastic boundary layer model for tracers (SBLM)
We developed a stochastic boundary layer model for tracers (SBLM) to estimate the synoptic rectifier effect, here defined as the component of CO 2 vertical concentration gradients due to the covariation of surface fluxes and atmospheric transport at synoptic timescales.We also used SBLM to quantify the effects of differences between observed and simulated synoptic rectifier effects on inverse surface flux estimates.SBLM consists of a tracer conservation equation for boundary layer CO 2 concentrations and a stochastic model for surface fluxes and meteorological variables entering the conservation equation.The stochastic component accounts for the effects of random but correlated noise in these variables on the synoptic rectifier effect.
The conservation equation is based on the well-mixed convective boundary layer approximation, and has appeared in similar forms as a diagnostic tool in previous studies of boundary layer CO 2 and other tracers (e.g., Betts, 1992;Styles et al., 2002;Williams et al., 2011).The boundary layer CO 2 mixing ratio is conserved according to where h, c f and c m , ρ, v, and ∇ are the boundary layer depth, free troposphere and boundary layer CO 2 mixing ratios, atmospheric molar density, horizontal wind velocity vector, and horizontal gradient operator, respectively; subscripts (f and m) denote quantities at the free-tropospheric level just above the boundary layer and averaged within the boundary layer, respectively.
The entrainment velocity (w e ) is given by where w is the resolved-scale vertical wind velocity evaluated at the mixed layer height (positive upward), and w cld is a sub-grid-scale subsidence velocity due to convective clouds (described further in Section 4.2).There is no entrainment into the boundary layer when the net movement of air is upward out of the boundary layer.The surface flux or net ecosystem exchange (F sfc , in Eq. 1, positive upward) balances the time rate of change of boundary layer CO 2 (first left-hand-side term), entrainment of free-tropospheric CO 2 due to the combination of boundary layer growth into the free troposphere and vertical advection by the subsiding flow (second term), and horizontal advection (third term).Rearranging Eq. ( 1) to obtain a prognostic equation for the vertical concentration gradient, where c = c m − c f (hereafter referred to as the vertical gradient) is evaluated at a single point in the horizontal, and and the entrainment rate is defined as The terms F , G, and E in Eq. ( 3) represent an external forcing similar to the forcing applied to single-column climate models for testing parameterization schemes.Those models are driven by time series of temperature and water vapor horizontal advective tendencies, as well as mass divergence (divergence is analogous to the entrainment term E in Eq. 3).We model F , E, and h as a linear system, which accounts for the covariance between these variables at each timescale (τ ), according to where is a matrix of impulse response functions relating the input N(t) to the output Y (t).The input is a vector of uncorrelated white-noise processes, randomly drawn from a Gaussian distribution with unitary variance, and the output   is a vector containing the time series we wish to model stochastically.
Since the impulse response matrix (q) completely describes any linear system in general (Jenkins and Watts, 1969), the problem of modeling the time series of F , E, and h reduces to finding the impulse response matrix that best approximates the system.The impulse response matrix is uniquely related to the cross-covariance matrix, given by with matrix conjugate transpose ( * ).The elements of R(τ ) are the cross-covariances of the variables in Y (t), which describe the covariances between each variable at each timescale (τ ).It is easier to transform the problem from the time to frequency domains using the Fourier transform, where the impulse response matrix becomes the frequency response matrix (Q).The Fourier transform of the crosscovariance matrix, known as the cross-spectral density matrix, is given by where f is the frequency (τ −1 ).The elements of S are the cross-spectral densities (abbreviated as CSD).Note that diagonal elements are the same as power spectral densities (PSD), and the matrix is Hermitian such that off-diagonal elements can be obtained from one another.Conceptually, the PSD distributes the variance as a function of frequency, whereas the CSD distributes the cross-covariance as a function of frequency.
The method of obtaining the impulse response matrix (q, or its Fourier transform, Q) from estimates of the crossspectral density (S) was derived by Ferraioli et al. (2010), and is outlined below (Eqs.5-7).Multiplying Eq. (4) by Y (t + τ ) and taking Fourier transforms yields where I is the unit matrix corresponding to the cross-spectral matrix of the input white-noise process (N).The eigendecomposition of the cross-spectral matrix is where V (f) and (f) are the eigenvector and eigenvalue matrices of S(f ), respectively.Since S is Hermitian and its eigenvector matrix is unitary, i.e., V (f )V * (f ) = I , Eqs. ( 5) and ( 6) can be combined to obtain the frequency response matrix where √ (f ) is a diagonal matrix whose elements are the square roots of those in (f ).
The impulse response matrix acts as a filter on the input white-noise time series (N ) to produce the stochastically modeled time series for F , E, and h according to Eq. ( 4).The numerical method to solve Eqs. ( 4)-( 7) (Ferraioli et al., 2010) requires only model cross-spectral densities, which can be estimated from transport model output or observations, as described in the following section.
We used SBLM to estimate the synoptic rectifier effect from observations and to compare modeled and observed rectifier effects.For example, deeper boundary layers often coincide with greater surface CO 2 uptake due to the effects of surface absorbed solar radiation on photosynthesis and surface sensible heating.The cross-spectral density of boundary layer heights and surface CO 2 fluxes can capture this covariation and its frequency dependence, and SBLM can be used to estimate the impact of this cross-spectral density on vertical concentration gradients.The ability to generate arbitrarily long synthetic time series with the same cross spectra as the observations or transport models enables us to detect statistically significant differences between observed and modeled synoptic rectifier effects.

Observations
Surface CO 2 fluxes were measured by eddy covariance on the 60 m tower of the US Southern Great Plains Atmospheric Radiation Measurement Climate Research Facility (36.61 • N, 97.49 • W; hereafter referred to as SGP) using a 3-D sonic anemometer and H 2 O and CO 2 densities from an infrared gas analyzer.H 2 O and CO 2 fluxes were calculated at 30 min intervals by block-averaging scalar quantities, rotating coordinates to zero mean vertical and meridional wind speed, and correcting for effects of covariance between air density and turbulent vertical wind on fluxes (Webb et al., 1980).H 2 O and CO 2 concentrations were obtained from a precision gas system (Bakwin et al., 1995) at 15 min intervals.More than 80% of the Southern Great Plains is managed for agriculture and grazing, with winter wheat from November through June over 40 % of all land to the southeast, and pasture (40%) and a mixture of C 3 and C 4 crops (20 %) from April through August (Cooley et al., 2005;Fischer et al., 2007;Riley et al., 2009).
Pressurized air samples have been collected approximately weekly in the free troposphere at SGP for subsequent analysis at NOAA/ESRL (Biraud et al., 2013).These flask-based measurements include, but are not limited to, CO 2 , CH 4 , CO, N 2 O, SF 6 , H 2 , and 13 C and 18 O in CO 2 .The sample collection procedures have been described in detail in Conway et al. (1994).Typically, air was pumped into glass flasks and slightly pressurized above ambient pressure.The freetropospheric concentrations used in this study were taken from the first available flask sample just above the boundary layer (samples were typically available at altitudes of 457.2, 609.6, 914, 1219.2, 1524, 1828.8, 2133.6, 2438.4, and 2743.2 m above ground, and higher, with an uncertainty of about 30 m).These measurements were used to construct a continuous monthly average of the c f tendency term in Eq. ( 3), onto which we added higher-frequency tendencies estimated from CarbonTracker (described in Sect.3.2 below).The effects of higher-frequency fluctuations in c f were assumed equal for observational and CarbonTracker estimates, as continuous daily observational estimates of c f were not available for the multi-year time period considered here.
Observational time series of boundary layer heights (h) were estimated from balloon sonde profiles of potential temperature measured four times daily (Heffter, 1980).Vertical velocities were obtained from the North American Regional Reanalysis (Mesinger et al., 2006) at 32 km horizontal resolution and 25 hPa vertical resolution every 3 h (NARR).Results using boundary layer heights from NARR were qualitatively similar to sonde-derived estimates (see Sect. 4.3, Fig. 6).

CarbonTracker (CT-TM5)
Measurement time series at SGP were compared to the 2010 release of CarbonTracker, a carbon data assimilation system (Peters et al., 2007).CarbonTracker assimilates surface CO 2 concentrations with modeled transport and surface fluxes in a global, two-way nested atmospheric transport model (Transport Model 5; Krol et al., 2005; hereafter TM5) driven by meteorological fields from the European Centre for Medium-Range Weather Forecasts (ECMWF).CarbonTracker data were averaged over the four horizontal grid points nearest to SGP, and include three-dimensional distributions of CO 2 concentrations and surface fluxes provided at 1 • × 1 • spatial resolution every 3 h.We used CT-TM5 concentrations averaged over the first two model levels above the boundary layer to define free-tropospheric concentrations.We hereafter refer to the CarbonTracker surface fluxes and ECMWF meteorological forcing collectively as the CT-TM5 data set.
Prior surface fluxes serve as surface CO 2 boundary conditions in TM5.These priors are provided by CASA, a biogeochemical model (Potter et al., 1993) that is driven by ECMWF meteorological fields, including temperature, precipitation, solar radiation, and vegetation specified by satellite normalized difference vegetation index (Randerson et al., 1997).TM5 integrates concentration fields forward in time with the surface boundary conditions specified above using three-dimensional winds from ECMWF to drive the transport.The effects of surface energy forcing on atmospheric boundary layer mixing are taken into account through the boundary layer diffusion scheme in TM5, which in turn depends on surface latent and sensible heat fluxes predicted by the land surface scheme underlying the ECMWF reanalysis data set (Viterbo and Beljaars, 1995).Transport and mixing by deep convection in TM5 is parameterized similarly to ECMWF.
The CarbonTracker products did not archive vertical velocities, so we used the ECMWF interim reanalysis vertical velocities (Dee et al., 2011) that are predicted based on the same general circulation model and parameterization schemes underlying transport in TM5.Differences between CT-TM5 and our recreated advective tendencies could arise due to interpolation of the ECMWF meteorological data to 1 • × 1 • resolution in the CT-TM5 system; however averaging CarbonTracker and ECMWF output over different numbers of grid cells centered over the study site did not significantly influence our results.
CarbonTracker uses a Kalman filter to optimize surface flux estimates with respect to available observations (Peters et al., 2007).The prior CASA fluxes are scaled by a linear factor updated each week to minimize differences between measured and modeled atmospheric CO 2 concentrations.Observational data assimilated in CT-TM5 currently include flask samples collected at surface sites as part of North American sampling networks and CO 2 time series from measurement towers.Modeled concentrations are compared to daytime-averaged measurements at SGP and other measurement towers, and these error estimates are temporally aggregated to a weekly timescale when estimating the CASA scaling factors.The scaled surface fluxes represent the assimilated surface fluxes in the final data product, often referred to as inverse surface flux estimates.

Horizontal advection
Horizontal advection was calculated from winds and spatial concentration gradients archived in CT-TM5 data products.The effect of horizontal advection was assumed to be identical for observations as in CT-TM5, since we did not have independent estimates of spatial concentration gradients from observations over the range of timescales considered here.We found that horizontal advection dampens fluctuations in vertical CO 2 gradients that are generated by synoptic covariation of surface fluxes and transport.Horizontal advection is therefore included in estimates of the net effects of synoptic covariation on vertical concentration gradients according to Eq. ( 3).Horizontal movement of the atmospheric boundary layer over the surface spatially averages the effects of heterogeneous surface fluxes on boundary layer concentrations (Gloor et al., 2001).More details on horizontal advection effects are provided at the end of Sect.4.2.

Model-data synthesis
Cross spectra were estimated for both observed and CT-TM5 time series using the Daniell spectral estimator, following standard methods (e.g., von Storch and Zwiers, 2004).Seasonality and trends were accounted for by separating the time series into spring (March-May), summer (June-August), autumn (September-November), and winter (December-February) seasons and subtracting a quadratic polynomial fitted to each season for each year.The spectra were tapered in the time domain using a cosine-bell weighting function.Our results were not sensitive to the choice of weighting function or the extent of smoothing using the Daniell window.Seasonal averages were taken over all years after estimating power and cross spectra for each season of each year separately.
Using the cross-spectral densities as input to the stochastic model described in Sect.2, we generated two sets of stochastically modeled time series for F , E, and h, one each for the CT-TM5 and observational data sets.We refer to these time series as the synthetic CT-TM5 and synthetic observations, and refer to the original time series as the original CT-TM5 and original observed time series.Note that the output time series of the stochastic model are Gaussian distributed due to the linearity assumption underlying Eq. ( 4).To account for possible effects of skewness or extreme values, we transformed the Gaussian distributions of the model output to match the original (univariate) distributions of F , E, and h using the inverse transform sampling method (see Supplement).Examples of the modeled and original time series are shown in Figs.S1, S2.
Equation (3) was integrated using finite differences by prescribing F , E, and h from the above linear system in order to obtain c.While the original time series are limited to the study period (8 yr), the synthetic time series can have arbitrary length (here spanning 300 yr).These long time series are useful in assessing statistically significant effects of weak to moderate covariation of F , E, and h on vertical concentration gradients.

Power and cross spectra of synthetic time series
Boundary layer height (h), entrainment (E), and CO 2 fluxes (F ) have enhanced power at lower frequencies, without pronounced spectral peaks (Fig. 1).There are broad peaks at about the 4-day timescale (frequencies between 0.2 and 0.3 cycles per day), particularly in power spectra of combined surface and horizontal CO 2 fluxes (F ) (red lines in Fig. 1), most notably during the autumn and winter months.Processes that can contribute to power in surface and horizontal CO 2 fluxes at these frequencies include synoptic weather variability such as changes in radiation, temperature, and humidity accompanying transient storm systems, in addition to horizontal transport of CO 2 .Lack of sharp synoptic spectral peaks is not surprising, as weather systems typically evolve from smaller to larger spatial scales over their life cycles, generating atmospheric variability over a broad range of timescales (e.g., Vallis et al., 2004).
The power (top row of Fig. 1) and coherence (middle row of Fig. 1) spectra of the synthetic time series closely match those of the original CT-TM5 time series (cf.dashed and solid lines).Each time series was normalized to have unit standard deviation (keeping the area under each PSD curve the same, preserving the shape of the curves) in order to compare spectra across all three variables.Spectra of the observed time series (Fig. 2, top row) have less power at low frequencies (i.e., shallower spectral slopes) than the CT-TM5 time series (Fig. 1, top row), particularly for surface and horizontal CO 2 fluxes (F ) and boundary layer height (h), during June-August, as discussed further in the following sections.
Cross spectra have both real (re) and imaginary (im) parts, shown as coherence and phase plots in the middle rows of Figs. 1 and 2. The (squared) coherence between two variables x and y is given by Coh which is analogous to the squared correlation coefficient if correlation were distributed as a function of frequency.We refer to squared coherence simply as coherence.The phase relation between the two variables x and y, Phase x,y (f ) = arctan re(CSD x,y (f )) im(CSD x,y (f )) , quantifies how much the crests of time series y lag the crests of time series x at each frequency in terms of a phase angle.For example, a phase angle of 0.5 cycles at a frequency of 0.1 cycles per day (period 10 days) indicates that crests of x are aligned with troughs of y and that crests in y follow troughs in x after 5 days.Note that phase angles of −0.5 and +0.5 cannot generally be distinguished.
The CT-TM5 time series are only weakly coherent at most frequencies (middle row of Fig. 1), with the exception of coherence between entrainment (E) and boundary layer height (h), which turns out to have a relatively small impact on vertical concentration gradients (shown in the following section).Both Coh F,h (red lines in Fig. 1, middle row) and Coh F,E (green lines in Fig. 1, middle row) have maxima at frequencies around 0.25 • 10 −1 cycles per day (period of 4 days) in autumn and winter, consistent with the synoptic timescale.Coherence is weaker in observations compared to CT-TM5 (cf.middle rows of Figs. 1, 2), suggesting that the dynamical models underlying CT-TM5 are overestimating the coupling of atmospheric dynamics and surface fluxes at this site.Spectra shown in Figs. 1 and 2 were estimated for time series extending from January 2002 to January 2010 for CT-TM5 and January 2003 to January 2010 for observations.We found that coherence varies from year to year, which we explore in Sect.4.3.Phase angles (lower panels of Figs. 1, 2) are around 0.5 cycles at frequencies where coherence and power are greatest (e.g., 0.1 cycles per day) and are similar between both observed and CT-TM5 data sets.CO 2 flux into the boundary layer (F ) is most positive when boundary layer height fluctuations are most negative, indicating that the shallowest daytime boundary layer heights coincide with reduced photosynthetic uptake (enhanced respiratory release).The strongest entrainment rates (most negative E) coincide with enhanced photosynthetic uptake and reduced respiratory release.These results are consistent with stronger solar heating at the surface driving enhanced photosynthesis, deeper mixing, and stronger vertical transport.

Synoptic rectifier effect
Long simulations with the stochastic boundary layer model (SBLM) were designed to detect responses in vertical concentration gradients to covariance between the external forcing variables F , E, and h.We performed 300 statistically independent integrations of SBLM for each season, using synthetic time series for F , E, and h as external forcing for each simulation (using Eq. 3).We tested the SBLM capability of predicting vertical concentration gradients by comparing power spectra of original CT-TM5 concentration gradients (i.e., those predicted by the CarbonTracker data assimilation system) to those predicted by SBLM forced with synthetic CT-TM5 time series of F , E, and h.SBLM (blue dashed lines in Fig. 3) and CT-TM5 (red lines in Fig. 3) power spectra were consistent after accounting for uncertainty in the entrainment rates diagnosed using Eq. ( 2), as shown by the range of gray shading in Fig. 3. Concentration gradient power spectra are shown in Fig. 3 in the form of log-log plots due to their larger range of values.
Without sub-grid-scale entrainment, entrainment of freetropospheric air into the boundary layer would occur only during resolved-scale subsidence or when boundary layer growth exceeds the rate of resolved-scale ascent.In models, however, resolved-scale ascent typically coincides with moist convection (Lawrence and Salzmann, 2008), which brings free-tropospheric air into the boundary layer through parameterized mixing intended to represent the sub-gridscale compensating subsidence surrounding cumulus clouds, and the effects of rainy, evaporatively cooled downdrafts.These parameterized mass fluxes are not archived in model or reanalysis data products and occur at smaller scales than represented by resolved vertical velocities in Eq. ( 2).

CT-TM5
SBLM with synthetic CT-TM5 forcing CO 2 Vertical Gradient Power Spectral Density Fig. 3. Normalized power spectral densities for boundary layer CO 2 concentrations predicted by the stochastic boundary layer model (SBLM) with synthetic CT-TM5 forcing (dashed blue), and spectral densities for the original boundary layer CO 2 concentrations from the CT-TM5 data assimilation system (red).The gray shading indicates the range of SBLM solutions for 0 < f < 1.25, where f is an estimated sub-grid-scale boundary layer entrainment expressed as a fraction of the resolved-scale vertical velocity (see text for details).
We accounted for the effects of boundary layer entrainment due to sub-grid-scale cloud mass fluxes by modifying the entrainment rate with a parameter f , expressed as a fraction of the resolved-scale vertical velocity.This parameter is a function of time and takes on the values We added this term to the entrainment velocity by defining the sub-grid-scale entrainment term in equation ( 2) as under the assumption that sub-grid-scale compensating subsidence scales (in magnitude) with the resolved-scale ascent (Lawrence and Salzmann, 2008).Our entrainment parameterization is similar in form to others (e.g., Stull, 1988;their Eq. 1.2.2.b), and reduces to the equilibrium boundary layer model in Bakwin et al. (2004) for f o = 1 (i.e., upward convective mass fluxes are exactly compensated for by subsidence).
We explored a range of solutions for entrainment parameters ranging from f o = 0 (no compensating subsidence) Fig. 4. Seasonally averaged vertical CO 2 gradients (c m − c f ) predicted by SBLM with synthetic CT-TM5 forcing (dashed blue, corresponding to dashed blue lines in Fig. 3), and the original vertical CO 2 gradient from the CT-TM5 data assimilation system (red).Gray shading indicates the range of solutions for a range of sub-grid-scale boundary layer entrainment (as in Fig. 3); gradients monotonically decrease in magnitude with stronger entrainment.Error bars indicate the maxima and minima of seasonally averaged gradients for the 8 yr period.
to f o = 1.25 (compensating subsidence 25 % greater than the resolved-scale ascent), as shown by the gray shading in Fig. 3. SBLM correctly predicted the original CT-TM5 spectra for f o ≈ 1 (dashed blue lines of Fig. 3).The distribution of power shifts toward higher frequencies (i.e., shallower spectral slopes) as f o increases (gray shading in top panels of Fig. 3), reflecting the tendency for resolved-scale ascent to occur over short time intervals associated with transient deep convection.
Seasonal concentration gradients in SBLM (dashed blue lines in Fig. 4) agree with the original CT-TM5 gradients (red line in Fig. 4) for the range of entrainment rates explored here.Boundary layer CO 2 becomes depleted in summer (June, July, and August), slightly enhanced in spring and autumn, and significantly enhanced in winter.The magnitude of the peak-to-trough seasonal cycle in vertical concentration gradients at SGP (3.4 ppm CO 2 ) is smaller than the average of the northern hemispheric sites reported at other continental sites across the Northern Hemisphere (4.8 ppm CO 2 , Stephens et al., 2007).The smaller seasonal cycle amplitude at SGP is not an error in CT-TM5 but rather a result of the observed seasonality of heterogeneous fluxes at the SGP site, where peak surface carbon uptake from winter wheat (typically in April) is not coincident with peak uptake from pasture (typically in June).
We experimented with different forcing time series to isolate the role of synoptic variability in generating vertical  (d, e, f) shows separate contributions to the rectified gradients by each bivariate cross spectrum: CSD F,E (green), CSD F,h (red), and CSD E,H (blue).Rectified gradients are defined relative to control gradient given by the solution to SBLM forced with synthetic time series that all have zero cross spectra.Error bars show the range of solutions for 25 % weaker and 25 % stronger sub-grid-scale entrainment (see text for details).
gradients.We quantified the effects of covariation of surface fluxes and boundary layer entrainment (as represented by CSD F,E ) by setting the other two cross-spectral densities (CSD F,h and CSD E,h ) to zero and generating synthetic forcing time series that only have covariation between surface fluxes and boundary layer entrainment (Fig. 5d, green lines).We similarly isolated the effects of CSD F,h (Fig. 5d, red lines) and CSD E,h (Fig. 5d, blue lines) for both the CT-TM5 (dashed lines) and observational time series (solid lines).The total synoptic rectifier effect due to all three covariance terms is shown in Fig. 5a.Uncertainty in cloud mass fluxes did not significantly change our conclusions about the relative strengths of each rectifier effect, as shown by varying (f ) between 0.75 and 1.25 (error bars in Fig. 5).Effects of distribution skewness (i.e., non-Gaussian distributions in forcing time series) were also tested but found to have little impact (results not shown), meaning that synoptic rectifier effects were not dominated by rare extreme events in CO 2 fluxes or boundary layer depths.
We expected entrainment and boundary layer depths to covary because the boundary layer grows by entrainment, such that deeper boundary layers tend to entrain more freetropospheric air.Also, the subsidence velocity is typically an increasing function of height in the lower troposphere (Williams et al., 2011), so a deeper boundary layer in equilibrium with the subsiding flow will tend to experience greater entrainment rates.Boundary layer concentrations were enhanced by covariation between surface fluxes and boundary layer height (0.2 ppm CO 2 in summer), but less so by covariation between entrainment and boundary layer height (0.1 ppm CO 2 in summer), shown in Fig. 5d (dashed red and blue lines, respectively).On the other hand, boundary layer concentrations were depleted by covariation of surface fluxes and entrainment (−0.1 ppm CO 2 in summer, dashed green lines in Fig. 5d).These opposing effects resulted in a net enhancement of about 0.13 ppm in boundary layer CO 2 in summer (Fig. 5a, JJA), or 9 % of the mean vertical gradient in summer at SGP, for the synthetic CT-TM5 forcing time series.Synoptic rectifier effects in other seasons were generally smaller than in summer.
Horizontal mixing and transport could suppress spatially localized synoptic rectifier effects over heterogeneous land surfaces such as SGP, because the covariation of atmospheric dynamics and surface CO 2 fluxes can depend on vegetation type.Since F includes both surface fluxes and horizontal advection, we ran separate simulations with horizontal advection removed and found stronger rectifier effects (0.4 ppm in summer, not shown), indicating that our results represent conservative estimates for homogeneous land surfaces where horizontal advection would likely play a smaller role.We leave analysis of spatial variability to future work and focus here on the important impact of time series non-stationarity due to interannual variability.

Interannual variability
Precipitation extremes can have a strong influence on ecosystem productivity and land-atmosphere coupling over the mostly non-irrigated Southern Great Plains.We repeated the above analyses for two separate years, 2006 and 2007, to quantify effects of non-stationarity resulting from unusually dry or wet summers.Synthetic forcing time series were fitted separately to cross spectra for 2006, the second driest year in Oklahoma, and 2007, the wettest year on record for central Oklahoma (Dong et al., 2011).The results indicate much larger rectifier effects in summer 2007, with synoptic variability generating vertical gradients of 0.3 ppm, or about 20 % of the average seasonal cycle amplitude in CT-TM5 (Fig. 5c, dashed lines), compared to negligible synoptic rectifier effects in summer 2006 (Fig. 5b, dashed lines).Covariation of surface fluxes and boundary layer heights explains most of the 0.3 ppm CO 2 synoptic rectifier effect during the wet year (Fig. 5f, dashed red lines).Observed rectifier effects were weaker than in CT-TM5, particularly in the summer (cf.black and dashed lines in Fig. 5c).Additionally, synoptic covariation depleted observed boundary layer concentrations in observations but enhanced CT-TM5 concentrations during spring of 2007.Results were similar when comparing CT-TM5 to observations using NARR boundary layer depths as opposed to sonde-based estimates (Fig. 6).These results suggest stronger rectifier effects in CT-TM5 than observed due to stronger covariation between PBL depth and surface fluxes, indicating that land-atmosphere coupling may be too strong in the land surface and atmospheric models used in CT-TM5.
Cross spectra explain why synoptic rectifier effects are stronger than average during the summer months of 2007 in CT-TM5 (Fig. 7).Coherence between surface fluxes and boundary layer depths is much larger in 2007 than over all years combined (cf.middle rows of Figs. 1 and 7, for Jun-Aug), particularly at lower frequencies corresponding to periods of about 10-20 days (thick lines in Fig. 7).The 2007 coherence spectra were separated into contributions from low (thick lines in Fig. 7) and high (dashed-dotted lines) frequencies with a separation frequency of 0.1 cycles per day (period of 10 days).We separated the effects of low-frequency variability by weighting the cross spectra by a window that varies from 100% for frequencies lower than 0.09 cycles per day (meaning that 100 % of the coherence in the original time series is retained) to 1 % for frequencies above 0.11, with a smooth taper in between (split cosine bell).We repeated this procedure for high frequencies (retaining 100 % of coherence greater than 0.11 cycles per day).Phase angles of 0.5 cycles at these lower frequencies mean that high PBL depth coincides with negative surface CO 2 flux (i.e.land uptake).Power spectra show enhanced variability in CO 2 flux at low frequencies, which also contributes to the larger rectifier effect in 2007 (cf.top panels of Figs. 1, 7), particularly for the spring and summer months.
We ran separate SBLM simulations with synthetic forcing time series that only have coherence at high or low frequencies, corresponding to the thick and dashed-dotted lines in Fig. 7 (middle row).Integrating SBLM for the high-and lowfrequency forcing separately shows that the 2007 summer synoptic rectifier effect in CT-TM5 resulted primarily from coherence between PBL depth and surface CO 2 fluxes at low frequencies having periods greater than 10 days (Fig. 8, solid line with downward-arrow markers).The higher frequencies account for one-third of the synoptic rectifier effect (Fig. 8, dashed-dotted line with upward-arrow markers).These results indicate the potential for improvement in surface flux estimates through assimilation of concentration data at higher frequencies, considering that periodicities of , CSD F,h (red), and CSD E,h (blue), separated into contributions from low frequencies (thick lines) and high frequencies (dasheddotted lines), with a separation timescale of 10 days (0.1 cycles per day).Bottom row: phase spectra corresponding to the cross-spectral densities CSD F,E (green), CSD F,h (red), and CSD E,h (blue).The horizontal black line (middle row) indicates the 90% confidence level for coherence.
10 days require sampling timescales of 5 days or shorter, and data assimilation systems currently update prior fluxes on weekly timescales.

Impact on surface flux inversions
Biased synoptic rectifier effects could lead to biased surface CO 2 flux estimates in carbon cycle data assimilation systems depending on how accurately the underlying prior flux models represent the true fluxes.We addressed this question by integrating SBLM in data-assimilation mode, using a simplified form of the Kalman filter used in transport model inversions, where we reduced the state vector to a scalar value for the spring surface flux at SGP. Effects of horizontal advection were included in SBLM as described in Sect.3.3.We used synthetic CT-TM5 entrainment rates, with synthetic Carbon-Tracker surface CO 2 fluxes being a stand-in for a prior flux model.The solution to SBLM consistent with this synthetic forcing would represent a measurement time series in the hypothetical case of perfect transport and prior fluxes, whereas here we subtract 0.2 ppm from this solution to simulate a +0.2 ppm bias in rectifier effects corresponding to the difference between CT-TM5 and observations from Fig. 5 for March-May (the observed period of peak net CO 2 uptake from winter wheat at SGP).Table 1 displays the results of the data assimilation tests for two different assumptions of prior flux uncertainty.A +0.2 ppm bias in the synoptic rectifier effect leads to a −0.15 µmol m −2 s −1 surface flux bias when the model is constrained to closely match observed concentrations by specifying a large prior flux uncertainty of one standard deviation in F (denoted P 1 in Table 1).This negative flux bias allows SBLM to reduce boundary layer concentrations relative to the free troposphere in order to better match the synthetic measurement time series.These biases represent 13 % of the peak-to-trough seasonal cycle amplitude at SGP.Though small relative to the seasonal cycle, the corresponding −0.15 µmol m −2 s −1 surface CO 2 flux biases compare in magnitude to interannual variability, as estimated from CT-TM5 standard deviations of springsurface CO 2 fluxes from 2003 to 2010 (0.45 µmol m −2 s −1 ).Reducing the prior flux uncertainties to 1 % of the standard deviation in F recovers the expected result that biases in concentration vertical gradients have no effect on surface flux estimates when constrained to match perfect prior fluxes (P 2 in Table 1).In that case, the 0.2 ppm CO 2 bias is left uncorrected in the assimilated data and the surface flux is perfectly estimated.This example illustrates why it can be difficult to separate the effects of biases in transport from biases in prior fluxes when comparing modeled and observed concentration profiles.

Discussion
The synoptic covariation examined here results in vertical concentration gradients, leaving the whole-column CO 2 concentration unchanged.Whole-column CO 2 retrievals from GOSAT and OCO-2 satellites (Crisp et al., 2004;Yokota et al., 2009) and ground-based spectrometers (Washenfelder et al., 2006) could therefore be used to discriminate between true surface fluxes and apparent surface fluxes due to vertical covariation, as true surface fluxes generate tendencies in the total mass of CO 2 in a column (Keppel-Aleks et al., 2011).Tendencies in whole-column CO 2 can also occur downwind of deep convection, since trace gases are transported upward through convective clouds and horizontally through convective outflows in the upper troposphere (Cotton et al., 1995).Aircraft vertical profile measurements in combination with global column CO 2 retrievals would therefore be useful in further constraining uncertainties in convective parameterizations and their effects on the three-dimensional convective transport in transport model inversions.
Convective parameterizations are not able to account for boundary layer thermals that penetrate the stable inversion layer (Hohenegger et al., 2009), nor are they able to account for the smooth transition between shallow and deep convection observed over diurnal to synoptic timescales (Zhang and Klein, 2010).These errors may be responsible for the drift of general circulation models into unrealistically per-Table 1. Spring (MAM) biases in inverse surface CO 2 flux estimates for two prior flux uncertainties (P 1 and P 2 ) with imposed 0.2 ppm bias in simulated vertical gradients.The P 1 case is tightly constrained to the concentration measurements, whereas the P 2 case is tightly constrained to the perfectly specified prior flux.Prior flux uncertainties are scaled according to σ , the standard deviation of F .sistent dry states (Klein et al., 2006).Errors in the simulation of cloud cover and depth are coupled to errors in simulated boundary layer heights, as deep convection can be triggered by higher surface sensible heat fluxes and deeper boundary layers (Pielke, 2001).Moist convection also sets in too early in the diurnal cycle in general circulation models, which could result in vertical transport of CO 2 that is more representative of nighttime respiration than afternoon photosynthesis.These errors are constrained by meteorological measurements assimilated in the reanalysis data products used to drive transport model inversions, but some biases in convection timing remain in these reanalysis products persist (Betts et al., 2009).Our results demonstrate that biases in the representation of diurnal to synoptic timescale boundary layer process can impact surface flux estimates on much longer, seasonal and interannual timescales.

Conclusions
The synoptic covariation of transport and surface fluxes generates vertical concentration gradients that atmospheric transport model inversions could falsely interpret as surface sources and sinks if not adequately represented in transport models.Comparisons of CarbonTracker to observations indicate that, based on offline simulations in a simplified transport model (SBLM), biases in synoptic covariation can result in surface flux biases of 13% of the seasonal cycle amplitude at SGP.Our results represent conservative estimates for homogeneous land surface types where horizontal advection would likely play a smaller role than at SGP. Two-thirds of the synoptic covariation of surface fluxes and transport occurs at periodicities greater than 10 days, suggesting that surface flux inversions would benefit from improved model representation of dynamics at the lower-frequency end of the synoptic spectrum, including soil moisture and leaf area index and their coupling to atmospheric transport.Interannual variability can provide useful analogs of climate change for observationally testing carbon cycle models, such as comparing drought effects on simulated and observed carbon fluxes.Our results demonstrate that the strength of synoptic rectifier effects varies interannually -with a drought year (2006) having negligible synoptic rectifier effects, and a wet year (2007) having large synoptic rectifier effects during the growing season -due only to differences in the covariation of surface fluxes and transport between these years in an atmospheric transport model inversion.This interannual variability in vertical gradients due to synoptic rectifier effects is of the same magnitude as that due to surface carbon sinks alone.Quantifying interannual variability in carbon sinks will therefore require transport models that accurately resolve the covariation of surface fluxes and transport at synoptic timescales.
Based on applications to the Southern Great Plains, we conclude that carbon data assimilation systems can be improved through better representation of synoptic-scale landatmosphere covariation.There remains potential for further observational constraints on regional CO 2 sources and sinks by optimizing inverse surface flux estimates on shorter timescales using daily CO 2 concentration measurements and by using remotely sensed soil moisture measurements to improve representation of high-frequency surface fluxes in the land surface model components of these inversions.These improvements will help meet the challenges of testing carbon cycle models with limited observational data at regional scales.

Fig. 1 .
Fig.1.Top row: power spectral density estimates for entrainment (E, green lines), surface and horizontal CO 2 flux (F , red lines), and boundary layer height (h, blue lines) for original CT-TM5 time series (solid) and synthetic CT-TM5 time series (dashed) for each season, normalized in order to have equal area under each curve (equal standard deviation).Middle row: coherence spectra corresponding to the cross-spectral densities CSD F,E (green), CSD F,h (red), and CSD E,h (blue).Bottom row: phase spectra corresponding to the cross-spectral densities CSD F,E (green), CSD F,h (red), and CSD E,h (blue).

Fig. 5 .
Fig. 5. Top row (a, b, c): rectified vertical CO 2 gradients (c m − c f ) predicted by SBLM forced with synthetic observations (solid), and synthetic CT-TM5 (dashed) with multivariate cross spectra fitted to the original time series, as shown in Figs.(1, 2).The stochastic model was fitted to cross spectra estimated over years 2002-2010 (a), and for years 2006 (b) and 2007 (c).The bottom row(d, e, f) shows separate contributions to the rectified gradients by each bivariate cross spectrum: CSD F,E (green), CSD F,h (red), and CSD E,H (blue).Rectified gradients are defined relative to control gradient given by the solution to SBLM forced with synthetic time series that all have zero cross spectra.Error bars show the range of solutions for 25 % weaker and 25 % stronger sub-grid-scale entrainment (see text for details).

Fig. 6 .
Fig. 6.Rectified vertical CO 2 gradients estimated from observations as in (c) of Fig. 5 but using NARR boundary layer depths (red lines) in addition to sonde-based boundary layer depths (solid black).CT-TM5 rectified vertical CO 2 gradients (dashed black) are shown as in Fig. 5c for comparison.

Fig. 7 .
Fig. 7.Top row: power spectral densities of E (green), F (red), and h (blue), estimated from CT-TM5 forcing time series for year 2007, normalized in order to have equal area under each curve.Middle row: coherence spectra corresponding to the cross-spectral densities CSD F,E (green), CSD F,h (red), and CSD E,h (blue), separated into contributions from low frequencies (thick lines) and high frequencies (dasheddotted lines), with a separation timescale of 10 days (0.1 cycles per day).Bottom row: phase spectra corresponding to the cross-spectral densities CSD F,E (green), CSD F,h (red), and CSD E,h (blue).The horizontal black line (middle row) indicates the 90% confidence level for coherence.

Fig. 8 .
Fig. 8. Rectified vertical CO 2 gradients predicted by SBLM with synthetic forcing fitted to modified cross-spectral densities having the coherence spectra shown in Fig. 7.The modified forcing time series have coherence spectra separated into low frequencies (downward triangles) or high frequencies (upward triangles).The unmodified forcing time series is shown for reference (dashed line, solid squares) and corresponds to the dashed line in Fig. 5c.Error bars indicate the standard error of the mean over 300 independent realizations of each season.