Improving the Model Performance of the Ecosystem Carbon Cycle by Integrating Soil Erosion–Related Processes

Abstract

Soil erosion is a key factor in soil quality degradation and carbon balance in arid ecosystems. However, many models ignore the soil erosion process in arid regions, which may lead to limits in our understanding of ecosystem processes in arid regions. In this study, we added the soil erosion process according to field observed data of soil hydrothermal regimes and carbon flux. We validated this coupling version of IBIS (Integrated Biosphere Simulator) and RUSLE (RU–IBIS) by examining four different vegetation types and the carbon budget in the arid region on the Loess Plateau (LP). Our results indicated that the coupling model (RU–IBIS) produced more reliable simulations of the soil water content (with the r from 0.23–0.90 to 0.71–0.97) and evaporation (ET) (the average r was 0.76) and significantly improved the simulation of the leaf area index (LAI) (the average r was 0.95) and net primary production (NPP) (the average r was 0.95). We also conducted sensitivity experiments to determine how soil texture and aerodynamic roughness (Z_0m) affect the soil water content. Moreover, it was revealed that specific leaf area (SLA) plays a key role in the simulation of NPP and NEE. Our study suggests that the coupled soil erosion process and parameterization can effectively improve the performance of IBIS in arid regions. These results need to be considered in future Earth system models.

1. Introduction

As a component of terrestrial ecosystems, arid regions account for about 46% of the global carbon pool [1]. Arid regions have low precipitation and high solar insolation, and their evaporation is much higher than their precipitation, making the ecosystem relatively fragile and more vulnerable to external climate impacts [2]. Insufficient water and human activities such as overgrazing and excessive land reclamation have also led to serious land degradation in arid regions [3]. Land degradation releases carbon into the atmosphere and may form a positive feedback loop within the carbon cycle [4]. Therefore, it is of great importance to study the processes of carbon cycling in arid region ecosystems in the context of climate change.

Simulating and predicting the arid ecosystem carbon cycle is inseparable from the dynamic global vegetation model (DGVM) [5]. Using the DGVM, many studies have been conducted to simulate the change in the carbon budget and its influencing factors in arid regions [6,7,8,9]. Under the influence of long–term soil erosion, the nutrient–rich topsoil layer is gradually denuded, and the soil nutrients are lost to runoff and sediment, which leads to a decrease in soil fertility and thus reduces the erosion resistance of the soil. This results in a vicious cycle of “soil erosion–nutrient loss–soil erosion” [10,11]. Soil erosion leads to the thinning of the soil layer and decreases soil aeration, permeability, and soil quality [12,13]. The soil erosion process strongly affects the energy and water exchange between the soil and the air and thus definitely affects the simulation accuracy of the carbon cycle in arid regions. Soil erosion must be included in the study of the carbon cycle in arid regions [14,15]. However, in many DVGMs, this process has only been expressed as simple soil erosion parameterization schemes [16,17,18].

IBIS was developed to illuminate the interactions between different biosphere factors on time scales relevant to the carbon cycle process [19]. In recent years, some studies have started to apply IBIS on a global scale [20]. Further refinements of the IBIS model include the creation of IBIS 2.0 and Agro–IBIS [21,22,23], which incorporate other models of atmospheric circulation, agricultural vegetation, and hydrological processes. In terms of model improvement, Liu et al. (2005) combined IBIS with the nitrogen cycle model to simulate the regional nitrogen cycle and to reveal the extent to which the nitrogen cycle affects the ecosystem [24]. While most ecosystems examined with IBIS are forests and farmland [25,26,27], recent updates to IBIS include adding more vegetation types and statistical models [28,29]. However, the IBIS model does not consider changes to the soil erosion regime. The absence of this process may produce model results that overestimate the moisture and carbon budget in arid regions.

It is important to use the DVGM to model the carbon cycle in arid ecosystems. Therefore, the main contents of this study include (1) incorporating the revised universal soil loss equation (RUSLE) into the soil erosion and verifying the accuracies of modeling soil water content, LAI, and NPP simulations in the Loess Plateau (LP) and (2) determining the sensitivity of the soil hydrothermal regime, ET (evapotranspiration), and the carbon flux to the changing model parameters. The modified IBIS model provides a useful tool for understanding the carbon cycle of arid regions in warming climate conditions.

2. Materials and Methods

2.1. Experimental Sites

The Loess Plateau is the largest loess accumulation area in the world [30], which is located in north–central of China, with soil erosion covering 71%, and the average thickness of loess in the region is 92.2 m. Due to its vegetation degradation, soil nutrient loss, land quality degradation, and ecological environment deterioration, the LP is very sensitive to climate change. We collected data from four regions in the LP (Figure 1). These four stations are Shapotou (SPT: cropland), Changwu (CW: cropland) Ansai (AS: cropland), and Eerduosi (ERDS: desert). Detailed information on the four stations is summarized in

Table 1. Site information for (SPT, CW, AS, and ERDS) in the LP.

Site		SPT	CW	AS	ERDS
Observation period		2001–2008	2004–2008	2004–2006	2005
Vegetation type		cropland	cropland	cropland	desert
Latitude	N°	37°27′	35°12′	36°51′	39°29′
Longitude	E°	104°57′	107°40′	109°19′	110°11′
Elevation	m	1250	1200	1068	1290
Mean precipitation	mm	186	580	500	310
Mean temperature	°C	9.6	9.1	8.8	6

.

2.2. Input and Observed Data

The meteorological data were collected at four stations operated under the auspices of the National Ecosystem Science Data Center, National Science and Technology Infrastructure of China (http://www.nesdc.org.cn, accessed on 1 March 2023). The period of single–point simulation is 1998–2008. The meteorological data include wind speed, relative humidity, pressure, radiation, daily maximum temperature, daily minimum temperature, daily average temperature, and precipitation. The basic processing procedures include anomaly rejection, blank interpolation, half–hourly synthesis of daily data, and unit conversion. The processed meteorological data are used as single–point simulated climate–driven data. Model validation data include soil water content, evapotranspiration, above–ground biomass, below–ground biomass, and leaf area index (LAI) [31,32,33,34]. In this study, the measured biomass values were converted to NPP values for model validation [35].

2.3. Model Description and Development

Foley et al. (1996) developed the biosphere–oriented, regional–scale Integrated Biosphere Simulator (IBIS) model [20]. In order to simulate the dynamic changes in the vegetation cover in Europe and Africa, the second iteration of the IBIS model was coupled with GENESIS (version 2.0), an atmospheric circulation model. The most recent version of IBIS is version 2.6 [35]. The IBIS model consists of five modules: land surface processes, vegetation dynamics, vegetation phenology, subterranean carbon and nitrogen cycling, and soil geochemistry (Figure 2). Our study focuses on improving and refining the land surface processes and vegetation dynamics modules.

(1) Land surface processes module

The land surface processes module is used to calculate the exchange of water, heat, and carbon dioxide. The thicknesses of the six soil layers, from top to bottom, are 0.10 m, 0.15 m, 0.30 m, 0.50 m, 1.00 m, and 2.00 m. In the soil module, the rate of change in the liquid soil water content is calculated using the Richard equation, and the vertical flux of water is determined using Darcy’s law. At the impermeable bottom boundary of the soil module, there is no diffusion of liquids or heat.

To simulate the soil hydrothermal regime, we calculated the soil temperature, the water content, and the ice content of a given soil layer with the Richard equation as follows [20]:

$$\frac{\partial \theta }{\partial t}=\frac{\partial }{\partial z}[D(\theta )\frac{\partial \theta }{\partial z}]+\frac{\partial K(\theta )}{\partial z}-ST(z,t)$$

(1)

where θ (m³m⁻³) is the actual soil water content of a soil layer at a certain time and θi (m³m⁻³) is the ice rate of the soil layer. D(θ) (mms⁻¹) represents the soil diffusivity of the soil volumetric water content. K(θ) (mms⁻¹) is the hydraulic conductivity; when the soil volumetric moisture content is K(θ), it needs to be calculated by the saturated hydraulic conductivity K_s. b is the index of the Campbell formula. ST(z,t) (s⁻¹) is a sink term resenting water uptake by plant roots as a function of ET. The Campbell formula is as follows:

$$K(\theta )=K_{\mathrm{S}}{\left(\frac{\theta }{{\theta }_{\mathrm{S}}}\right)}^{2b+3}$$

(2)

Using the heat diffusion equation [20], we can quantify the heat transfer in the soil:

$$C\frac{\partial T(z,t)}{\partial t}=\frac{\partial }{\partial z_{}}\left[\frac{\lambda \partial T(z,t)}{\partial z}\right]+H_{fm}$$

(3)

$$H_{\mathrm{fm}}={\rho }_{ice}{L}_f\frac{\partial {\theta }_{ice}}{\partial t}$$

(4)

where C (Jkg⁻¹k⁻¹) and λ (Wm⁻¹k⁻¹) are the volumetric heat capacity and the thermal conductivity, respectively, t (s) is the time, and z (mm) is the vertical depth. T (k) is the temperature. H_fm (Jkg⁻¹) is the latent heat released or consumed during the phase change of soil water. is the density of ice, The density of ice is ρ_ice (917 kgm⁻³), θ_ice (m³m⁻³) is the partial volume of ice content, and L_f (0.3336 × 10⁶ Jkg⁻¹) is the latent heat of fusion.

(2) Vegetation dynamic module

After registering a change in the biomass for each plant type, this module calculates the gross primary productivity (GPP), the net primary productivity (NPP), and the net ecosystem exchange (NEE) in the vegetation coverage and type. The vegetation productivity calculation includes both the total photosynthesis and the breathing per hour. Vegetation biomass is the sum of the root system biomass and the plant leaves. The GPP calculation [20] is:

$$GPP={\displaystyle \int A_g}dt$$

(5)

and the NPP calculation is:

$$NPP=(1-\eta ){\displaystyle \int (A_g}-R_{\mathrm{leaf}}-R_{\mathrm{stem}}-R_{root})dt$$

(6)

where A_g (CO₂/m²·s) is the net photosynthetic rate of the canopy; R_leaf (CO₂/m²·s) is the leaf respiration rate; R_stem (CO₂/m²·s) is the respiration rate; and R_root (CO₂/m²·s) is the root respiration rate. The rate of photosynthesis is simulated according to the Farquhar equations [36]. In accord with Amthor [37], η is 0.33 or 0.3, which is the proportion of carbon loss in tissue construction.

(3) Vegetation Phenology Module

In the phenological module, it is assumed that deciduous plants begin to germinate in winter when the growth of a particular plant functional type exceeds a certain threshold. This threshold is defined by the day when two conditions are met: (1) the average temperature (Avt10) is greater than 5 °C for 10 consecutive days and (2) at least 100 days have experienced temperature increases greater than 0 °C (GDD₀). When Avt10 is less than 5 °C, the vegetation enters the withered yellow period [38]. This threshold day forms the basis of the leaf area index (LAI) calculations [20]:

$$LAI=LA{I}_{\max }\times \mathrm{f}(T)$$

(7)

$$\begin{array}{c}\hfill \mathrm{f}\left(T\right)=\left\{\begin{array}{c}\left(GD{D}_0-T_{\mathrm{b}}\right)/b\hspace{1em}GD{D}_0>T_b \mathrm{and} avt10>5\hfill \\ \hspace{1em}\hspace{1em} \hspace{0.17em}1\hfill \\ \mathrm{f}\left(T\right)-1/dfac\hspace{1em}avt10<5\hfill \end{array}\right.\end{array}$$

(8)

where T_b is the threshold of GDD₀ in the vegetation green–back period, b is the coefficient of controlling the leaf spreading speed, and dfac is the coefficient of controlling the defoliation speed.

(4) Modified soil erosion parameterization scheme

Numerous studies have pointed out that soil erosion and the erosion of poorly vegetated land are accompanied by changes in soil nutrients, structure, moisture, texture, etc. The soil microbial community will also be transformed during the restoration of soil erosion and erosion–damaged land [39,40]. As the main driver of soil organic carbon transformation, microbial community changes in erosion and revegetation systems largely determine the rate of organic carbon mineralization and fixation and the microbial contribution to soil organic carbon. The microbial communities of erosion and revegetation systems largely determine the rate of organic carbon mineralization and fixation and the apparent functional redundancy of microorganisms in the soil. Therefore, we used RUSLE to improve the soil erosion process module. The expression of general soil loss equation is:

$$A=R·K·LS·C·P$$

(9)

where A is the soil loss (kg km⁻²), R is the rainfall erosivity factor (J m⁻²), K is the soil erodibility factor (dimensionless), LS is the slope length and gradient factor (dimensionless), C is the vegetation cover and management factor (dimensionless), and P is the soil and water conservation measures factor (dimensionless).

(5) Model parameterization

Previous modifications of the IBIS model [41] suggested that accurate model results were dependent upon the accurate calibration of stomatal conductance coefficients (m, b), leaf respiration parameter (γ), and the maximum rubisco capacity (Vmax). For the photosynthesis of vegetation in arid regions, these parameters were set to γ = 0.025, m = 4.0, b = 0.04, and Vmax = 25 molCO₂ m⁻² s⁻¹ [42]. We set these four parameters in the vegetation dynamic module.

Aerodynamic roughness is a parameter that is important for describing the aerodynamic characteristics in surface heat flux parameterization schemes, which results in large uncertainties in soil hydrothermal condition calculations at smaller spatial scales [43]. Therefore, we changed the initial Z_0m (0.005 m) to be equal to the two annual average Z_0m for the two vegetation types along LP based on the measured data, the cropland (0.29), and the desert (0.0003) (Table 2).

The specific leaf area (SLA) can be used not only to infer the processes of photosynthesis, evaporation, and ET but also to estimate the net productivity of terrestrial ecosystems. In an individual plant or group of plants, the amount of light received by the vegetation is inversely proportional to the SLA value. The SLA value can also reflect the ability of plants to obtain light; plants with low SLA values can better adapt to poor treatment and drought environments, while plants with high SLA values can more easily preserve vital nutrients [44,45]. This study used SLA that correspond to the vegetation conditions in arid regions to improve the quality of our model results (Table 2) [46,47].

Table 2. Parameter differences between the O–IBIS and RU–IBIS models (M–IBIS).

Site		SPT		CW		AS		ERDS
		O–IBIS	RU–IBIS	O–IBIS	RU–IBIS	O–IBIS	RU–IBIS	O–IBIS	RU–IBIS
Z0m	(m)	0.005	0.29	0.005	0.29	0.005	0.29	0.005	0.0003
SLA	(m2/kg)	20	16.04	20	11.69	20	10.9	20	11.69

Note: Z_0m: thermodynamic roughness. SLA: specific leaf area. O–IBIS: original IBIS. RU–IBIS: modified IBIS. References: [46,47].

Table 2. Parameter differences between the O–IBIS and RU–IBIS models (M–IBIS).

Site		SPT		CW		AS		ERDS
		O–IBIS	RU–IBIS	O–IBIS	RU–IBIS	O–IBIS	RU–IBIS	O–IBIS	RU–IBIS
Z0m	(m)	0.005	0.29	0.005	0.29	0.005	0.29	0.005	0.0003
SLA	(m2/kg)	20	16.04	20	11.69	20	10.9	20	11.69

Note: Z_0m: thermodynamic roughness. SLA: specific leaf area. O–IBIS: original IBIS. RU–IBIS: modified IBIS. References: [46,47].

2.4. Model Validation

2.4.1. Correlation Analysis

We used many statistical methods to quantitatively evaluate the applicability of our models, including the relative error (RE), the correlation coefficient r, the root mean square error (RMSE), and the coefficient of determination R² between the simulated and measured values. Detailed descriptions of these statistical tests and formulas are found in Willmott et al. (1982) [48].

2.4.2. Parameter Sensitivity Analysis

To investigate how parameters within IBIS could influence the simulation results, we ran 16 experiments. Individual model parameters were individually increased or decreased by 20% in each run. Overall model sensitivity to NPP, ET, SW, and ST (multilayer soil depth averages) was evaluated by comparing the simulated annual averages from the modified parameter runs to the control runs.

2.4.3. Redundancy Analyses

We conducted redundancy analyses (RDAs), both with and without variance partitioning to estimate both the climate effects and the nature of the relationship between the site and the climate effects [49]. Here, RDAs were used to analyze the relationships between carbon flux and climate variables. These analyses were performed using the vegan package in R–4.2.2.

3. Results

3.1. Soil Water Content Conditions

A comparison of the simulated and measured soil water content values at the four sites is shown in Figure 3, where the simulated daily temperature trends are in good agreement with the observations. RU–IBIS generated improved soil water content estimates for the four stations, with a higher range of r values (from 0.23–0.90 to 0.71–0.97) and RMSE, with RE values of 0.15–5.69 m³/m³ and −0.12–0.27 m³/m³, respectively (

Table 3. Statistical measures of modeled soil water content at a monthly time scale.

Stations		10 cm		30 cm		50 cm		100 cm
		RU–IBIS	O–IBIS	RU–IBIS	O–IBIS	RU–IBIS	O–IBIS	RU–IBIS	O–IBIS
SPT	r	0.71	0.45	0.88	0.78	0.91	0.82	0.92	0.90
	R2	0.51	0.20	0.78	0.62	0.84	0.67	0.84	0.82
	RMSE	3.23	4.42	0.71	1.42	0.60	1.43	1.05	3.58
	RE	−0.08	0.09	−0.12	0.17	−0.12	−3.02	0.17	−0.04
AS	r	0.87	0.71	0.88	0.58	0.86	0.60	0.82	0.59
	R2	0.75	0.51	0.77	0.33	0.73	0.36	0.68	0.34
	RMSE	2.67	3.62	1.77	5.60	1.81	5.25	1.89	4.81
	RE	0.07	−0.01	−0.03	−0.46	0.08	−0.31	−0.07	−0.43
CW	r	0.94	0.88	0.97	0.46	0.90	0.65	0.91	0.63
	R2	0.88	0.77	0.94	0.21	0.82	0.42	0.82	0.40
	RMSE	5.69	7.26	1.44	5.29	0.64	3.78	1.06	3.69
	RE	0.27	0.37	−0.05	0.09	−0.11	−0.42	−0.09	−0.13
ERDS	r	0.84	0.34	0.87	0.74	0.91	0.60	0.76	0.23
	R2	0.70	0.11	0.76	0.55	0.83	0.35	0.58	0.06
	RMSE	1.12	1.67	0.29	0.32	0.15	0.26	0.34	0.73
	RE	−0.08	−0.19	−0.10	−0.39	0.07	−2.91	0.19	−0.35

). The underestimation of soil water content was greatly improved. In terms of different depths, the correlation between the simulated results and the measured values of soil water content at 10 cm is the weakest. With the increase in soil depth, the deviation value decreases sequentially with the increase in soil depth, which is related to the fact that the deeper the soil depth, the smaller the actual amplitude. The RU–IBIS model has the highest simulation accuracy in the cropland (CW station), and the simulation effect is relatively poor in the desert region (ERDS). Our statistical analyses for these parameters are shown in Table 3.

3.2. ET Simulation Results

Due to the limitation of the monitoring conditions, only monthly average daily ET data of SPT (2003–2004), AS (2004–2006), and CW (2004) stations were observed. In Figure 4, it can be seen that the simulation result of evapotranspiration is higher than the measured data. The r values of ET were 0.61 (SPT), 0.83 (AS), and 0.86 (CW). The RMSE and RE were 2.06 mm, −0.05 mm (SPT), 2.07 mm and −0.20 mm (AS), and 0.75 mm and 0.17 mm (CW) (

Table 4. Statistical measures of the modeled ET.

	SPT				AS				CW
	r	R2	RMSE	RE	r	R2	RMSE	RE	r	R2	RMSE	RE
RU–IBIS	0.61	0.38	2.06	−0.05	0.83	0.69	2.07	−0.20	0.86	0.75	0.75	0.17
O–IBIS	0.23	0.05	2.97	−23.93	0.41	0.17	3.49	0.76	0.59	0.35	1.18	0.95

).

3.3. LAI and NPP Simulation

The original IBIS underestimated the LAI values; the RU–IBIS values, which had smaller uncertainties than the O–IBIS values, generally agreed with the observed LAI values (Figure 5). The RU–IBIS model generated LAI correlation coefficient (0.95) and R² (0.90) results that were higher than those of the O–IBIS model (0.70 and 0.50, respectively). Furthermore, the RMSE and RE values simulated by the RU–IBIS results revealed improving O–IBIS results on LAI and NPP bases (

Table 5. Statistical analyses of modeled monthly LAI and NPP.

		SPT				AS
		r	R2	RMSE	RE	r	R2	RMSE	RE
LAI	RU–IBIS	0.94	0.88	0.60	0.02	0.96	0.92	0.47	−0.06
	O–IBIS	0.77	0.60	1.39	0.56	0.63	0.40	1.46	−9.67
NPP	RU–IBIS	0.96	0.92	474.11	0.29	0.93	0.87	318.18	−0.38
	O–IBIS	0.23	0.05	1188.11	−0.43	0.81	0.65	467.82	0.38

). Clearly, the RU–IBIS dynamic vegetation results are more accurate than those of the O–IBIS model.

Figure 5c,d shows the results between the observed O–IBIS and RU–IBIS NPP at SPT and AS stations, which demonstrated that the RU–IBIS simulations maintained consistency with the observed NPP. The NPP of O–IBIS was much lower than the observed value. After modifying IBIS, the underestimation phenomenon improved significantly. The correlation coefficients were 0.52 (O–IBIS) and 0.95 (RU–IBIS), respectively. The statistical metrics for accuracy and error improved when we used the RU–IBIS results instead of the O–IBIS results.

3.4. Parameter Sensitivity

As shown in

Table 6. The results of the model parameter sensitivity experiments.

Sites	Parameters	Value	Sensitivity (+20%)						Sensitivity (−20%)
			ST	SW	ET (%)	NPP	NEE	Re	ST	SW	ET (%)	NPP	NEE	Re
SPT	Z0m (m)	85	−26.68 ※	86.43	42.71	1.02	−18.25	38.07	−9.62	38.02	16.84	−0.31	−23.09	46.37
	SLA (m2/kg)	6.85	0.17	0.00	0.01	−0.01	−0.15	0.01	−1.20	−0.01	−0.02	0.03	0.28	−0.02
AS	Z0m (m)	66.4	−0.42	−28.13	−15.64	0.16	−10.70	−32.72	8.35	0.39	0.13	−0.17	0.56	0.07
	SLA (m2/kg)	3.55	0.00	0.00	0.00	0.00	0.00	0.00	−0.31	0.01	0.00	0.01	−0.01	0.00
CW	Z0m (m)	80	−1.30	−15.81	−5.27	0.49	−62.27	−23.44	5.69	44.15	22.12	−0.27	10.93	41.47
	SLA (m2/kg)	4.34	0.96	0.00	0.01	−0.02	0.02	0.00	−0.37	0.00	−0.02	0.01	−0.05	−0.01
ERDS	Z0m (m)	40.2	13.15	−22.28	9.20	0.60	−53.31	30.55	3.17	0.10	0.10	−0.11	−4.06	0.12
	SLA (m2/kg)	54.05	22.19	−22.16	9.46	0.34	−54.24	30.81	−4.83	−0.17	−0.09	0.15	7.20	−0.25

Note: +20%: model parameters increased by 20%; −20%: model parameters decreased by 20%; −26.68 ^※: when the Z_0m increases by 20%, the soil temperature decreases by 26.68%.

, the simulation results are most sensitive to changes in the Z_0m. For the cropland, a 20% increase in the Z_0m decreased ST and NEE by 26.68% and 18.25%, respectively, and increased SW, ET, and Re by 86.43%, 42.71%, and 38.07%, respectively. Therefore, the accuracy of the input parameters plays an important role in the simulation results. While the +20% change in the SLA value resulted in a moderate NPP increase of 5–8%, the sensitivity of NEE to SLA is much higher than that of NPP.

3.5. RDA Analysis

Figure 6 and Figure 7 show that the arid region ecosystem is highly sensitive to changes in temperature. While warming temperatures are minimally impacted by precipitation, NPP tends to rise along with the temperature. The two consecutive RDA axes (Figure 6) suggest the pattern of carbon fluxes with climate change and illustrate that the first RDA axis (66.84%) is absolutely dominant. The second RDA axis accounts for 17% of the carbon flux variation. Carbon flux variation on the second axis is indicative of warmer, higher radiation ecosystems with higher NPP, NEE, and Re values. Cloud amount, temperature, and soil temperature have the greatest influence on the carbon flux. In the LP, vegetation growth and soil respiration are very sensitive to temperature and radiation changes.

The pie chart in Figure 7 summarizes the contributions of climate and site to the different ecosystem carbon fluxes. Climate is shown nested within research sites because differences in climate are determined by the site location. The three temperature parameters jointly account for 31% of the NEE variation (Figure 7), the most important controlling factor of NPP is the cloud amount (29.2%), and soil respiration is most sensitive to changes in the site (63.2%). Different sites mean different soil textures, vegetation types, and ecosystems. Different ecosystems and ecosystem locations lead to different growth season onsets, growth season durations, and carbon emission intensities.

4. Discussion

4.1. Improvement of the IBIS Model

Our results revealed that soil erosion in O–IBIS is not realistic, as it largely underestimates soil water content (Figure 3). Also, the simulated ET, LAI, and NPP in O–IBIS are significantly lower than our observations (Figure 4 and Figure 5). We tentatively attribute this deficiency to the original IBIS, as most of the soil hydrothermal regime improvements and the soil erosion process of arid regions in IBIS require more validation and improvement. Previous works have suggested that the reduction in soil fertility and the deterioration of soil structure under the influence of long–term hydraulic erosion structural deterioration will further affect the growth and development of plants. In our study, we not only added the soil erosion process to the IBIS model, but we also recalibrated the physiological parameters according to the actual climate conditions in the study area. These modifications resulted in better fits with the observed data and reduced model uncertainties. We suggest further collection of such detailed ecological factor data at multiple locations to improve the arid vegetation modeling and reduce errors due to multi–site heterogeneity.

Further development always requires field–measured data for calibration. The LP is dominated by the C3 plants. The physiological parameters of different vegetation types differ greatly between the C3 plants. The RU–IBIS model, with its more complete soil erosion process, may provide better opportunities for parameter calibration in arid regions.

4.2. Parameters Affecting the Simulation Results

Soil water content and ET are greatly affected by dynamic roughness, which reflects the weakening effect of different ground surface characteristics on the wind speed [43]. The vegetation characteristics and the wind speed cannot be ignored in sparsely vegetated areas. Z_0m is negatively correlated with wind speed, and wind speed is positively correlated with ET [50]. An increase in the Z_0m parameter represents a decrease in the wind speed, which in turn leads to a reduction in the ET value.

SLA values can directly affect the efficiency of vegetation in capturing light and CO₂ [49]. Higher SLA values mean that plants usually have higher productivity values and that plants with lower specific leaf area values experience more efficient photosynthesis in limited–resource conditions [51]. SLA is involved in the area–based biochemical parameterization of photosynthetic enzyme concentrations linking LAI to NPP and provides structural parameters linking leaf carbon distribution to LAI [44,52].

4.3. Potential Effects of Climate Change on the Carbon Cycle

With the concentrated rainfall, heavy storms, and human interference in the Loess Plateau, soil erosion of the LP is inevitable. Improving the activity and quantity of soil microorganisms results in accelerated soil respiration and CO₂ emissions. Additionally, the organic carbon in the soil, which was originally sealed in the arid region, was released into the ground surface and eventually made its way into the carbon cycle. However, the temperature increase can improve the photosynthetic rate of vegetation and prolong the growth period of vegetation, resulting in a longer period over which vegetation can draw CO₂ out of the atmosphere. Similarly, arid region erosion results in decreased soil water content, which inhibits vegetation growth and CO₂ absorption. Increased soil erosion can reduce the soil water content volumetric under the same suction, which is unfavorable to soil water storage and reduces the drought–resistant performance of the soil. The infiltration rate of soil in the 0–10 cm layer gradually decreased with the increase in soil erosion, and the difference in infiltration rate of each soil layer in different vertical depths gradually decreased with the increase in soil erosion [49]. Therefore, the potential impact of arid region erosion on ecosystem carbon emissions mainly depends on the arid region balance between the increase in the soil respiration potential caused by warming and the diminution of vegetation photosynthesis in conjunction with the increase in the soil water content [52,53], which is supported by additional field data from the LP [54]. NEE is numerically equal to the difference between NPP and soil heterotrophic respiration. Therefore, NPP could be largely canceled out by a higher ER value, which tends to respond positively to the ecosystem’s annual carbon emissions [52].

5. Conclusions

The IBIS model has been validated in a variety of different ecosystems, regions, and temperature zones worldwide. We used RU–IBIS to determine the appropriateness of applying this model to arid regions in the LP. The results showed that:

• The O–IBIS model generally underestimates the soil water content in all stations. After adding the soil erosion process scheme and adjusting the dynamic roughness, the ability of the IBIS model to simulate soil water content was greatly improved, though still with underestimated bias (RE = −0.013 m³/m³). These results indicate that soil erosion plays an important role in controlling soil hydrothermal regimes in arid regions. We should pay attention to the reliability of the Z_0m value in the future simulation of the barren region with IBIS.
• The underestimation of NPP and LAI was reduced by 68.7% and 88.5%, respectively, through the adjustment of the SLA. The modified model (RU–IBIS) accurately reproduces the vegetation growth conditions in the cropland region and does a poor job of constraining the conditions in the desert. The results show that the complexities of applying to heterogeneous sites remain, and it needs further improvement in a wider range area.
• Sensitivity tests showed that the soil hydrothermal conditions and ET were sensitive to changes in the Z_0m values. SLA largely affects the simulated NPP and NEE at all stations. An increase in the SLA produces corresponding increases in the NPP and NEE values.
• We infer that the temperature, radiation, and soil water content dictate the inter–annual variability of the carbon emission (absorption) strength in this simulation. Moreover, any soil water content lost during soil erosion may exacerbate water deficits and thus reduce NPP. However, it is unknown whether the arid region erosion will result in a net increase or decrease in the total carbon emissions in this region.

The RU–IBIS model improves the accuracy of the soil and vegetation modules. By incorporating the soil erosion process into IBIS, we hope to extend our single–site simulation to the whole arid region of the LP and explore the feedback of arid regions on global climate change in the future.