**1. Introduction**

The carbon cycle is a popular topic in ecological research [1]. Soil organic carbon (SOC) is the largest carbon pool on the planet excluding the ocean's and rock's sediments; thus, small changes in SOC have a great impact on the atmosphere [2]. The carbon pool of the agro-ecosystem is one of the most active parts of the global carbon cycle, in which soil organic carbon storage in farmland accounts for 8–10% of that in all types of land [3]. The SOC in farmland is vulnerable to disturbances from human activities [4], but this SOC can be artificially regulated on a short-time scale [5]. In addition, China has a total paddy soil area of 45.7 Mha, accounting for approximately one-fifth of the total cultivated land area in the world [6]. Thus, paddy fields have a considerable carbon sequestration potential. At the same time, the dynamics of SOC in paddy fields are affected by many factors, such as

**Citation:** Jiang, Z.; Yang, S.; Ding, J.; Sun, X.; Chen, X.; Liu, X.; Xu, J. Modeling Climate Change Effects on Rice Yield and Soil Carbon under Variable Water and Nutrient Management. *Sustainability* **2021**, *13*, 568. https://doi.org/10.3390/ su13020568

Received: 20 November 2020 Accepted: 29 December 2020 Published: 8 January 2021

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2021 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

temperature, precipitation, irrigation, and fertilization [7]. However, few studies on SOC changes in paddy fields have focused on the impact of coupling water-saving irrigation and fertilizer management. In recent years, water-saving irrigation technology has been widely used in China and has changed the soil moisture status and organic carbon content. Thus, evaluating the impact of water and carbon management measures on the dynamic changes in SOC is important to maintaining agricultural productivity.

Moreover, our understanding of climate change as an important factor affecting SOC and rice yield remains limited. Thus, improving our understanding of the impact of environmental change and field management on nutrient cycling and crop growth is necessary. Despite the growing importance of industry, agricultural production, as one of the most sensitive sectors to climate change [8], plays an important role in ensuring food security throughout the world, especially in China [9]. Rice paddies are an important source of both global food production and greenhouse gas (GHG) [10,11], and rice yield is extremely sensitive to agricultural measures, such as irrigation and fertilization [12]. At present, China's sustainable agricultural development is facing challenges in maintaining optimal yields while mitigating environmental impacts [13,14]. Therefore, addressing climate change and optimizing management measures for paddy fields are problems that should be urgently resolved.

The combination of process-based modeling and various experimental data provides opportunities for quantifying the impacts of different management practices and future climate change on soil C dynamics [15]. In fact, comprehensively and accurately evaluating SOC change is difficult due to the low speed of SOC dynamics and time-consuming and laborious on-site sampling; thus, a calibration model is necessary. Agronomists and scientists have worked diligently in the past to devise and promote the use of agricultural practices that can maintain or increase SOC levels. With the continuous development of agricultural research methods in addition to physical sampling and analysis of soil profiles for SOC, dynamic modelling of SOC can be used to effectively monitor soil organic carbon storage under different agricultural management. Among the relatively mature related models, including CENTURY, denitrification and decomposition (DNDC), NCSOIL, and RothC, the DNDC model is widely used due to its simple parameter inputs and accurate result simulation. The DNDC model can satisfactorily simulate SOC conversion in paddy fields and crop growth under climate change [4].

Kunshan is located in the Tai-Lake region in the middle and low reaches of the Yangtze River paddy soil region of China, which is a typical rice production area in the country [16]. Many recent studies have revealed that the paddy soils in this area have high SOC sequestration potential [17,18]. Hence, combining the experimental data from Kunshan with that from the DNDC model is feasible. Although we are encouraged by the tests of and improvements in DNDC for crop yields and environmental impact estimation in the past two decades, the widespread application of this tool in China has several limitations [19]. For instance, the constant 50-cm soil depth leads to overestimation of soil water content [20]. In addition, some soil properties, such as bulk density, porosity, and hydraulic parameters are assumed to be constant across all layers (down to a depth of 50 cm). However, most soil properties vary inherently between layers. Additionally, the traditional flood irrigation mode is the only irrigation mode for paddy fields, which makes it difficult to simulate the increasingly popularized water-saving irrigation mode. These factors may decrease the accuracy of irrigation simulation.

Although the DNDC model has been improved and applied in China through a two-decade effort, only four models exist for paddy fields under flood irrigation, namely, continuous flooding (the field water level is maintained at 10 cm), alternative wet and dry flooding (water level fluctuates between −5 to 5 cm), and rain-fed and empirical parameters. These four existing modes are inconsistent with the water-saving irrigation model in China. In a previous study [21], DNDC was used to simulate methane emissions from paddy fields under medium-term drainage, intermittent irrigation, and continuous flooding. In contrast to the above irrigation methods, under the condition of controlled

irrigation, which is widely applied in China, a shallow water layer is reserved on the surface of the field after transplanting seedlings until the regreening stage and the soil remains not flooded on the surface of the irrigation field in each subsequent growth stage, usually 60–80% of the time [22]. The irrigation time and irrigation amount were determined with the root-layer soil moisture as the control index. The existing model cannot simulate the controlled irrigation conditions. Thus, urgently modifying the DNDC model for paddy fields under water-saving irrigation is necessary to decrease site-specific suitability [23]. Given these problems, this research modified the 50-cm soil layer in the model to the depth of the root layer and controlled the upper and lower limits of paddy irrigation with soil moisture content. Additionally, the limits were modified in accordance with the needs of different growth stages of rice to adapt to the local water-saving irrigation mode. We hypothesized that crop growth and SOC dynamics could be simulated by improving the soil moisture module of this model. On this basis, the effects of different water and carbon management on SOC and rice yield in future climate conditions were studied.

Interest is growing in terms of finding ways to simulate climate change by using General Circulation Models (GCMs), which is the main current approach to predict future climate change and its responses. Substantial progress in global and regional modeling at medium to high resolutions and in downscaling methods has provided the basis for an increasing number of studies that attempt to simulate the effect of future climate change [24]. Predicting the dynamic changes in SOC and rice yield in paddy fields in the future is important for formulating agricultural management measures to save water, to increase yield, and to promote sustainable development. We carried out this study on the basis of the modified DNDC model and four climate scenarios under four GCMs weighted by Bayesian model averaging (BMA). The objectives of the study are (1) to validate the relevant parameters and to simulate changes in SOC and rice yield in Kunshan for the next 80 years, and (2) to extend the paddy field irrigation module in the DNDC model to provide a theoretical basis for optimizing field management measures to cope with climate change.

#### **2. Materials and Methods**

#### *2.1. Experimental Site*

The experiment was conducted in 2015 and 2016 at the State Key Laboratory of Hydrology-Water Resources and Hydraulic Engineering of Hohai University, Kunshan Irrigation and Drainage Experiment Station (31◦15015" N, 120◦57043" E), Jiangsu Province, China (Figure 1). The study area has a subtropical monsoon climate with a mean annual precipitation of 1097 mm, an average annual air temperature of 15.5 ◦C, a sunshine duration of 2086 h, and a frost-free period of 234 days·y −1 . The locals are accustomed to a rotation of rice and wheat planting. The paddy soil is classified as a hydragric anthrosol, which has a heavy loam texture, with a bulk density of 1.32 g cm−<sup>3</sup> at 0–30 cm and an initial pH of 7.4 at 0–18 cm. The organic matter is 21.71 g kg−<sup>1</sup> for the top 0–18 cm layer, and total K, total P, and total N are 20.86, 1.40, and 1.79 g kg−<sup>1</sup> for the 0–30 cm layer, respectively.

#### *2.2. Field Management*

The experiment was laid out (plot size 150 m<sup>2</sup> ) in a randomized block design with six treatments in triplicate. The six treatments were a combination of irrigation and fertilizer management systems: the two irrigation managements regimes were flood irrigation (FI) and controlled irrigation (CI), and the three fertilizer managements were wheat straw returning (S), organic fertilizer management (O), and farmer fertilizer practices (FFP). The six treatments were FS (FI and S), FO (FI and O), FF (FI and FFP), CS (CI and S), CO (CI and O), and CF (CI and FFP), with a total of 18 cells. Rain-fed wheat was grown in the plots during the non-rice planting season.

*Sustainability* **2021**, *13*, x FOR PEER REVIEW 4 of 26

**Figure 1.** Location of the experimental station. **Figure 1.** Location of the experimental station.

*2.2. Field Management*  The experiment was laid out (plot size 150 m2) in a randomized block design with six treatments in triplicate. The six treatments were a combination of irrigation and fertilizer management systems: the two irrigation managements regimes were flood irrigation (FI) and controlled irrigation (CI), and the three fertilizer managements were wheat straw returning (S), organic fertilizer management (O), and farmer fertilizer practices (FFP). The six treatments were FS (FI and S), FO (FI and O), FF (FI and FFP), CS (CI and S), CO (CI and O), and CF (CI and FFP), with a total of 18 cells. Rain-fed wheat was grown in the plots during the non-rice planting season. The rice variety in the experiment was Japonica Rice Nanging 46. Three or four seedlings per hill were transplanted in late June, with a plant spacing of 13.0 cm × 25.0 cm, and The rice variety in the experiment was Japonica Rice Nanging 46. Three or four seedlings per hill were transplanted in late June, with a plant spacing of 13.0 cm × 25.0 cm, and were harvested in late October. Local nitrogen (N) fertilizer was adopted in FFP (Table 1). The chemical fertilizer management of the S treatment was similar to that of the FFP treatment, and 3000 kg ha−<sup>1</sup> of straw from the previous wheat crop (the organic carbon content of wheat straw was 441 g kg−<sup>1</sup> , while the C/N ratio was 50:1 and the organic carbon input through wheat straw was 1322 g kg−<sup>1</sup> ) was returned to the S paddy fields both years. Additionally, 7500 kg ha−<sup>1</sup> of well-decomposed chicken manure (23% moisture content, 16.3 g kg−<sup>1</sup> N, 261 g kg−<sup>1</sup> organic carbon, 15.4 g kg−<sup>1</sup> P2O5, and 20.7 g kg−<sup>1</sup> K2O (Shijiazhuang Jitian Biotechnology Co., Ltd., China) was applied to the O paddy fields in 2015 and 2016. The base fertilizer and wheat straw were mixed into the muddy soil during tillage, and surface application was adopted for all other fertilizers.

were harvested in late October. Local nitrogen (N) fertilizer was adopted in FFP (Table 1). The chemical fertilizer management of the S treatment was similar to that of the FFP treatment, and 3000 kg ha−1 of straw from the previous wheat crop (the organic carbon content **Table 1.** Date and rate of nitrogen fertilization during the rice-growing season in farmer fertilizer practice (FFP) (kg N ha−<sup>1</sup> ).


The base fertilizer and wheat straw were mixed into the muddy soil during tillage, and surface application was adopted for all other fertilizers. Dates in brackets are when the fertilizer was applied in 2015 and 2016, respectively. CF: compound fertilizer (N, P2O5, and K2O contents were 16.0%, 12.0%, and 17.0% in 2015 and 2016), AB: ammonium bicarbonate (N content was 17.1%), U: urea (N content was 46.2%).

**Table 1.** Date and rate of nitrogen fertilization during the rice-growing season in farmer fertilizer practice (FFP) (kg N ha<sup>−</sup>1). **Activity 2015 2016**  Base fertilizer (29 and 28 June) 155.2 (72.0CF+83.2AB) 72.0 (72.0CF) Tillering fertilizer (16 Jul) 69.3 (U) 97.0 (U) Panicle fertilizer (9 and 11 Aug) 58.9 (U) 104.0 (U) Total nitrogen 283.4 273.0 Dates in brackets are when the fertilizer was applied in 2015 and 2016, respectively. CF: compound fertilizer (N, P2O5, and K2O contents were 16.0%, 12.0%, and 17.0% in 2015 and 2016), AB: The irrigation water layer of the CI paddy fields was maintained at 5–25 mm in the regreening stage. Irrigation was applied only to keep the soil moist, and standing water was avoided in the other stages except during periods of pesticide and fertilizer application. In accordance with local rice planting habits, a 30–50 mm shallow water layer was retained in the FI paddy fields after transplantation except during the midseason drainage period and the yellow maturity stage of rice. Rainfall was deflected with a canopy to accurately control soil moisture. The root zone soil water content criteria in different rice growth stages for CI are shown in Table 2.

The irrigation water layer of the CI paddy fields was maintained at 5–25 mm in the regreening stage. Irrigation was applied only to keep the soil moist, and standing water

ammonium bicarbonate (N content was 17.1%), U: urea (N content was 46.2%).


**Table 2.** Limits for irrigation in different stages of rice under controlled irrigation.

*θ*s1, *θ*s2, and *θ*s3 represent saturated volumetric soil moisture for the 0–20, 0–30, and 0–40 cm layers, respectively. <sup>a</sup> In the case of pesticide, fertilizer application, and rainfall, standing irrigation water at a depth of up to 5 cm was maintained for less than 5 days. <sup>b</sup> The data show the water depth during the re-greening stage.

#### *2.3. Yield Measurement, Soil Sampling, and Analysis*

Rice yield was estimated by artificially harvesting crops per unit area of each plot. Three hills of rice were randomly chosen to evaluate the filled grain number, setting percentage, thousand kernel weight, and panicle number of each treatment.

A total of 108 soil samples were collected from each plot following an S-shaped pattern at 0–10, 10–20, and 20–40 cm depths during the whole growth stage of rice in 2015 (23 June, 12 July, 20 August, 23 August, 21 September, and 25 October) and 2016 (29 June, 27 July, 21 August, 4 September, 21 September, and 25 October). After harvesting with a spiral drill (diameter, 38 mm; length, 1 m), three samples of 0–40 cm soil were randomly collected in each plot and fully stirred. Then, samples from the same depth were homogenized by mixing, separated from visible debris and crop residues, and divided into two parts. One part of the fresh samples was stored at 4 ◦C, and the other was air-dried, ground, and screened with a sieve of 0.149 mm; 12.5 g of fresh soil samples for dissolved organic carbon (DOC) was placed in a conical flask, immersed into 50 mL of 0.5 mol L−<sup>1</sup> K2SO<sup>4</sup> solution, and shaken for 30 min before the extracts were separated with a 0.45-µm filter. SOC was measured by the potassium dichromate external heating method, and the oxidation correction coefficient was considered. Besides, soil water content was recorded by a time domain reflectometer (Soil Moisture Equipment, Ltd., Corp. USA), and vertical rulers were used to monitor water layer at 8 a.m. everyday. The amount of irrigation water for each plot was calculated by using the water meter.

#### *2.4. DNDC Model*

#### 2.4.1. Overview of the DNDC Model

The DNDC model is a process-cased biogeochemical model written in Visual C++ 6.0 language for C and N dynamics in agro-ecosystems. This model has evolved over decades of development since it was developed by Li et al. [25]. Various soil hydrological processes were included in the present model. The DNDC model has been used worldwide because of its simple parameter input and accurate simulation results. It was designated as the preferred biogeochemical model in Asia by the International Symposium on Global Change in the Asia-Pacific region in 2000 [26]. The DNDC model has good adaptability in China [27,28], but studies on predicting SOC dynamics in paddy soil under water-saving irrigation and water-carbon coupling based on future climatic conditions are few. Therefore, the present study improved the irrigation module of DNDC95, which is the latest version of the DNDC model, to realize simulation of paddy fields under controlled irrigation and to optimize the irrigation module of paddy fields in the model on the basis of experimental data. More detail can be found in the Supplementary Materials.

#### 2.4.2. Input Data

Daily meteorological data, soil properties, and agricultural management measures were collected to support DNDC simulation. Soil physical and chemical properties, including initial soil C and N content, texture, and field capacity, were obtained through field sampling and laboratory analysis. The value of SOC at surface soil (0–10 cm) used as an input to the model was based on a measured total SOC value (11.1 g C kg−<sup>1</sup> ). The

1

contents of TN, NH<sup>4</sup> + -N, and NO<sup>3</sup> −-N served as a pre-fertilization input value of DNDC. Agricultural management measures were obtained on the basis of field records and local farmers' habits. The meteorological data used in this paper were as follows: historical meteorological observation data and GCMs from the Meteorological Information Center of China Meteorological Administration (http://data.cma.gov.cn/). Data included the daily maximum temperature, minimum temperature, radiation, wind speed, and precipitation. The future climate projections were acquired from four GCMs participating in the Coupled Model Intercomparison Project (CMIP5) experiment, including BCC-CSM1.1 (m), MIROC-ESM-CHEM, GFDL-ESM2M, and HadGEM2-ES (Table 3) [24,29]. In accordance with the new emissions scenarios proposed by CMIP5, representative concentration pathways (RCPs) and four climate scenarios, RCP 2.6, RCP 4.5, RCP 6.0, and RCP 8.5, were selected. RCP 2.6 is a low peak-and-decay scenario (the radiation force reaches its maximum near the middle of the 21st century before falling to 2.6 W m−<sup>2</sup> ), RCP 8.5 is a high-emissions scenario (the radiation force rises to 8.5 W m−<sup>2</sup> by 2100), and RCP 6.0 and RCP 4.5 are two intermediate scenarios (with a radiation force stability in 6.0 W m−<sup>2</sup> and 2.6 W m−<sup>2</sup> , respectively, by 2100).

**Table 3.** Four general climate models used in this study.


: Atmosphere and Ocean Research Institute (The University of Tokyo) and National Institute of Environmental Studies.

#### 2.4.3. BMA Method

As an advanced statistical method based on Bayesian theory and in consideration of model uncertainty, BMA has been proposed to combine multiple climate models to provide good performance models with high weights. BMA has been widely used in multimodel ensemble predictions of future climate. Therefore, four future climate models were predicted using the BMA weighted set in the present study and estimated on the basis of two statistical downscaling methods: back-propagation neural network and Statistical Downscaling Model (SDSM) developed by Wilby et al. [30,31]. Their mathematical expressions are as follows [24]:

Assume that *y* is the prediction variable, and its posterior distribution is as follows:

$$p(y|f\_{1'}f\_{2'},...,f\_{k'}D) = \sum\_{k=1}^{K} p(y|f\_{k'}D)p(f\_k|D) \tag{1}$$

On the premise of satisfying the minimum mean squared error, the combined prediction formula on the basis of the basic principle of Bayesian theorem is as follows:

$$E\_{BMA}(y|D) = \sum\_{k=1}^{K} p(f\_k|D) \operatorname{E}[p\_k(y|f\_{k'}D)] = \sum\_{k=1}^{K} \omega\_{k'} f\_k \tag{2}$$

where *p*(*fk*|*D*) denotes the posterior probability that model *f<sup>k</sup>* is correct given the training data and is calculated with Bayes' theory; *p*(*y*|*f<sup>k</sup>* , *D*), estimated from the training data, is the predictive probability density function based on model *y*|*f<sup>k</sup>* alone; and *k* is the number of models being combined, which is equal to four in this study. This formula uses the posterior probability *p*(*fk*|*D*) of the model as the weight for all possible model predictions *E*(*D*|*f<sup>k</sup>* , *D*) and obtains the combined predicted value.

Based on the field experimental data, we modified and verified the DNDC model to simulate soil organic carbon in paddy fields under different water and carbon management systems. The controlled irrigation module was added to the irrigation module of DNDC to realize the simulation of paddy fields under controlled irrigation. Then, combined with the climate model and climate scenarios after the BMA-weighted average, the simulation of SOC and rice yield under the corresponding water and carbon management systems in the next 80 years was conducted.

#### *2.5. Data Analysis*

Validation of the model results in the current study mainly included the average deviation method, correlation coefficient method, relative error method, and root mean squared method [32]. The absolute root mean squared error (*RMSEa*), normalized root mean squared error (*RMSEn*), coefficient of model efficiency (*EF*), and coefficient of determination (*R* 2 ) were used to quantitatively assess the goodness-of-fit between the simulated results and measured (observed) results. Their mathematical expressions are as follows:

$$EF = 1 - \frac{\sum\_{i=1}^{n} \left( SM\_i - OBS\_i \right)^2}{\sum\_{i=1}^{n} \left( OBS\_i - \overline{OBS} \right)^2} \tag{3}$$

$$EF = 1 - \frac{\sum\_{i=1}^{n} \left( SM\_i - OBS\_i \right)^2}{\sum\_{i=1}^{n} \left( OBS\_i - \overline{OBS} \right)^2} \tag{4}$$

$$RMSE\_n = \frac{100 \times RMSE\_a}{OBS\_{avg}} \tag{5}$$

$$R^2 = (\frac{\sum\_{i=1}^{n} \left(OBS\_i - OBS\_{avg}\right) \left(SM\_i - SM\_{avg}\right)}{\sqrt{\sum\_{i=1}^{n} \left(OBS\_i - OBS\_{avg}\right)^2 \sum\_{i=1}^{n} \left(SM\_i - SM\_{avg}\right)^2}}\tag{6}$$

where *OBS<sup>i</sup>* is the observed value, *OBSavg* is the average observed value, *SM<sup>i</sup>* is the simulated value, *SMavg* is the average simulated value, and *n* is the sample size. Higher *R* 2 and lower *RMSE<sup>n</sup>* indicated a good fit between the simulated and observed data. The smaller the *RMSE<sup>n</sup>* value is, the higher the fitting degree between the simulated value and the observed value. A value less than 10% indicates good consistency between the simulated value and the observed value. The results between 10% and 20% indicate an ordinary simulation effect, and a value higher than 30% indicates an unsatisfactory simulation effect [33,34]. The Taylor diagram is a polar-style graph, which summarizes the three statistical indices, i.e., the correlation coefficient between simulations and observations (*R*), the root mean squared error (*RMSE*), and the standard deviation (*STD*) using a single point. Given its comprehensiveness and visibility, Taylor diagrams are particularly beneficial in evaluating the relative accuracy of the different models. The radial distance from the origin reflects *STD*, the cosine of the azimuth angle denotes *R*, and the radial distance from the observed points is proportional to the *RMSE* difference. A main criterion can usually be summarized: the closer a point is to the observed data, the better the fit between the observed and simulated data [35].

Origin 9.1 software (OriginLab Corporation, Northampton, MA, USA) and MATLAB 2017 (MathWorks Corporation, USA) were used to calculate data and construct the relevant charts. Statistical analysis was carried out using standard procedures on a randomized plot design (SPSS 22.0). Significance was calculated on the basis of a Least significant difference (LSD) test at the 0.05 probability level.

The Mann–Kendall trend test, which we used in this study based on MATLAB 2017, is one of the widely used distribution-free tests of trend in time series. A standard normal variate *Z* is calculated as follows:

$$Z = \left\{ \begin{array}{c} \frac{S-1}{\sqrt{Var(S)}}, S > 0 \\ 0, S = 0 \\ \frac{S+1}{\sqrt{Var(S)}}, S < 0 \end{array} \right\} \tag{7}$$

$$\text{LIF}\_k = \frac{\text{S}\_k - E(\text{S}\_k)}{\sqrt{Var(\text{S}\_k)}}, k = 1, 2, \dots, n \tag{8}$$

$$
\Delta IB\_k = \left\{ \begin{array}{c} -\mathsf{U}F\_{k'}k = n\_\prime n - 1, \dots, 1 \\ 0, k = 1 \end{array} \right\} \tag{9}
$$

In a two-sided test for the trend, the null hypothesis of no trend is rejected if |*Z*| > *Zα*/2 where α is the significance. The calculation method of *Var*(*S*) and *S* can be found in the literature [36], where *Z* > 0 indicates an upward trend and *Z* < 0 indicates a downward trend. In addition, *UF* is the standardized result of *S*, which is a statistical sequence calculated in time sequence and obeys normal distribution, while *UB* is repeatedly calculated in reverse chronological order.

#### **3. Results**

#### *3.1. Model Modification and Validation*

#### 3.1.1. Model Modification

On the basis of the source code of DNDC95, this study improved the module on paddy field flooding in the farmland management menu. The two methods for the original water management module are the following: continuous flooding (water level is maintained at 10 cm) and alternative irrigation (water level fluctuates between −5 to 5 cm). The problems in the model were solved by improving the following three aspects: (1) the 50-cm constant soil layer assumed in the original DNDC model was adjusted to a value that varied with the depth of the rice root layer; (2) the fluctuation range of the water level was adjusted in accordance with the upper and lower limits of irrigation water controlled by soil moisture content; and (3) the upper and lower limits of irrigation with controlled irrigation were changed with the rice growth period, controlled irrigation with rice growth period was implemented, and the corresponding parameters were adjusted. Controlled irrigation was monitored in accordance with the soil moisture and water layer indicators in Table 2. The amount of irrigation water simulated by DNDC under controlled irrigation and traditional flooding irrigation after the modification was consistent with the observed irrigation water amount (Table 4). Additionally, crop parameters were calibrated in this study. The maximum crop yield, biomass allocation, and C/N ratio of the crops were modified on the basis of the observed results, and some internal parameters were modified to simulate actual conditions in the field. For example, the chromic acid wet oxidation method [37] and the Kjeldahl method [38] were used to estimate the total carbon nitrogen ratio of stems, leaves, and grains at the heading and maturing stages of Nanging 46. The total C/N ratios used for model correction were 55 for the root, 75 for the stem and leaf, and 75 for the grain. The maximum biomass production of grain was modified to 4700 kg C ha−<sup>1</sup> to stay consistent with our observed data.


**Table 4.** Comparison of observed and simulated irrigation values of the Denitrification Decomposition (DNDC) model simulation. sition (DNDC) model simulation. **Year Treatments Observed/mm Simulated/mm** *RMSEn*

**Table 4.** Comparison of observed and simulated irrigation values of the Denitrification Decompo-

Notes: Observed and simulated denote the observed irrigation amount and the simulated irrigation amount, respectively. 3.1.2. Model Calibration and Validation

#### 3.1.2. Model Calibration and Validation Model Calibration

#### Model Calibration The comparisons of DOC and SOC measured values and simulated values in the test

The comparisons of DOC and SOC measured values and simulated values in the test area in 2015 are shown in Figures 2 and 3. The dynamic changes in SOC and DOC in paddy soil under different water and carbon management systems in one year were well fitted through the modified DNDC model. The simulated values were consistent with the observed values. Tables 5 and 6 reflect the evaluation results of the SOC and DOC simulation values, respectively. The *RMSE<sup>a</sup>* values of the SOC and DOC simulations were 0.35–1.62 g kg−<sup>1</sup> and 23.63–38.49 mg kg−<sup>1</sup> , respectively. The *RMSE<sup>n</sup>* values of the SOC and DOC simulations were 3.54–17.59% and 8.79–13.93%, respectively. The regression coefficient *R* <sup>2</sup> of DOC was 0.80–0.99, and the EF values of SOC and DOC were close to 1. The SOC regression coefficients of the partial treatments (FS and FO) were closer to 1, which indicated that the modified DNDC model can accurately simulate the effects of different water and carbon management systems on SOC and DOC dynamics in paddy soil. area in 2015 are shown in Figures 2 and 3. The dynamic changes in SOC and DOC in paddy soil under different water and carbon management systems in one year were well fitted through the modified DNDC model. The simulated values were consistent with the observed values. Tables 5 and 6 reflect the evaluation results of the SOC and DOC simulation values, respectively. The *RMSEa* values of the SOC and DOC simulations were 0.35– 1.62 g kg−1 and 23.63–38.49 mg kg−1, respectively. The *RMSEn* values of the SOC and DOC simulations were 3.54–17.59% and 8.79–13.93%, respectively. The regression coefficient *R2* of DOC was 0.80–0.99, and the EF values of SOC and DOC were close to 1. The SOC regression coefficients of the partial treatments (FS and FO) were closer to 1, which indicated that the modified DNDC model can accurately simulate the effects of different water and carbon management systems on SOC and DOC dynamics in paddy soil.

**Figure 2.** Simulation of dissolved organic carbon (DOC) (0–10 cm soil) change in each treatment during the calibration period (2015), where (**a**–**f**) present the CS, FS, CO, FO, CF, and FF treatments, respectively. **Figure 2.** Simulation of dissolved organic carbon (DOC) (0–10 cm soil) change in each treatment during the calibration period (2015), where (**a**–**f**) present the CS, FS, CO, FO, CF, and FF treatments, respectively.

*Sustainability* **2021**, *13*, x FOR PEER REVIEW 10 of 26

**Figure 3.** Simulation of soil organic carbon (SOC) change in each treatment during the calibration period (2015), where (**a**–**f**) present the CS, FS, CO, FO, CF, and FF treatments, respectively. **Figure 3.** Simulation of soil organic carbon (SOC) change in each treatment during the calibration period (2015), where (**a**–**f**) present the CS, FS, CO, FO, CF, and FF treatments, respectively.

**Table 5.** Estimation of SOC results for each treatment by using the modified DNDC model during the calibration period (units of SOC: g kg<sup>−</sup>1). **Table 5.** Estimation of SOC results for each treatment by using the modified DNDC model during the calibration period (units of SOC: g kg−<sup>1</sup> ).


 CO 6 7.60(0.52) 7.27(0.04) 0.24 \* −0.06 7.67 0.59 0.65 8.88 0.99 FF 6 7.68(0.68) 7.78(0.03) 0.76 \* −0.01 7.80 0.81 0.70 8.95 0.99 FS 6 10.12(1.31) 9.20(0.03) 0.19 \* −0.02 9.41 0.62 1.62 17.59 0.97 FO 6 9.32(0.31) 9.20(0.03) 0.46 \* −0.05 8.39 0.79 0.35 3.82 1.00 Notes: *N* is the number of samples; *X*obs is the average observed value; *X*sim is the average simulated value; SD is standard deviation; *P*(*t\**) is *t*-test significance; *α* and *β* are the slope and intercept of the linear correlation between simulated values and observed values, respectively; and *R* 2 is the coefficient of determination between the simulated value and the observed value. In *P*(*t\**), \* means that the difference between the simulated value and the observed value is not significant and that the credibility is 95%.

> Notes: *N* is the number of samples; *Xobs* is the average observed value; *Xsim* is the average simulated value; *SD* is standard deviation; *P*(*t\**) is t-test significance; *α* and *β* are the slope and intercept of the linear correlation between simulated values and observed values, respectively; and *R*2 is the coefficient of determination between the simulated value and the observed value. In *P*(*t\**), \* means that the difference between the simulated value and the observed value is not significant and that

the credibility is 95%.


**Table 6.** Evaluation of DOC simulation results of each treatment by using a modified DNDC model during the calibration period and verification period (units of DOC: mg kg−<sup>1</sup> ).

Notes: *N* is the number of samples; *X*obs is the average observed value; *X*sim is the average simulated value; SD is standard deviation; *P*(*t\**) is *t*-test significance; *α* and *β* are the slope and intercept of the linear correlation between simulated values and observed values, respectively; and *R* 2 is the coefficient of determination between the simulated value and the observed value. In *P*(*t\**), \* means that the difference between the simulated value and the observed value is not significant and that the credibility is 95%.

#### Validation of Model Parameters

This study validated the modified DNDC model with 2016 data. The comparison between the simulated and observed values of DOC and SOC with different treatments during the verification period is shown in Figures 4 and 5. In most cases, the modified DNDC model with calibration parameters can simulate the dynamics of DOC and SOC in paddy fields under different water and carbon management systems. On the time scale of one year, DOC in paddy fields clearly changed with time, showing an increasing first and then decreasing trend, whereas the SOC content had a negligible change. In addition, the vertical distribution of SOC in paddy fields under different water and carbon management systems was relatively consistent. The SOC in the paddy field decreased as the soil depth increased, and the SOC fluctuation of 0–10 cm was larger than the SOC fluctuations of 10–20 cm and 20–40 cm. These results were essentially consistent with those of previous studies [39]. The results (Figure 6) showed that the simulated values of rice yield under different water and carbon treatments in the calibration and verification periods were close to the observed data, that is, to the line 1:1.

#### Comparison of Observed and Simulated Values

The parameter evaluation results for DOC (Table 6) and SOC (Table 7) in paddy fields simulated by the modified DNDC model showed the relationship between the simulated and observed values. *RMSE<sup>a</sup>* and *RMSE<sup>n</sup>* were small, indicating that the simulation was good. The model verification results indicated that irrigation and fertilization management had a great impact on SOC and DOC in paddy fields. Irrigation affected the dynamics of SOC and DOC. SOC under controlled irrigation was lower than that under flooding irrigation, but DOC was higher. Controlled irrigation is beneficial to the oxidative decomposition of paddy soil, which may be the cause of this phenomenon. In addition, the SOC contents of the organic fertilizer and straw returning treatments were significantly higher than the SOC content of the conventional fertilizer treatment, indicating that the appropriate fertilization method was beneficial to SOC accumulation in paddy fields.

**DOC (mg kg-1)**

**DOC (mg kg-1)**

**DOC (mg kg-1)**

*Sustainability* **2021**, *13*, x FOR PEER REVIEW 12 of 26

*Sustainability* **2021**, *13*, x FOR PEER REVIEW 12 of 26

**Figure 4.** Simulation of DOC (0–10 cm soil) dynamics in each treatment during the verification period (2016), where (**a**–**f**) present the CS, FS, CO, FO, CF, and FF treatments, respectively. **Figure 4.** Simulation of DOC (0–10 cm soil) dynamics in each treatment during the verification period (2016), where (**a**–**f**) present the CS, FS, CO, FO, CF, and FF treatments, respectively. **Figure 4.** Simulation of DOC (0–10 cm soil) dynamics in each treatment during the verification period (2016), where (**a**–**f**) present the CS, FS, CO, FO, CF, and FF treatments, respectively.

**Figure 5.** Simulation of SOC changes in each treatment during the verification period (2016), where (**a**–**f**) present the CS, FS, CO, FO, CF, and FF treatments, respectively. **Figure 5.** Simulation of SOC changes in each treatment during the verification period (2016), where (**a**–**f**) present the CS, FS, CO, FO, CF, and FF treatments, respectively. **Figure 5.** Simulation of SOC changes in each treatment during the verification period (2016), where (**a**–**f**) present the CS, FS, CO, FO, CF, and FF treatments, respectively.

*Sustainability* **2021**, *13*, x FOR PEER REVIEW 13 of 26

**Figure 6.** Simulation of yield changes in each treatment during the validation period (2015) and calibration period (2016): the solid line is a 1:1 relationship. **Figure 6.** Simulation of yield changes in each treatment during the validation period (2015) and calibration period (2016): the solid line is a 1:1 relationship.

Comparison of Observed and Simulated Values The parameter evaluation results for DOC (Table 6) and SOC (Table 7) in paddy fields **Table 7.** Evaluation of SOC simulation results of each treatment by using modified DNDC model during the verification period (units of SOC: g kg−<sup>1</sup> ).


0–10 cm CS 6 11.02(0.59) 11.58(0.54) 0.05 \* 0.57 10.00 0.76 0.75 6.46 0.84 CO 6 13.69(0.78) 12.35(0.38) 0.01 0.35 10.83 0.62 1.45 11.77 0.83 FF 6 10.88(0.15) 11.11(0.25) 0.09 \* 0.52 5.58 0.85 0.34 3.09 0.92 FS 6 12.96(0.81) 12.35(0.24) 0.16 \* 0.82 13.37 0.84 1.02 8.26 0.56 Notes: *N* is the number of samples; *X*obs is the average observed value; *X*sim is the average simulated value; SD is standard deviation; *P*(*t\**) is *t*-test significance; *α* and *β* are the slope and intercept of the linear correlation between simulated values and observed values, respectively; and *R* 2 is the coefficient of determination between the simulated value and the observed value. In *P*(*t\**), \* means that the difference between the simulated value and the observed value is not significant and that the credibility is 95%.

FO 6 13.18(0.60) 12.56(0.19) 0.10 \* 0.09 11.38 0.87 0.93 7.37 0.87

SOC CF 6 8.76(0.08) 9.01(0.05) 0.01 −0.25 11.75 0.80 0.28 3.10 0.89 10–20 cm CS 6 9.28(0.45) 9.61(0.39) 0.01 0.82 2.00 0.87 0.37 3.83 0.93 CO 6 10.56(1.23) 10.41(0.05) 0.79 \* 0.02 10.18 0.82 1.21 11.67 0.83 *3.2. Projection of SOC and Rice Yield in Paddy Fields Based on BMA and Modified DNDC* 3.2.1. BMA Method Evaluation of Predicted Values of Meteorological Parameters Required by DNDC

 FF 6 10.10(0.50) 10.24(0.37) 0.69 \* −0.29 17.17 0.73 0.75 7.30 0.83 FS 6 11.45(0.33) 11.13(0.03) 0.09 \* −0.02 11.41 0.87 0.46 4.16 0.93 FO 6 10.78(0.36) 11.22(0.03) 0.04 0.01 11.21 0.81 0.57 5.09 0.87 SOC CF 6 7.48(0.34) 7.66(0.10) 0.33 \* −0.08 8.40 0.73 0.42 5.43 1.00 20–40 cm CS 6 8.04(0.52) 7.91(0.02) 0.61 \* 0.02 7.78 0.84 0.53 6.66 1.00 CO 6 8.83(0.63) 7.92(0.02) 0.02 0.02 7.76 0.67 1.10 13.83 0.98 FF 6 8.55(0.77) 8.70(0.62) 0.25 \* 0.77 2.25 0.91 0.30 3.42 −1.16 FS 6 9.70(0.07) 9.42(0.03) 0.00 −0.12 10.83 0.99 0.29 3.03 1.00 FO 6 9.08(4.07) 9.14(0.03) 0.80 \* −0.03 9.40 0.84 0.49 5.37 1.00 Notes: *N* is the number of samples; *Xobs* is the average observed value; *Xsim* is the average simulated value; *SD* is standard deviation; *P*(*t\**) is *t*-test significance; *α* and *β* are the slope and intercept of the linear correlation between simulated values and observed values, respectively; and *R*2 is the coefficient of determination between the simulated value and the observed value. In *P*(*t\**), \* means that the difference between the simulated value and the observed value is not significant and that Different GCMs should be combined to provide detailed and accurate climate data in the context of climate change. In the present study, four GCMs processed by BMA were used to obtain four climate variables as required by the modified DNDC model: maximum temperature, minimum temperature, wind speed, and radiation (Figure 7). The performance of the BMA ensemble multi-model to predict future climate variations was evaluated with a Taylor chart (Figure 8). Numerous studies have shown that the prediction effect of BMA parameters is improved by extending the model training time [40,41]. This study used 40 years (1961–2000) to train BMA weights, and current and future climate parameters were generated in the remaining stages (2001–2099). The comparison between simulated and observed precipitation values in 2015 and 2016 treated by BMA is shown in Figure 9. In the calibration and verification period of the model, the simulated and the observed rainfall values treated by BMA had a good fitting effect. The simulated precipitation value and the observed value were relatively close except for the slightest

the credibility is 95%.

quired by DNDC

occurrence of a peak value. In Figure 7, the meteorological parameters generated by BMA were more consistent on the daily scale than at other scales measured by any single model. Figure 8 shows the relative accuracy of the model with a Taylor diagram. The results of the BMA method (point E) were closer to the points marked "observed" than to the data measured by any single model (points A, B, C, and D). Thus, BMA exhibited a good correlation and small RMSE. Except for the analog value matching the effect of wind speed, which was slightly poor (even if R of the BMA method was also approximately 0.7), the prediction of the other meteorological factors was good. rence of a peak value. In Figure 7, the meteorological parameters generated by BMA were more consistent on the daily scale than at other scales measured by any single model. Figure 8 shows the relative accuracy of the model with a Taylor diagram. The results of the BMA method (point E) were closer to the points marked "observed" than to the data measured by any single model (points A, B, C, and D). Thus, BMA exhibited a good correlation and small RMSE. Except for the analog value matching the effect of wind speed, which was slightly poor (even if R of the BMA method was also approximately 0.7), the prediction of the other meteorological factors was good.

*Sustainability* **2021**, *13*, x FOR PEER REVIEW 14 of 26

*3.2. Projection of SOC and Rice Yield in Paddy Fields Based on BMA and Modified DNDC*  3.2.1. BMA Method Evaluation of Predicted Values of Meteorological Parameters Re-

Different GCMs should be combined to provide detailed and accurate climate data in the context of climate change. In the present study, four GCMs processed by BMA were used to obtain four climate variables as required by the modified DNDC model: maximum temperature, minimum temperature, wind speed, and radiation (Figure 7). The performance of the BMA ensemble multi-model to predict future climate variations was evaluated with a Taylor chart (Figure 8). Numerous studies have shown that the prediction effect of BMA parameters is improved by extending the model training time [40,41]. This study used 40 years (1961–2000) to train BMA weights, and current and future climate parameters were generated in the remaining stages (2001–2099). The comparison between

observed rainfall values treated by BMA had a good fitting effect. The simulated precipitation value and the observed value were relatively close except for the slightest occur-

**Figure 7.** Time series of daily mean maximum temperature (**a**), minimum temperature (**b**), wind speed (**c**), and radiation (**d**) from 2012 to 2016: observed is the measured value, and BCC-CSM1.1 (m), GFDL-ESM2M, HadGEM2-ES, and MIROC2SM-CHEM represent the four climate models in Table 2, respectively. BMA (Bayesian Model Averaging) represents the value after BMA-weighted average.

average.

average.

**Figure 7.** Time series of daily mean maximum temperature (**a**), minimum temperature (**b**), wind speed (**c**), and radiation (**d**) from 2012 to 2016: observed is the measured value, and BCC-CSM1.1

**Figure 7.** Time series of daily mean maximum temperature (**a**), minimum temperature (**b**), wind speed (**c**), and radiation (**d**) from 2012 to 2016: observed is the measured value, and BCC-CSM1.1

*Sustainability* **2021**, *13*, x FOR PEER REVIEW 15 of 26

**Figure 8.** Taylor diagrams for meteorological factors in Kunshan, 2012–2016: this diagram is a comparison between the projected and measured values of four meteorological parameters required by a modified DNDC model. The four figures are as follows: (**a**) maximum temperature, (**b**) minimum temperature, (**c**) wind speed, and (**d**) radiation. Observed is the observed value, A is BCC-CSM1.1 (m), B is GFDL-ESM2M, C is HadGEM2-ES, D is MIROC-3SM-CHEM, and E is the BMA-weighted value. **Figure 8.** Taylor diagrams for meteorological factors in Kunshan, 2012–2016: this diagram is a comparison between the projected and measured values of four meteorological parameters required by a modified DNDC model. The four figures are as follows: (**a**) maximum temperature, (**b**) minimum temperature, (**c**) wind speed, and (**d**) radiation. Observed is the observed value, A is BCC-CSM1.1 (m), B is GFDL-ESM2M, C is HadGEM2-ES, D is MIROC-3SM-CHEM, and E is the BMA-weighted value. **Figure 8.** Taylor diagrams for meteorological factors in Kunshan, 2012–2016: this diagram is a comparison between the projected and measured values of four meteorological parameters required by a modified DNDC model. The four figures are as follows: (**a**) maximum temperature, (**b**) minimum temperature, (**c**) wind speed, and (**d**) radiation. Observed is the observed value, A is BCC-CSM1.1 (m), B is GFDL-ESM2M, C is HadGEM2-ES, D is MIROC-3SM-CHEM, and E is the BMA-weighted value.

**Figure 9.** Comparison of simulated and actual precipitation values in 2015 (**a**) and 2016 (**b**) treated by BMA.

3.2.2. SOC Dynamics Prediction in Paddy Fields under Water and Carbon Regulation in Future Climate Conditions

On the basis of the modified DNDC model and the BMA method, this study predicted the SOC changes (0–10 cm) in paddy fields under four climate scenarios (i.e., RCP 2.6, RCP 4.5, RCP 6.0, and RCP 8.5) over the next 80 years (2020–2099), as shown in Figure 10. The average predicted SOC under different climate scenarios consistently occurred in the following order FO > CO > FS > CS > FF > CF. The trend lines of the SOC change in paddy fields under the four climate scenarios were estimated via linear square fitting (Figure 10). This trend indicated that the effect of fertilizer management on SOC changes in paddy fields over the long term was very large in the four scenarios. To some extent, this phenomenon explained the similar results found for the different climate scenarios, i.e., the SOC of the CF and FF treatments decreased with prolonged time, while the CS, CO, FS, and FO treatments showed an increasing trend with an extended time. Fertilizer management obviously affected the long-term trend of SOC in paddy fields under the same irrigation treatment. Irrigation had a certain impact on SOC in paddy fields over a short time, but only a negligible difference was observed over the long term. The overall trend in the SOC changes in paddy fields under flooding irrigation and controlled irrigation treatments was consistent and showed that SOC decreased in the conventional fertilizer treatment and increased in the treatment with organic fertilizer and straw application. In comparison with that in the 2020s, in the 2090s, the average values of the CF and FF treatments decreased by 4.98%, 5.86%, 6.07%, and 7.49% in the RCP 2.6, RCP 4.5, RCP 6.0, and RCP 8.5 scenarios, respectively, while the average values of the other treatments in the 2090s increased by 102.97%, 99.68%, 99.57%, and 97.54%, respectively. In addition, in the first 5 years, the CS and CO treatments showed an unexpected downward trend and then increased rapidly, which was different from the results of the model verification period. This may have been due to the frequent alternation of drying and wetting under controlled irrigation conditions, which promoted soil respiration. Therefore, the SOC of paddy fields decreased in the short term, while the long-term application of organic fertilizer and straw application can offset this carbon loss effect. However, the SOC of the organic fertilizer treatment under the RCP 4.5 and RCP 6.0 scenarios increased in 2100, which were because both the low peak attenuation and high emissions scenarios were not conducive to the *Sustainability* **2021**, *13*, x FOR PEER REVIEW accumulation of SOC in paddy fields. 17 of 26

**Figure 10.** Prediction of SOC change in paddy fields with different treatments in the next 80 years under different climate scenarios (0–10 cm): the dashed lines in different colors in the figure correspond to the corresponding trend lines, and each trend line was derived from a series of annual values. The annual SOC is the final content at the end of the growth period of each treatment in the next 80 years. **Figure 10.** Prediction of SOC change in paddy fields with different treatments in the next 80 years under different climate scenarios (0–10 cm): the dashed lines in different colors in the figure correspond to the corresponding trend lines, and each trend line was derived from a series of annual values. The annual SOC is the final content at the end of the growth period of each treatment in the next 80 years.

SOC of the conventional fertilizer treatment decreased rapidly in the first 10 years but gradually flattened. The soil organic carbon levels in the CF and FF treatments decreased by 14.18% and 13.50%, respectively. The SOC of the CS treatment abnormally decreased by 8.13% and increased rapidly. The effect of climate scenario on the SOC in paddy fields was not obvious (Figure 11). The SOC of the organic fertilizer treatment under the various climate scenarios increased with time. Compared with that under baseline conditions (2020), the SOC in the CO treatment under RCP 2.6 increased from 45.89% in 2040 to 149.39% in 2080 and the SOC in the CS treatment under RCP 4.5 increased from 3.07% in 2040 to 41.05% in 2080. In addition, the decline in the SOC in the CF and FF treatments

**Table 8.** Changes in the SOC of paddy fields with different treatment in the next 80 years under

**Period CF CO CS FF FO FS**  2020–2029 −14.18% 13.97% −8.13% −13.50% 18.86% 6.26% 2030–2039 −0.74% 24.33% 10.15% −3.18% 19.04% 8.74% 2040–2049 4.22% 19.66% 12.47% 5.17% 17.96% 8.36% 2050–2059 0.79% 13.70% 7.67% 0.41% 12.21% 5.36% 2060–2069 0.57% 10.85% 6.61% 0.21% 9.68% 4.66% 2070–2079 0.46% 8.88% 5.07% −0.04% 7.81% 3.67% 2080–2089 0.66% 6.99% 5.63% 2.04% 6.96% 4.40% 2090–2099 0.71% 6.07% 4.26% 0.11% 5.44% 3.02% Notes: The values above denote simulated SOC change every 10 years (compared with the baseline 10 years ago) of the CF, CO, CS, FF, FO, and FS treatments in the 2020s, 2030s, 2040s, 2050s,

was the largest in the first 20 years and remained unchanged.

the RCP 2.6 scenario.

2060s, 2070s, 2080s, and 2090s.

In Table 8, the dynamics of SOC every 10 years under different treatments in the next 80 years is reflected by the RCP 2.6 scenario as an example. The results showed that the SOC of the conventional fertilizer treatment decreased rapidly in the first 10 years but gradually flattened. The soil organic carbon levels in the CF and FF treatments decreased by 14.18% and 13.50%, respectively. The SOC of the CS treatment abnormally decreased by 8.13% and increased rapidly. The effect of climate scenario on the SOC in paddy fields was not obvious (Figure 11). The SOC of the organic fertilizer treatment under the various climate scenarios increased with time. Compared with that under baseline conditions (2020), the SOC in the CO treatment under RCP 2.6 increased from 45.89% in 2040 to 149.39% in 2080 and the SOC in the CS treatment under RCP 4.5 increased from 3.07% in 2040 to 41.05% in 2080. In addition, the decline in the SOC in the CF and FF treatments was the largest in the first 20 years and remained unchanged.

**Table 8.** Changes in the SOC of paddy fields with different treatment in the next 80 years under the RCP 2.6 scenario.


Notes: The values above denote simulated SOC change every 10 years (compared with the baseline 10 years ago) of the CF, CO, CS, FF, FO, and FS treatments in the 2020s, 2030s, 2040s, 2050s, 2060s, 2070s, 2080s, and 2090s.

#### 3.2.3. Projection of Rice Yield Changes

On the basis of the modified DNDC model and BMA method, we predicted rice yield changes under different water and carbon management systems over the next 80 years under four climate scenarios (Figure 12). The relationship between the different treatments was essentially the same under various climate scenarios, which showed that the rice yield of the CS treatment was the highest and that of the CF treatment was the lowest. Thus, the long-term return of straw can significantly promote an increase in rice yield. Similar to the regulation of water and carbon regulation of SOC dynamics in paddy fields, irrigation and carbon management both affected the yield under the same climate conditions while the combination of appropriate fertilization and controlled irrigation evidently increased rice yield. The rice yield in the CS and CO treatments in most cases was higher than that in the FS and FO treatments. This study provides a trend line of each rice yield with time (Figure 12). Overall, the rice yields of the different treatments have good synchronization and almost simultaneously changed at different stages of the 21st century. In comparison with that in the 2020s, the average rice yield of each treatment in the 2090s decreased by 18.41%, 38.59%, 65.11%, and 65.62% in RCP 2.6, RCP 4.5, RCP 6.0, and RCP 8.5, respectively. In addition, the climate scenarios resulted in clear effects on rice yields under the same water and carbon management mode. The rice yield tended to increase in the first 20 years as the radiative force increased. However, under the high emissions scenario of RCP 8.5, the rice yield of the CS treatment initially remained unchanged but declined rapidly with increased time. Taking RCP 2.6 as an example, the results of the Mann–Kendall trend test [42] are shown in Figure 13. The yields of the CF and FF treatments increased in 2020–2023 and 2087 (UF > 0), while the UF values of the CO, CS, FO, and FS treatments were less than zero within the 95% confidence interval, except for the increase in 2020–2023, which indicated that maintaining rice yield via excessive carbon input might be difficult to sustain.

*Sustainability* **2021**, *13*, x FOR PEER REVIEW 18 of 26

**Figure 11.** SOC in different treatments in four climate scenarios, where (**a**–**f**) present the CS, FS, CO, FO, CF, and FF treatments, respectively: the red, blue, and black lines represent the changes of SOC in paddy soil in 2040, 2060, and 2080, respectively, compared with the baseline (2020). The horizontal and vertical coordinates are the percentage values of the changes. **Figure 11.** SOC in different treatments in four climate scenarios, where (**a**–**f**) present the CS, FS, CO, FO, CF, and FF treatments, respectively: the red, blue, and black lines represent the changes of SOC in paddy soil in 2040, 2060, and 2080, respectively, compared with the baseline (2020). The horizontal and vertical coordinates are the percentage values of the changes.

On the basis of the modified DNDC model and BMA method, we predicted rice yield changes under different water and carbon management systems over the next 80 years under four climate scenarios (Figure 12). The relationship between the different treatments was essentially the same under various climate scenarios, which showed that the rice yield of the CS treatment was the highest and that of the CF treatment was the lowest. Thus, the long-term return of straw can significantly promote an increase in rice yield.

3.2.3. Projection of Rice Yield Changes

be difficult to sustain.

irrigation and carbon management both affected the yield under the same climate conditions while the combination of appropriate fertilization and controlled irrigation evidently increased rice yield. The rice yield in the CS and CO treatments in most cases was higher than that in the FS and FO treatments. This study provides a trend line of each rice yield with time (Figure 12). Overall, the rice yields of the different treatments have good synchronization and almost simultaneously changed at different stages of the 21st century. In comparison with that in the 2020s, the average rice yield of each treatment in the 2090s decreased by 18.41%, 38.59%, 65.11%, and 65.62% in RCP 2.6, RCP 4.5, RCP 6.0, and RCP 8.5, respectively. In addition, the climate scenarios resulted in clear effects on rice yields under the same water and carbon management mode. The rice yield tended to increase in the first 20 years as the radiative force increased. However, under the high emissions scenario of RCP 8.5, the rice yield of the CS treatment initially remained unchanged but declined rapidly with increased time. Taking RCP 2.6 as an example, the results of the Mann– Kendall trend test [42] are shown in Figure 13. The yields of the CF and FF treatments increased in 2020–2023 and 2087 (UF > 0), while the UF values of the CO, CS, FO, and FS treatments were less than zero within the 95% confidence interval, except for the increase in 2020–2023, which indicated that maintaining rice yield via excessive carbon input might

**Figure 12.** Prediction of rice yield change under different climate scenarios and treatments in the next 80 years: the trend lines in black, red, and blue in the figure represent conventional fertilizer treatment, organic fertilizer treatment, and straw returning treatment, respectively, while the solid and dotted lines represent conventional irrigation and controlled irrigation. **Figure 12.** Prediction of rice yield change under different climate scenarios and treatments in the next 80 years: the trend lines in black, red, and blue in the figure represent conventional fertilizer treatment, organic fertilizer treatment, and straw returning treatment, respectively, while the solid and dotted lines represent conventional irrigation and controlled irrigation.

*Sustainability* **2021**, *13*, x FOR PEER REVIEW 20 of 26

**Figure 13.** Mann–Kendall test charts of rice yield changes under different water and carbon treatments (with RCP 2.6 as an example), where (**a**–**f**) represent the CF, CO, CS, FF, FO, and FS treatments, respectively. The ordinate axis represents the values of UF and UB. UF > 0 indicates an upward trend, and UF < 0 indicates a downward trend. UB exceeded the upper and lower straight lines, indicating a significant upward or downward trend (*p* < 0.05). A sudden change point is **Figure 13.** Mann–Kendall test charts of rice yield changes under different water and carbon treatments (with RCP 2.6 as an example), where (**a**–**f**) represent the CF, CO, CS, FF, FO, and FS treatments, respectively. The ordinate axis represents the values of UF and UB. UF > 0 indicates an upward trend, and UF < 0 indicates a downward trend. UB exceeded the upper and lower straight lines, indicating a significant upward or downward trend (*p* < 0.05). A sudden change point is indicated when UF intersects UB and is between the upper and lower lines.

#### **4. Discussion**

**4. Discussion** 

indicated when UF intersects UB and is between the upper and lower lines.

#### *4.1. Performance of the Modified DNDC Model and Limitations 4.1. Performance of the Modified DNDC Model and Limitations*

The default parameters of the DNDC model did not meet the needs of simulating SOC dynamic changes [23], and the model should be calibrated to reduce uncertainties in new systems or environments [20]. The results of this study showed that the modified DNDC model had good adaptability to SOC and yield simulation of paddy fields in the Kunshan area. The modified DNDC model successfully predicted the irrigation situation under water-saving irrigation and flood irrigation, and the effects of different irrigation and fertilization conditions on the SOC, DOC, and rice yield in paddy fields can be simulated. In addition, current research has mainly focused on water consumption and water use efficiency [43] and less on the effect of climate change on SOC in rice fields, and climate factors, such as temperature and precipitation, are important driving forces in SOC The default parameters of the DNDC model did not meet the needs of simulating SOC dynamic changes [23], and the model should be calibrated to reduce uncertainties in new systems or environments [20]. The results of this study showed that the modified DNDC model had good adaptability to SOC and yield simulation of paddy fields in the Kunshan area. The modified DNDC model successfully predicted the irrigation situation under water-saving irrigation and flood irrigation, and the effects of different irrigation and fertilization conditions on the SOC, DOC, and rice yield in paddy fields can be simulated. In addition, current research has mainly focused on water consumption and water use efficiency [43] and less on the effect of climate change on SOC in rice fields, and climate factors, such as temperature and precipitation, are important driving forces in SOC

change [44], which have a far-reaching impact on agricultural production [9]. In this study, SOC prediction and rice yield were based on the modified DNDC model, local irrigation, fertilization management, and four GCMs integrated with BMA. The results weighted by BMA were closer to the observed points than to any single model in the Taylor diagram; thus, integrating multiple climate models with BMA is reliable, which is consistent with the results of Wang et al. [24]. Interpretation based on the single model was one of the limitations of this study. The uncertainty could be reduced by the method of multi-model ensemble [45]. In addition, it is desirable to calibrate the model results with data from more sites and long-term series of observed data under different water and carbon management.

#### *4.2. Effects of Water and Carbon Management Systems on SOC in Paddy Fields and Rice Yield*

The present study found that the combination of irrigation and fertilization patterns can markedly increase SOC and rice yield, which was consistent with the findings of Kamoni et al. [46]. This result may be due to irrigation improving the availability of soil N, thereby increasing productivity. The mechanism of the effects of irrigation on organic carbon remains unclear. Some studies have found that irrigation affects SOC mineralization and transfer [47], while others found that waterlogging affects rice residue input and the decomposition rate of SOC under anaerobic conditions, thus affecting SOC accumulation [48]. For example, Kelliher et al. [49] found that irrigation reduced SOC by 61%, while Houlbrooke et al. [50] found that irrigation had little effect on SOC, which may be related to environmental conditions, soil development stages and types, irrigation water quality, and years. This study found that the SOC of controlled irrigation paddy fields was lower than that of fields with conventional irrigation, which may be due to the frequent dry–wet alternation of controlled irrigation promoting microbial activities, increased soil fertility, and soil respiration, thus increasing soil carbon loss [51]; this finding is different from the results of Zhao et al. [52]. Zhao et al. found that optimized irrigation and fertilization treatments increased SOC in the North China Plain, which may be related to the retention of residue in the experiment every year. In addition, the present study found that controlled irrigation reduced the SOC of paddy fields while reducing irrigation water; the SOC content evidently increased after the combination of irrigation with straw returning or application of organic fertilizer. Thus, applying organic fertilizer or straw returning under controlled irrigation conditions can reduce the water footprint while addressing SOC. Combining controlled irrigation with organic fertilizer and returning straw to the fields, which is a feasible alternative water and carbon management mode, saved a large amount of water resources and increased rice yield and SOC content.

The dynamics of SOC in paddy fields are the net result of organic matter input and output. Irrigation schedules and fertilization affect soil organic carbon in paddy fields by changing the input of energy or material [53]. SOC dynamics are difficult to measure in the short term. This process-based model is a good tool for predicting future trends. The results of long-term simulation of the SOC changes in paddy fields under different water and carbon management systems (Figure 10) showed that the combination of controlled irrigation and suitable organic fertilizer application is a satisfactory water and carbon regulation mode. SOC growth was rapid, and yield was maintained at a high level with prolonged time. In addition, fertilizer management has a considerable impact on the longterm evolution of SOC on farmlands, which was consistent with the results of previous studies. For example, Wan et al. and Wang et al. [40,54] found through model research that an SOC of 0–30 cm on farmland in China would decrease to 7.8–8.2 t ha−<sup>1</sup> in 2080 without fertilizer management but would increase markedly if organic fertilizer or straw was applied to the field. This study found a synergistic relationship between SOC content and rice yield, and rice yield was high in the treatments with high SOC content, such as the CS and CO treatments, which was similar to the conclusion of Qiu et al. [55].

#### *4.3. Effects of Climate Scenarios on SOC and Rice Yield in Paddy Fields and Possible Countermeasures*

Impacts in climate scenarios have a considerable impact on rice yield, but their effect on SOC is less than that of agricultural management measures, which may be because climate change affects the decomposition of SOC, while agricultural management measures affect the soil carbon input quantity [56]; excessive carbon input may mask the impact of SOC decomposition. Additionally, the change in SOC was negatively correlated with initial SOC concentrations [57], and a high carbon input and low initial SOC would increase the pool of soil carbon. Conversely, the conversion of excessive carbon input into soil may offset the carbon loss caused by soil respiration, which explains to some extent why climate scenario impacts have a negligible effect on SOC changes in paddy fields. Unlike the current conclusion that fertilization can maintain high rice yields over the long term, although excessive fertilization can maintain high rice yields in the short term under future climate conditions, rice yields may still decrease in the long term (Figure 12). This phenomenon is attributed to the decline in rice yield caused by high temperatures and water stress that may have exceeded the impacts of promotion by fertilizer. The SOC of the controlled irrigation treatment increased rapidly in the late period but decreased in 2025, 2040, and 2083 in all treatments. This result may have been caused by the impact of climate conditions in such years. The average rice yields in all the treatments after 80 years decreased under the RCP 2.6, RCP 4.5, RCP 6.0, and RCP 8.5 scenarios by 18.41%, 38.59%, 65.11%, and 65.62%, respectively, compared with that in the baseline treatment (2020). In addition, previous studies [58] found that a variety of improvements can offset the decline in rice yield caused by climate warming, which might be a possible strategy to address climate change in the future.

Overall, paddy fields play a significant role in mitigating climate change through carbon sequestration, but the impact of different climate scenarios on SOC changes in paddy fields is less obvious than that of water and carbon management measures. Yu et al. [56] found that maintaining existing farmland management measures can maintain China's paddy soil carbon sequestration potential over the next 20–40 years; however, this result depends on long-term continuation of the current excessive carbon input management, which is closely related to the current policy of vigorously promoting and subsidizing straw returning and organic fertilizer application in China [59]. In accordance with the report released by the agricultural sector, most crop residues were removed from farmlands before the 1980s and used as fuel and animal feed in rural areas. This trend was reversed by the government through a policy in the 1990s to encourage farmers to recycle crop straw as much as possible, and the policy achieved considerable results [60]. At the same time, farmers stopped using crop straw as fuel due to improvements in living standards, which have caused serious environmental pollution in the past [61]. In addition, unreasonable fertilization leads to soil degradation, water pollution, soil acidification, and serious agricultural nonpoint source pollution [62]. Thus, how to promote straw returning in many developing countries across the world and to reduce its pollution is the direction of further study.

In addition, notably, in the future climate model, although water and carbon management will increase production and carbon sequestration, whether it will increase GHGs still needs further study. For example, excessive carbon input may increase greenhouse gases, such as CO<sup>2</sup> and CH4, while SOC changes are sensitive to CO<sup>2</sup> concentrations. Thus, the benefits of carbon sequestration may be offset. The predicted results showed that the rice yields of all the treatments will decrease in the future, especially after the middle of the 21st century. Although the rice yield decreased under the coupling of controlled irrigation with straw returning and organic fertilizer, the rice yield was always higher than that in conventional fertilizer treatments. Thus, finding an appropriate amount of organic fertilizer or straw application to balance carbon sequestration is necessary to increase production and to reduce greenhouse gas emissions.

#### **5. Conclusions**

This study modified the DNDC model to adapt to the common water-saving irrigation mode in China, especially in the middle and lower reaches of the Yangtze River. The parameters related to SOC and rice yield were calibrated. In addition, the dynamics of SOC and rice yield in Kunshan over the next 80 years under different water and carbon management were predicted on the basis of the four climate scenarios synthesized via the BMA method. The results showed that the modified DNDC model had good adaptability to the simulation of SOC and rice yield under different water and carbon management. The *RMSE<sup>n</sup>* values of the SOC and DOC simulations were 3.45% to 17.59% and 8.79% to 13.93%, respectively. The *R* <sup>2</sup> of DOC was between 0.80 and 0.99, and the model efficiency coefficient *EF* values of SOC and DOC were all close to 1. In comparison to the single model, the BMA method can better simulate the changes in climate factors. Climate scenarios significantly affect rice yield, but their impact on SOC is less than agricultural management measures. Unfavorable climate will reduce yields in the future climate in spite of long-term over fertilization. Compared with traditional water and carbon management systems, the combination of controlled irrigation and organic fertilizer application or straw returning can obviously increase the SOC content and rice yield in the long-term simulation under the four climate scenarios, and the yield of the straw-returning treatment was higher. The SOC of controlled irrigation paddy fields was lower than that of conventional irrigation, but only a negligible difference was observed over the long term. Therefore, combining controlled irrigation and appropriate organic fertilizer can balance water conservation, can maintain SOC and a stable rice yield in paddy fields, and is the recommended water and carbon management system for paddy fields.

**Supplementary Materials:** The following are available online at https://www.mdpi.com/2071-105 0/13/2/568/s1, Figure S1: Structure of the DNDC model, Table S1: Input parameters required for regional simulation with DNDC.

**Author Contributions:** Conceptualization, Z.J. and S.Y.; methodology, Z.J.; software, Z.J.; validation, X.S., J.D. and X.C.; formal analysis, X.L.; investigation, X.S.; resources, J.D.; data curation, X.C.; writing—original draft preparation, Z.J.; writing—review and editing, Z.J. and S.Y.; visualization, Z.J.; supervision, J.X.; project administration, X.L.; funding acquisition, S.Y. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research was funded by the National Natural Science Foundation of China (no. 51879076 and 51579070) and by the Fundamental Research Funds for the Central Universities (no. 2019B67814).

**Institutional Review Board Statement:** Not applicable.

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** Data is contained within the article or supplementary material.

**Acknowledgments:** Thanks to Li Changsheng, one of the main founders of the DNDC95 model. Thanks to the National Meteorological Information Center, China Meteorological Administration for offering the meteorological data and the Working Group of the World Climate Research Program on Coupled Modeling, which is responsible for CMIP5.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **Abbreviations**



## **References**

