Soil-related developments of the Biome-BGCMuSo v6.2 terrestrial ecosystem model

Abstract

Terrestrial biogeochemical models are essential tools to quantify climate–carbon cycle feedback and plant–soil relations from local to global scale. In this study, a theoretical basis is provided for the latest version of the Biome-BGCMuSo biogeochemical model (version 6.2). Biome-BGCMuSo is a branch of the original Biome-BGC model with a large number of developments and structural changes. Earlier model versions performed poorly in terms of soil water content (SWC) dynamics in different environments. Moreover, lack of detailed nitrogen cycle representation was a major limitation of the model. Since problems associated with these internal drivers might influence the final results and parameter estimation, additional structural improvements were necessary. In this paper the improved soil hydrology as well as the soil carbon and nitrogen cycle calculation methods are described in detail. Capabilities of the Biome-BGCMuSo v6.2 model are demonstrated via case studies focusing on soil hydrology, soil nitrogen cycle, and soil organic carbon content estimation. Soil-hydrology-related results are compared to observation data from an experimental lysimeter station. The results indicate improved performance for Biome-BGCMuSo v6.2 compared to v4.0 (explained variance increased from 0.121 to 0.8 for SWC and from 0.084 to 0.46 for soil evaporation; bias changed from $-$0.047 to$-$0.007 ${\mathrm{m}}^{\mathrm{3}}{\mathrm{m}}^{-\mathrm{3}}$ for SWC and from $-$0.68 to $-$0.2 $\mathrm{mm}{\mathrm{d}}^{-\mathrm{1}}$ for soil evaporation). Simulations related to nitrogen balance and soil ${\mathrm{CO}}_{\mathrm{2}}$ efflux were evaluated based on observations made in a long-term field experiment under crop rotation. The results indicated that the model is able to provide realistic nitrate content estimation for the topsoil. Soil nitrous oxide (${\mathrm{N}}_{\mathrm{2}}\mathrm{O}$) efflux and soil respiration simulations were also realistic, with overall correspondence with the observations (for the ${\mathrm{N}}_{\mathrm{2}}\mathrm{O}$ efflux simulation bias was between $-$0.13 and$-$0.1 $\mathrm{mg}\mathrm{N}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{d}}^{-\mathrm{1}}$, and normalized root mean squared error (NRMSE) was 32.4 %–37.6 %; for${\mathrm{CO}}_{\mathrm{2}}$ efflux simulations bias was 0.04–0.17 $\mathrm{g}\mathrm{C}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{d}}^{-\mathrm{1}}$, while NRMSE was 34.1 %–40.1 %). Sensitivity analysis and optimization of the decomposition scheme are presented to support practical application of the model. The improved version of Biome-BGCMuSo has the ability to provide more realistic soil hydrology representation as well as nitrification and denitrification process estimation, which represents a major milestone.

1 Introduction

The construction and development of biogeochemical models (BGMs) represent the response of the scientific community to address challenges related to climate change and human-induced global environmental change. BGMs can be used to quantify future climate–vegetation interaction including climate–carbon cycle feedback, and as they simulate plant production, they can be used to study a variety of ecosystem services that are related to human nutrition and resource availability (Asseng et al., 2013; Bassu et al., 2014; Huntzinger et al., 2013). Similarly to models describing various and complex environmental processes, the structure of biogeochemical models reflects our current knowledge about a complex system with many internal processes and interactions.

Processes of the atmosphere–plant–soil system take place on different temporal (sub-daily to centennial) scales and are driven by markedly different mechanisms that are quantified by a large diversity of modeling tools (Schwalm et al., 2019). Plant photosynthesis is an enzyme-driven biochemical process that has its own mathematical equation set and related parameters (and a large body of literature; e.g., Farquhar et al., 1980; Medlyn et al., 2002; Smith and Dukes, 2013; Dietze, 2013). Allocation of carbohydrates in the different plant compartments is studied extensively and also has a large body of literature and mathematical tool set (Friedlingstein et al., 1999; Olin et al., 2015; Merganičová et al., 2019). Plant phenology is quantified by specific algorithms that are rather uncertain components of the models (Richardson et al., 2013; Hufkens et al., 2018; Peaucelle et al., 2019). Soil biogeochemistry is driven by microbial and fungal activity and also has its own methodology and a vast body of literature (Zimmermann et al., 2007; Kuzyakov, 2011; Koven et al., 2013; Berardi et al., 2020). Emerging scientific areas like the quantification of the dynamics of non-structural carbohydrates (NSCs) in plants have a separate methodology that requires mathematical representation in models (Martínez-Vilalta et al., 2016). Simulation of land surface hydrology including evapotranspiration is typically handled by some variant of the Penman–Monteith equation that is widely studied and thus represents a separate scientific field (McMahon et al., 2013; Doležal et al., 2018).

Putting it all together, if we are about to construct and further improve a biogeochemical model to consider novel findings and track global changes, we need comprehensive knowledge that integrates many, almost disjunct scientific fields. Clearly, transparent and well-documented development of a biogeochemical model is of high priority but is challenging from the very beginning and requires cooperation of researchers from various scientific fields.

Continuous model development is inevitable, but it has to be supported by extensive comparison with observations and some kind of implementation of the model–data fusion approach (Keenan et al., 2011). It is well documented that structural problems might trigger incorrect parameter estimation that might be associated with distorted internal processes (Sándor et al., 2017; Martre et al., 2015). In other words, one major issue with BGMs (and in fact with all models using many parameters) is the possibility to get good simulation results for wrong reasons (which means incorrect parameterization) due to compensation of errors (Martre et al., 2015). In order to avoid this issue, any model developer team has to make an effort to also focus on internal ecosystem conditions (e.g., soil volumetric water content – SWC), nutrient availability, stresses) and other processes (e.g., decomposition) rather than the main simulated processes (e.g., photosynthesis, evapotranspiration).

Historically, biogeochemical models have been developed to simulate the processes of undisturbed ecosystems with simple representation of the vegetation (Levis, 2010). As the focus was on the carbon cycle, water and nitrogen cycles as well as related soil processes were not well represented. Incorrect representation of SWC dynamics is still an issue with models, especially in drought-prone ecosystems (Sándor et al., 2017). Additionally, human intervention representation (management) is still incomplete in some state-of-the-art BGMs; e.g., thinning, grass mowing, grazing, tillage, or irrigation is missing in some models (see Table A1 in Friedlingstein et al., 2020).

In contrast, crop models with different complexity were used for about 50 years or so to simulate the processes of managed vegetation (Jones et al., 2017; Franke et al., 2020). As the focus of crop models is on final yield due to economic reasons, the carbon balance or the full greenhouse gas balance was not addressed or was just partially addressed originally. Crop models typically have a sophisticated representation of soil water balance with a multilayer soil module that usually calculates plant response to water stress as well. Nutrient stress, soil conditions during planting, consideration of multiple phenological phases, heat stress during anthesis, vernalization, manure application, fertilization, harvest, and many other processes have been implemented over the decades (Ewert et al., 2015). Therefore, it seems to be straightforward to exploit the benefits of crop models and implement sound and well-tested algorithms into BGMs.

Starting from the well-known Biome-BGC model originally developed to simulate undisturbed forests and grasslands using a simple single-layer soil submodel (Running and Hunt, 1993; Thornton and Rosenbloom, 2005), we developed a complex, more sophisticated model (Hidy et al., 2012, 2016). Biome-BGCMuSo v4.0 (Biome-BGC with Multilayer Soil module) uses a seven-layer soil module and is capable of simulating different ecosystems from natural grassland to cropland including several management options (mowing, grazing, thinning, planting, and harvest), taking into account many environmental effects (Hidy et al., 2016). The developments included improvements regarding both soil and plant processes. In a nutshell, the most important soil-related developments were the improvement of the soil water balance module by implementing routines for estimating percolation, diffusion, pond water formation, and runoff as well as the introduction of multilayer simulation for belowground processes in a simplified way. The most important plant-related developments involved the implementation of a routine for estimating the effect of drought on vegetation growth and senescence, the improvement of stomatal conductance calculation considering atmospheric ${\mathrm{CO}}_{\mathrm{2}}$ concentration, the integration of selected management modules, the implementation of new plant compartments (e.g., yield), the implementation of a C${}_{\mathrm{4}}$ photosynthesis routine, the implementation of photosynthesis and respiration acclimation of plants and temperature-dependent $Q_{\mathrm{10}}$, and empirical estimation of methane and nitrous oxide soil efflux.

Problems found with the Biome-BGCMuSo v4.0 simulation result (namely the poor representation of soil water content, the lack of sophisticated layer-specific soil nitrogen dynamics representation, or model-structure-related problems, such as the lubber parameterization of the model) marked the path for further developments.

The aim of the present study is to provide detailed documentation of the current, improved version of Biome-BGCMuSo v6.2, which has many new features and facilitates various in-depth investigations of ecosystem functioning. Due to large number of developments, this paper focuses only on the soil-related model improvements. Case studies are also presented to demonstrate the capabilities of the new model version and to provide guidance for the model user community.

2 The original Biome-BGC model

Biome-BGC was developed from the Forest-BGC mechanistic model family in order to simulate vegetation types other than forests. Biome-BGC was one of the earliest biogeochemical models that included explicit carbon and nitrogen cycle modules. Biome-BGC simulates the storage and fluxes of water, carbon, and nitrogen within and between the vegetation, litter, and soil components of terrestrial ecosystems. It uses a daily time step and is driven by daily values of maximum and minimum temperatures, precipitation, solar radiation, and vapor pressure deficit (Running and Hunt, 1993). The model calculations apply to a unit ground area that is considered to be homogeneous.

The three most important components of the model are the phenological, carbon uptake and release, and soil flux modules. The core logic that is described below in this section remained intact during the developments. The phenological module calculates foliage development that affects the accumulation of C and N in leaf, stem (if present), and root and consequently the amount of litter. In the carbon flux module gross primary production (GPP) of the biome is calculated using Farquhar's photosynthesis routine (Farquhar et al., 1980) and the enzyme kinetics model based on Woodrow and Berry (2003). Autotrophic respiration is separated into maintenance and growth respiration. Maintenance respiration is calculated as the function of the N content of living plant pools, while growth respiration is an adjustable but fixed proportion of the daily GPP. The single-layer soil module simulates the decomposition of dead plant material (litter) and soil organic matter, N mineralization, and N balance in general (Running and Gower, 1991). The soil module uses the so-called converging cascade method (Thornton and Rosenbloom, 2005) to simulate decomposition, carbon and nitrogen turnover, and related soil ${\mathrm{CO}}_{\mathrm{2}}$ efflux.

The simulation has two basic steps. During the first (optional) spin-up simulation the available climate data series is repeated as many times as required to reach a dynamic equilibrium in the soil organic matter content to estimate the initial values of the carbon and nitrogen pools. The second, normal simulation uses the results of the spin-up simulation as initial conditions and runs for a given, predefined time period (Running and Gower, 1991). The so-called transient simulation option (which is the extension of the spin-up routine) is a novel feature in Biome-BGCMuSo v6.2 relative to the previous versions in order to ensure smooth transition between the spin-up and normal phase (Hidy et al., 2021).

In Biome-BGC, the main parts of the simulated ecosystem are defined as plant, litter, and soil. The most important pools include leaf (C, N, and intercepted water), root (C, N), stem (C, N), soil (C, N, and water), and litter (C, N). Plant C and N pools have sub-pools (actual pools, storage pools, and transfer pools). The actual sub-pools store C and N for the current year of growth. The storage sub-pools (essentially the non-structural carbohydrate pool, the source for the cores or buds) contain the amount of C and N that will be active during the next growing season. The transfer sub-pools inherit the entire content of the storage pools at the end of every simulation year. Soil C also has sub-pools representing various organic matter forms characterized by considerably different decomposition rates.

In spite of its popularity and proven applicability, the development of Biome-BGC was temporarily stopped (the latest official NTSG version is Biome-BGC 4.2; https://www.ntsg.umt.edu, last access: November 2021). One major drawback of the model was its relatively poor performance in modeling managed ecosystems and the simplistic soil water balance submodel using a single soil layer only.

Our team started to develop the Biome-BGC model further in 2006. According to the logic of the team, the new model branch was planned to be the continuation of the Biome-BGC model with regard to the original concept of the developers (keeping the model code open-source, providing detailed documentation, and providing support for users).

The starting point of our model development was Biome-BGC v4.1.1 that was a result of the model improvement activities of the Max Planck Institute (Vetter et al., 2008). Development of the Biome-BGCMuSo model branch has a long history by now. Previous model developments were documented in Hidy et al. (2012, 2016). Below, we provide a detailed description of the new developments that are included in Biome-BGCMuSo v6.2, which is the latest version released in September 2021. A comprehensive review of the input data requirement of the model, together with an explanation of the input data structure, is available in the user guide (Hidy et al., 2021). In this paper we refer to some input files (e.g., soil file, plant file) that are described in the user guide in detail.

One of the most important novelties and advantages of the new model version (Biome-BGCMuSo v6.2) compared to any previous versions due to the extensive and detailed soil parameter set (the current version has 79, MuSo 4.0 has 39, and the original model version has only 6 adjustable soil-related parameters) is that the parameterization of the model is much more flexible. But it might, of course, be a challenging task to define all of the input parameters. In order to support practical application of the model, the user guide contains proposed values for most of the new parameters (Hidy et al., 2021).

3 Soil-hydrology-related developments

In Biome-BGCMuSo v6.2 a 10-layer soil submodel was implemented. Previous model versions included a seven-layer submodel, which turned out to be insufficient to capture hydrological events like drying of the topsoil layers with sufficient accuracy. The thicknesses of the layers from the surface to the bottom are 3, 7, 20, 30, 30, 30, 30, 50, 200, and 600 cm. The center of the given layer represents the depth of each soil layer. Soil texture can be defined by the percentage of sand and silt for each layer separately along with the most important physical and chemical parameters (pH, bulk density, characteristic SWC values, drainage coefficient, hydraulic conductivity) in the soil input file (Hidy et al., 2021).

The water balance module of Biome-BGCMuSo has five major components to describe soil-water-related processes at daily resolution (listed here following the order of calculation): pond water accumulation and runoff, infiltration and downward gravitational flow (percolation), water movement within the soil (diffusion) driven by water potential gradient, evaporation and transpiration (root water uptake), and the downward and upward fluxes to and from groundwater. In the following subsections these five major components are described.

3.1 Pond water accumulation and runoff

Precipitation can reach the surface as rain or snow (below 0 ${}^{\circ }\mathrm{C}$ snow accumulation is assumed). Snow water melts from the snowpack as a function of temperature and radiation and is added to the precipitation input.

The canopy can intercept rain. The intercepted volume goes into the canopy water pool, which can evaporate. No canopy interception of snow is assumed. The throughfall (complemented with the amount of melted snow) gives the potential infiltration.

A new development in Biome-BGCMuSo v6.2 is that maximum infiltration is calculated based on the saturated hydraulic conductivity and the SWC of the topsoil layers. If the potential infiltration exceeds the maximum infiltration, pond water can be formed. If the sum of the precipitation and the actual pond height minus the maximum infiltration rate is greater than the maximum pond height, the excess water is added to surface runoff detailed below (Balsamo et al., 2009). The maximum pond height is an input parameter. Water from the pond can infiltrate into the soil at a rate at which the topsoil layer can absorb it. Evaporation of pond water is assumed to be equal to the potential evaporation.

Surface runoff is the water flow occurring on the surface when a portion of the precipitation cannot infiltrate into the soil. Two types of surface runoff processes can be distinguished: Hortonian and Dunne. Hortonian runoff is unsaturated overland flow that occurs when the rate of precipitation exceeds the rate at which water can infiltrate. The other type of surface runoff is the Dunne runoff (also known as the saturation overland flow), which occurs when the entire soil is saturated but the rain continues to fall. In this case the rainfall immediately triggers pond water formation and (above the maximum pond water height) surface runoff. The handling of these processes is presented in the soil hydrological module of Biome-BGCMuSo v6.2.

Calculation of Hortonian runoff ($\mathrm{kg}{\mathrm{H}}_{\mathrm{2}}\mathrm{O}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{d}}^{-\mathrm{1}}$) is based on a semi-empirical method and uses the precipitation amount ($\mathrm{cm}{\mathrm{d}}^{-\mathrm{1}}$), the unitless runoff curve number (RCN), and the actual moisture content status of the topsoil (Rawls et al., 1980; this method is known as the SCS runoff curve number method; SCS: soil conservation service). This type of runoff simulation can be turned off by setting RCN to zero. A detailed description can be found in Sect. S1 in the Supplement. The amount of runoff as a function of the soil type and the actual SWC is presented in Fig. 1.

Figure 1

Hortonian runoff as a function of rainfall intensity, soil type, and actual soil water content of the topsoil layer. Sand soil means 92 % sand, 4 % silt, and 4 % clay; silt soil means 8 % sand, 86 % silt, and 6 % clay; clay soil means 20 % sand, 20 % silt, and 60 % clay. SWC and SAT denote soil water content and saturation, respectively.

[Figure omitted. See PDF]

3.2 Infiltration, percolation, and diffusion

There are two optional methods in Biome-BGCMuSo v6.2 to calculate soil water movement between soil layers and actual SWC layer by layer. The first one is a cascade method (also known as tipping bucket method), and the second is a physical method based on the Richards equation. The tipping bucket method has a long history in crop modeling and is considered to be a successful, well-evaluated algorithm that can accurately simulate downward water flow in the soil.

The cascade method uses a semi-empirical input parameter (DC: drainage coefficient, ${\mathrm{d}}^{-\mathrm{1}}$) to calculate the downward water flow rate. When the SWC of a soil layer exceeds field capacity (FC), a fraction (equal to DC) of the water amount above FC goes to the layer next below. If DC is not set in the soil input file, it is estimated from the saturated hydraulic conductivity: $\text{DC}=\mathrm{0.1122}\cdot K_{\text{SAT}}^{\mathrm{0.339}}$ ($K_{\text{SAT}}$: saturated hydraulic conductivity, $\mathrm{cm}{\mathrm{d}}^{-\mathrm{1}}$; the user can set its value or the model based on soil texture estimates it internally; see Hidy et al., 2016). A detailed description of the method can be found in Sect. S2 in the Supplement. Drainage from the bottom layer is a net loss for the soil profile.

Water diffusion that is the capillary water flow between the soil layers is calculated to account for the relatively slow movement of water. The flow rate is a function of the water content difference of two adjacent layers and the soil water diffusivity at the boundary of the layers, which is determined based on the average water content of the two layers. A detailed mathematical description of the method can be found in Sect. S3 in the Supplement.

A detailed description of the Richards method can be found in Hidy et al. (2012). To support efficient and robust calculations of soil water fluxes a dynamically changing time step was introduced in version 4.0 (Hidy et al., 2016). The implementation of the more sophisticated Richards method is still in an experimental phase requiring rigorous testing and validation in the future.

3.3 Evapotranspiration

Biome-BGCMuSo, like its predecessor Biome-BGC, estimates evaporation of leaf-intercepted water, bare soil evaporation, and transpiration to estimate the total evapotranspiration at a daily level. The potential rates of all three processes are calculated based on the Penman–Monteith (PM) method. PM equation requires net radiation (minus soil heat flux) and conductance values by definition using different parameterization for the different processes. The model calculates leaf- and canopy-level conductances of water vapor and sensible heat fluxes to be used in Penman–Monteith calculations of canopy evaporation and canopy transpiration. Note that in the Biome-BGC model family the direct wind effect is ignored but can be considered indirectly by adjusting boundary layer conductance to site-specific conditions. A possible future direction might be the extension of the model logic to consider the wind effect directly.

3.3.1 Canopy evaporation

If there is intercepted water, this portion of evaporation is calculated using the canopy resistance (reciprocal of conductance) to evaporated water and the resistance to sensible heat. The time required for the water to evaporate based on the average daily conditions is calculated and subtracted from the day length to get the effective day length for evapotranspiration. Combined resistance to convective and radiative heat transfer is calculated based on canopy conductance of vapor and leaf conductance of sensible heat, both of which are assumed to be equal to the boundary layer conductance. Besides the conductance and resistance parameters the canopy-absorbed shortwave radiation drives the calculation. Note that the canopy evaporation routine was not modified significantly in Biome-BGCMuSo.

3.3.2 Soil evaporation

In order to estimate soil evaporation, first the potential evaporation is calculated, assuming that the resistance to vapor is equal to the resistance to sensible heat and assuming no additional resistance component. Both resistances are assumed to be equal to the actual aerodynamic resistance. Actual aerodynamic resistance is a function of the actual air pressure and air temperature as well as the potential aerodynamic resistance (${\text{potR}}_{\text{air}}$ in $\mathrm{s}{\mathrm{m}}^{-\mathrm{1}}$). ${\text{potR}}_{\text{air}}$ was a fixed value in the previous model versions (107 $\mathrm{s}{\mathrm{m}}^{-\mathrm{1}}$). Its value was derived from observations over bare soil in tiger bush in southwestern Niger (Wallace and Holwill, 1997). In Biome-BGCMuSo v6.2, the ${\text{potR}}_{\text{air}}$ is an input parameter that can be adjusted by the user (Hidy et al., 2021). Another new development in Biome-BGCMuSo v6.2 is the introduction of an upper limit for daily potential evaporation (${\text{evap}}_{\text{limit}}$) that is determined by the available energy (incident shortwave flux that reaches the soil surface):

1$${\text{evap}}_{\text{limit}}={\displaystyle \frac{\text{irad}\cdot \text{dayl}}{{\text{LH}}_{\text{vap}}}},$$ where irad is the incident shortwave flux density ($\mathrm{W}{\mathrm{m}}^{-\mathrm{2}}$), dayl is the length of the day in seconds, and ${\text{LH}}_{\text{vap}}$ is the latent heat of vaporization (the amount of energy that must be added to liquid to transform into gas; $\mathrm{J}{\mathrm{kg}}^{-\mathrm{1}}$). This feature was missing from previous model versions, resulting in considerable overestimation of evaporation on certain days that was caused by the missing energy limitation on evaporation.

A new feature in Biome-BGCMuSo v6.2 is the calculation of the actual evaporation from the potential evaporation and the square root of time elapsed since the last precipitation (expressed by days; Ritchie, 1998). This is another method that has been used by the crop modeler community for many years. A detailed description of the algorithm can be found in Sect. S4 in the Supplement.

One major new development in Biome-BGCMuSo v6.2 is the simulation of the reducing effect of surface residue or mulch cover on bare soil evaporation. Here we use the term “mulch” to quantify surface residue cover in general, keeping in mind that mulch is typically a human-induced coverage. Surface residue includes aboveground litter and coarse woody debris as well.

The evaporation reduction effect (evapREDmulch; unitless) is a variable between 0 and 1 (0 means full limitation, and 1 means no limitation) estimated based on a power function of the surface coverage (mulchCOV in %) and a soil-specific constant set by the user (pREDmulch; see Hidy et al., 2021). If variable mulchCOV reaches 100 % it means that the surface is completely covered. If mulchCOV is greater than 100 % it means the surface is covered by more than one layer. Surface coverage is a power function of the amount of mulch (mu, $\mathrm{kg}\mathrm{C}{\mathrm{m}}^{-\mathrm{2}}$) with parameters ${\text{p1}}_{\text{mulch}}$, ${\text{p2}}_{\text{mulch}}$, and ${\text{p3}}_{\text{mulch}}$ (soil parameters) based on the method of Rawls et al. (1980). $$\begin{array}{}\text{2}& {\displaystyle }\text{mulchCOV}={\text{p1}}_{\text{mulch}}\cdot (\text{mu}/{\text{p2}}_{\text{mulch}}{)}^{{\text{p3}}_{\text{mulch}}}\text{3}& {\displaystyle }\text{evapREdmulch}={\text{pREDmulch}}^{\frac{\text{mulchCOV}}{\mathrm{100}}}\end{array}$$ Another simulated effect of surface residue cover is the homogenization of soil temperature between 0 and 30 cm depth (layers 1, 2, and 3). The functional forms of surface coverage and the evaporation reduction factor are presented in Fig. 2.

Figure 2

Surface coverage as a function of the amount of surface residue or mulch (a) and the evaporation reduction factor (evapREDmulch) as a function of mulch coverage (b) using different mulch-specific soil parameters (pREDmulch). See text for details.

[Figure omitted. See PDF]

3.3.3 Transpiration

In order to simulate transpiration, first transpiration demand (TD, $\mathrm{kg}{\mathrm{H}}_{\mathrm{2}}\mathrm{O}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{d}}^{-\mathrm{1}}$) is calculated using the Penman–Monteith equation separately for sunlit and shaded leaves. TD is a function of leaf-scale conductance to water vapor, which is derived from stomatal, cuticular, and leaf boundary layer conductances. A novelty in Biome-BGCMuSo v6.2 is that potential evapotranspiration is also calculated using the maximal stomatal conductance instead of the actual stomatal conductance, which means that stomatal aperture is not affected by the soil moisture status (in contrast to the actual one).

TD is distributed across the soil layers according to the actual root distribution using an improved method (the logic was changed since Biome-BGCMuSo v4.0). From the plant-specific root parameters and the actual root weight Biome-BGCMuSo calculates the number of layers in which roots can be found together with the root mass distribution across the layers (Jarvis, 1989; Hidy et al., 2016). If there is not enough water in a given soil layer to fulfill the transpiration demand, the transpiration flux from that layer is limited, and below the wilting point (WP) it is set to zero. The sum of layer-specific transpiration fluxes across the root zone gives the actual transpiration flux. A detailed description of the algorithm can be found in Sect. S5 in the Supplement.

3.4 Effect of groundwater

Simulation of the groundwater effect was introduced in Biome-BGCMuSo v4.0 (Hidy et al., 2016), but the method has been significantly improved, and the new algorithm it is now available in Biome-BGCMuSo v6.2. In the recent model version there is an option to provide an additional input file with the daily values of the groundwater table depth (GWdepth in meters).

Groundwater may affect soil hydrological and plant physiological processes if the water table is closer to the root zone than the thickness of the capillary fringe (that is the region saturated from groundwater via the capillary effect). The thickness of the capillary fringe (CF in meters) is estimated using literature data and depends on the soil type (Johnson and Ettinger model; Tillman and Weaver, 2006). Groundwater table distance (GWdist in m) for a given layer is defined as the difference between GWdepth and the depth of the midpoint of the layer.

The layers completely below the groundwater table are assumed to be fully saturated. In the case of layers within the capillary fringe ($\text{GWdist}<\text{CF}$), the calculation of water balance changes: the FC rises, and thus the difference between saturation (SAT) and FC decreases, and the layer charges gradually until the increased FC value is reached. The FC rising effect of groundwater for the layers above the water table is calculated based on the ratio of the groundwater distance and the capillary fringe thickness, but only after the water contents of the layers below have reached their modified FC values. A detailed description of the groundwater effect can be found in Sect. S6 in the Supplement.

3.5 Soil moisture stress

In the original Biome-BGC model the effect of changing soil water content on photosynthesis and decomposition of soil organic matter is expressed in terms of soil water potential ($\mathrm{\Psi }$). Instead of $\mathrm{\Psi }$, the SWC is also widely used to calculate the limitation of stomatal conductance and decomposition. A practical advantage of using SWC as a factor in the stress function is that it is easier to measure in the field and the changes in the driving function are much smoother than in the case of $\mathrm{\Psi }$. The disadvantage is that SWC is not comparable among different soil types (in contrast to $\mathrm{\Psi }$).

The maximum of SWC is the saturation value; the minimum is the wilting point or the hygroscopic water depending on the type of simulated process. The novelty of Biome-BGCMuSo v6.2 is that the hygroscopic water, the wilting point, the field capacity, and the saturation values are calculated internally by the model based on the soil texture data, or they can be defined in the input file layer by layer.

In Biome-BGCMuSo v6.2 the so-called soil moisture stress index (SMSI) is calculated to represent overall soil stress conditions. SMSI is affected by the length of the drought event (SMSE: extent of soil stress) and the severity of the drought event (SMSL: length of soil stress), and it is aggravated by extreme temperature (extremT: effect of extreme heat). SMSI is equal to zero if no soil moisture limitation occurs and equal to 1 in the case of full soil moisture limitation. SMSI is used by the model for plant senescence calculations (presentation of plant-related processes is the subject of a forthcoming publication). The components of SMSI are explained in detail below.

4$$\text{SMSI}=\mathrm{1}-\text{SMSE}\cdot \text{SMSL}\cdot \text{extremT}$$

The magnitude of soil moisture stress (SMSE) is calculated layer by layer based on SWC. Regarding soil moisture stress two different processes are distinguished: drought (i.e., low SWC close to or below WP) and anoxic conditions (i.e., after large precipitation events or in the presence of a high groundwater table; Bond-Lamberty et al., 2007). An important novelty of Biome-BGCMuSo v6.2 is the soil curvature parameter ($q$), which is introduced to provide a mechanism for soil-texture-dependent drought stress as it can affect the shape of the soil stress function (which means a possibility for a nonlinear ramp function): 5$$\begin{array}{rl}{\text{SMSE}}^i& =\mathrm{0};\\ & {\text{SWC}}^i<{\text{SWC}}_{\text{WP}}^i\\ {\text{SMSE}}^i& ={\left({\displaystyle \frac{{\text{SWC}}^i-{\text{SWC}}_{\text{WP}}^i}{{\text{SWC}}_{\text{drought}}^i-{\text{SWC}}_{\text{WP}}^i}}\right)}^q;\\ & {\text{SWC}}_{\text{WP}}^i\le \text{SWC}<{\text{SWC}}_{\text{drought}}^i\\ {\text{SMSE}}^i& =\mathrm{1};\\ & {\text{SWC}}_{\text{drought}}^i\le \text{SWC}\le {\text{SWC}}_{\text{anoxic}}^i\\ {\text{SMSE}}^i& ={\displaystyle \frac{{\text{SWC}}_{\text{SAT}}^i-{\text{SWC}}^i}{{\text{SWC}}_{\text{SAT}}^i-{\text{SWC}}_{\text{anoxic}}^i}};\\ & {\text{SWC}}^i>{\text{SWC}}_{\text{anoxic}}^i\end{array}$$ where $q$ is the curvature of the soil stress function, and ${\text{SWC}}_{\text{drought}}^i$ and ${\text{SWC}}_{\text{anoxic}}^i$ are critical SWC values for calculating soil stress.

In order to make the SWC values comparable between different soil types, ${\text{SWC}}_{\text{drought}}^i$ and ${\text{SWC}}_{\text{anoxic}}^i$ can be set in normalized form (such as in Eq. 4) as part of the ecophysiological parameterization of the model. More details about the adjustment of the critical SWC values can be found in Hidy et al. (2021).

The layer-specific soil moisture stress extent values are summed across the root zone using the relative quantity of roots in the layers (${\text{RP}}^i$) as weighting factors to obtain the overall soil moisture stress extent (SMSE): $$\begin{array}{}\text{6}& {\displaystyle }\text{SMSE}=\sum _{i=\mathrm{0}}^{i=\text{nr}}{\text{SMSE}}^i\cdot {\text{RP}}^i,\text{7}& {\displaystyle }{\text{RP}}^i=\text{RD}{\displaystyle \frac{z^i}{\text{RL}}}\cdot e^{-\text{RD}\cdot ({\text{mid}}^i/\text{RL})},\end{array}$$ where nr is the number of soil layers in which roots can be found, RL is the actual length of roots, and RD is the rooting distribution parameter (ecophysiological parameter; see details in the user guide; Hidy et al., 2021). In the current model version SMSE can also affect the entire photosynthetic machinery by the introduction of an empirical parameter. This mechanism is responsible for accounting for the non-stomatal effect of drought on photosynthesis (details about this algorithm will be published in a separate paper). Since there is no mechanistic representation behind this empirical downregulation of photosynthesis, further tests are needed for the correct setting of this parameter preferentially using eddy covariance data.

The factor (SMSL) related to soil moisture stress length is the ratio of the critical soil moisture stress length (ecophysiological parameter) and the sum of the daily ($\mathrm{1}-\text{SMSE}$) values. This cumulated value restarts if SMSE is equal to 1 (no stress). Extreme heat (extremT) is also considered and is taken into account in the final stress function (see above) by using a ramp function. Its parameterization thus requires the setting of two critical temperature limits that define the ramp function (set by the ecophysiological parameterization; see Hidy et al., 2021). Its characteristic temperature values can be set by parameterization (ecophysiological input file).

4 Soil carbon and nitrogen cycles

4.1 Soil–litter module

We made substantial changes in the soil biogeochemistry module of the Biome-BGC model. Previous model versions already offered solutions for multilayer simulations (Hidy et al., 2012, 2016), but some pools still inherited the single-layer logic of the original model. In the new model version all relevant soil processes are separated layer by layer, which is a major step forward.

Instead of defining a single litter, soil organic carbon (SOC), and nitrogen pool, we implemented separate carbon and nitrogen pools for each soil layer in the form of soil organic matter (SOM) and litter in Biome-BGCMuSo v6.2. The changes in the mass of the carbon and nitrogen pools are calculated layer by layer. Mortality fluxes (whole plant mortality, senescence, litterfall) of aboveground plant material are transferred into the litter pools of the topsoil layers (0–10 cm, layers 1–2). Mortality fluxes of belowground plant material are transferred into the corresponding soil layers based on their location within the root zone. Due to ploughing and leaching, carbon and nitrogen can also be relocated to deeper layers. The plant material turning into the litter compartment is divided between the different types of litter pools (labile, unshielded cellulose, shielded cellulose, and lignin) according to the parameterization. Litter and soil decomposition fluxes (carbon and nitrogen fluxes from litter to soil pools) are calculated layer by layer, depending on the actual temperature and SWC of the corresponding layers. Vertical mixing of soil organic matter between the soil layers (e.g., bioturbation) is not implemented in the current model version.

Figure 3 shows the most important simulated soil and litter processes. N fixation (Nf) is the N input from the atmosphere to soil layers in the root zone by microorganisms. The user can set its annual value as an input parameter. N deposition (Nd) is the N input from the atmosphere to the topsoil layers (see below). The user can set its annual value as a site-specific parameter in the initialization input file. Nitrogen deposition can be provided by annually varying values as well. Plant uptake (PU) is the absorption of mineral N by plants from the soil layers in the root zone. Mineralization (MI) is the release of plant-available nitrogen (flux from soil organic matter to mineralized nitrogen). Immobilization (IM) is the consumption of inorganic nitrogen by microorganisms (flux from mineralized nitrogen to soil organic matter). Nitrification (NI) is the biological oxidation of ammonium to nitrate through nitrifying bacteria. Denitrification (DN) is a microbial process whereby nitrate (${{\mathrm{NO}}_{\mathrm{3}}}^{-}$) is reduced and converted to nitrogen gas (${\mathrm{N}}_{\mathrm{2}}$) through intermediate nitrogen oxide gases. Leaching ($L$) is the loss of water-soluble mineral nitrogen from the soil layers. If leaching occurs in the lowermost soil layer that means loss of N from the simulated system. Litterfall (LI) is the plant material transfer from plant compartments to litter. Decomposition is the C and N transfer from litter to soil pools and between soil pools. In the case of woody vegetation coarse woody debris (CWD) contains the woody plant material after litterfall before physical fragmentation. Litter has also four sub-pools based on composition: labile (L1), unshielded and shielded cellulose (L2, L3), and lignin (L4). Soil organic matter also has four sub-pools based on turnover rate: labile (S1), medium (S2), slow (S3), and passive (recalcitrant; S4). The soil-mineralized nitrogen pool contains the inorganic N forms of the soil: ammonium and nitrate.

Figure 3

Soil- and litter-related simulated carbon and nitrogen fluxes (arrows) as well as pools (rectangles) in Biome-BGCMuSo v6.2. HR: heterotrophic respiration, IM: immobilization, MI: mineralization, PU: plant uptake, LI: litterfall, NI: nitrification, $D$: decomposition ($D_{\mathrm{L}}$: decomposition of litter, $D_{\mathrm{S}}$: decomposition of SOM, $D_{\mathrm{C}}$: fragmentation of coarse woody debris), $L$: leaching, Nf: nitrogen fixation, Nd: nitrogen deposition, DN: denitrification. $L$ represents loss of C and N from the simulated system.

[Figure omitted. See PDF]

4.2 Decomposition

In the decomposition module (i.e., converging cascade scheme; Thornton, 1998) the fluxes between litter and soil pools are calculated layer by layer. The potential fluxes are modified in the case of N limitation when the potential gross immobilization is greater than the potential gross mineralization.

To explain the decomposition processes implemented in Biome-BGCMuSo v6.2 the main carbon and nitrogen pools as well as fluxes between litter and soil organic and inorganic (mineralized) matter are presented in Fig. 4.

Figure 4

Overview of the converging cascade model of litter and soil organic matter decomposition that is implemented in Biome-BGCMuSo v6.2. The notation rf represents the respiration fraction of the different transformation fluxes, and $\mathit{\tau }$ is the residence time (reciprocal of the rate constants, which is the turnover rate). IM/MI: immobilization and mineralization fluxes, HR: heterotrophic respiration. Note that both the respiration fraction and the turnover rate parameters can be adjusted through parameterization.

[Figure omitted. See PDF]

In the original Biome-BGC and in previous Biome-BGCMuSo versions the $\mathrm{C}:\mathrm{N}$ ratio (CN) of the soil pools was fixed in the model code. One of the new features in Biome-BGCMuSo v6.2 is that the CN of the passive soil pool (S4 in Fig. 4; recalcitrant soil organic matter) can be set by the user as a soil parameter. The CN of the other soil pools (labile, medium, and slow; S1, S2, and S3) is calculated based on the proportion of fixed CN values of the original Biome-BGC (${\text{CN}}_{\text{labile}}/{\text{CN}}_{\text{passive}}=\mathrm{1.2}$, ${\text{CN}}_{\text{medium}}/{\text{CN}}_{\text{passive}}=\mathrm{1.2}$, ${\text{CN}}_{\text{slow}}/{\text{CN}}_{\text{passive}}=\mathrm{1}$). Note that the CN of the donor and acceptor pools is used in decomposition calculations (see details in Sect. S7 in the Supplement), and as a result these parameters set the $\mathrm{C}:\mathrm{N}$ ratio of the soil pools. The donor and acceptor pools can be seen in Figs. 3 and 4.

For the calculation of nitrogen mineralization, first respiration cost (respiration fraction) is estimated. Mineralization is then a function of the remaining part of the pool and its $\mathrm{C}:\mathrm{N}$ ratio. The nitrogen mineralization fluxes of the SOM pools are functions of the potential rate constant (reciprocal of residence time) and the integrated response function that accounts for the impact of multiple environmental factors. The integrated response function of decomposition is a product of the response functions of depth, soil temperature, and SWC ($F_{\mathrm{r}}(d{)}_{\mathrm{D}}$, $F_{\mathrm{r}}(T{)}_{\mathrm{D}}$, $F_{\mathrm{r}}{\text{(SWC)}}_{\mathrm{D}}$; Fig. 5). Its detailed description can be found in Sect. S7. The dependence of the three different factors on depth, temperature, and SWC with default parameters is presented in Fig. 5.

Figure 5

The dependence of the individual factors that form the complex environmental response function of decomposition on depth ($F_{\mathrm{r}}(d{)}_{\mathrm{D}}$), temperature ($F_{\mathrm{r}}(T{)}_{\mathrm{D}}$), and SWC in the case of different soil types ($F_{\mathrm{r}}(\text{SWC}{)}_{\mathrm{D}}$). ED is the $e$-folding depth, which is one of the adjustable soil parameters of the model. For the definition of sand, silt, and clay see Fig. 1.

[Figure omitted. See PDF]

4.3 Soil nitrogen processes

In Biome-BGCMuSo v6.2 separate ammonium (sNH4) and nitrate (sNO3) soil pools are implemented instead of a general mineralized nitrogen pool. This was a necessary step for the realistic representation of many internal processes like plant nitrogen uptake, nitrification, denitrification, consideration of the effect of different mineral and organic fertilizers, and ${\mathrm{N}}_{\mathrm{2}}\mathrm{O}$ emission.

It is important to introduce the availability concept that Biome-BGCMuSo uses and is associated with the ammonium and nitrate pools. We use the logic proposed by Thomas et al. (2013), which means that the plant has access only to a part of the given inorganic nitrogen pool. The unavailable part is buffered as it is associated with soil aggregates and is unavailable for plant uptake. The available part of ammonium is calculated based on the ${\mathrm{NH}}_{\mathrm{4}}$ mobile proportion (that is a soil parameter set to 10 % according to Thomas et al., 2013; Hidy et al., 2021) and the actual pool. The available part of nitrate is assumed to be 100 %.

The amount of ammonium and nitrate is determined layer by layer controlled by input and output fluxes ($F$, $\mathrm{kg}\mathrm{N}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{d}}^{-\mathrm{1}}$) listed below: $$\begin{array}{cc}{\displaystyle }{F}_{\text{sNH4}}^i=& {\displaystyle }{\text{IN}}_{\text{sNH4}}^i-L_{\text{sNH4}}^i+L_{\text{sNH4}}^{i-\mathrm{1}}-{\text{PU}}_{\text{sNH4}}^i-{\text{IM}}_{\text{sNH4}}^i\\ \text{8}& {\displaystyle }& {\displaystyle }+{\text{MI}}_{\text{sNH4}}^i-{\text{NI}}_{\text{sNH4}}^i,{\displaystyle }{F}_{\text{sNO3}}^i=& {\displaystyle }{\text{IN}}_{\text{sNO3}}^i-L_{\text{sNO3}}^i+L_{\text{sNO3}}^{i-\mathrm{1}}-{\text{PU}}_{\text{sNO3}}^i-{\text{IM}}_{\text{sNO3}}^i\\ \text{9}& {\displaystyle }& {\displaystyle }+{\text{MI}}_{\text{sNO3}}^i-{\text{DN}}_{\text{sNO3}}^i,\end{array}$$ where ${\text{IN}}_{\text{sNH4}}^i$ and ${\text{IN}}_{\text{sNO3}}^i$ are the input fluxes to the ammonium and nitrate pools, respectively; $L_{\text{sNH4}}^i$, $L_{\text{sNH4}}^{i-\mathrm{1}}$, $L_{\text{sNO3}}^i$, and $L_{\text{sNO3}}^{i-\mathrm{1}}$ represent the amount of leached mineralized ammonium and nitrate from a layer ($i$) or from the upper layer ($i-\mathrm{1}$), respectively; ${\text{PU}}_{\text{sNH4}}^i$ and ${\text{PU}}_{\text{sNO3}}^i$ are the plant uptake fluxes of ammonium and nitrate, respectively; ${\text{IM}}_{\text{sNH4}}^i$ and ${\text{IM}}_{\text{sNO3}}^i$ are the immobilization fluxes of ammonium and nitrate, respectively; ${\text{MI}}_{\text{sNH4}}^i$ and ${\text{MI}}_{\text{sNO3}}^i$ are the mineralization fluxes of ammonium and nitrate, respectively; ${\text{NI}}_{\text{sNH4}}^i$ is the nitrification flux of ammonium; and ${\text{DN}}_{\text{sNO3}}^i$ is the denitrification flux of nitrate.

In the following subsections the different terms of the equations are described in detail.

Input to the sNH4 and sNO3 pools (IN in Eqs. 6 and 7)

According to the model logic N fixation occurs in the root zone layers. Its distribution between sNH4 and sNO3 pools is calculated based on their actual available proportion in the actual layer (${\text{NH4prop}}^i$):

10$${\text{NH4prop}}^i={\text{sNH4avail}}^i+{\text{sNO3avail}}^i,$$ where ${\text{sNH4avail}}^i$ and ${\text{sNO3avail}}^i$ are the available part of the sNH4 and sNO3 pools in the actual layer.

N-deposition-related nitrogen input is associated with the 0–10 cm soil layers assuming uniform distribution across layers 1–2 in the model, and the distribution between sNH4 and sNO3 pools is calculated based on the proportion of the ${\mathrm{NH}}_{\mathrm{4}}$ flux of the N deposition soil parameter (Hidy et al., 2021).

Organic and inorganic fertilization is also an optional nitrogen input. The amount and composition (${{\mathrm{NH}}_{\mathrm{4}}}^{+}$ and ${{\mathrm{NO}}_{\mathrm{3}}}^{-}$ content) can be set in the fertilization input file.

Leaching – downward movement of mineralized N ($L$ in Eqs. 6 and 7)

The amount of leached mineralized N (mobile part of the given N pool) from a layer is directly proportional to the amount of drainage and the available part of the sNH4 and sNO3 pools. Leaching from the layer above is a net gain, while leaching from actual layer is a net loss for the actual layer. Leaching is described in Sect. 4.5.

Plant uptake by roots (PU in Eqs. 6 and 7)

N uptake required for plant growth is estimated in the photosynthesis calculations, and the amount is distributed across the layers in the root zone. The partition of the N uptake between sNH4 and sNO3 pools is calculated based on their actual available proportion in each layer.

Mineralization and immobilization (MI and IM Eqs. 6 and 7)

Mineralization and immobilization calculations are detailed in Sect. 4.2. The distribution of these N fluxes between sNH4 and sNO3 pools is calculated based on their actual available proportion in each layer.

Nitrification (NI Eqs. 6 and 7)

Nitrification is a function of the soil ammonium content, the net mineralization, and the response functions of temperature, soil pH, and SWC ($F_{\mathrm{r}}(\text{pH}{)}_{\text{NI}}$, $F_{\mathrm{r}}(T{)}_{\text{NI}}$, and $F_{\mathrm{r}}(\text{SWC}{)}_{\text{NI}}$, respectively) based on the method of Parton et al. (2001) and Thomas et al. (2013). Its detailed mathematical description can be found in Sect. S8 in the Supplement. The response functions with proposed parameters are shown in Fig. 6.

Figure 6

The dependence of the individual factors of the environmental response function of nitrification on soil pH ($F_{\mathrm{r}}(\text{pH}{)}_{\text{NI}}$), temperature ($F_{\mathrm{r}}(T{)}_{\text{NI}}$), and SWC $F_{\mathrm{r}}(\text{SWC}{)}_{\text{NI}}$ in the case of different soil types. The pH and temperature response functions are independent of the soil texture.

[Figure omitted. See PDF]

Denitrification (DN Eqs. 6 and 7)

Denitrification flux is estimated with a simple formula (Thomas et al., 2013): 11$${\text{DN}}^i=\text{DNcoeff}\cdot {\text{SOMresp}}^i\cdot {\text{sNO3avail}}^i\cdot {\text{WFPS}}^i,$$ where DN of the actual layer is the product of the available nitrate content (sNO3avail, $\mathrm{kg}\mathrm{N}{\mathrm{m}}^{-\mathrm{2}}$), ${\text{SOMresp}}^i$ ($\mathrm{g}\mathrm{C}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{d}}^{-\mathrm{1}}$) is the SOM-decomposition-related respiration cost, ${\text{WFPS}}^i$ is the water-filled pore space, and DNcoeff is the soil-respiration-related denitrification rate ($\mathrm{g}{\mathrm{C}}^{-\mathrm{1}}$), which is an input soil parameter (Hidy et al., 2021). The unitless water-filled pore space is the ratio of the actual and saturated SWC. SOM-decomposition-associated respiration is the sum of the heterotrophic respiration fluxes of the four soil compartments (S1–S4, Fig. 4).

4.4

${\mathrm{N}}_{\mathrm{2}}\mathrm{O}$ emission and N emission

During both nitrification and denitrification ${\mathrm{N}}_{\mathrm{2}}\mathrm{O}$ emission occurs, which (added to the ${\mathrm{N}}_{\mathrm{2}}\mathrm{O}$ flux originating from grazing processes if applicable) contributes to the total ${\mathrm{N}}_{\mathrm{2}}\mathrm{O}$ emission of the examined ecosystem.

In Biome-BGCMuSo v6.2 a fixed part (set by the coefficient of the ${\mathrm{N}}_{\mathrm{2}}\mathrm{O}$ emission of the nitrification input soil parameter; Hidy et al., 2021) of nitrification flux is lost as ${\mathrm{N}}_{\mathrm{2}}\mathrm{O}$ and not converted to ${\mathrm{NO}}_{\mathrm{3}}$.

During denitrification, nitrate is transformed into ${\mathrm{N}}_{\mathrm{2}}$ and ${\mathrm{N}}_{\mathrm{2}}\mathrm{O}$ gas depending on the environmental conditions: ${\mathrm{NO}}_{\mathrm{3}}$ availability, total soil respiration (proxy for microbial activity), SWC and pH. The denitrification-related ${\mathrm{N}}_{\mathrm{2}}/{\mathrm{N}}_{\mathrm{2}}\mathrm{O}$ ratio input soil parameter is used to represent the effect of the soil type on the ${\mathrm{N}}_{\mathrm{2}}/{\mathrm{N}}_{\mathrm{2}}\mathrm{O}$ ratio (del Grosso et al., 2000; Hidy et al., 2021). Detailed mathematical description of the algorithm can be found in Sect. S9 in the Supplement.

4.5 Leaching of dissolved matter

Leaching of nitrate, ammonium, and dissolved organic carbon and nitrogen (DOC and DON) content from the actual layer is calculated as the product of the concentration of the dissolved component in the soil water and the amount of water (drainage plus diffusion) leaving the given layer either downward or upward. The dissolved component (concentration) of organic carbon is calculated from the SOC pool contents and the corresponding fraction of the dissolved part of SOC soil parameters. The dissolved component of the organic nitrogen content of the given soil pool is calculated from the carbon content and the corresponding $\mathrm{C}:\mathrm{N}$ ratio. The downward leaching is net loss from the actual layer and net gain for the next layer below; the upward flux is net loss for the actual layer and net gain for the next layer up. The downward leaching of the bottom active layer (ninth) is net loss for the system. The upward movement of dissolved substance from the passive (10th) layer is net gain for the system.

5 Case studies

5.1 Evaluation of soil hydrological simulation

In order to evaluate the functioning of the new model version (and to compare simulation results made by the current and previously published model version), a case study is presented regarding soil water content and soil evaporation simulations. The results of a bare soil simulation (i.e., no plant is assumed to be present) are compared to observation data from a weighing lysimeter station installed at Martonvásár, Hungary ($\mathrm{47}{}^{\circ }{\mathrm{18}}^{\prime }{\mathrm{57.6}}^{\prime \prime }$ N, $\mathrm{18}{}^{\circ }{\mathrm{47}}^{\prime }{\mathrm{25.6}}^{\prime \prime }$ E), in 2017. The climate of the area is continental with a 30-year average temperature of 11.0 ${}^{\circ }\mathrm{C}$ ($-$1 ${}^{\circ }\mathrm{C}$ in January and 21.2 ${}^{\circ }\mathrm{C}$ in July) and annual rainfall of 548 mm based on data from the on-site weather station.

The station consists of 12 scientific lysimeter columns 2 m deep with 1 m diameter (Meter Group Inc., USA) and with soil temperature, SWC, and soil water potential sensors installed at 5, 10, 30, 50, 70, 100, and 150 cm depth. Observation data for 2020 from six columns without vegetation cover (i.e., bare soil) were used to validate the model.

Raw lysimeter observation data were processed using standard methods. Bare soil evaporation values were derived based on changes in the mass of the soil columns also considering the mass change in the drainage water. Additionally, experience has shown that wind speed is related to the high-frequency mass change in the soil column mass. To reduce noise, five-point (5 min) moving averages were used based on Marek et al. (2014). After quality control of the data, the corrected and smoothed lysimeter mass values were used for the calculations. SWC observations were averaged to daily resolution to match the time step of the model.

Observed local meteorology was used to drive the models for the year 2020. Soil physical model input parameters (field capacity, wilting point, bulk density, etc.) were determined in the laboratory using 100 ${\mathrm{cm}}^{\mathrm{3}}$ undisturbed soil samples taken from various depths during the installation of the lysimeter station. Regarding other soil parameters the proposed values were used. A detailed description of the input soil parameters and their proposed values is presented in the user guide (Hidy et al., 2021).

In Fig. 7 the simulated and observed time series of soil evaporation are presented for Martonvásár for 2020. The figure shows that the soil evaporation simulation by v6.2 is more realistic than by v4.0. Biome-BGCMuSo v4.0 provides very low values during summer on some days, which is not in accordance with the observations. Biome-BGCMuSo v6.2 provides more realistic values during this time period.

Figure 7

The simulated (blue line: v4.0; red line: v6.2) and observed (grey dots) daily soil evaporation values at Martonvásár during 2020. Vertical grey lines associated with the observations represent the standard deviation of the observations from six lysimeter columns. The improved model clearly outperforms the earlier version.

[Figure omitted. See PDF]

In Fig. 8 the simulated and observed SWCs at 10 cm depth are presented with the daily sum of precipitation representing the bare soil simulation in Martonvásár for 2020. The soil water balance simulation seems to be realistic using v6.2, since the annual course captures the low and high end of the observed values. In contrast, Biome-BGCMuSo v4.0 underestimates the range of SWC and provides overestimations during the growing season (from spring to autumn). With a couple of exceptions, the simulated values using v6.2 fall into the uncertainty range of the measured values defined by the standard deviation of the six parallel measurements. This is not the case for the simulations with the 4.0 version.

Figure 8

The simulated (blue line: v4.0; red line: v6.2) and observed (grey dots) soil water content values at 10 cm depth (right $y$ axis) with the daily sums of precipitation (left axis; black columns) during 2020 at the Martonvásár lysimeter station. Vertical grey lines associated with the observations represent $\pm$ 1 standard deviation around the observations. The improved model clearly outperforms the earlier version in simulating soil water balance.

[Figure omitted. See PDF]

Model performance was evaluated by quantitative measures such as the coefficient of determination ($R^{\mathrm{2}}$), mean absolute error (MAE), and mean signed error (MSE). In Fig. 9 the comparison of the simulated and observed daily evaporation is presented. Based on the performance indicators it is obvious that the simulation with the new model version (v6.2) is much closer to observations than the old version (v4.0). Biome-BGCMuSo v6.2 slightly underestimated the observations.

Figure 9

Comparison of the simulated (a: v4.0; b: v6.2) and observed daily soil evaporation representing the means of measured data obtained from six weighing lysimeter columns with bare soil at Martonvásár in 2020. $R^{\mathrm{2}}$, MAE, and MSE denote the coefficient of determination, mean absolute error, and mean signed error (bias) of the simulated values, respectively.

[Figure omitted. See PDF]

In Fig. 10 the comparison of the simulated and observed daily SWC from the lysimeter experiment is presented. Based on the model evaluation it seems that the simulation with the new model version is much closer to observations than with old version (4.0). The results obtained from v4.0 are consistent with earlier findings about the incorrect representation of the annual SWC cycle (Hidy et al., 2016; Sándor et al., 2017).

Figure 10

Comparison of the simulated (a: v4.0; b: v6.2) and observed daily SWC representing the means of measured data obtained from six weighing lysimeter columns with bare soil at Martonvásár in 2020. $R^{\mathrm{2}}$, MAE, and MSE denote the coefficient of determination, mean absolute error, and mean signed error (bias) of the simulated values, respectively.

[Figure omitted. See PDF]

A thorough validation of the improved model based on observed SWC and ET datasets from eddy covariance sites is planned to be published in an upcoming paper about the plant-related improvements.

5.2 Evaluation of the soil nitrogen balance module and the simulated soil respiration

Soil-related developments were evaluated with a case study focusing on topsoil nitrate content, soil ${\mathrm{N}}_{\mathrm{2}}\mathrm{O}$ efflux, and soil respiration.

Experimental data were collected in a long-term fertilization experiment that was set up in 1959 at Martonvásár, Hungary ($\mathrm{47}{}^{\circ }{\mathrm{18}}^{\prime }{\mathrm{41}}^{\prime \prime }$ N, $\mathrm{18}{}^{\circ }{\mathrm{46}}^{\prime }{\mathrm{50}}^{\prime \prime }$ E). According to the FAO-WRB classification system (IUSS Working Group, 2015), the soil is a Haplic Chernozem, with 51.4 % sand, 34 % silt, and 14.6 % clay content. Bulk density is 1.47 $\mathrm{g}{\mathrm{cm}}^{-\mathrm{3}}$, pH(${\mathrm{H}}_{\mathrm{2}}\mathrm{O}$) is 7.3, ${\mathrm{CaCO}}_{\mathrm{3}}$ content is 0 %–1 %, and the mean soil organic matter content in the topsoil is 3.2 %. The plant-available macronutrient supply in the soil was poor for phosphorus and medium to good for potassium based on the ProPlanta plant nutrition advisory system (Fodor et al., 2011). In the long-term fertilization experiment the treatments were arranged in a random block design with $\mathrm{6}\mathrm{m}\times \mathrm{8}$ m plots in four replicates. Eight different treatments were set up: control (zero artificial fertilizer applied), only N, only P, NPK – with farmyard manure, absolute control (zero nutrient supply), only N, only P, and NPK – without farmyard manure. The crop in the 4-year fertilizer cycles was maize in the first and second years and winter wheat in the third and fourth years. Here we used data from the absolute control and from the farmyard manure (FYM) treatments only. FYM was applied once every 4 years at a rate of 35 $\mathrm{t}{\mathrm{ha}}^{-\mathrm{1}}$ in autumn.

Topsoil nitrate content was measured during 2017, 2018, and 2020 on a few occasions by wet chemical reactions using a stream distillation method after KCl extraction of soil samples (Hungarian Standards Institution MSZ 20135:1999; Akhtar et al., 2011).

Dynamic-chamber-based soil ${\mathrm{N}}_{\mathrm{2}}\mathrm{O}$ efflux observations were available from 2020 and 2021. The ${\mathrm{N}}_{\mathrm{2}}\mathrm{O}$ efflux measurements with a gas incubation time of 10 min were performed by using a Picarro G2508 (Picarro, USA) cavity ring-down spectrometer (Christiansen et al., 2015; Zhen et al., 2021). The cylinder-shaped transparent gas incubation chamber was 16.5 cm in diameter, and its height was 30 cm. ${\mathrm{N}}_{\mathrm{2}}\mathrm{O}$ flux measurements were executed in six replicates per treatment on a biweekly (2020) and precipitation-event-related (2021) basis. Soil respiration was measured with the same Picarro gas analyzer. Sampling numbers and points were identical to those of the ${\mathrm{N}}_{\mathrm{2}}\mathrm{O}$ efflux measurements. ${\mathrm{CO}}_{\mathrm{2}}$ and ${\mathrm{N}}_{\mathrm{2}}\mathrm{O}$ effluxes were calculated by a linear equation (Widen and Lindroth, 2003) based on gas concentration data.

For the simulations we used a maize parameterization from previous studies (Fodor et al., 2021). A winter wheat parameterization was constructed based on a country-scale optimization using the AgroMo software package (https://github.com/hollorol/AgroMo, last access: November 2021) and the NUTS 3 level long-term (1991–2020) yield database of the Hungarian Central Statistical Office. For nitrogen-cycle-related parameters we mainly used the values presented in the user guide (Hidy et al., 2021). Two soil parameters were adjusted (coefficient of ${\mathrm{N}}_{\mathrm{2}}\mathrm{O}$ emission for nitrification and ${\mathrm{N}}_{\mathrm{2}}/{\mathrm{N}}_{\mathrm{2}}\mathrm{O}$ ratio multiplier for denitrification-related N gas flux; Del Grosso et al., 2000; Parton et al., 2001; Thomas et al., 2013; Hidy et al., 2021) to match the simulated ${\mathrm{N}}_{\mathrm{2}}\mathrm{O}$ efflux to the observations.

Figure 11 shows the comparison of the simulated and observed ${\mathrm{NO}}_{\mathrm{3}}$ content of the topsoil for the two selected treatments. The results indicate that the model underestimates the topsoil ${\mathrm{NO}}_{\mathrm{3}}$ content in the case of both C and FYM (bias is $-$2.3 and $-$2.4 ppm, respectively) treatments, but the simulation error is in an acceptable range (NRMSE is 45.5 % and 37.6 % for C and FYM, respectively).

Figure 11

Comparison of the simulated and observed ${\mathrm{NO}}_{\mathrm{3}}$ content of the topsoil for the absolute control (C; a) and for the farmyard manure (FYM; b) treatment between May 2017 and November 2021 at Martonvásár.

[Figure omitted. See PDF]

Figure 12 shows the comparison of the observed and simulated ${\mathrm{N}}_{\mathrm{2}}\mathrm{O}$ efflux for the 2020–2021 time period. Measurement uncertainties are also indicated. Note that the uncertainty of the observations (e.g., due to spatial heterogeneity and sample number, soil disturbance, improper chamber design, methods of sample analysis) is remarkable due to known features of the chamber technique (Chadwick et al., 2014; Pavelka et al., 2018). The model captured more of the magnitude of ${\mathrm{N}}_{\mathrm{2}}\mathrm{O}$ efflux peaks and less of their timing. Overall the model underestimated the observed values in both cases (bias is $-$0.13 and $-$0.1 $\mathrm{mg}\mathrm{N}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{d}}^{-\mathrm{1}}$ for C and FYM, respectively), with NRMSE of 32.4 % and 37.6 % for C and FYM, respectively.

Figure 12

Comparison of the simulated and observed soil ${\mathrm{N}}_{\mathrm{2}}\mathrm{O}$ efflux for two treatments: absolute control (C; a) and application of farmyard manure (FYM; b) between January 2020 and December 2021 at Martonvásár. Whiskers indicate the uncertainty ($\pm$ 1 standard deviation) of the measurements.

[Figure omitted. See PDF]

Figure 13 presents the comparison of the observed and simulated soil respiration for the same time period as for the soil ${\mathrm{N}}_{\mathrm{2}}\mathrm{O}$ efflux. Observation uncertainty is indicated that represents 1 standard deviation of the replicates. The model mostly captured the magnitude and variability of soil respiration flux. The model overestimated the observed values in both cases with bias of 0.17 and 0.04 $\mathrm{g}\mathrm{C}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{d}}^{-\mathrm{1}}$ for C and FYM, respectively. The NRMSE is 34.1 % and 40.1 % for C and FYM, respectively. It is interesting to note that the observations and the simulations are particularly different after harvest time in both years (i.e., beginning of October). The simulated respiration has peaks corresponding to harvest when the amount of litter sharply increases due to the by-products left behind (decomposition of residues left on the site after harvest is accounted for in the model). The chamber-based ${\mathrm{CO}}_{\mathrm{2}}$ efflux data did not really show similar peaks, likely because of methodological issues (litter is removed from the soil surface before placing of the chambers).

Figure 13

Comparison of the simulated and observed soil respiration flux for two treatments: absolute control (C; a) and application of farmyard manure (FYM; b) between January 2020 and December 2021 at Martonvásár. Whiskers indicate the uncertainty ($\pm$ 1 standard deviation) of the measurements.

[Figure omitted. See PDF]

Overall, the model provided nitrate content, ${\mathrm{N}}_{\mathrm{2}}\mathrm{O}$ emission, and soil respiration simulation results that are consistent with the observations. The model was capable of estimating the observed values with comparable efficiency reported in similar studies (Gabrielle et al., 2002; Andrews et al., 2020).

5.3 Sensitivity analysis and optimization of the soil biogeochemistry scheme

Here we present another case study that provides insight into the functioning of the converging cascade (decomposition) scheme that is implemented in Biome-BGCMuSo v6.2. A large-scale in silico experiment is also presented, the main aim of which was to perform model self-initialization (i.e., spin-up) at country scale (for the entire area of Hungary) with the resulting soil organic matter pools expected to be consistent with the observations.

The observation-based, gridded, multilayer SOC database of Hungary (DOSoReMI database; Pásztor et al., 2020; see Figs. S3–S4 in the Supplement) and the FORESEE meteorological database (Kern et al., 2016) were used for the sensitivity analysis of the soil scheme as well as for optimizing the most important soil parameters when the model was calibrated to the observation-based SOC values. As a first step, the area of the country was divided into 1104 grid cells (regular grid with 0.1${}^{\circ }$ by 0.1${}^{\circ }$ resolution, corresponding to approximately 10 km resolution). The 1104 grid cells of the DOSoReMI database were grouped based on their dominant land use type (cropland, grassland, or forest based on the CORINE-2018 database; EEA, 2021; Figs. S1–S2 in the Supplement) as well as the soil texture class (12 classes according to the USDA system; USDA, 1987) and SOC content (high and low; high is greater than the group mean, while low is less than the mean) of the topsoil (0–30 cm layer). As some of the theoretically possible 72 groups had no members (e.g., there is no soil in Hungary with sandy–clay texture) soils of the 1104 grid cells were categorized into 51 groups. For each group one single cell (so-called representative cell) was selected based on the topsoil SOC content. The representative cell was the one with the smallest absolute deviation from the group mean SOC content (land use maps for Hungary are presented in Sect. S10 in the Supplement: Figs. S1–S2).

The grassland ecophysiological parameterization without management was used in the spin-up phase to initialize SOC pools for croplands. For the transient phase, the cropland parameterization was used with fertilization, ploughing, planting, and harvest settings. In the case of grasslands, a grassland parameterization was used during both the spin-up and transient phases, and in the transient phase mowing was assumed once a year. In the case of forests a generic deciduous broadleaf forest parameterization was used for both the spin-up and transient phases with thinning in the latter phase. For our parameterization presented in the MS the generic, plant-functional-type-specific ecophysiological parameter sets published by White et al. (2000) served as starting points. These Biome-BGCMuSo-specific parameter sets are available at the website of the model

http://nimbus.elte.hu/bbgc/files/generic_EPC_set_6.2.zip (last access: November 2021).

.

Soil parameters in Biome-BGCMuSo v6.2 were classified into six groups: (1) 4 generic soil parameters; (2) 24 parameters related to decomposition, nitrification, and denitrification; (3) 14 rate scalars for the converging (decomposition) cascade scheme; (4) 19 soil-moisture-related parameters; (5) 7 methane-related parameters; and (6) 11 soil composition and characteristic values (can be set layer by layer). A detailed description and proposed value of each soil parameters can be found in the user guide (Hidy et al., 2021).

As methane simulation was not the subject of the present case study we neglected the related parameters. Regarding the soil composition and characteristic values we used the DOSoReMI database (Pásztor et al., 2020). From the remaining 61 parameters, soil depth, runoff curve number, and the three soil-moisture-related parameters (tipping bucket method) were not included in the analysis. The groundwater module was switched off in this case (no groundwater is assumed) and the related parameters were not studied. The remaining 53 parameters were used in the sensitivity analysis and are listed in Table 1.

Table 1

Soil parameters of Biome-BGCMuSo v6.2 (referring to SOC simulation) that were used during the sensitivity analysis. The “Value” column shows the originally proposed values (Hidy et al., 2021). See Fig. 4 for an explanation of the compartment names. The parameters that were included in the second phase of the sensitivity analysis are marked with italics (see text).

Group	Parameter name	Abbreviation	Value
Generic soil	$\mathrm{C}:\mathrm{N}$ ratio of stable soil pool (soil4)	soil4CN	12
parameters	${\mathrm{NH}}_{\mathrm{4}}$ mobile proportion	amMP	0.1
	aerodynamic resistance	potRair	107
Decomposition,	parameter 1 for temperature response function of decomp.	Tp1decomp	1.75
nitrification,	parameter 2 for temperature response function of decomp.	Tp2decomp	17
denitrification	parameter 3 for temperature response function of decomp.	Tp3decomp	2.6
parameters	parameter 4 for temperature response function of decomp.	Tp4decomp	40
	minimum T for decomposition and nitrification	Tp5decomp	$-\mathit{5}$
	$e$-folding depth of decomposition rate's depth scalar	EFD	10
	net mineralization proportion of nitrification	NITRnetMINER	0.2
	maximum nitrification rate	NITRmaxRATE	0.1
	coefficient of ${\mathrm{N}}_{\mathrm{2}}\mathrm{O}$ emission of nitrification	NITRratioN2O	0.02
	parameter 1 for pH response function of nitrification	pHp1nitrif	0.15
	parameter 2 for pH response function of nitrification	pHp2nitrif	1
	parameter 3 for pH response function of nitrification	pHp3nitrif	5.2
	parameter 4 for pH response function of nitrification	pHp4nitrif	0.55
	parameter 1 for Tsoil response function of nitrification	Tp1nitrif	1
	parameter 2 for Tsoil response function of nitrification	Tp2nitrif	12
	parameter 3 for Tsoil response function of nitrification	Tp3nitrif	2.6
	parameter 4 for Tsoil response function of nitrification	Tp4nitrif	2.6
	minimum WFPS for scalar of nitrification calculation	minWFPS	0.1
	lower optimum WFPS for scalar of nitrification	opt1WFPS	0.45
	higher optimum WFPS for scalar of nitrification	opt2WFPS	0.55
	minimum value for saturated WFPS scalar of nitrification	minWFPSscalar	0.2
	soil-respiration-related denitrification rate	DENITcoeff	0.05
	denitrification-related ${\mathrm{N}}_{\mathrm{2}}/{\mathrm{N}}_{\mathrm{2}}\mathrm{O}$ ratio multiplier	DNratioN2O	2
	critical WFPS value for denitrification	critWFPSdenitr	0.50
Rate scalars	respiration fractions for fluxes between compartments (l1s1)	RFl1s1	0.39
	respiration fractions for fluxes between compartments (l2s2)	RFl2s2	0.55
	respiration fractions for fluxes between compartments (l4s3)	RFl4s3	0.29
	respiration fractions for fluxes between compartments (s1s2)	RFs1s2	0.28
	respiration fractions for fluxes between compartments (s2s3)	RFs2s3	0.46
	respiration fractions for fluxes between compartments (s3s4)	RFs3s4	0.55
	potential rate constant of labile litter pool	RCS1	0.7
	potential rate constant of cellulose litter pool	RCS2	0.07
	potential rate constant of lignin litter pool	RCS3	0.014
	potential rate constant of fast microbial recycling pool	RCS4	0.07
	potential rate constant of medium microbial recycling pool	RCS5	0.014
	potential rate constant of slow microbial recycling pool	RCS6	0.0014
	potential rate constant of recalcitrant SOM (humus) pool	RCS7	0.0001
	potential rate constant of physical fragmentation of wood	RCS8	0.001
	maximum height of pond water	MP	5
	curvature of soil stress function	q	1
	fraction of dissolved part of S1 organic matter	fD1	0.005
	fraction of dissolved part of S2 organic matter	fD2	0.004
	fraction of dissolved part of S3 organic matter	fD3	0.003
	fraction of dissolved part of S4 organic matter	fD4	0.002
	mulch parameter: critical amount	CAmulch	1
	parameter 1 for mulch function	p1mulch	100
	parameter 2 for mulch function	p2mulch	0.75
	parameter 3 for mulch function	p3mulch	0.75
	mulch parameter: evaporation reduction	ERmulch	0.5

As a first step sensitivity analysis was carried out for the selected 53 soil parameters by running the Biome-BGCMuSo v6.2 model in spin-up mode until a quasi-equilibrium in the total SOC was reached (that is the usual logic of the spin-up run). The model was run for each representative cell 2000 times with varying model parameters using the Monte Carlo method. All model parameters were varied randomly within the $\pm \mathrm{10}\mathit{\%}$ range of their initial values that were inherited from the Biome-BGC model or were set according to the literature. The least square linearization (LSL) method (Verbeeck et al., 2006) was used for dividing output uncertainty into its input-parameter-related variability. As a result of the LSL method, the total variance of the model output and the sensitivity coefficient of each parameter can be determined. Sensitivity coefficients show the percent of total variance for which the given parameter is responsible.

Figure 14

The sensitivity coefficients of the soil parameters as the result of the sensitivity analysis. Black columns refer to the crop, light grey to the grass, and dark grey to the forest simulations. The sensitivity coefficients are calculated as the mean pixel-level sensitivity coefficient for the given land use type. The horizontal line indicates the 5 % threshold that was used to select the final parameter set for optimization.

[Figure omitted. See PDF]

In order to simplify the workflow and decrease the degree of freedom another sensitivity analysis was performed. In this second step, the sensitive parameters (sensitivity coefficient $>\mathrm{1}$ % for at least one land use type; a total of 18 parameters) were used in the following sensitivity analysis with 6000 iteration steps. These 18 parameters are marked with italics in Table 1.

Figure 14 shows the summary of the second sensitivity analysis wherein the overall importance of the parameters is calculated as the mean of all selected pixels in a given land use category. It can be seen in Fig. 4 that from the 18 parameters (selected during the first phase), the soil carbon ratio of the recalcitrant pool (soil4CN), the temperature dependence parameters of the decomposition function (Tp1decomp, Tp2decomp, Tp3decomp, Tp4_decomp), the respiration fraction of the S2–S3 and S3–S4 decomposition process (RFs2s3 and RFs3s4), the curvature of the soil stress function ($q_{\text{soilstress}}$), and the fraction of the dissolved part of S4 organic matter (fD4) are the most important for all land use types. Among the other parameters the critical WFPS of denitrification (critWFPSdentir) for grasslands has a remarkably high sensitivity (greater than 35 %). It means that in the case of grasslands the nitrogen availability seems to be an important limitation on primary production, probably because there are only natural sources of nitrogen (no fertilization is assumed here) and the rooting zone is shallower than in the case of forest, which involves limited mineralized N access. Thus, in the case of higher values of critical WFPS of denitrification, the simulated production of grassland (and therefore the final SOC) seems to be significantly underestimated.

The 10 selected soil-biogeochemistry-related parameters were optimized for each of the 51 groups separately using maximum likelihood estimation. For each group, the parameter set providing the smallest deviation between the simulated and observed values of the weighted average SOC content (weight factor of 5 for 0–30 cm and weight factor of 1 for 30–60 cm soil layers) was considered to be the final (optimized) model parameter set.

The differences of the simulated and observed SOC content for the 0–30 cm layer (SOC0-30) using the initial (Table 1) and final soil parameters (not shown here) are presented in Fig. 15. In the upper panel the signed relative error of SOC0-30 simulation before optimization can be seen, while in the lower panel the signed relative error of SOC0-30 simulation after optimization can be seen. It is clearly visible that because of optimization the overestimation of the SOC0-30 simulation significantly decreased.

Figure 15

Differences (expressed as signed relative error, %) between the simulated and observed SOC data for the 0–30 layer (SOC0-30) using the initial (a) and optimized (b) soil parameters. Visual comparison of the maps reveals the success of the optimization in terms of capturing the overall SOC for the whole country area.

[Figure omitted. See PDF]

Table 2

Comparison of the internal processes simulated in Biome-BGC 4.1.1, Biome-BGCMuSo v4.0, and Biome-BGCMuSo v6.2.

Routine	original Biome-BGC	Biome-BGCMuSo v4.0	Biome-BGCMuSo v6.2
Runoff	no	based on simple, empirical formulation	distinguishing Hortonian and Dunne runoff
Pond water	no	simple solution	development of pond water formation (based on infiltration capacity)
Soil evaporation	Based on Penman–Monteith equation Calculation of the actual evaporation from the potential evaporation and the square root of time elapsed since the last precipitation	Based on Penman–Monteith equation Calculation of the actual evaporation from the potential evaporation and the square root of time elapsed since the last precipitation	Based on Penman–Monteith equation Parameterization possibility of actual aerodynamic resistance Introduction of an upper limit for daily potential evaporation that is determined by the available energy Calculation of the actual evaporation is based on the method of Ritchie (1981) Simulation of the reducing effect of surface residue or mulch cover on bare soil evaporation
Transpiration	Transpiration from one-layer bucket soil	Transpiration from seven-layer soil based on soil stress	Transpiration from 10-layer soil based on available water
Groundwater	no	Simple groundwater simulation	Improvement of the simulation of groundwater effect (using capillary fringe) Introduction of two different methods
Soil moisture stress	no	Relative SWC data are used to calculate soil water stress The hygroscopic water, the wilting point, the field capacity, and the saturation values of the soil layers can be defined in the input file layer by layer The soil moisture stress index is affected by the length and the day since the drought event	The hygroscopic water, the wilting point, the field capacity, and the saturation values of the soil layers can be defined in the input file layer by layer The soil moisture stress index is affected by the length and the severity of the drought event, aggravated by the extreme temperature Introduction of the soil curvature parameters to provide mechanism for soil-texture-dependent drought stress since it can affect the shape of the soil stress function Normalized SWC data are used to calculate soil moisture stress index
Organic carbon and nitrogen	One-layer soil module with one organic carbon and nitrogen pool	Multilayered soil module without soil carbon and nitrogen profile	Instead of defining a single litter, soil organic carbon, and nitrogen pool, separate carbon and nitrogen pools for each soil layer in the form of soil organic matter and litter were implemented Separation of aboveground and belowground litter pools Litter and soil decomposition fluxes (carbon and nitrogen fluxes from litter to soil pools) are calculated layer by layer, depending on the actual temperature and SWC of the corresponding layers Leaching of dissolved organic carbon and nitrogen
Inorganic nitrogen	One-layer soil module with one mineralized N pool	Multilayer soil module with an empirical inorganic N profile (no layer-by-layer calculations, only estimation of the sub-pools in the different soil layer based on the root length proportion)	Separation of ammonium (sNH4) and nitrate (sNO3) soil pools instead of a general mineralized nitrogen pool Nitrification fluxes are calculated layer by layer, depending on the actual pH, temperature, and SWC of the given layers Denitrification fluxes are calculated layer by layer, depending on the depth, actual temperature, and SWC of the given layers

We do not claim, of course, that the optimized parameters have universal value. Site history is neglected during the spin-up simulations, and we use many simplifications like disregarding land use change and present-day ecophysiological parameterization. In this sense, the optimized parameter set can best be considered a pragmatic solution to provide initial conditions (equilibrium SOC pools) for the model at country scale consistent with the observations.

6 Concluding remarks

In this paper, we presented a detailed description of Biome-BGCMuSo v6.2 terrestrial ecosystem model developments related to soil hydrology as well as the carbon and nitrogen budget. We mostly focused on changes relative to the previously published Biome-BGCMuSo v4.0 (Hidy et al., 2016), but our intention was also to provide a complete, stand-alone reference for the modeling community with mathematical equations (detailed in the Supplement). Table 2 summarizes the structural changes that we made during the developments starting from Biome-BGC v4.1.1 also including the previously published Biome-BGCMuSo v4.0 (Hidy et al., 2016).

Earlier model versions used a soil hydrology scheme based on the Richards equation, but the results were not satisfactory. Sándor et al. (2017) presented results from the first major grassland model intercomparison project (executed within the framework of FACCE MACSUR) wherein Biome-BGCMuSo v2.2 was used. That study demonstrated the problems associated with proper representation of soil water content, which was a common shortcoming of all included models. In the Hidy et al. (2016) paper, wherein the focus was on Biome-BGCMuSo v4.0, the SWC-related figures clearly indicated problems with the simulations compared to observations. The SWC amplitude was not captured well, which clearly influences drought stress, decomposition, and other SWC-driven processes like nitrification and denitrification. For the latter two processes this is especially critical as they are associated with contrasting SWC regimes (nitrification is aerobic, while denitrification is an anaerobic process). This is a good example of an erroneous internal process representation that may lead to improper results. Note that the functions currently used for nitrification and denitrification are also subject to uncertainty that needs to be addressed in the future (Heinen, 2006). Nevertheless, the presented model developments might contribute to more realistic soil process simulations and improved results.

The algorithm ensemble approach is already implemented in Biome-BGCMuSo. Algorithm ensemble means that the user has more than one option for the representation of some processes. Biome-BGCMuSo v6.2 has alternative phenology routines (Hidy et al., 2012) and two alternative methods for soil temperature (Hidy et al., 2016), soil hydrology (described in this study), photosynthesis, and soil moisture stress calculation. We plan to extend the algorithm ensemble by providing alternative decomposition schemes to the model. One possibility is the implementation of a CENTURY-like structure (Koven et al., 2013) that is a promising direction and might improve the quality of the equilibrium (spin-up) simulations as well as the simulated N mineralization related to SOM decomposition. Reported problems related to the rapid decomposition of litter in the current model structure (Bonan et al., 2013) need to be addressed in future model versions as well.

Plant growth- and allocation-related developments were not addressed in this study but of course have many inferences with the presented model logic (i.e., parameterization and related primary production define the amount and quality of litter). A forthcoming publication will provide a comprehensive overview of the plant growth- and senescence-related model modifications with elements from crop models also included.

Biome-BGCMuSo is still an open-source model that can be freely downloaded from its website with a detailed user guide and other supplementary files. We also encourage users to test the so-called RBBGCMuso package (available at GitHub) that has many advanced features to support model application and optimization. A graphical environment, called AgroMo (also available at GitHub: https://github.com/hollorol/AgroMo, last access: November 2021), was also developed around Biome-BGCMuSo to help users carry out simulations with either site-specific plot-scale data or with gridded databases representing large regions.

Code and data availability

The current version of Biome-BGCMuSo, together with sample input files and a detailed user guide, is available from the website of the model at http://nimbus.elte.hu/bbgc/download.html (last access: November 2021) under the GPL-2 license. Biome-BGCMuSo v6 is also available at GitHub: https://github.com/bpbond/Biome-BGC/tree/Biome-BGCMuSo_v6 (last access: November 2021). The exact version of the model (v6.2 alpha) used to produce the results in this paper is archived on Zenodo (10.5281/zenodo.5761202; Hidy and Zoltán, 2021). Experimental data and model parameterization used in the study are available from the corresponding author upon request.

The supplement related to this article is available online at: https://doi.org/10.5194/gmd-15-2157-2022-supplement.

Author contributions

DH developed Biome-BGCMuSo, maintained the source code, and executed the sample simulations. The study was conceived and designed by DH, ZB, and NF, with assistance from TÁ, LD, and RH. It was directed by DH and ZB. TÁ and LD contributed with model benchmarking. RH participated with the construction of a modeling framework for Biome-BGCMuSo. TF, DI, DZ, LP, MD, and ET contributed with experimental data. DH, ZB, NF, and KM prepared the paper and the Supplement with contributions from all co-authors. All authors reviewed and approved the present article and the Supplement.

Competing interests

The contact author has declared that neither they nor their co-authors have any competing interests.

Disclaimer

Publisher's note: Copernicus Publications remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Acknowledgements

We are grateful to Galina Churkina for reviewing the paper.

Financial support

The research was funded by the Széchenyi 2020 program, the European Regional Development Fund, and the Hungarian Government (grant no. GINOP-2.3.2-15-2016-00028). This research was supported by the NRDI Fund FK 20 (grant no. 134547) as well and also supported by the grant “Advanced research supporting the forestry and wood-processing sector's adaptation to global change and the 4th industrial revolution” (grant no. CZ.02.1.01/0.0/0.0/16_019/0000803) financed by OP RDE. Katarína Merganičová was also financed by the project “Scientific support of climate change adaptation in agriculture and mitigation of soil degradation” (grant no. ITMS2014+ 313011W580) supported by the Integrated Infrastructure Operational Programme funded by the ERDF.

Review statement

This paper was edited by Sam Rabin and reviewed by two anonymous referees.