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

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 m³ m⁻³ for SWC and from −0.68 to −0.2 mm d⁻¹ for soil evaporation). Simulations related to nitrogen balance and soil CO₂ 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 (N₂O) efflux and soil respiration simulations were also realistic, with overall correspondence with the observations (for the N₂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 CO₂ 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 CO₂ concentration, the integration of selected management modules, the implementation of new plant compartments (e.g., yield), the implementation of a C₄ photosynthesis routine, the implementation of photosynthesis and respiration acclimation of plants and temperature-dependent Q₁₀, 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 CO₂ 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 ^∘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 (cm d⁻¹), 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.

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, d⁻¹) 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_SAT: saturated hydraulic conductivity, cm d⁻¹; 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 (potR_air in s m⁻¹). potR_air was a fixed value in the previous model versions (107 s m⁻¹). 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 potR_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 (evap_limit) that is determined by the available energy (incident shortwave flux that reaches the soil surface):

$$\begin{array}{}\text{(1)}& {\text{evap}}_{\text{limit}}={\displaystyle \frac{\text{irad}\cdot \text{dayl}}{{\text{LH}}_{\text{vap}}}},\end{array}$$

where irad is the incident shortwave flux density (W m⁻²), dayl is the length of the day in seconds, and LH_vap is the latent heat of vaporization (the amount of energy that must be added to liquid to transform into gas; J kg⁻¹). 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, kg C m⁻²) with parameters p1_mulch, p2_mulch, and p3_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.

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 (GWdist<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 (Ψ). Instead of Ψ, 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 Ψ. The disadvantage is that SWC is not comparable among different soil types (in contrast to Ψ).

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.

$$\begin{array}{}\text{(4)}& \text{SMSI}=\mathrm{1}-\text{SMSE}\cdot \text{SMSL}\cdot \text{extremT}\end{array}$$

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):

$$\begin{array}{}\text{(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}\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 (RPⁱ) 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 (1−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 (N₂) 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_L: decomposition of litter, D_S: decomposition of SOM, D_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.

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 τ 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.

In the original Biome-BGC and in previous Biome-BGCMuSo versions the C: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 C: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 C: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_r(d)_D, F_r(T)_D, F_r(SWC)_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_r(d)_D), temperature (F_r(T)_D), and SWC in the case of different soil types (F_r(SWC)_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.

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 N₂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 NH₄ 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}{ll}{\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−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 (NH4propⁱ):

$$\begin{array}{}\text{(10)}& {\text{NH4prop}}^i={\text{sNH4avail}}^i+{\text{sNO3avail}}^i,\end{array}$$

where sNH4availⁱ and sNO3availⁱ 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 NH₄ 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_r(pH)_NI, F_r(T)_NI, and F_r(SWC)_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_r(pH)_NI), temperature (F_r(T)_NI), and SWC F_r(SWC)_NI in the case of different soil types. The pH and temperature response functions are independent of the soil texture.

Denitrification (DN Eqs. 6 and 7)

Denitrification flux is estimated with a simple formula (Thomas et al., 2013):

$$\begin{array}{}\text{(11)}& {\text{DN}}^i=\text{DNcoeff}\cdot {\text{SOMresp}}^i\cdot {\text{sNO3avail}}^i\cdot {\text{WFPS}}^i,\end{array}$$

where DN of the actual layer is the product of the available nitrate content (sNO3avail, kg N m⁻²), SOMrespⁱ ($\mathrm{g}\mathrm{C}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{d}}^{-\mathrm{1}}$) is the SOM-decomposition-related respiration cost, WFPSⁱ is the water-filled pore space, and DNcoeff is the soil-respiration-related denitrification rate (g C⁻¹), 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 N₂O emission and N emission

During both nitrification and denitrification N₂O emission occurs, which (added to the N₂O flux originating from grazing processes if applicable) contributes to the total N₂O emission of the examined ecosystem.

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

During denitrification, nitrate is transformed into N₂ and N₂O gas depending on the environmental conditions: NO₃ 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 C: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 ^∘C (−1 ^∘C in January and 21.2 ^∘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 cm³ 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.

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 ± 1 standard deviation around the observations. The improved model clearly outperforms the earlier version in simulating soil water balance.

Model performance was evaluated by quantitative measures such as the coefficient of determination (R²), 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², MAE, and MSE denote the coefficient of determination, mean absolute error, and mean signed error (bias) of the simulated values, respectively.

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², MAE, and MSE denote the coefficient of determination, mean absolute error, and mean signed error (bias) of the simulated values, respectively.

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 N₂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 g cm⁻³, pH(H₂O) is 7.3, CaCO₃ 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 6 m×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 t ha⁻¹ 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 N₂O efflux observations were available from 2020 and 2021. The N₂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. N₂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 N₂O efflux measurements. CO₂ and N₂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 N₂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 N₂O efflux to the observations.

Figure 11 shows the comparison of the simulated and observed NO₃ content of the topsoil for the two selected treatments. The results indicate that the model underestimates the topsoil NO₃ 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 NO₃ 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 12 shows the comparison of the observed and simulated N₂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 N₂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 N₂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 (± 1 standard deviation) of the measurements.

Figure 13 presents the comparison of the observed and simulated soil respiration for the same time period as for the soil N₂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 CO₂ 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 (± 1 standard deviation) of the measurements.

Overall, the model provided nitrate content, N₂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^∘ by 0.1^∘ 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¹.

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).

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 ±10 % 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.

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 >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_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.

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

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.