Dynamics of Carbon and Water Fluxes over Cropland and Agroforest Ecosystems on the Southern Chinese Loess Plateau

Abstract

Studies on the spatiotemporal dynamics in ecosystem carbon and water exchanges are essential in predicting the effects of climate change on regional carbon and energy budgets. Using the eddy covariance technique, carbon and water fluxes were observed in a typical winter wheat ecosystem (WWE) and an agroforest ecosystem (AFE) in the southern Loess Plateau from 2004 to 2010. The seasonal and inter-annual variability in gross primary productivity (GPP), net ecosystem exchange (NEE), evapotranspiration (ET), and water use efficiency (WUE) were examined and the main influencing factors were identified using the Pearson correlation. The results indicate that the seasonal GPP and NEE showed a bimodal distribution in WWE, while this was unimodal in AFE. The sinusoidal function did well in the characterization of seasonal ET dynamics for both ecosystems, with the determination coefficients being 0.85 and 0.94, respectively. In WWE and AFE, the annual mean GPP were 724.33 and 723.08 g C m⁻² a⁻¹, respectively, and the corresponding ET were 392.22 and 410.02 mm a⁻¹. However, the difference in NEE between the two ecosystems was obvious, NEE were −446.28 and −549.08 g C m⁻² a⁻¹, respectively, showing a stronger carbon sink in AFE. There were strong coupling relationships between the GPP and ET of both ecosystems; the overall slopes were 1.71 and 1.69, respectively. The seasonal trend of WUE was bimodal in WWE, with peak values of 3.94 and 3.65 g C kg⁻¹ H₂O, occurring in November and April, respectively. However, the monthly WUE in AFE had one single peak of 4.07 g C kg⁻¹ H₂O in January. Photosynthetically active radiation (PAR) and soil temperature (T_s) were most positively correlated with GPP, net radiation (R_n) and T_s were the major factors influencing ET, while vapor pressure deficit (VPD) and soil water content (SWC) were the major influencing factors for WUE. These results provide observational support for regional carbon neutrality simulations.

1. Introduction

Evapotranspiration (ET) and photosynthesis are coupled processes that mediate water, carbon, and energy fluxes between terrestrial ecosystems and atmosphere, which are important components of energy–water balances and global biogeochemical cycling [1,2,3]. In the context of climate change, observations of ecosystem water and carbon exchange, in order to deeply understand the trade-off between ET and the amount of terrestrial carbon assimilated, are currently one of the hottest issues in ecological hydrology research [4,5]. Studies on the characteristics of water and carbon fluxes and their coupling relationships in agroecosystems are essential to improving carbon and ET models [6,7], as well as predicting the effects of climate change on regional carbon and energy budgets [8]. Meanwhile, such study is also of great significance for maintaining and improving agricultural productivity and further understanding the material cycle of terrestrial ecosystems.

The eddy covariance technique, which provides direct and dependable flux measurements, has been used to characterize the temporal and spatial variability of ecosystem carbon dioxide (CO₂) fluxes (net ecosystem carbon exchange, NEE), ET, and climatic variables [9]. Based on this technique, measurements of high-frequency NEE allow us to investigate the dynamics of gross primary productivity (GPP), as well as estimations of ecosystem water use efficiency (WUE, the ratio of GPP to ET) [10,11]. In recent decades, the Global Flux Observation Network (FLUXNET) has provided a large amount of data for the study of ecosystem water–carbon flux and its response to climate change [12], covering almost all major terrestrial ecosystems, including forests [13,14,15,16], grasslands [17,18,19], and farmlands [20,21,22,23]. Most previous studies have focused only on the water flux [20] or carbon budget [21,22,24] of single types of ecosystem separately, while ignoring the possible effects of parts outside the observational range on the ecosystem.

The responses of GPP, ET, and their parameters related to environmental factors differed among ecosystems, resulting in uncertainties in WUE estimations [10,25], which makes it urgent to find biophysical controls for inter-annual variations in GPP, NEE, ET, and WUE in the agroecosystems. Jia et al. (2016) found that low soil moisture conditions in spring coincided with lower seasonal and annual carbon and ET, as well as lower WUE [26]. Wagle et al. (2018) studied the relationship between maize NEE and photosynthetically active radiation (PAR) by using the eddy covariance technique, and the results showed that NEE was almost linear with PAR under favorable climatic conditions [22]. Studies have shown that vapor pressure deficit (VPD) is an important environmental factor for WUE [27,28]. The results of Tong et al. (2009) indicated that, when LAI and solar radiation were large enough, WUE showed negative correlation with VPD in irrigated wheat and maize fields on the North China Plain [4]. Other factors also contributed to variation in WUE, and Tang et al. estimated the WUE of different ecosystem types at 32 FLUXNET stations in the northern hemisphere, showing that the WUE had obvious latitude zonality [29]. However, water–carbon couplings in cropland and other agroecosystems had their own particularities and the influencing factors were more complex. Few studies carry out comprehensive analysis and evaluation of water–carbon-energy processes based on years of continuous observation, making it difficult to fully reveal regional ecohydrological processes and their characteristics.

The Changwu tableland is an important grain-producing area located in the warm temperate subhumid continental monsoon climate zone, south of the Loess Plateau, China. In recent decades, climate change and continuous improvement in agricultural productivity have significantly changed regional water–carbon processes. In the present study, we summarized nearly six years of consecutive measurement data of NEE, ET, and environmental factors such as air temperature (T_a), soil temperature (T_s), relative humidity (RH), wind speed (W_s), precipitation (P), soil water content (SWC), VPD, photosynthetically active radiation (PAR), and net radiation (R_n), covering the period from September 2004 to June 2010, which included six winter wheat seasons. The main objectives were to: (1) analyze seasonal and interannual variations in GPP, NEE, ET, and WUE in a winter wheat ecosystem (with a single crop from September to the next June, observed at a 2 m height, denoted as WWE), and an agroforest ecosystem (with both cropland and apple orchard, observed at a 32 m height, denoted as AFE); (2) to explore the applicability of sinusoidal function in the simulation of interannual ET; and (3) to quantify the correlation between environmental factors and the GPP, NEE, and WUE of the two different agroecosystems using the Pearson correlation.

2. Materials and Methods

2.1. Site Description

The study site is located at Changwu Agro-ecological Experimental Station, Chinese Academy of Sciences (CAS), located in Changwu County, Shaanxi Province, China (35°14′ N, 107°41′ E, 1220 m above sea level, Figure 1a). The climate is controlled by the sub-humid warm continental monsoon, with a long-term mean annual temperature of 9.1 °C, precipitation of 584 mm, and potential evapotranspiration of 949.3 mm (1957–2013) [30]. The major soil type is dark loessial soil, with a soil organic carbon content of 6.74 g kg⁻¹ (0–20 cm profile), and a saturated soil water content of 0.29 cm³ cm⁻³. The main land use types include wheat and maize cultivation and apple orchards. The growing season for winter wheat is from September to the next June, with a summer fallow period in July to early September. A rotary cultivator was used to a depth of nearly 30 cm during sowing, and the nitrogen, phosphorus, potassium applied as base fertilizer were 120, 90, and 60 kg ha⁻¹, respectively. The apple orchards covered 160.6 km² in 2018, accounting for 70% of the cultivated land. The tableland of the study site has a flat surface with little runoff, the groundwater depth is 50–80 m below the surface [31], and agricultural production mainly relies on precipitation with no irrigation.

2.2. Field Measurements

Flux and meteorological instruments were mounted on a 32 m high tower erected at the center of the observation field in Changwu station, the tower is about 500, 1400, 1000, and 1800 m away from the surrounding steep ravines in the directions of east, south, west, and northwest, respectively [32,33]. Two sets of eddy covariance (EC) systems were used for measuring latent heat flux (LE), sensible heat flux (H_s), CO₂ flux (F_c) at 2 m and 32 m above ground. The area around the experimental site is quite variable, but presumably statistically homogeneous at the meso-gamma scale which ranges roughly from 2 to 20 km, as the variabilities are constant in space, as described by Brutsaert (1998) [34] and Brutsaert et al. (2017) [33]. The details of the topography around the experimental site are given in Li et al. (2007) [32] and Brutsaert et al. (2017) [33]. According to the footprint study of Chu et al. (2009) [35] at this site, when the measurement height changed from 2 m to 12 m above ground, the largest observation distance on prevailing wind direction (southeast) changed from 172.8 m to 1555.2 m, and it could be larger when the measurement height is increased [36]. Thus, the footprint of the EC system placed at 2 m covered the winter wheat continuous cropping field around the tower, while the footprint of the EC system at 32 m could represent the complex of wheat field, apple orchard, and other ecosystems (the AFE) within a distance of more than 2000 m in the direction of the prevailing wind (Figure 1b).

An open-path infrared gas analyzer (LI7500, LI-COR, Inc., Lincoln, NE, USA) was used to measure fluctuations of CO₂ and H₂O density. Wind velocity and virtual temperature fluctuations were measured using a three-dimensional sonic anemometer (CSAT3, Campbell Scientific, Inc., Logan, UT, USA). The sampling frequency of the eddy covariance system was 10 Hz. T_a and RH (HMP45C, Vaisala Inc., Vantaa, Finland), W_s (3-cup anemometer, Vector Instruments, Hudson, NH, USA) were all measured at 2 m and 32 m above the ground. P was measured using a tipping bucket rain gauge (RA-1, Ota Keiki Seisakusho Co., Ltd., Tokyo, Japan) at 0.4 m above the ground. PAR (LI190SB, LI-COR, Inc., Lincoln, NE, USA) was measured at 2 m above the ground. T_s and SWC were measured at 2, 10, 20, 40, and 80 cm below the ground using a thermometer (Campbell Scientific, Inc., Logan, UT, USA) and TDR sensor (IMKO, Ettlingen, Germany). Climatological and soil parameters were recorded at 30 min intervals with a CR10X data logger (Campbell Scientific, Inc., Logan, UT, USA). VPD and R_n were calculated based on the above parameters.

2.3. EC Data Processing and Gap Filling

Raw data from the EC systems was processed offline using the improved EdiRe software (V1.5.0.32) [37]. The output of half-hourly LE and H_s was calculated from the mean covariance of water and temperature fluctuations with fluctuating vertical velocity; half-hourly F_c was computed by the eddy covariance method as the mean covariance between fluctuations in vertical wind velocity (w′, m s⁻¹) and the CO₂ density (ρ_c′, mg m⁻³). For data quality control, spectral loss correction [38], 3D coordinate rotation [39], sonic virtual temperature conversion [40], Webb–Pearman–Leuning (WPL) correction [41], and spike detection [42] were included in the program. In the spike detection, the values that exceeded the mean value by ±5.5 times the standard deviations in a window of 10 values were labelled as spike. The unfavorable data were deleted and gap-filled to calculate the daily mean values. Short gaps (<2 h) were filled with linear interpolation, and larger data gaps were filled using methods based on the empirical relationships between LE, H_s, F_c, and meteorological data such as air temperature, PAR, etc. [43]. Where meteorological data were missing, the mean diurnal variation method was applied using a 7-day data window [12].

2.4. Calculations of the Main Parameters

ET (mm) was calculated as:

$$ET=LE/\lambda$$

(1)

where LE (w m⁻²) is the output latent heat flux, λ (MJ kg⁻¹) is the latent heat of vaporization of water, which can be calculated from air temperature (T, °C):

λ = (2.501 − 0.00237 × T) × 10⁶

(2)

The NEE was linked to the eddy covariance CO₂ flux (F_c, μmol m⁻² s⁻¹) and the rate of canopy CO₂ storage change below the sensor:

$$NEE=F_{\mathrm{c}}+F_{\mathrm{st}}$$

(3)

where F_st is the storage flux, which reflects the accumulation and depletion of CO₂ in the canopy volume, and was estimated based on the assumption that half hour changes in CO₂ concentration at the 2 m and 32 m levels are representative of the whole layer below the eddy covariance systems [44], as follows:

$$F_{\mathrm{st}}={\rho }_{\mathrm{a}}h∆c/∆t$$

(4)

where ρ_a is the molar density of air (mol m⁻³); h is the height of the eddy covariance system (m), Δc is the change in CO₂ concentration (μmol mol⁻¹) during the time interval, Δt (s).

Nighttime NEE is equal to R_eco, and daytime R_eco was estimated from a regression model using the relationship between nighttime NEE and soil temperature (T_s) at a depth of 10 cm. This method is based on the Vant Hoff equation as described by Lloyd and Taylor (1994) [45]:

$$R_{\mathrm{eco}}=R_{\mathrm{ref}}exp\left(b{T}_{\mathrm{s}}\right)$$

(5)

where R_ref (the reference respiration at T_s of 0 °C) and b are the regression parameters.

GPP was calculated as the difference between R_eco and NEE, which represent the CO₂ assimilation and release of the ecosystem, respectively.

$$GPP=R_{\mathrm{eco}}-NEE$$

(6)

Based on gas exchange approaches, water use efficiency can be calculated at ecosystem, canopy, and leaf scales. In this study, the daily ecosystem WUE (g C kg⁻¹ H₂O) was calculated through dividing GPP (g C m⁻² d⁻¹) by ET (kg H₂O m⁻² d⁻¹) [8]:

$$WUE=GPP/ET$$

(7)

To distinguish these parameters between the two ecosystems, the GPP, NEE, ET, and WUE were denoted as GPP_w, NEE_w, ET_w, and WUE_w in WWE, while they were denoted as GPP_a, NEE_a, ET_a, and WUE_a in AFE.

2.5. Model Description and Application

To quantify the seasonal variations in ET and make the parametric representations, the sinusoidal function considering phase transition was used [46,47,48]:

$$ET (t)=\overline{ET} (1 + {\delta }_{ET}\sin (\frac{2\pi }{\tau } \frac{ t-S_{ET}}{12}))$$

(8)

where ET(t) is the monthly ET with the annual mean value of $\overline{ET}$; t is the time (in months); τ is the cycle of seasonality, with 6 months in the tropics and 1 year outside the tropics; δ_ET is the dimensionless seasonal amplitude of monthly ET, representing the extent of the change during the year; S_ET is the phase shift of monthly ET, in this study, S_ET takes the values of 1.0, 1.2, 1.4, …, 5.8, and 6.0, for trial calculations.

2.6. Statistical Analyses

The root mean square error (RMSE) and determination coefficient (R²) were calculated to verify the simulation performance of the sinusoidal functions, including the extent of overestimation or underestimation in the simulation, as well as the proportion of variation in dependent variables explained by the independent variable. These indexes could be calculated as follows:

$$\mathrm{RMSE}=\sqrt{\frac{{\sum }_{i=1}^{\mathrm{n}}{\left(S_i{-Q}_i\right)}^2}{n}}$$

(9)

$$R^2=1-\frac{{\sum }_{i=1}^n{(Q_i-{\widehat{Q}}_i)}^2}{{\sum }_{i=1}^n{(Q_i-\overline{Q})}^2}$$

(10)

where S_i is the simulated value, Q_i is the observed data, ${\widehat{Q}}_i$ is the predicted value according to the regression coefficients, $\overline{Q}$ is the average value of the observed data, and the variable n represents the number of pairs of simulated and observed data.

One-way analysis of variance (ANOVA) was performed to evaluate the significance level of the differences in ET, GPP, NEE, and WUE between the WWE and AFE. Multiple means comparisons were made using least significant difference (LSD) with statistical significance denoted at the 5% level. Pearson correlation analysis was conducted using Origin 2021 Pro to examine the correlation between T_a, T_s, W_s, RH, e, R_n, PAR, and SWC on GPP, ET, and WUE (OriginLab Corporation, Northampton, MA, USA). The statistical analyses were conducted using SPSS 22 for Windows (SPSS Inc., Chicago, IL, USA).

3. Results

3.1. Seasonal and Interannual Variations in Ecosystem Carbon Fluxes and Environmental Factors

The results showed strong seasonal patterns in GPP and biophysical factors (Figure 2). Daily mean T_a varied seasonally from low values of about −14.6 °C in winter to high values exceeding 27.2 °C in summer (Figure 2a). T_s showed a similar pattern, but had a higher annual mean T_s (11.2 °C) than annual mean T_a (9.4 °C) during the study period. SWC fluctuated following large rainfall events (e.g., >20 mm day⁻¹, Figure 2b).

The seasonal dynamics of GPP_w and NEE_w were closely related to crop development and phenology (Figure 2c and Figure 3). The GPP_w increased from 33.30 g C m⁻² month⁻¹ in September to an initial peak of 58.76 g C m⁻² month⁻¹ in November (Figure 3a), while NEE_w decreased from −4.53 g C m⁻² month⁻¹ to −45.04 g C m⁻² month⁻¹ (Figure 3b). During the wintering stage, GPP_w decreased to the minimum value at 13.19 g C m⁻² month⁻¹ in January, and increased to the maximum value at 158.62 g C m⁻² month⁻¹ in May; the corresponding NEE_w decreased from −5.03 g C m⁻² month⁻¹ to −125.82 g C m⁻² month⁻¹. After the ripening and harvesting stage, GPP_w decreased to 37.57 g C m⁻² month⁻¹ in July, while NEE_w increased to 4.80 g C m⁻² month⁻¹, the ecosystem changed from carbon sink to carbon source (Figure 3).

The seasonal trends of GPP_a and NEE_a in AFE showed significant differences compared to the WWE from March to September (Figure 3). Monthly mean GPP_a showed a unimodal trend, with the minimum value of 14.42 g C m⁻² month⁻¹ in January and the maximum of 121.30 g C m⁻² month⁻¹ in July (Figure 3a). Conversely, NEE_a reached its maximum in January and its minimum in July, with values of −8.56 g C m⁻² month⁻¹ and −92.66 g C m⁻² month⁻¹ (Figure 3b).

The annual sums of GPP_w ranged from 677.49 g C m⁻² a⁻¹ to 784.96 g C m⁻² a⁻¹, with the minimum and maximum occurred in the 2004–2005 and 2008–2009 seasons, respectively. The maximum GPP_a value occurred in the 2004–2005 season, which was 798.94 g C m⁻² a⁻¹ (

Table 1. Annual summaries of GPP, NEE, ET, and WUE of the winter wheat and agroforest ecosystems.

Year	GPP (g C m−2 a−1)		NEE (g C m−2 a−1)		ET (mm)		WUE (g C kg−1 H2O)
	GPPw	GPPa	NEEw	NEEa	ETw	ETa	WUE w	WUE a
2004–2005	677.49	798.94	−450.73	−637.25	413.24	453.44	2.31	2.68
2005–2006	732.80	750.49	−437.35	−545.42	374.23	376.64	2.34	2.48
2006–2007	699.50	624.16	−357.70	−448.59	338.32	382.82	2.9	2.07
2007–2008	726.91	757.11	−464.88	−595.65	403.21	395.86	2.37	2.63
2008–2009	784.96	684.68	−520.72	−518.52	432.09	441.35	2.24	2.12
2009–2010 *	542.90	385.90	−327.46	−256.69	310.63	298.22	2.10	1.79
Average **	724.33	723.08	−446.28	−549.08	392.22	410.02	2.43	2.40

* Data statistics were from September to the next August during the 2004–2009 seasons, and from September to June in the 2009–2010 season. ** The average data were calculated for the 2004–2009 seasons.

). The difference between the annual average GPP was not significant, while the NEE values showed a larger difference, with the mean annual NEE_w and NEE_a being −446.28 and −549.08 g C m⁻² a⁻¹, respectively. This indicates that both of the WWE and the AFE are strong carbon sinks during the wheat growth seasons.

3.2. Characteristics of the Temporal Dynamics in ET

The dynamics of ET_w and ET_a showed obvious seasonal trends (Figure 2d), which were in the form of sinusoidal curves (Figure 4). The minimum values of monthly ET_w and ET_a occurred in January, while the maximum values occurred in May for ET_w (61.47 mm month⁻¹), and in July for ET_a (71.27 mm month⁻¹), respectively (Figure 4).

During the 2004–2010 growth seasons, the phase shifts (S_ET) of monthly ET_w varied from 2.8 to 4.2, while it was 4.0 in each year for ET_a. The amplitudes (δ_ET) in the WWE and AFE ranged from 0.53 to 0.86, and 0.75 to 1.01, respectively (

Table 2. Characteristic parameters calculated by the sinusoidal functions for ET in different growth seasons.

Seasons/ Parameters	ETw				ETa
	SET	δET	R2	RMSE	SET	δET	R2	RMSE
2004–2005	4.0	0.86	0.92	5.98	4.0	1.01	0.97	6.06
2005–2006	2.8	0.72	0.76	8.89	4.0	0.81	0.88	6.34
2006–2007	4.2	0.53	0.86	5.62	4.0	0.75	0.95	6.28
2007–2008	4.0	0.79	0.85	8.12	4.0	0.91	0.90	6.59
2008–2009	3.0	0.75	0.78	9.58	4.0	0.89	0.95	5.00
2009–2010 *	3.0	0.79	0.88	6.98	4.0	0.87	0.95	4.68

* Data statistics were from September to the next August during the 2004–2009 seasons, and from September to June in the 2009–2010 season. SET is the phase shift of monthly ET, δ_ET is the dimensionless seasonal amplitude, R² is the determination coefficient, RMSE is the root mean square error.

), indicating that the inner-annual variability and fluctuation of ET_a was larger. The minimum and maximum amplitudes for both ecosystems occurred in the 2006–2007 and 2004–2005 growth seasons. On an annual scale, the annual sums of ET_w and ET_a were 392.22 and 410.02 mm, respectively, the ranges were from 338.32 mm to 432.09 mm and from 376.64 mm to 453.44 mm (Table 1).

3.3. Credibility of the Sinusoidal Function in ET Simulation

The model testing indicated that the simulated values from the sinusoidal function were close to the measured monthly ET in both ecosystems, with the determination coefficients (R²) being 0.85 and 0.94, respectively (Figure 5). During the 2004–2010 growth seasons, the comparisons between the simulated and measured ET values showed that the R² and RMSE ranged from 0.76~0.92 and 5.62~9.58 mm month⁻¹ for ET_w; the ranges were 0.88~0.97 and 4.68~6.59 mm month⁻¹ for ET_a, respectively. Moreover, the dynamics of the simulated ET using the sinusoidal function were in good agreement with the measured values for both the WWE and AFE (Figure 5).

3.4. Seasonal and Interannual Patterns of WUE

Figure 6 showed that there was a strong positive correlation between daily GPP and ET, the overall slopes for WWE and AFE were 1.71 and 1.69, respectively. The AFE had much better linear fit according to the determination coefficient; however, the determination coefficient of the WWE was enhanced if only when considering the wheat growing stage (Figure 6). The seasonal trend of WUE_w was bimodal, with peak values of 3.94 and 3.65 g C kg⁻¹ H₂O occurring in November and April, respectively, and the lowest values of WUE_w occurred during the summer fallow stage (Figure 7). The monthly WUE_a had one single peak of 4.07 g C kg⁻¹ H₂O in January. The peak values of WUE in both ecosystems occurred in the winter and spring months, when the GPP and ET were quite low, and did not coincide with the GPP and ET peaks. The WUE_w was significantly larger than WUE_a in November, April, and May, but was significantly smaller from June to September (Figure 7).

On an annual scale, the averaged WUE_w ranged from 2.24 to 2.90 g C kg⁻¹ H₂O, and from 2.07 to 2.68 g C kg⁻¹ H₂O for WUE_a, with the average values being 2.43 and 2.40 g C kg⁻¹ H₂O, respectively (Table 1).

3.5. Effects of Environmental Factors on Ecosystem Water–Carbon Fluxes

For both the WWE and AFE, there were positive correlations between GPP, ET, and most environmental factors, while there were negative correlations between WUE and the environmental factors. Specifically, the correlation coefficients between GPP_w and T_a, R_n, PAR, W_s, and T_s ranged from 0.27 to 0.52, and from 0.34 to 0.89 for GPP_a, PAR, and T_s, exhibiting a highly significant positive correlation with GPP_w and GPP_a (Figure 8), with the correlation coefficients being 0.52 and 0.89, respectively. The correlation coefficients between T_a, R_n, PAR, W_s, VPD, T_s, SWC, and ET_w ranged from 0.32 to 0.89, and from 0.31 to 0.95 for ET_a, where R_n and T_s were the major factors, with the correlation coefficients being 0.89 and 0.95, respectively (Figure 8). Additionally, these factors demonstrated significant negative correlations with WUE_w and WUE_a; the correlation coefficients ranged between −0.37 to −0.68 (Figure 8a) and from −0.16 to −0.56 (Figure 8b), respectively. VPD and SWC were the major influencing factors for WUE_w and WUE_a, with the correlation coefficients being −0.68 and −0.56, respectively.

4. Discussion

4.1. Possible Impacts of Crop Growth Stage on Ecosystem WUE

In the present study, the seasonal dynamic trends of GPP in the WWE were bimodal, while they were unimodal in the AFE; such obvious differences might be caused by different water and carbon uptake capacities of plants at different growth stages (Figure 2 and Figure 3). As a result, the ecosystem WUE varied with different seasons within the year (Figure 7). During the wintering stage (December to the next February), ET and GPP were very low because of the dormancy of the crops [20]. Afterwards, with the increasing temperature and vigorous crop growth, ET and GPP increased quickly to the maximum of the year. A dual relationship between GPP and ET appeared in Figure 6, which is due to the decreased ecosystem productivity during the fallow period after the wheat harvest, but the ET is still high under high temperatures. However, the overall trend of the AFE was less affected by the growth and harvest of the winter wheat, indicating that the source area was much larger when the observation height increased from 2 m to 32 m [35,36]. The maximum WUE and GPP appeared at different times, and the maximum WUE occurred at lower GPP stages when the temperature was nearly 0 °C. This kind of inconsistency might occur because the trade-off between carbon assimilation and transpiration of the WWE and AFE was influenced by soil water limitation or drought, which is consistent with the results of Gao et al. (2017) [49]. The relationship between crop biomass formation and water consumption should be like this, but it may also be related to crop type.

4.2. Complex Responses of Water and Carbon Fluxes to Different Ecosystem Types

Although we cannot predict future ET variations using the sinusoidal model, it provided a new insight by describing ET variations under different ecosystems. Although the simulated ET in both ecosystems, composed of different vegetation types, had similar sinusoidal trends, their characteristic parameters showed large differences. In the sinusoidal simulations, the ET of the AFE showed a better fit between the observed and simulated values than that of the WWE, with a larger R² and constant phase shifts (Figure 5, Table 2). Previous studies reported large variations in ET, GPP, and NEE due to differences in geographical climate, soil type, crop rotation system, field management, and data processing methods. Most of the ecosystems behaved as carbon sinks, and the main difference between agroecosystems and natural ecosystems in terms of regional-scale carbon budget is that crops are harvested and transformed to CO₂, which is again released to the atmosphere, many ecosystems changed from carbon sink to source once this element of carbon cycling was taken into consideration [50,51]. The results of Veeck et al. (2022) showed that, taking into account the carbon release in fallow periods, a soybean–wheat succession ecosystem in Brazil was a net C sink with NEE at −50 g C m⁻² a⁻¹. Lei and Yang (2010b) [21] obtained similar results in a wheat and maize rotation system in the North China Plain. However, this study was mainly concerned with ET and GPP correlations for the two kinds of ecosystems; the grain carbon transformation needs further calculation and comparison.

4.3. Implications and Further Scopes of This Study

Understanding the dynamics of ET, GPP, NEE, and their related factors in agroecosystems is important for evaluating the impact of climate change on regional carbon sequestration and energy balance [52,53]. However, the mechanisms of the water and carbon exchanges for the WWE and AFE are still poorly understood and require further study.

First, the water and carbon exchanges between agroecosystem and the atmosphere were affected by many factors, such as geographic location [29], vegetation ages [13], ad crop growing and rotation [20,21], as well as climate factors like summer drought [48] and rainfall events [54]. This study discussed the influence of some factors on WUE; however, the contributions of each factor need to be furtherly quantified, and the dominant factor of WUE for the WWE and AFE identified in order to assess the effect of climate change on the water and carbon fluxes of agroecosystems.

Second, the energy balance closure has a significant influence on the calibration and validation of models which are made for weather prediction or climate models [53]. This study uses two sets of equipment installed at different heights covering different scales; it is important to make comparisons for the energy balance closure of the WWE and the AFE in this region in order to provide parameters for model validation.

Finally, at the annual scale, the WUE calculated from crop yield (or total biomass production) and water consumption, as well as the carbon budget from the total production, can be seen as the validation of the results of eddy covariance. Further studies should incorporate crop productivity, empirical models, and increasing replicates to minimize flux loss and uncertainty [55,56], which will also provide the basis for the development of local crop moisture–yield management and water-saving agriculture.

5. Conclusions

The interactions between carbon–water fluxes and their influencing factors vary significantly among different ecosystems, and these differences reflect various responses to climate change. In this study, the inter-annual variation characteristics of GPP, NEE, ET, and WUE were analyzed for a typical rain fed winter wheat ecosystem and an agroforest ecosystem in the southern Loess Plateau. Quantified correlation using the Pearson correlation analysis method and the sinusoidal function were used to simulate the interannual variation in ET. We found that the seasonal dynamics of GPP_w and NEE_w were closely related to crop development and phenology, and the ecosystem changed from a carbon sink to a carbon source after the wheat harvest. On an annual scale, both the WWE and the AFE are strong carbon sinks, with the AFE fixing more carbon. The sinusoidal function can be used in seasonal ET simulations for both WWE and AFE, and the inner-annual variability and fluctuation of ET in the AFE was larger with greater amplitudes (δ_ET). There was strong coupling relationship between daily GPP and ET, and WUE_w was significantly larger than WUE_a in November, April, and May, when the winter wheat grows vigorously, but the apple trees were not yet in full bloom. PAR and T_s were most positively correlated with GPP_w and GPP_a, R_n and T_s were the major factors influencing ET_w and ET_a, while VPD and SWC were the major influencing factors for WUE_w and WUE_a, indicating the key effects of SWC on yield formation and WUE in the water-deficient areas of the Loess Plateau. The present study provides new insights into the estimation of carbon and water fluxes in different ecosystem types, which may have important implications for regional carbon neutralization and water cycle modeling under climate change.