Decadal variation in CO₂ fluxes and its budget in a wheat and maize rotation cropland over the North China Plain

Carbon sequestration in agroecosystems has great potential to mitigate global greenhouse gas emissions. To assess the decadal trend of CO₂ fluxes of an irrigated wheat–maize rotation cropland over the North China Plain, the net ecosystem exchange (NEE) with the atmosphere was measured by using an eddy covariance system from 2005 to 2016. To evaluate the detailed CO₂ budget components of this representative cropland, a comprehensive experiment was conducted in the full 2010–2011 wheat–maize rotation cycle by combining the eddy covariance NEE measurements, plant carbon storage samples, and a soil respiration experiment that differentiated between heterotrophic and below-ground autotrophic respirations. Over the past decade (from 2005 to 2016), the cropland exhibited a statistically nonsignificant decreasing carbon sequestration capacity; the average of total NEE, gross primary productivity (GPP), and ecosystem respiration (ER), respectively, were −364, 1174, and 810 gC m⁻² for wheat and −136, 1008, and 872 gC m⁻² for maize. The multiple regression revealed that air temperature and groundwater depth showed pronounced correlations with the CO₂ fluxes for wheat. However, in the maize season, incoming shortwave radiation and groundwater depth showed pronounced correlations with CO₂ fluxes. For the full 2010–2011 agricultural cycle, the CO₂ fluxes for wheat and maize were as follows: for NEE they were −438 and −239 gC m⁻², for GPP 1078 and 780 gC m⁻², for ER 640 and 541 gC m⁻², for soil heterotrophic respiration 377 and 292 gC m⁻², for below-ground autotrophic respiration 136 and 115 gC m⁻², and for above-ground autotrophic respiration 128 and 133 gC m⁻². The net biome productivity was 59 gC m⁻² for wheat and 5 gC m⁻² for maize, indicating that wheat was a weak CO₂ sink and maize was close to CO₂ neutral to the atmosphere for this agricultural cycle. However, when considering the total CO₂ loss in the fallow period, the net biome productivity was −40 gC m⁻² yr⁻¹ for the full 2010–2011 cycle, implying that the cropland was a weak CO₂ source. The investigations of this study showed that taking cropland as a climate change mitigation tool is challenging and that further studies are required for the CO₂ sequestration potential of croplands.

1 Introduction

The widely used eddy covariance technique (Aubinet et al., 2000; Baldocchi et al., 2001; Falge et al., 2002a, b) has enabled us to better understand the terrestrial CO₂ exchange with the atmosphere, and has thereby fostered our understanding of the mechanisms through which terrestrial ecosystems contribute to mitigating ongoing climate change (Falkowski et al., 2000; Gray et al., 2014; Poulter et al., 2014; Forkel et al., 2016). Agroecosystems play an important role in regulating the global carbon balance (Lal, 2001; Bondeau et al., 2007; Özdoğan, 2011; Taylor et al., 2013; Gray et al., 2014) and are believed to have great potential to mitigate global carbon emissions through cropland management (Sauerbeck, 2001; Freibauer et al., 2004; Smith, 2004; Hutchinson et al., 2007; van Wesemael et al., 2010; Ciais et al., 2011; Schmidt et al., 2012). Furthermore, some studies proposed using agroecosystems as “natural climate solutions” to mitigate global carbon emissions (e.g., Griscom et al., 2017; Fargione et al., 2018). Field management practices (e.g., irrigation, fertilization and residue removal, etc.) impact the cropland CO₂ fluxes (Baker and Griffis, 2005; Béziat et al., 2009; Ceschia et al., 2010; Eugster et al., 2010; Drewniak et al., 2015; de la Motte et al., 2016; Hunt et al., 2016; Vick et al., 2016), but their relative importance in determining the cropland CO₂ budget remain unclear because of limited field observations (Kutsch et al., 2010), motivating comprehensive CO₂ budget assessments across different cropland management styles.

Over the past 2 decades, CO₂ investigations of agroecosystems have mainly focused on the variations in the net ecosystem exchange with the atmosphere (i.e., net ecosystem exchange, NEE) or its two derived components (i.e., gross primary productivity, GPP, and ecosystem respiration, ER) using the eddy covariance method. To date, these evaluations have been widely conducted for wheat (Gilmanov et al., 2003; Anthoni et al., 2004a; Moureaux et al., 2008; Béziat et al., 2009; Vick et al., 2016), maize (Verma et al., 2005), sugar beet (Aubinet et al., 2000; Moureaux et al., 2006), potato (Anthoni et al., 2004b; Fleisher et al., 2008), soybean–maize rotation cropland (Gilmanov et al., 2003; Hollinger et al., 2005; Suyker et al., 2005; Verma et al., 2005; Grant et al., 2007) and winter wheat–summer maize cropland (Zhang et al., 2008; Lei and Yang, 2010). However, the long-term variations in the cropland CO₂ fluxes remain limited, leaving our knowledge of cropland potential as a future climate change mitigation tool incomplete.

The widely used eddy covariance technique has fostered our understanding of the integrated fluxes of NEE, GPP and ER but cannot provide detailed CO₂ budget components, which consist of carbon assimilation (i.e., GPP), soil heterotrophic respiration (R_H), above-ground autotrophic respiration (R_AA), below-ground autotrophic respiration (R_AB), lateral carbon export at harvest, and import at sowing or through organic fertilization (Ceschia et al., 2010). These different CO₂ components result from different biological and biophysical processes (Moureaux et al., 2008) that may respond differently to climatic conditions, environmental factors and management strategies (Ekblad et al., 2005; Zhang et al., 2013). Differentiating among these components is a prerequisite for understanding the response of terrestrial ecosystems to changing environment (Heimann and Reichstein, 2008), thus the carbon budget evaluations have been reported for a few croplands (e.g., Moureaux et al., 2008; Ceschia et al., 2010; Wang et al., 2015; Demyan et al., 2016; Gao et al., 2017). In particular, to account for the literal carbon export, the net biome productivity (NBP) is often estimated by combining the eddy covariance technique and field carbon measurements associated with harvests and residue treatments (Ceschia et al., 2010; Kutsch et al., 2010). As a detailed CO₂ budget might facilitate better predictions of agroecosystems' responses to climate change, CO₂ budget evaluations in different croplands remain necessary.

The North China Plain (NCP) is one of the most important food production regions in China, and it guarantees national food security by providing more than 50 % and 33 % of the nation's wheat and maize, respectively (Kendy et al., 2003). Irrigation by diverting water from the Yellow River is common to alleviate water stress during spring in the NCP, resulting in a very shallow groundwater depth (usually range from 2 to 4 m) along the Yellow River (Cao et al., 2016) (Fig. 1). Wang et al. (2015) suggested that groundwater-fed cropland in the NCP had been losing carbon, and other studies also reported croplands in this region as carbon sources (e.g., Li et al., 2006; Luo et al., 2008). However, the long-term variations (e.g., >10 years) of the CO₂ fluxes over the NCP remain lacking, leaving the trend of carbon sequestration capacity of this region unknown.

Figure 1 Location of the experimental site. The background is the shallow groundwater depth in early September of 2011 provided by the © Water Information Center in the Ministry of Water Resources, P.R. China.

To this end, this study is designed to assess the long-term variation in CO₂ fluxes and its budget of the representative wheat–maize rotation cropland in the NCP. The eddy covariance system was used to measure the CO₂ exchange from 2005 to 2016. For the full 2010–2011 agricultural cycle, we measured soil respiration and sampled crops to quantify the detailed CO₂ budget components. These measurements allow us to (1) investigate the decadal CO₂ flux (NEE, GPP and ER) trend over this cropland; (2) provide detailed CO₂ budget components; and (3) estimate the net primary productivity (NPP), net ecosystem productivity (NEP), and NBP.

2 Materials and methods

2.1 Site description and field management

The experiment was conducted in a rectangular-shaped (460 m × 280 m) field of the representative cropland over the NCP (36^∘39^′ N, 116^∘03^′ E, Weishan site of Tsinghua University, Fig. 1). The soil is silt loam with a field capacity of 0.33 m³ m⁻³ and saturation point of 0.45 m³ m⁻³ for the top 5 cm of the soil. The mean annual precipitation is 532 mm and the mean air temperature is +13.3 ^∘C. The winter wheat–summer maize rotation system is the representative cropping style in this region. On average, the winter wheat is sown around 17 October and harvested around 16 June of the following year with crop residues left on the field; summer maize is sown following the wheat harvest around 17 June and harvested around 16 October. Prior to sowing wheat of the next season, the field is thoroughly plowed to fully incorporate maize residues into the top 20 cm of the soil. The canopies of both wheat and maize are very uniform across the whole season. Nitrogen fertilizer is commonly applied at this site with the amount being 35 gN m⁻² for wheat and 20 gN m⁻² for maize. The crop density is 775 plants per square meter for wheat with a ridge spacing of 0.26 m and 4.9 plants per square meter for maize with a ridge spacing of 0.63 m on average. Wheat is commonly irrigated with water diverted from the Yellow River and the irrigation is about 150 mm every year; maize is rarely irrigated because of the high precipitation in the summer. During the 2010–2011 agricultural cycle, when CO₂ budget components were evaluated, winter wheat was sown on 23 October 2010 and subsequently harvested on 10 June 2011 and summer maize was sown on 23 June 2011 and harvested on 30 September 2011. The entire year from 23 October 2010 to 22 October 2011 was studied for the annual CO₂ budget evaluation.

2.2 Eddy covariance measurements

A flux tower was set up at the center of the experiment field in 2005 (Lei and Yang, 2010; Zhang et al., 2013). The NEE was measured at 3.7 m above ground with an eddy covariance system consisting of an infrared gas analyzer (LI-7500, LI-COR Inc., Lincoln, NE, USA) and a three-dimensional sonic anemometer (CSAT3, Campbell Scientific Inc., Logan, UT, USA). The 30 min averaged NEE was calculated from the 10 Hz raw measurements with TK2 (Mauder and Foken, 2004) from 2005 to 2012 and TK3 software package (Mauder and Foken, 2011) from 2013 to 2016. The storage flux was calculated by assuming a constant CO₂ concentration profile. Nighttime measurements under stable atmospheric conditions with a friction velocity lower than 0.1 m s⁻¹ were removed from the analysis (Lei and Yang, 2010). In the gap-filling procedure, gaps less than 2 h were filled using linear regression, while other short gaps were filled using the Mean Diurnal Variation (MDV) method (Falge et al., 2001); gaps longer than 4 weeks were not filled. NEE was further partitioned to derive GPP and ER using the nighttime method (Reichstein et al., 2005; Lei and Yang, 2010), which assumes that daytime and nighttime ER follow the same temperature response, which thereby estimates the daytime ER using the regression model derived from the nighttime measurements. In particular, this study adopted the method proposed by Reichstein et al. (2005) to quantify the short-term temperature sensitivity of ER from nighttime measurements as described by the van 't Hoff equation,

$$\begin{array}{}\text{(1)}& \mathrm{ER}={\mathrm{ER}}_{\mathrm{ref}}\exp \left(b{T}_{\mathrm{S}}\right),\end{array}$$

where T_S is soil temperature, ER_ref is the reference respiration at 0 ^∘C and b is a parameter associated with the commonly used temperature sensitivity coefficient Q₁₀,

$$\begin{array}{}\text{(2)}& Q_{\mathrm{10}}=\exp \left(\mathrm{10}b\right).\end{array}$$

The long-term temperature sensitivity b of the season (either wheat or maize) was determined by averaging all the estimated short-term b in each of the 4 d windows with the inverse of the standard error as a weighing factor. The long-term temperature sensitivity b was then used to estimate the ER_ref parameter in each of the 4 d windows by fitting Eq. (1). Following this, ER_ref of each day was estimated by using the least-squares spline approximation (Lei and Yang, 2010).

To quantify the contribution of source areas to the CO₂ flux measurement of the eddy covariance, we used an analytical footprint model (Hsieh et al., 2000),

$$\begin{array}{}\text{(3)}& f(\mathit{\chi },z_{\mathrm{m}})={\displaystyle \frac{\mathrm{1}}{{\mathit{\kappa }}^{\mathrm{2}}{\mathit{\chi }}^{\mathrm{2}}}}D{z}_u^P{\left|L\right|}^{\mathrm{1}-P}\exp \left({\displaystyle \frac{-\mathrm{1}}{{\mathit{\kappa }}^{\mathrm{2}}\mathit{\chi }}}D{z}_u^P{\left|L\right|}^{\mathrm{1}-P}\right),\end{array}$$

where D=0.28 and P=0.59 are similarity constants for unstable condition (Hsieh et al., 2000), κ=0.4 is von Karman constant, χ represents the horizontal coordinate, L represents the Obukhov length, z_m represents the measurement height, and z_u represents the length scale expressed as follows:

$$\begin{array}{}\text{(4)}& z_u=z_{\mathrm{m}}\left[\ln \left({\displaystyle \frac{z_{\mathrm{m}}}{z_{\mathrm{0}}}}\right)-\mathrm{1}+{\displaystyle \frac{z_{\mathrm{0}}}{z_{\mathrm{m}}}}\right],\end{array}$$

where z₀ represents the roughness height set to be 0.1 H_c (canopy height).

Note that the eddy covariance system failed from 23 October 2010 to 1 April 2011 during the wheat dormant season. To evaluate the seasonal CO₂ budget of this rotation cycle, the flux gap of this period was filled by using the machine learning Support Vector Regression (SVR) algorithm (Cristianini and Shave-Taylor, 2000), which has been proved to be an appropriate tool for flux gap filling (e.g., Kang et al., 2019; Kim et al., 2019) (see Appendix A).

2.3 Meteorological and environmental condition measurements

The meteorological variables were measured at 30 min intervals by a standard meteorological station on the tower. Among these variables were the air temperature (T_a) and relative humidity (RH) (HMP45C, Vaisala Inc, Helsinki, Finland) at a height of 1.6 m and precipitation (P) (TE525MM, Campbell Scientific Inc), incoming shortwave radiation (R_si) (CRN1, Kipp & Zonen, Delft, Netherlands), and photosynthetic photon flux density (PPFD) (LI-190SA, LI-COR Inc) at a height of 3.7 m. The 30 min interval edaphic measurements included soil temperature (T_S) (109-L, Campbell Scientific Inc.) and volumetric soil moisture (θ) (CS616-L, Campbell Scientific Inc.) for the top 5 cm of the soil; soil matric potential (ψ) (257-L, Campbell Scientific Inc.) has been measured since 2010 at the same depth. The groundwater depth (WD) (CS420-L, Campbell Scientific Inc.) was measured at a location close to flux tower in 30 min intervals.

2.4 Biometric measurements and crop samples

To trace crop development and carbon storage, we measured canopy height (H_c), leaf area index (LAI), crop dry matter (DM) and carbon content of crop organs at an interval of 7–10 d in the footprint of eddy covariance. Due to inclement weather, measurement intervals were occasionally extended to 2 weeks or longer. The H_c was measured with a ruler, and LAI was measured with LAI-2000 (LI-COR Inc.) at 10 locations randomly distributed in the field. For crop samples, four locations were randomly selected at the start of the growing season, and crop samples were then collected close to these four locations throughout the experimental period. At each location, 10 crop samples were collected for wheat and 3 crop samples were collected from maize. To reduce the sample uncertainty at harvest, 200 crops and 5 crops were collected in each location for wheat and maize, respectively. The crop organs were separated and oven-dried at 105 ^∘C for kill-enzyme torrefaction for 30 min and then oven-dried at 75 ^∘C until a constant weight. The crop samples were used to estimate the average field biomass (Dry Matter). The carbon content was analyzed using the combustion–oxidation–titration method (National Standards of Environmental Protection of the People's Republic of China, 2013) to estimate carbon storage. The crop samples provided a direct estimate of the NPP.

2.5 Soil respiration measurements

Soil respiration was measured every day in the footprint of the eddy covariance between 13:00 and 15:00 UTC+8 from March to September 2011 using a portable soil respiration system LI-8100 (LI-COR Inc.). Below-ground autotrophic respiration and heterotrophic respiration were differentiated using the root exclusion method (Zhang et al., 2013). The total soil respiration (R_S) and R_H were measured at treatments with and without roots, respectively, and the corresponding difference is R_AB. To reduce the uncertainty associated with spatial variability, we set three replicate pairs of comparative treatments (i.e., with root and without root) randomly in the field. The uniform field condition contributes to reducing the measurement uncertainty associated with the spatial variability (see Zhang et al., 2013). To assess the seasonal variations and total amount of soil respirations, the seasonal continuous R_H was constructed using the Q₁₀ model by incorporating soil moisture as follows (Zhang et al., 2013):

$$\begin{array}{}\text{(5)}& {\displaystyle }{R}_{\mathrm{H}}=A\exp \left(B{T}_{\mathrm{S}}\right)\cdot f\left(\mathit{\theta }\right),\text{(6)}& {\displaystyle }f\left(\mathit{\theta }\right)=\left\{\begin{array}{l}\mathrm{1},\mathit{\theta }\le {\mathit{\theta }}_{\mathrm{f}}\\ a(\mathit{\theta }-{\mathit{\theta }}_{\mathrm{f}}{)}^{\mathrm{2}}+\mathrm{1},\mathit{\theta }>{\mathit{\theta }}_{\mathrm{f}}\end{array}\right.,\end{array}$$

where θ_f is the field capacity. The parameters were inferred by fitting the R_H and T_S measurements by using the least-squares method (see Zhang et al., 2013), where A=1.16 µmol m⁻² s⁻¹, B=0.0503 ^∘C⁻¹ and $a=-\mathrm{44.9}$ (unitless) (see Zhang et al., 2013). Note that the plant biomass was negligible before 14 March, during which R_H was set to equal to the ecosystem respiration and the R_AB was assumed to be 0. R_AB of other periods was estimated based on the R_H measurement and the ratio of R_AB to R_S estimated previously (Zhang et al., 2013), and the continuous R_AB∕R_S ratio was interpolated from the daily records (Fig. 2). This estimation method is robust because the R_AB∕R_S ratio is nearly constant around its diurnal average (Zhang et al., 2015b).

Figure 2 Seasonal variations in the ratio of below-ground autotrophic respiration (R_AB) to total soil respiration (R_S). Two vertical dashed lines represent the date of harvesting wheat and sowing maize, respectively; this is also used in Figs. 5, 6, 9 and 10.

2.6 Synthesis of the CO₂ budget components

The CO₂ budget components were derived by combining the eddy covariance measurements, soil respiration experiments and crop samples. Eddy-covariance-measured NEE is the difference between carbon assimilation (i.e., GPP) and carbon release (i.e., ER). The ER consists of R_H, R_AB (i.e., root respiration) and above-ground autotrophic respiration (R_AA). The total soil respiration is the sum of R_H and R_AB,

$$\begin{array}{}\text{(7)}& R_{\mathrm{S}}=R_{\mathrm{H}}+R_{\mathrm{AB}}.\end{array}$$

The total autotrophic respiration (R_A) is the difference between the eddy-covariance-derived ER and R_H,

$$\begin{array}{}\text{(8)}& R_{\mathrm{A}}=\mathrm{ER}-R_{\mathrm{H}}.\end{array}$$

The above-ground autotrophic respiration (R_AA) is the difference between the eddy-covariance-derived ER and R_S in Eq. (6),

$$\begin{array}{}\text{(9)}& R_{\mathrm{AA}}=\mathrm{ER}-R_{\mathrm{S}}.\end{array}$$

NPP is plant biomass carbon storage and can be quantified as the difference between GPP and R_A,

$$\begin{array}{}\text{(10)}& {\mathrm{NPP}}_{\mathrm{EC}}=\mathrm{GPP}-R_{\mathrm{A}},\end{array}$$

where the subscript “EC” represents that the NPP is estimated from the eddy-covariance-derived GPP. In parallel, NPP can also be directly inferred from biomass samples as follows:

$$\begin{array}{}\text{(11)}& {\mathrm{NPP}}_{\mathrm{CS}}=C_{\mathrm{cro}},\end{array}$$

where the subscript “CS” indicates that NPP is based on crop samples and C_cro is the plant biomass carbon storage at harvest. We used the average of the two independent NPPs as the measurement for this site.

NEP is commonly estimated by the NEE measurement (NEP_EC = −NEE). In this study, the crop samples and soil respiration measurements also provided an independent estimate as follows:

$$\begin{array}{}\text{(12)}& {\mathrm{NEP}}_{\mathrm{CS}}={\mathrm{NPP}}_{\mathrm{CS}}-R_{\mathrm{H}}.\end{array}$$

We used the average of the two NEPs as the measurement for this site.

At this site, there were no fire and insect disturbances and no manure fertilizer application. The carbon input from seeds was negligible, and all crop residues were returned to the field. Thus, NBP can be quantified as the difference between NEP and grain export carbon loss (C_gra),

$$\begin{array}{}\text{(13)}& \mathrm{NBP}=\mathrm{NEP}-C_{\mathrm{gra}}.\end{array}$$

3 Results

3.1 Meteorological conditions and crop development

The interannual variations in major meteorological variables are shown in Fig. 3, and they showed no clear trend for both wheat and maize seasons. For the full 2010–2011 cycle with comprehensive experiments, the average R_si and T_a were very close to other years; however, the P during maize season was a little higher than other years (Fig. 3c), leading to a shallow WD in the maize season (Fig. 3d). The intra-annual variations in field microclimates for the full 2010–2011 cycle are shown in Fig. 4. The seasonal maximum and minimum T_a occurred in July and January, respectively, and the variations in vapor pressure deficit (VPD) followed the T_a well. The WD mainly followed the irrigation events in winter and spring but followed P in summer and autumn. In particular, the WD varied from 0 to 3 m throughout the year. The wet soil conditions prohibited the field from experiencing water stress (Fig. 4d) because even the lowest soil matric potential (−187.6 kPa) remained a lot higher than the permanent wilting point of crops (around −1500.0 kPa).

Figure 3 The seasonal (a) total incoming shortwave radiation (R_si), (b) average air temperature (T_a), (c) total precipitation (P) and (d) average groundwater depth (D) for both wheat and maize evaluated for the period from 2005 to 2016. Note that incoming shortwave radiation in the 2013 season is missing due to equipment malfunction.

Figure 4 Seasonal variations in the environmental variables of (a) air temperature (T_a) and vapor pressure deficit (VPD); (b) photosynthetic photon flux density (PPFD); (c) precipitation (P), irrigation (I), and groundwater depth (WD); and (d) volumetric soil moisture (θ) and soil matric potential (ψ_m).

Figure 5 shows the seasonal variations in H_c and LAI, reflecting the crop development for the full 2010–2011 cycle. The maximum LAI was 4.2 m² m⁻² for wheat and 3.6 m² m⁻² for maize. The variations in H_c and LAI distinguished the different stages of crop development. During the wheat season, the stages of regreening, jointing, booting, heading and maturity started approximately on 1 March, 20 April, 1 May, 7 May and 5 June, respectively. The seasonal variations in DM agreed well with the crop stages (Fig. 6), and the wheat biomass mainly accumulated in April and May, while maize biomass mainly accumulated in July and August. The total DM was 1718 g m⁻² for wheat and 1262 g m⁻² for maize at harvest. Upon harvest, the wheat DM was distributed as 3 % root, 43 % stem, 9 % leaf and 45 % grain, while the maize DM was distributed as 2 % root, 29 % stem, 7 % green leaf, 5 % dead leaf, 4 % bracket, 7 % cob and 46 % grain. The seasonal average carbon contents of the root, stem, green leaf, dead leaf, and grain were 410, 439, 486, 452, and 457 gC kg⁻¹ DM for wheat and 408, 438, 477, 457, and 456 gC kg⁻¹ DM for maize (see Table 1 for the seasonal variation).

Figure 5 Seasonal variations in canopy height (H_c) and leaf area index (LAI). The error bars denote 1 standard deviation of the 10 points.

Figure 6 Seasonal variations in the total dry biomass (DM) and its major components of root, stem, green leaf and grain. The error bars of total biomass denote 1 standard deviation of the four sample points.

Table 1 Carbon content of different parts of each crop (gC kg⁻¹ DM).

3.2 The interannual variations in the NEE, GPP and ER

For the period from 2005 to 2016, if grain export was not considered, wheat was a consistent CO₂ sink, as the seasonal total NEEs were consistently negative, and maize was a CO₂ sink in most years, except for 2012 and 2013 when NEE was positive (Fig. 7a). NEEs of both wheat and maize fields became less negative during the past decade (though not in a statistically significant way), implying a progressive decline of the carbon sequestration potential of this cropland. The GPPs of both wheat and maize showed an increasing trend, though they were not statistically significant (Fig. 7b). The ERs of both wheat and maize also showed an increasing trend in these years, but only the trend of maize was significant (Fig. 7c). The decadal average of NEE, GPP, and ER were −364 (SD ± 98), 1174 (SD ± 189), and 810 (SD ± 161) gC m⁻² for wheat and −136 (SD ± 168), 1008 (SD ± 297), and 872 (SD ± 284) gC m⁻² for maize.

Figure 7 The temporal trend of annual (a) net ecosystem exchange (NEE), (b) gross primary productivity (GPP) and (c) ecosystem respiration (ER) for both wheat and maize from 2005 to 2016. Note that though most gaps of carbon fluxes were filled, the wheat of 2007 was excluded as it had a large gap accounting for 26 % of annual records that we were unable to fill. Maize was not planted in the growing season of 2010. Note that the solid line shows where the temporal trend passes F test at p<0.05 significance level, while the dashed line shows where the temporal trend does not pass the F test at p<0.05 level.

The NEE, GPP and ER for both wheat and maize were correlated with the three main environmental variables of R_si, T_a and WD using the multiple regression (see Appendix B for details). In the wheat season, T_a showed its relatively great importance (compared to R_si and WD) to all three of the CO₂ fluxes with a higher T_a increasing both GPP and ER and also enhancing NEE (more negative) (Fig. 8a). WD correlated negatively with GPP, thereby reducing net carbon uptake (less negative NEE). WD exhibited almost no effect on ER. R_si exhibited almost no effect on all three CO₂ fluxes. Therefore, T_a explained most of the interannual variations in NEE, GPP and ER, followed by WD. In the maize season, WD had good correlations with all three fluxes of GPP, ER and NEE, where a deeper WD contributed to lower both GPP and ER and also drove higher net carbon uptake (more negative NEE). T_a showed almost no effect on all three CO₂ fluxes. R_si had a positive correlation with ER but almost no correlation with GPP (Fig. 8b). Ultimately, higher R_si in the maize season lowered the net carbon uptake (more positive NEE). Overall, R_si and WD showed their great importance in influencing the interannual variation in maize NEE, with R_si having a positive correlation and WD having a comparable negative correlation (Fig. 8b).

Figure 8 The result of multiple regression for NEE, GPP and ER with incoming shortwave radiation (R_si), air temperature (T_a) and groundwater depth (WD) for both (a) wheat and (b) maize. Note that ^* denotes that the regression passes p<0.05 significance level and that NS indicates nonsignificant.

3.3 Intra-annual variations in the NEE, GPP and ER

The intra-annual variations in NEE, GPP and ER exhibited a bimodal curve corresponding with the two crop seasons (Fig. 9). All three CO₂ fluxes were almost in phase, with peaks appearing at the start of May during the wheat season and in the middle of August during the maize season. During some of the winter season, the field still sequestered a small amount of CO₂ because of the weak photosynthesis, which was confirmed by leaf level gas exchange measurement (data not shown). Net carbon emission happened during the fallow periods, in addition to the start of the maize season when the plant was small and high temperatures enhanced heterotrophic respiration. During the wheat season, two evident spikes appeared on 21 April and 8 May with positive NEE values (i.e., net carbon release). These spikes resulted from the radiation decline during the inclement weather (Fig. 4b), which suppressed the photosynthesis rate; similar phenomena also appeared during the maize season.

Figure 9 Seasonal variations in gross primary productivity (GPP), net ecosystem exchange (NEE) and ecosystem respiration (ER) (those before 2 April were calculated with the SVR method).

Figure 10 shows the variations in ER and its components. During the wheat season, the variation in ER closely followed crop development and temperature, but there were two evident declines at the end of April and the start of May due to low temperatures associated with the inclement weather. During the early growing stage of maize, R_H was the main component of ER. When waterlogging conditions occurred in late August and early September, both R_H and R_AB were suppressed to zero.

Figure 10 Seasonal variations in the components of ecosystem respiration (ER), total soil respiration (R_S) and soil heterotrophic respiration (R_H). The difference between ER and R_S denotes above-ground autotrophic respiration (R_AA), and the difference between R_S and R_H denotes below-ground autotrophic respiration (R_AB).

3.4 CO₂ budget synthesis in the 2010–2011 agricultural cycle

CO₂ budget analysis showed that this wheat–maize rotation cropland has the potential to uptake carbon from the atmosphere (Fig. 11). In the full 2010–2011 cycle, the total NEE, GPP, and ER values were −438, 1078, and 640 gC m⁻² for wheat and −239, 780, and 541 gC m⁻² for maize. The NPP values were 750 and 815 gC m⁻² for wheat based on crop samples and the eddy covariance and that complemented with soil respiration measurements, respectively, and were 592 and 532 gC m⁻² for maize based on the two methods. We used the average of these two methods for NPP measurements, which were 783 (SD ± 46) gC m⁻² for wheat and 562 (SD ± 43) gC m⁻² for maize. We also used the average of NEP from the two independent methods for the measurement, and the NEP was 406 gC m⁻² for wheat and 269 gC m⁻² for maize. Furthermore, when considering the carbon loss associated with the grain export, the NBP values were 59 gC m⁻² for wheat and 5 gC m⁻² for maize, respectively. Considering the net CO₂ loss of −104 gC m⁻² during the two fallow periods, NBP of the whole wheat–maize crop cycle was −40 gC m⁻² yr⁻¹, suggesting that the cropland was a weak carbon source to the atmosphere under these specific climatic conditions and field management practices.

Figure 11 Carbon budget of wheat (a), maize (b) and the full wheat–maize rotation cycle, with fallow periods included (c). Note that the absolute value of NEE is shown here; NBPs of wheat and maize are the average of two independent methods (i.e, an eddy covariance-based method and a crop-sample-based method).

4 Discussion

This study investigated the decadal variations in the NEE, GPP and ER for an irrigated wheat–maize rotation cropland over the North China Plain, and the results exhibited a decreasing trend of the CO₂ sink capacity during the past decade. The interannual variations in the carbon fluxes of wheat showed close dependence on temperature and groundwater depth, while those of maize were mostly regulated by solar radiation and groundwater depth. Furthermore, the detailed CO₂ budget components were quantified for a full wheat–maize agricultural cycle. Investigating the decadal trend of the CO₂ fluxes and quantifying the detailed CO₂ budget components for this representative cropland will provide useful knowledge for regional greenhouse gas emission evaluation over the North China Plain.

4.1 Comparison with other croplands

The cropland has been reported as carbon neutral to the atmosphere (e.g., Ciais et al., 2010), as a carbon source (e.g., Anthoni et al., 2004a; Verma et al., 2005; Kutsch et al., 2010; Wang et al., 2015; Eichelmann et al., 2016) and also as a carbon sink (e.g., Kutsch et al., 2010). Such inconsistency probably results from the different crop types and management practices (residue removal, the use of organic manure, etc.), in addition to variations in the climatic conditions (Béziat et al., 2009; Smith et al., 2014) and fallow period length (Dold et al., 2017). Our results show that the fully irrigated wheat–maize rotation cropland with a shallow groundwater depth was a weak CO₂ sink during both the wheat and maize seasons in the full 2010–2011 cycle, but the CO₂ loss during the fallow period reversed the cropland from a sink into a weak source with an NBP of −40 gC m⁻² yr⁻¹. These results are consistent with previous studies that reported the wheat–maize rotation cropland as a carbon source (Li et al., 2006; Wang et al., 2015). However, the net CO₂ loss was much lower at our site, most likely due to the shallow groundwater depth.

Field measurements of the long-term CO₂ fluxes over croplands remain lacking, and we found the carbon sequestration capacity of this cropland has been progressively decreasing, though it was not statistically significant. The cropland has been widely suggested as a climate change mitigation tool (e.g., Lal, 2001), but the potential in the future is challenging. However, since cropland management greatly impacts the carbon balance of cropland (Béziat et al., 2009; Ceschia et al., 2010), it remains required investigating if the management adjustment can foster the cropland carbon sink capacity over the long term.

The annual total NPP of 1345 gC m⁻² yr⁻¹ at our site is approximately twice the average of the model-estimated NPP for Chinese croplands (714 gC m⁻² yr⁻¹) with a rotation index of 2 (i.e., two crop cycles within 1 year) (Huang et al., 2007), more than 3 times the value estimated by MODIS (400 gC m⁻² yr⁻¹) (Zhao et al., 2005) and slightly higher than the value of the same crop rotation at the Luancheng site (1144 gC m⁻² yr⁻¹) (Wang et al., 2015). The higher NPP at our site may partially result from the sufficient irrigation and fertilization (Huang et al., 2007; Smith et al., 2014).

The contrasting respiration partitioning of the same crop in different regions (Table 2) indicate that the respiration processes may also be subject to climatic conditions and management practices. Though the ratio of ecosystem respiration to GPP at our site is comparable to other studies, the ratio of autotrophic respiration to GPP is much lower at our site, and while the ratio of heterotrophic respiration to ecosystem respiration is greater at our site, these findings are different from those at the other sites with similar crop variety (Moureaux et al., 2008; Aubinet et al., 2009; Suleau et al., 2011; Wang et al., 2015; Demyan et al., 2016), as they showed that ecosystem respiration is usually dominated by below-ground and above-ground autotrophic respirations. The higher soil heterotrophic respiration at our site probably results from the full irrigation and shallow groundwater, which both alleviate soil water stress.

Table 2 Various ratios associated with carbon fluxes in croplands.

^a The values in parentheses indicate that the value is calculated by the equation R_A ∕ GPP = 1 − NPP∕GPP. ^b The data were from 2012, and the estimation is based on the average of the static and dynamic methods. ^c R_A and R_H are the averaged values of the two corresponding methods.

4.2 The effects of groundwater on carbon fluxes

The groundwater table at our site is much closer to the surface because of the irrigation by water diverted from the Yellow River. In contrast, the nearby Luancheng site (Wang et al., 2015) is groundwater-fed with a very deep groundwater depth (approximately 42 m) (Shen et al., 2013), and their CO₂ budget components had some differences with our study. Comparing the net CO₂ exchange of wheat, the GPP at our site is a little higher than the Luancheng site, implying the irrigation at our site may better sustain the photosynthesis rate for wheat; ER at our site is also a little higher than the Luancheng site. For maize, both sites are not irrigated due to the high summer precipitation. GPP and ER at our site were comparable to Luancheng site, implying that the irrigation method prior to the maize season had no discernible effect on the integrated CO₂ fluxes for maize. However, the three components of ER in our study showed pronounced differences from the Luancheng site, where they reported the R_AA was 411 gC m⁻² for wheat and 428 gC m⁻² for maize, 3 times the results of our study (128 gC m⁻² for wheat and 133 gC m⁻² for maize). However, their R_AB for wheat (36 gC m⁻²) and maize (16 gC m⁻²) were less than a quarter of our results (136 gC m⁻² for wheat and 115 gC m⁻² for maize). Their R_H of wheat (245 gC m⁻²) was less than our estimate (377 gC m⁻²), but R_H of maize (397 gC m⁻²) was greater than our result (292 gC m⁻²). In general, the above-ground crop parts in our site respired more carbon than the Luancheng site, possibly because the shallow groundwater depth at our site increased the above-ground biomass allocation but lowered the root biomass allocation (Poorter et al., 2012). These independent cross-site comparisons demonstrate that carbon budget components may be subject to the specific groundwater depth influenced by the irrigation type, and even the same crop under similar climatic conditions can behave differently in carbon consumption.

Our site experienced a short period of waterlogging during the 2010–2011 cycle due to the combined effects of full irrigation and the high precipitation during the summer. This distinct field condition reduced soil carbon losses in the maize season, potentially maintaining the CO₂ captured by the cropland. Waterlogging events were occasionally reported in upland croplands. For example, Terazawa et al. (1992) and Iwasaki et al. (2010) suggested that waterlogging causes damage to plants, resulting in a decline in GPP as reported by Dold et al. (2017) and our study. Our study further shows that waterlogging reduces ER to a greater degree than GPP, possibly because of the low soil oxygen conditions, and thereby reduces the overall cropland CO₂ loss. However, the CH₄ released over the short term may be pronounced in waterlogged soils. As CH₄ emission in this kind of cropping system over the North China Plain cropland remains lacking, additional field experiments are required to understand how irrigation and water saturation field condition impact the overall carbon budget.

4.3 Uncertainty in the estimation and limitation of this study

In the comprehensive experiment period for the full 2010–2011 agricultural cycle, the NEE of the wheat season from 23 October 2010 to 1 April 2011 was calculated using a calibrated SVR model. The SVR model performs well for predicting GPP and ER with very high R² of 0.95 and 0.97 and an acceptable uncertainty level of 22.9 % and 15.2 % for GPP and ER, respectively. Hence, these estimates should have a negligible effect on the seasonal total carbon evaluation. The footprint analysis showed that 90 % of the measured eddy flux comes from the nearest 420 and 166 m in wheat and maize crops under unstable conditions, respectively, confirming that both soil respiration experiments and crop samples paired well with the EC measurements .

Root biomass was difficult to measure, but the uncertainty should be low because the root ratio (the ratio of the root weight to the total biomass weight) accounts for 15 %–16 % of the crop for wheat and maize (Wolf et al., 2015), and our measurements are very close to these values; i.e., the averaged seasonal root ratio was 15 % for wheat and 10 % for maize at our site. However, the relatively low root ratios (3 % for wheat and 2 % for maize) at harvest probably result from the root decay associated with plant senescence. The estimates of annual soil respiration are based on the Q₁₀ model validated by the field measurements that may generate some uncertainty in the soil respiration budget due to the hysteresis response of soil respiration to temperature (Phillips et al., 2011; Zhang et al., 2015a, 2018). However, the Q₁₀ model remains robust in soil respiration estimations if it is well validated (Tian et al., 1999; Zhang et al., 2013; Latimer and Risk, 2016), allowing for confidence in the estimates.

During the wheat season, the cumulative curves of NPP_EC and NPP_CS were not perfectly consistent in the main growing season, as clear differences emerged during the dormant season of wheat from 15 December 2010 to 8 March 2011 (Fig. 12). These differences may result from the small wheat sample number. However, the sample number at harvest was sufficiently big, and no discernible difference was found between the two NPPs at harvest. These two independent estimates of NPP were similar throughout the maize season (Fig. 12).

Figure 12 Seasonal variations in the cumulative net primary productivity (NPP) with two independent methods of crop sampling (NPP_CS) and eddy covariance (NPP_EC) complemented with soil respiration measurements.

This study provides a comprehensive quantification of the CO₂ budget components of the cropland, but it remains limited to a relatively wet year (see Fig. 3c and d). The integrated carbon fluxes (NEE, GPP and ER) have pronounced interannual variations, also suggesting further investigations are required on the interannual variations in the carbon budget components.

5 Conclusion

Based on the decadal measurements of CO₂ fluxes over an irrigated wheat–maize rotation cropland over the North China Plain, we found the cropland was a strong CO₂ sink if grain export was not considered. When considering the grain export, the cropland was a weak CO₂ source, with an NBP of −40 gC m⁻² yr⁻¹ in the full 2010–2011 agricultural cycle. The net CO₂ exchange during the past decade from 2005 to 2016 showed a statistically nonsignificant decreasing trend, implying a decreasing carbon sequestration capacity of this cropland, discouraging the potential of taking agroecosystems as the mitigation tool of climate change. In the wheat season, air temperature showed the best correlation with the CO₂ fluxes followed by the groundwater depth, whereas in the maize season both shortwave radiation and groundwater depth showed good correlation with the CO₂ fluxes. The comprehensive investigation showed most of the carbon sequestration occurred during the wheat season, while maize was close to being CO₂ neutral. Soil heterotrophic respiration in this cropland contributes substantially to CO₂ loss in both the wheat and maize seasons. This study provides detailed knowledge for estimating regional carbon emissions over the North China Plain.

Appendix A: Flux calculation of the period with equipment failure

A1 Support Vector Regression method

The Support Vector Regression (SVR) method is a machine learning technique-based regression, which transforms regression from nonlinear into linear by mapping the original low-dimensional input space to higher-dimensional space (Cristianini and Shave-Taylor, 2000). The SVR method has two advantages: (1) the model training always converges to global optimal solution, with only a few free parameters to adjust, and no experimentation is needed to determine the architecture of SVR. (2) The SVR method is robust to small errors in the training data (Ueyama et al., 2013). The Support Vector Machine (SVM) software package obtained from LIBSVM (Chang and Lin, 2005) is used in this study.

A2 Data processing and selection of explanatory variables

Gross primary productivity (GPP) is influenced by several edaphic, atmospheric and physiological variables, among which air temperature (T_a), relative humidity (RH), leaf area index (LAI), net photosynthetically active radiation (PAR) and soil moisture (θ) are the dominant factors. Hence, we select T_a, RH, LAI, PAR and θ as explanatory variables of GPP. Ecosystem respiration (ER) consists of total soil respiration and above-ground autotrophic respiration. The total soil respiration is largely influenced by soil temperature and soil moisture, while above-ground autotrophic respiration is largely influenced by air temperature and above-ground biomass. Therefore, we select T_a, soil temperature at 5 cm (T_S), θ and LAI as explanatory variables of ER. LAI is estimated from the Wide Dynamic Range Vegetation Index derived from the MOD09Q1 reflectance data (250 m, 8 d average, https://ladsweb.modaps.eosdis.nasa.gov/missions-and-measurements/products/MOD09Q1/, last access: 15 May 2016; see Lei et al., 2013).

The three wheat seasons of 2005–2006, 2009–2010 and 2010–2011 are selected for model training, and the original half-hourly measurements of GPP and ER, together with the explanatory variables, are averaged to the daily scale, but we remove days missing more than 25 % of half-hourly data. We have GPP available from 466 d and ER from 483 d for model training. The explanatory variables for the equipment failure are also averaged into daily scale, which will be used to calculate GPP and ER with the trained model described in the following section.

A3 SVR model training and flux calculation

In order to eliminate the impact of variables with different absolute magnitudes, we rescale all the variables in the training data set to the [0, 1] range prior to SVR model training. In the training process, the radial basis function (RBF, a kernel function of SVR) is used and the width of insensitive error band is set as 0.01. The SVR model training follows these steps:

1. All training data samples are randomly divided into five nonoverlapping subsets, and four of them are selected as the training sets (also calibration set); the remaining subset is treated as the test set (or validation set). This process is repeated five times to ensure that every subset has a chance to be the test set.

2. For the selected training set, the SVR parameters (cost of errors c and kernel parameter σ) are determined using a grid search with a five-fold cross-validation training process. In this approach, the training set is further randomly divided into five nonoverlapping subsets. Training is performed on each of the four subsets within this training set, with the remaining subset reserved for calculating the root-mean-square error (RMSE), and model parameters (c and σ) yielding the minimum RMSE value are selected.

3. The SVR model is trained based on the training set from step (1) and initialized by the parameters (c and σ) derived from step (2).

4. The test set from the step (1) is used to evaluate the model obtained from the step (3) by using the coefficient of determination (R²) and RMSE.

5. The model is trained with all of the available samples that achieved good performance, as R² are 0.95 and 0.97 for GPP and ER, respectively, and the mean RMSE is 1.28 and 0.44 gC m⁻² d⁻¹. The RMSE can be further used as a metric quantifying uncertainty, which accounts for 22.9 % and 15.2 % for the averaged GPP and ER, respectively. GPP and ER during the equipment failure period are then calculated with the trained model complemented with the observed explanatory variables, and NEE is derived as the difference of GPP and ER.

Appendix B: Multiple regression for NEE, GPP and ER with microclimate variables

The flux of NEE, GPP or ER is correlated with incoming shortwave radiation (R_si), air temperature (T_a) and groundwater depth (WD) as flux = $a{R}_{\mathrm{si}}+b{T}_{\mathrm{a}}+c\mathrm{WD}+d$, where flux is NEE, GPP or ER; a, b, c and d are regression parameters. All the variables are normalized to derive their z score before the regression, where the z score is calculated by subtracting the mean from the data and dividing the result by the standard deviation. The coefficient of each variable represents the relative importance of the corresponding variable in contributing to the dependent variable.