A Sensitivity Analysis of the SPACSYS Model

Abstract

A sensitivity analysis is critical for determining the relative importance of model parameters to their influence on the simulated outputs from a process-based model. In this study, a sensitivity analysis for the SPACSYS model, first published in Ecological Modelling (Wu, et al., 2007), was conducted with respect to changes in 61 input parameters and their influence on 27 output variables. Parameter sensitivity was conducted in a ‘one at a time’ manner and objectively assessed through a single statistical diagnostic (normalized root mean square deviation) which ranked parameters according to their influence of each output variable in turn. A winter wheat field experiment provided the case study data. Two sets of weather elements to represent different climatic conditions and four different soil types were specified, where results indicated little influence on these specifications for the identification of the most sensitive parameters. Soil conditions and management were found to affect the ranking of parameter sensitivities more strongly than weather conditions for the selected outputs. Parameters related to drainage were strongly influential for simulations of soil water dynamics, yield and biomass of wheat, runoff, and leaching from soil during individual and consecutive growing years. Wheat yield and biomass simulations were sensitive to the ‘ammonium immobilised fraction’ parameter that related to soil mineralization and immobilisation. Simulations of CO₂ release from the soil and soil nutrient pool changes were most sensitive to external nutrient inputs and the process of denitrification, mineralization, and decomposition. This study provides important evidence of which SPACSYS parameters require the most care in their specification. Moving forward, this evidence can help direct efficient sampling and lab analyses for increased accuracy of such parameters. Results provide a useful reference for model users on which parameters are most influential for different simulation goals, which in turn provides better informed decision making for farmers and government policy alike.

1. Introduction

Process-based models for agricultural systems provide a widely used and efficient tool for understanding the complex interactions between soil water, carbon (C), nitrogen (N), phosphorus (P), and plant growth [1,2], where both production under various environmental conditions [3] and nutrient cycling [4] can be simulated. To run a simulation successfully, models need to be informed with parameters that accurately quantify individual processes, together with local soil conditions, weather, and management practices. However, parameters that might vary with environmental conditions and have varying, possibly complex distributions themselves [5] are often difficult to precisely characterise, especially if they are costly to measure, requiring long-term monitoring [6]. In turn, this parameter variability, uncertainty, or quality directly affects the reliability of the model simulations [1,7]. In general, agricultural simulation models require a large amount of input data and parameters, many of which are difficult to source or collect. At the same time, it is common for only a few parameters to strongly influence (and maximize) the variability of model outputs, while the majority of parameters only provide a weak influence in this respect [8]. Therefore, understanding the likely influence of each parameter on a model’s outputs is crucial, as this allows the number of parameters to be safely reduced without a worrying loss of model accuracy, while at the same time, focus can be placed on data collection for parameters that most strongly influence model performance.

To identify how parameters impact model outputs, sensitivity diagnostics can be calculated through running multiple model simulations while varying the parameters across specific ranges [9]. Sensitivity analysis for model parameterization and outputs have been applied to various cropping systems, soil types and climate conditions using a range of models. For example, Specka et al. [10] applied a parameter sensitivity analysis for modelling above-ground biomass with regard to a future model calibration and an improved understanding of model response patterns; Jabloun et al. [11] assessed the sensitivity of outputs (crop yield and N leaching) from a crop simulation model, where 128 parameters were assessed. A sensitivity analysis quantifies those parameters which are most influential on model outputs, and in this respect can guide the efforts towards improving their accuracy as well as the model’s output accuracy [12].

The SPACSYS model is implemented in a modular approach and provides a field-scale and weather-driven dynamic process-based simulation model of water, C and N cycling between plants, soils, and microbes [13]. The model has been widely used to assess the impact of climate change [14,15,16], tillage [17], fertilizer application [16,18], and different cultivars [19] on agricultural systems in terms of crop yields, C and N budgets, soil physical properties, and soil water redistribution. However, the model requires over 200 parameters for the simulation of various processes and plant growth and development. It would therefore be informative and beneficial to assess which parameters the SPACSYS outputs are the most and the least sensitive to, so that there can be guidelines for its use when all 200 parameters are not readily available.

A sensitivity analysis has also been used in building and understanding the structure of agro-ecosystem models [20,21]. Distinctions are made between a local and global sensitivity analysis (LSA and GSA, respectively), where the former is also known as a ‘one at a time’ (OAT) approach, while the latter considers multiple parameters at the same time. Thus, an LSA does not consider any interactions between parameters, while a GSA does. Both GSA [22] and LSA [23,24,25,26] have merit, where in ideal situations both should be applied and their results objectively assessed and compared [27,28]. In studies where uncertainty and interactions in model parameters are minimal, then an LSA may suffice [27] over the inherently more complex, GSA. Specifically, Link et al. [27] concludes that an LSA may suffice for situations where the ranking of model parameters is of importance, while GSA should be used for a precise attribution of variance in model outputs. Of course, it should be noted that results and insights from literature are always dependent on the characteristics of the process-based models employed.

Given this study is the first to identify parameter input sensitivity for the SPACSYS model, we chose to present results in two stages: first, an LSA as presented here, and second, a GSA (which is prep.) for presentation elsewhere. This two-stage reporting approach is followed for computational reasons, together with differences in associated interpretations, visualisations, and comparisons of the LSA and GSA outputs. Given SPACSYS requires over 200 parameters, assessing and reporting the sensitivity of them all would be problematic, for both an LSA and a GSA alike. Further, different subsets of parameters would be chosen according to their likely importance in simulating different processes under different background conditions. Because of computational and associated interpretation constraints, an LSA can investigate the sensitivity of outputs to a greater number of input parameters and investigate the distribution of each input parameter more intensively, than that feasible with a GSA. A GSA for, say, only four parameters entails the added investigation of six 2-parameter interactions, four 3-parameter interactions, and one 4-parameter interaction, which themselves vary on how the parameter input distributions are described (likely done simply, with say low, medium, and high values only). Given these thorny design and implementation issues for a GSA, for this study, we only consider the effects of a single parameter change on outputs via an LSA. This allows for a greater number of parameters to be investigated coupled with more detailed descriptions of the parameter distributions than would be viable within a GSA.

For this study, we carried out an LSA on the parameters that control soil water redistribution, and C and N cycling with the SPACSYS model, using a winter wheat field experiment conducted at Rothamsted Research in Harpenden, UK as the case study. To investigate potential influences of climatic and soil conditions on the analysis, simulations with the same configurations were run with different climatic and soil conditions. The study objectives were: (1) to identify those parameters that have the maximum influence on the simulation outputs with respect to wheat production and environment risks under UK soils and climatic conditions, and (2) to examine and understand the relationships between parameter sensitivity and model outputs.

2. Materials and Methods

2.1. Study Site

Data from a winter wheat field experiment was used for this study. The experiment was conducted for three growing seasons from 2011 to 2014 in the experimental plots of “Exhaustion Land”, one of the classical long-term experiments at Rothamsted Research in Harpenden (51°49′ N, 0°21′ W and 128 m a.s.l.). The winter wheat cultivar was Xi-19. The soil is classified as a Chromic Luvisol (FAO classification) with a silty clay loam texture topsoil. Daily maximum and minimum air temperatures and precipitation during the growing years are shown in Figure 1. Compared to the mean climatic conditions over the growing season between 1981 and 2010 (mean temperature is 8.7 °C, precipitation 611mm), the first growing season and the third growing season were warmer and wetter (9.4 °C and 770 mm, 10.1 °C and 814 mm, respectively), but the second season was cooler and dryer (7.9 °C and 672 mm, respectively). The total sunshine hours of the three growing seasons (1314, 1268, and 1363 h, respectively) were all higher than the 1981 to 2010 mean (1243 h). Data for the study period were downloaded from the electronic Rothamsted Archive (e-RA, http://www.era.rothamsted.ac.uk/, accessed on 30 June 2021).

2.2. Parameters and Simulated Outputs

A total of 61 model required parameters, describing both water redistribution and C and N processes, were chosen to investigate the responses of the system. To reduce complexity in the analysis, parameters that relate to crop growth and development were set as the usual default values, and thus not part of the 61 parameters investigated. The ranges (maximum to minimum) of the chosen parameters were ± 50% of the default values used in previous studies [29,30,31]. Each parameter was set to one of 100 different values to ensure sampling points are uniformly distributed across the probability distribution. A total number of 6100 SPACSYS simulations were run, each with a 3-years length, meaning that only one of 61 parameters (Table A1) was investigated in turn, whilst all others were set at their default. Details on the linkage of each parameter to each process can be found in various applied SPACSYS model studies [13,32] and also in subsequent, developmental studies [31].

For the simulation outputs, a total of 27 variables were considered. Variables include those for grain yield and above-ground dry matter biomass of winter wheat, soil water dynamics, losses (runoff, leaching, and release) from soil and nutrient pool sizes (Table A2).

2.3. Sensitivity Analysis and Diagnostics

As indicated above, only the effect of a single parameter change on the simulated outputs was considered—i.e., an LSA was followed. We analysed sensitivity to the 61 parameters by ranking the values of the root mean square deviation (RMSD) of simulated outputs as follows:

(1)$\mathrm{RMSD}=\sqrt{\frac{1}{n}{\displaystyle \sum }_{i=1}^n{\left(D_{}-S_i\right)}^2}$

where D represents the simulated value using the default parameter values (Wu, et al., 2019) and S_i represents the simulated values as defined above, at the ith step of parameter changes and n is equal to 100 (the number of simulations). In order to ensure that all results from the SPACSYS runs, with each of 61 parameters varying in turn, can be assessed and objectively compered to each other, we used the normalized RMSD (NRMSD) as follows:

(2)$\mathrm{NRMSD}=\frac{RMSD}{S_{max}}\times 100\%$

where S_max is the maximum value of all the simulated outputs for the whole set (100 simulations) of a single parameter over the entire simulated period. The NRMSD diagnostic combines RMSD values across multiple output variables to produce an overall measure of model sensitivity, while normalising for the different scales of the outputs. It can be interpreted as a fraction of the maximum values. The higher the NRMSD, the more sensitive the output variables are to the given parameter specification. A parameter that yields an NRMSD equal to 0% means it contributes very little on the simulated outputs and can be taken as a fixed parameter [33]. Simulations are considered sensitive to parameters yielding NRMSD values > 1%.

In order to capture the effects of different weather conditions on parameter sensitivity, we calculated the NRMSD in each individual growing year (2011/12, 12/2013, and 13/2014) and consecutive years (2011/13, 12/2014, and 11/2014). Further clarity was provided with averaged NRMSDs for the individual growing years and consecutive years:

(3)$NRMS{D}_{IN}=\frac{NRMS{D}_1+NRMS{D}_2+NRMS{D}_3}{3}$

(4)$NRMS{D}_{CO}=\frac{NRMS{D}_{12}+NRMS{D}_{23}+NRMS{D}_{123}}{3}$

where NRMSD_IN and NRMSD_CD represent the averaged NRMSD for the individual and the consecutive years, respectively. NRMSD₁, NRMSD₂, NRMSD₃, NRMSD₁₂, NRMSD₂₃, and NRMSD₁₂₃ represent the NRMSD for years 2011/12, 2012/13, 2013/14, 2011/13, 2012/14, and 2011/14, respectively. All analyses were carried out within the R statistical programming environment (version 3.5.0).

2.4. Simulated Climate and Soil Data

In order to assess the impact of climate and soil type on the results of the sensitivity analysis, four additional simulation runs were built with the same configuration as the main simulation experiment described—replacing the weather over the simulation period with another set collected for South West England (50°46′10′′ N, 3°54′05′′ W) and replacing soil properties with (a) a sandy soil with 10% clay, 70% sand, and 20% silt content in the top 10 cm layer, with 12% clay, 72% sand, and 16% silt content in the 10–20 cm layer and with 15% clay, 75% sand, and 10% silt content in the 20–30 cm layer; (b) a loam soil with 20% clay, 40% sand, and 40% silt content in the top 10 cm layer, with 22% clay, 42% sand, and 36% silt content in the 10–20 cm layer and with 25% clay, 45% sand, and 30% silt content in the 20–30 cm layer; and (c) a clay soil with 80% clay, 10% sand, and 10% silt content in the top 10 cm layer, with 76% clay, 12% sand, and 12% silt content in the 10–20 cm layer and with 70% clay, 15% sand, and 15% silt content in the 20–30 cm layer. Additional simulations were configured as described in Section 2.2, where the use of default parameter values was considered reasonable.

3. Results

Figure 2 maps the influence of each of the 61 parameters with respect to their effect on the simulations of 27 output variables both for the individual growing years and the consecutive years. The values of NRMSD range between 0 to 16%, with 87.83% between 0 and 1%. Output variables were insensitive to 11 of the 61 parameters, as their NRMSDs were all 0. The 11 parameters were: maximum nitrifier growth rate (MNG), maximum growth yield on NO (MGY-NO), maximum growth yield on N₂O (MGY-N₂O), half ammonium concentration (HAC), half NOx concentration (HNC), water content interval to unity (WCI), unsaturated conductivity decrease (UCD), half saturation global radiation intensity (HSG), drain pipe diameter (DPD), and finally, half closure vapour pressure deficit (HCV).

The influences of the parameters on the chosen output variable under alternative weather conditions are shown in Figure 3 and with the three different soil types in Figure 4 (sandy), Figure 5 (loam), and Figure 6 (clay). The results from these additional climate and soil simulations were hardly different to those from the main simulation experiment (Figure 2), where the sensitive parameters with NRMSD greater than 1% were almost identical. However, the number of fixed parameters reduced, especially in the sandy soil where none of the 61 parameters could be taken as fixed.

As the most sensitive parameters had similar trends across all the simulations (main and additional), regardless of climate or soil type, we now focus our attention only on the relationships between the parameters and the output variables from the main simulation experiment. Here the greatest sensitivity to parameter specification was observed for simulated outputs grouped as: (1) soil water dynamics, (2) crop dry matter accumulation, and (3) N and C losses and pools.

3.1. Soil Water Dynamics

The drainpipe level (DPL) parameter had the greatest influence on the simulations for ground water flow (GWF) and water content change (WSC) with NRMSD values > 7% (Figure 2), especially when the specified value of DPL was from the surface to the 2 m soil depth (Figure 7b,e,h). Here GWF increased and then plateaued at around 2 m, while WSC decreased moving from positive to negative values. For DPL values between 0 and −2 m, surface runoff (SRO) decreased from about 570 mm/year to values of zero. Likewise, other parameters relating to the drainage system, such as distance between drainpipes (DBD) and minimum roughness length (MRL) were similarly influential with NRMSD values > 1% and > 3% (Figure 2) for GWF, WSC, and SRO outputs (Figure 7). Here GWF, WSC, and SRO outputs all decreased as the MRL parameter increased from its minimum (0.001 m) to maximum specification (1 m) (Table A1 in Appendix A), while only GWF and SRO outputs decreased as the DBD parameter increased from its minimum (2 m) to maximum specification (100 m) (Table A1). The influence of minimum hydraulic conductivity (MHC), runoff first order rate coefficient (RFO) and maximum surface storage (MSS) parameters on GWF, WSC, and SRO outputs are given in Figure A1 in Appendix A.

3.2. Dry Matter of Winter Wheat

The simulation of dry matter of leaves, stems, and grains for winter wheat (LDM, SDM, and GDM) are particularly sensitive to parameters DBD, MRL, DPL, and AIF (ammonium immobilisation fraction) following their NRMSD values (Figure 2). DBD was the most sensitive parameter (NRMSD > 3%) for LDM, SDM, and GDM in the individual years, with increases in dry matter simulation only for DBD < 40 m. However, SDM and GDM were most sensitive to AIF with NRMSD > 4%, where dry matter simulation increased slightly across the full range of the AIF values specified (i.e., 0 to 1). Increases in dry matter were found for DPL > −2 m, while a stepped but small increase in dry matter was observed for MRL across its minimum to maximum values specified (Figure A2 in Appendix A).

3.3. Nitrogen and Carbon Losses from Soil

3.3.1. Losses with Surface Runoff

The simulations for nitrate (NO3) loss with surface runoff (NOR) were also highly sensitive to DPL, with NRMSD values of 1.38% and 1.73% in the individual and pairs of consecutive years, respectively (Figure 2). For DPL > −2 m, NOR displayed an increasing trend (Figure 8a). As would be expected, simulations for NOR were also sensitive to the fertiliser dissolution rate (SFD), where NOR simulations displayed a non-linear increase for SFD < 0.2. (Figure 8b). With NRMSD values > 0%, the simulations of dissolved N loss (NDR) and C loss (CDR) with surface runoff were sensitive to the loss of litter parameter (LLF) (Figure 8c), and the simulations of residue N loss (NRR) and residue C loss (CRR) were sensitive to the loss of residue parameter (RLF) (Figure 8d).

3.3.2. Nitrogen and Carbon Leaching

With the specification of parameters DBD, MRL, DPL, MHC, and ESP (the empirical scale in pore shape), each had a clear influence on the simulations of leaching for NO₃ (NOL), NH₄ (NHL), N (NDL), and C (CDL) in the individual and pairs of consecutive years (Figure 2). The same variable simulations were also sensitive to parameters AIF, LLF, RLF, DFF (a transferring fraction of decomposed fresh litter to dissolved organic matter), DFL (a transferring fraction of litter to dissolved organic matter), CWF (coefficient in the water function for decomposition), and DPR (potential decomposition rate of dissolved organic matter). DBD was a highly influential parameter for NOL and NHL with NRMSDs > 3% whilst DPL was highly influential for NDL and CDL with NRMSDs > 4%. Relationships between output simulations of NOL, NHL, NDL, and CDL and parameters DBD, MRL, DPL, AIF, MHC, ESP, DFF, DFL, CWF, DPR, LLF, and RLF are given in Figure A3 and Figure A4 in Appendix A.

3.3.3. Gas Emissions

Simulations of the N₂O emission rate (N2O) were most sensitive to changes in the maximum autotrophic nitrification rate (MAN), secondly, the NO production fraction from nitrification (PFF-NO), and thirdly, the Q10 temperature coefficient for denitrification (Q1D) (Figure 2). Specification of the latter two parameters also strongly influenced the NO emission rate (NOE) and N₂ emission rate (N2E) (Figure 2). In addition, N2O simulations were found to be relatively sensitive to DBD, MRL, and DPL (Figure 2). Figure 9 displays the relationships between the nitrogenous gas emissions outputs (NOE, N2O, and N2E) and a range of values for parameters PFF-NO and Q1D. For the range of PFF-NO specified, the simulation of NOE, N2O, and N2E were erratic with no clear pattern. Further plots for NOE only are given in Figure A5 in Appendix A.

Parameters CWF, Q1M (Q10 temperature coefficient for mineralization), and ASF (assimilation factor), in general, provided the highest NRMSDs (Figure 2) with respect to influencing CO₂ emissions from various soil organic C pools (DRE: dissolved release; FLR: fresh litter release; HRE: humus release; MRE: microbial release). However, the simulation of HRE was most sensitive to the humus potential decomposition rate (HPD) (NRMSD > 15%) (Figure 2). Figure 10 displays the relationships between CO₂ emissions DRE, FLR, HRE, MRE with respect to the top-ranking sensitive parameters CWF, Q1M, ASF over the simulation period. Almost all CO₂ emission outputs deceased as CWF, Q1M, and ASF parameters were increased through their minimum/maximums set (Table A1); the exception being MRE with ASF. Further CO₂ emission plots depicting their sensitivity to key parameters are given in Figure A6, Figure A7 and Figure A8 in Appendix A.

3.4. Changes of Soil C and N Pools

Although simulations of soil N and C pools in humus (NHP and CHP, respectively) were most sensitive to HPD (Figure 2), no discernible trend in this sensitivity was observed (Figure A9a in Appendix A). Simulations of N and C dissolved (NDP and CDP, respectively) were most sensitive to CWF, DFL (dissolved fraction in litter), HFL (humus fraction from fresh litter), and DBD (Figure 2), where in this case, clear increasing trends in this sensitivity were observed, some plateauing (Figure A9b,e for CWF and DBD, respectively), some not (Figure A9c,d for DFL and HFL, respectively). Further soil N and C pools plots depicting their sensitivity to key parameters are given in Figure A10 and Figure A11 in Appendix A.

4. Discussion

Parameters describing drainage implementation (DBD and DPL) play a critical role in the performance of the agriculture system reflecting changes in the weather and soil condition. For this study, simulations for soil water dynamics, crop dry matter accumulation and N and C losses, were all found to be highly sensitive to the specification of DBD and DPL. Our results corroborate with Ballantine and Tanner [34], and Ritzema [35], where both studies demonstrated how the drainage system was an important externality for agricultural production to prevent waterlogging and salinization of the soil in arid and semi-arid regions and strategically control drainage of excess water from the soil profile. It has also been reported that water use efficiency for wheat with a controlled drainage with a varying depth was 40% higher than that with the conventional drainage, depending on the crop stage [36]. Tomic et al. [37] showed that a DBD value of 25 m provided the highest wheat yield among seven different drainage systems, which is in agreement with our result in that DBD is one of the most sensitive parameters for the dry matter of winter wheat. However, this might not be transferrable to other crops, e.g., drainage control was found to have no significant influence on soybean yield [38]. Designing an appropriate drainage system in humid areas can reduce N loss from an agricultural field [39]. By raising the water level in the soil profile, which is equivalent to increasing DPL in our study here, controlled drainage has the potential to increase denitrification in the anaerobic zone [40]. Design of subsurface drainpipes involves the determination of depth (DPL), spacing (DBD) and pipe diameter (DPD) [41]. As a parameter of a drainage system, however, DPD had no influence on any of the 27 selected output variables (Figure 2), which implied the effects of vertical water flow were stronger than horizontal water flow towards the pipes in our study’s wheat field.

Our analysis suggests that simulations were sensitive to changes in MRL (minimum roughness length), where it is known to strongly affect water and heat fluxes on the soil surface [42,43]. Wang et al. [44] have shown that an accurate assessment of soil moisture might be problematic if MRL is inappropriately parameterised. However, our sensitivity analysis showed that SRI, HSG, and HCV (all related to evaporation and transpiration processes, see Table A1) had no influence on any of the 27 selected outputs.

4.1. Sensitive Parameters for Water Dynamics

Soil water dynamics is a key component for nutrient cycling and crop growth and development. Our study showed that soil water dynamics were sensitive to MHC (minimum hydraulic conductivity), DPL, DBD, and MRL (Figure 2 and Figure A1). This is not surprising given that: (1) MHC is a key parameter to control the soil water infiltration process where water moves downward from the surface [45] and (2) soil hydraulic conductivity determines surface energy and water fluxes and then soil water content [46,47,48].

4.2. Sensitive Parameters for Yield and Biomass

For simplicity, we did not include in our sensitivity analysis, the plant biological parameters. Therefore, the identified sensitivities influence biomass accumulation and the grain yield of winter wheat indirectly via other processes. Here DBD, MRL, DPL, and AIF are all important parameters for yield and aboveground biomass simulation (where DBD, MRL, and DPL also influence water dynamics in the soil, as discussed above). Unarguably, soil water content itself plays an important role in yield and aboveground dry matter accumulation under non-irrigation conditions [49,50]. It regulates N cycling and its availability to crops is important, as this in turn, affects crop growth and yield formation [51,52,53]. AIF that can be immobilised when the immobilisation process occurs, impacts crop growth and yield through changing mineralisation, and in sequence, the status of soil N and C content [1,54].

4.3. Sensitive Parameters for Losses from Soil

The drainage parameters DPL and DBD also had significant influence on C and N losses through surface runoff and leaching from soils. In addition, SFD (fertiliser dissolution rate), LLF (loss of litter), and RLF (loss of residue) that are each related to external nutrient inputs also played a critical role in the simulations of surface runoff losses. Results showed that SFD was an influential parameter for the simulations of NHL (NH4 leach); LLF was influential for the simulations of NDR and CDR (dissolved N and C loss); and RLF was influential for the simulations of NRR and CRR (residue N and C loss). External nutrient inputs can increase soil nutrient pools, and in turn, nutrient losses when surface runoff occurs.

Clearly, parameters related to nitrification and denitrification processes have a great influence on nitrogenous gas emissions. A higher nitrification rate increases N availability for the denitrification process that is enhanced by Q1D (Q10 temperature coefficient for denitrification). Specification of MAN (maximum autotrophic nitrification rate) had the biggest influence on NOE (NO emission rate), and both PFF-NO (NO production fraction from nitrification) and Q1D influenced NOE, N2O (N₂O emission rate), and N2E (N₂ emission rate). Conversely, simulations for soil CO2 release were sensitive to a different set of parameters that control nutrient cycling. Here, the specification for Q1M (Q10 temperature coefficient for mineralization), CWF (coefficient in the water function for decomposition), and ASF (assimilation factor) all had a significant influence on DRE (dissolved release), FLR (fresh litter release), HRE (humus release), and MRE (microbial release) (i.e., outputs variables that generate from decomposition from various soil C pools). Our results are supported by previous studies [55,56,57]. The soil organic matter decomposition rate depends on substrates [58,59], a potential decomposition rate and other abiotic and biotic factors including soil temperature (related to Q1M), soil moisture (i.e., CWF), and microbial activities (i.e., ASF) [60].

5. Conclusions

The local sensitivity analysis for the SPACSYS model using data from a winter wheat field experiment found that parameters related to drainage not only affected drainage fluxes but also grain yield and aboveground dry matter accumulation of winter wheat, and losses of water, C and N from soil under given weather and soil conditions. In addition, parameters related to mineralization and immobilisation processes influenced crop yield and biomass accumulation significantly. Further, parameters that control the nitrification and denitrification processes had a great influence on nitrogenous gas emissions. All results were found to be largely insensitive to different climatic and soil conditions, where identified parameter sensitivity remained relatively the same.

Overall, this study has provided evidence of which parameters require the most care in their specification and in turn, which parameters need to be determined accurately. This evidence can lead to efficient and cost-effective sampling and lab analyses for such parameters. Results provide a reference for model users on which parameters are key for different simulation goals. However, although this study has provided some important advances in model understanding and use, we have only considered the influence of single parameter change at a time on the model outputs. Interaction effects through a global sensitivity analysis were not assessed, and as such, are a focus of current and complementary research.

Author Contributions

Conceptualization, L.W.; methodology, L.W., P.H. and Y.S.; software, L.W.; validation, Y.S.; formal analysis, Y.S.; writing—original draft preparation, Y.S., M.H., P.H. and L.W.; writing—review and editing, Y.S., M.H., P.H. and L.W.; visualization, Y.S. and P.H.; supervision, L.W. All authors have read and agreed to the published version of the manuscript.

Funding

Y.S. was funded by the Chinese Scholarship Council (grant number 201706300103). The research was funded by the BBSRC Institute Strategic Programme (ISP) grant, “Soils to Nutrition” work package 2 grant number BBS/E/C/000I0320.

Institutional Review Board Statement

Not applicable.

Acknowledgments

We thank Andrew Mead for the design of statistical analysis.

Conflicts of Interest

The authors declare no conflict of interest.

Appendix A

Table A1

List of the parameters chosen for the sensitivity analysis, together with their maximum and minimum units specified.

No.	Abbreviations	Variable Description	Unit	Min	Max
1	SFD	Specific fertiliser dissolution rate	1/day	0.0001	0.5
2	Q1N	Coenzyme Q10 temperature coefficient for nitrification	-	0.1	5
3	Q1M	Coenzyme Q10 temperature coefficient for mineralization	-	0.1	5
4	Q1D	Coenzyme Q10 temperature coefficient for denitrification	-	0.1	5
5	PFF-NO	NO production fraction from nitrification	-	0.00001	0.1
6	PFF-N2O	N2O production fraction from nitrification	-	0.000001	0.01
7	MHC	Minimum hydraulic conductivity	mm/day	0	10
8	MMR	Microbial maintenance respiration rate	1/day	0.01	2
9	MDG-NO3	Maximum NO3 denitrifier growth rate	1/day	1	40
10	MDG-N2O	Maximum N2O denitrifier growth rate	1/day	1	20
11	MDG-NO2	Maximum NO2 denitrifier growth rate	1/day	1	40
12	MDG-NO	Maximum NO denitrifier growth rate	1/day	1	20
13	MNG	Maximum nitrifier growth rate	1/day	0.1	15
14	MND	Maximum nitrifier death rate	1/day	0.01	5
15	MGY- NO3	Maximum growth yield on NO3	gC/gN	0.01	1
16	MGY- NO2	Maximum growth yield on NO2	gC/gN	0.01	1
17	MGY- NO	Maximum growth yield on NO	gC/gN	0.01	3
18	MGY- N2O	Maximum growth yield on N2O	gC/gN	0.01	4
19	MAN	Maximum autotrophic nitrification rate	1/day	0.001	0.1
20	MCC	Maintenance coefficient on carbon in denitrification	1/day	0.0001	0.02
21	HPD	Humus potential decomposition rate	1/day	0.0000001	0.001
22	HFL	Partitioning fraction to humus from decomposed fresh litter	-	0	1
23	HFD	Partitioning fraction to humus from decomposed dissolved organic matter	-	0	1
24	HDC	Michaelis constant on dissolved organic carbon concentration	g/m3	1	40
25	HAC	Michaelis constant on ammonium concentration	gN/m3	1	30
26	HNO	Michaelis constant on NOx concentration	gN/m3	50	200
27	FLD	Fresh litter potential decomposition rate	1/day	0.0001	0.1
28	DBD	distance between drainpipes	m	2	100
29	DPR	Potential decomposition rate of dissolved organic matter	1/day	0.0001	0.1
30	UCD	Unsaturated conductivity decrease	-	0	10
31	SCF	Soil cover fraction to prevent infiltration	-	0	1
32	RFO	Runoff first order rate coefficient	1/day	0.01	1
33	MRL	Minimum roughness length	m	0.001	1
34	MSS	Maximum surface storage (no runoff)	mm	0.01	10
35	HSG	Half saturation global radiation intensity	J/m2/day	0	10,000,000
36	ESP	Empirical scale in pore shape	-	0.01	10
37	DPL	Drain pipe level, negative downwards	m	-10	0
38	DPD	Drain pipe diameter	m	1	10
39	CAC	Corresponding water amount that the ground is fully covered	mm	1	200
40	RAP	Relative activity at porosity	-	0	1
41	WCI	Water content interval to unity	vol%	1	15
42	OAW	Optimal water content at which there is no adverse effect from soil moisture	vol%	3	40
43	CWF	Coefficient in the water function for decomposition	-	0	100
44	ASF	Assimilation factor	-	0	1
45	LLF	Unit loss fraction of litter with surface runoff	-	0	1
46	RLF	Unit loss fraction of residue with surface runoff	-	0	1
47	BTM	Reference temperature at which the reaction function is unity for mineralization	°C	0	35
48	RLT	A transferring fraction of residue to the litter pool	1/day	0	0.5
49	DFH	A transferring fraction of decomposed humus pool to dissolved organic matter	-	0	1
50	DFF	A transferring fraction of decomposed fresh litter to dissolved organic matter	-	0	1
51	DFL	A transferring fraction of litter to dissolved organic matter	-	0	1
52	HNC	Michaelis constant on nitrate concentration	gN/m3	5	15
53	BTN	Reference temperature at which the reaction function is unity for nitrification	°C	0	35
54	BTD	Reference temperature at which the reaction function is unity for denitrification	°C	0	35
55	WCA	Water content interval activity	vol%	0	100
56	FWD- NH4	NH4 fraction in wet deposition	-	0	1
57	FDD- NH4	NH4 fraction in dry deposition	-	0	1
58	AFI	Fraction of ammonium that cannot move freely with water	-	0	1
59	MCN	C:N ratio in microbial biomass	-	5	15
60	AIF	A fraction of ammonium that can be used for immobilisation	-	0	1
61	HCV	Vapour pressure deficit at which leaf stomata half closed

Table A2

List of simulated output variables from the SPACSYS model.

No.	Abbreviations	Output	Unit
1	GWF	Groundwater water flux	mm/year
2	WSC	Soil water content change in the soil profile	%
3	GDM	Grain dry matter	gDM/m2
4	LDM	Leaf dry matter	gDM/m2
5	SDM	Stem dry matter	gDM/m2
6	NOR	NO3 loss with surface runoff	gN/m2/year
7	NDR	Dissolved N loss with surface runoff	gN/m2/year
8	NRR	Residue N loss with surface runoff	gN/m2/year
9	CDR	Dissolved C loss with surface runoff	gC/m2/year
10	CRR	Residue C loss with surface runoff	gC/m2/year
11	SRO	Surface runoff	Mm/year
12	NHL	NH4 leaching	gN/m2/year
13	NOL	NO3 leaching	gN/m2/year
14	NDL	N dissolved leach	gN/m2/year
15	CDL	C dissolved leach	gN/m2/year
16	NOE	NO emission rate	gN/m2/year
17	N2O	N2O emission rate	gN/m2/year
18	N2E	N2 emission rate	gN/m2/year
19	DRE	Dissolved release	gC/m2/year
20	FLR	Fresh litter release	gC/m2/year
21	HRE	Humus release	gC/m2/year
22	MRE	Microbial release	gC/m2/year
23	NHP	N in humus	gN/m2
24	NDP	Dissolved N	gN/m2
25	NMR	Mineralization rate	gN/m2/year
26	CDP	Dissolved C	gC/m2
27	CHP	C in humus	gC/m2

Figure A1

The relationships between ground water flow (GWF), water content change (WSC) and surface runoff (SRO) with changes in parameters with NRMSE > 1% (minimum hydraulic conductivity (MHC), ammonium immobilised fraction (AIF), specific fertiliser dissolution rate (SFD) and maximum autotrophic nitrification rate (MAN)).

[Figure omitted. See PDF]

Figure A2

Response of grains (GDM), leaves (LDM) and stems (SDM) dry matter to the parameters—distance between drainpipes (DBD), minimum roughness length (MRL), drain pipe level (DPL) and ammonium immobilised fraction (AIF). All with NRMSD > 1% over the simulation period.

[Figure omitted. See PDF]

Figure A3

Relationships between simulated leached losses of N and C to their most influential parameters (by NRMSD). Simulated variables: NHL: NH4 leach; NOL: NO3 leach; NDL: N dissolve leach; CDL: C dissolve leach. Model parameters: DBD: distance between drainpipes; MRL: minimum roughness length; DPL: drain pipe level; AIF: ammonium immobilised fraction; MHC: minimum hydraulic conductivity; and ESP: empirical scale in pore shape.

[Figure omitted. See PDF]

Figure A4

The relationships between leached losses of N and C and the parameters with NRMSEs higher than 1% (CDL: C dissolve leach; NDL: N dissolve leach; DFF: dissolved fraction from fresh litter; DFL: dissolved fraction in litter; CWF: coefficient in water function in the process of organic matter decomposition; DPR: potential decomposition rate of dissolved organic matter; LLF: litter unit loss fraction; RLF: residue unit loss fraction).

[Figure omitted. See PDF]

Figure A5

The relationships between NOE and sensitive parameters (DBD: distance between drainpipes; MRL: minimum roughness length; DPL: drain pipe level; AIF: ammonium immobilised fraction; MAN: maximum autotrophic nitrification rate; SCF: soil cover fraction prevent infiltration; CWF: coefficient in water function in the process of organic matter decomposition; BTD: base temp to unity denitrification; BTN: base temp to unity nitrification).

[Figure omitted. See PDF]

Figure A6

The relationships between carbon dissolved release (DRE) and sensitive parameters (DBD: distance between drainpipes; MRL: minimum roughness length; DPL: drain pipe level; MHC: minimum hydraulic conductivity; HPD: humus potential decomposition rate; LLF: litter unit loss fraction; OAW: optimal available water content; HFD: humus fraction from fresh litter; DFF: dissolved fraction from fresh litter; DFL: dissolved fraction in litter; RLT: residue to litter transfer fraction; FLD: fresh litter potential decomposition rate; RLF: residue unit loss fraction; DPR: dissolved potential decomposition rate; BTM: base temp to unity mineralization).

[Figure omitted. See PDF]

Figure A7

The relationships between fresh litter release (FLR) and sensitive parameters (MRL: minimum roughness length; DPL: drain pipe level; MHC: minimum hydraulic conductivity; LLF: litter unit loss fraction; OAW: optimal available water content; RLT: residue to litter transfer fraction; FLD: potential decomposition rate of fresh litter; RLF: residue unit loss fraction; BTM: base temp to unity mineralization; DFL: dissolved fraction in litter; DFF: dissolved fraction from fresh litter; HFL: partitioning fraction to humus from decomposed fresh litter).

[Figure omitted. See PDF]

Figure A8

The relationships between microbial release (MRE) and sensitive parameters (DBD: distance between drainpipes; MRL: minimum roughness length; DPL: drain pipe level; MHC: minimum hydraulic conductivity; HPD: humus potential decomposition rate; HFL: partitioning fraction to humus from decomposed fresh litter; LLF: litter unit loss fraction; OAW: optimal available water content; HFD: partitioning fraction to humus from decomposed dissolved organic matter; DFF: dissolved fraction from fresh litter; DFL: dissolved fraction in litter; RLT: residue to litter transfer fraction; FLD: potential decomposition rate of fresh litter; RLF: Residue unit loss fraction; DPR: potential decomposition rate of dissolved organic matter; BTM: base temp to unity mineralization; MMR: microbial maintenance respiration rate).

[Figure omitted. See PDF]

Figure A9

Relationships between changes in soil organic C and N pools and the parameters they are most sensitive to over the simulation period. NHP: N humus pool; CHP: C humus pool; NDP: N dissolve pool; CDP: C dissolve pool; HPD: humus potential decomposition rate; CWF: coefficient in water function; DFL: dissolved fraction in litter; HFL: humus fraction from fresh litter; DBD: distance between drainpipes.

[Figure omitted. See PDF]

Figure A10

The relationships between nitrogen dissolved pool (NDP) and sensitive parameters (MRL: minimum roughness length; DPL: drain pipe level; MHC: minimum hydraulic conductivity; HPD: humus potential decomposition rate; Q1M: Q10 value for mineralization; ASF: assimilation factor; LLF: litter unit loss fraction; OAW: optimal available water content; DFF: dissolved fraction from fresh litter; FLD: potential decomposition rate of fresh litter; RLF: residue unit loss fraction; DPR: potential decomposition rate of dissolved organic matter; BTM: base temp to unity mineralization).

[Figure omitted. See PDF]

Figure A11

The relationships between carbon dissolved pool (CDP) simulations and sensitive parameters (MRL: minimum roughness length; DPL: drain pipe level; MHC: minimum hydraulic conductivity; CWF: coefficient in water function in the process of organic matter decomposition; Q1M: Q10 value for mineralization; ASF: assimilation factor; LLF: litter unit loss fraction; DFF: dissolved fraction from fresh litter; RLT: residue to litter transfer fraction; FLD: potential decomposition rate of fresh litter; RLF: residue unit loss fraction; DPR: potential decomposition rate of dissolved organic matter; BTM: base temp to unity mineralization).

[Figure omitted. See PDF]

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

Figures

Figure 1. Daily rainfall and mean air temperature over the observational period (2011 to 2014).

Figure 2. Heat maps of averaged NRMSDs for the 61 parameters (rows) with respect to the 27 output variables (columns) under the experimental condition from individual years (left) and the pairs of consecutive years (right). A yellow pixel represents no interaction.

Figure 3. Heat maps of averaged NRMSDs for the 61 parameters (rows) with respect to the 27 output variables (columns) with the climate condition in South West England from individual years (left) and the pairs of consecutive years (right). A yellow pixel represents no interaction.

Figure 4. Heat maps of averaged NRMSDs for the 61 parameters (rows) with respect to the 27 output variables (columns) in the sandy soil from consecutive years. A yellow pixel represents no interaction.

Figure 5. Heat maps of averaged NRMSDs for the 61 parameters (rows) with respect to the 27 output variables (columns) in the loam soil from consecutive years. A yellow pixel represents no interaction.

Figure 6. Heat maps of averaged NRMSDs for the 61 parameters (rows) with respect to the 27 output variables (columns) in the clay soil from consecutive years. A yellow pixel represents no interaction.

Figure 7. Response of averaged drainage water fluxes and the change in soil water content change (i.e., groundwater flow (GWF, a–c), water content change (WSC, d–f) and surface runoff (SRO, g–i)) to parameter change for distance between drainpipe (DBD), drain pipe level (DPL), and minimum roughness length (MRL) over the simulation period.

Figure 8. Relationships for simulations of: (a) NO3 runoff (NOR) and parameter drain pipe level (DPL), (b) NOR and parameter specific fertiliser dissolution rate (SFD), (c) dissolved N loss (NDR) and C dissolved loss (CDR) and parameter litter loss (LLF), and (d) residue N loss (NRR) and residue C loss (CRR) and parameter residue loss (RLF) over the simulation period.

Figure 9. Relationships between simulations of nitrogenous gas emissions (NOE: NO emission rate (a,b), N2O: N2O emission rate (c,d), and N2E: N2 emission rate (e,f)) and the parameters of NO production fraction from nitrification (PFF-NO) and Q10 value for denitrification (Q1D).

Figure 10. Relationships between simulations of CO2 emissions from soil organic C pools (DRE: dissolved release; FLR: fresh litter release; HRE: humus release; MRE: microbial release) and the parameters they are most sensitive to (CWF: coefficient in water function; Q1M: Q10 value for mineralization; ASF: assimilation factor) over the simulation period. DRE and FLR vs CWF (a), Q1M (b) and ASF (c). HRE and MLR vs CWF (d), Q1M (e) and ASF (f).