Soil related developments of the Biome-BGCMuSo terrestrial ecosystem model by integrating crop model components

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, theoretical basis is provided for the latest version of 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. During the developments we took advantage of experiences from the crop modeller community where internal process representation has a long history. In this paper the improved soil hydrology and soil carbon/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 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 m3 m−3 for SWC, and from −0.68 mm day−1 to −0.2 mm day−1 for soil evaporation). Sensitivity analysis and optimization of the decomposition scheme is presented to support practical application of the model. The improved version of Biome-BGCMuSo has the ability to provide more realistic soil hydrology representation and nitrification/denitrification process estimation which represents a major milestone.



Runoff
Calculation of RUNOFF [kg H 2 O m -2 day -1 ] is based on a semi-empirical method that uses the precipitation amount (PRCP), the dimensionless runoff curve number (RCN), and the actual moisture content status of the topsoil (Williams, 1991). If PRCP is less than a critical value (product of a soil type and a moisture content dependent factor; k type and k moist ), runoff is equal to zero. Otherwise: , if 1 where: , and are the saturated, the hygroscopic and the actual volumetric water content of the top soil layer.
The amount of runoff is affected also by the setting of pond water calculation: first the pond water accumulation is calculated (see details in Section 3.1), only the surplus above the maximum pond depth leaves as runoff.

Downward (gravitational) water flow
The percolation calculations (tipping bucket method) are carried out layer by layer starting from the 1 st soil layer down to the 9 th soil layer. The bottom (10 th ) soil layer is special: the percolating water (and dissolved materials) from the bottom layer is a net loss from the system. It is assumed that the soil can hold moisture against gravity up to its field capacity.
One of the most important variables of percolation calculation is the drainage coefficient (DC) which determines the amount of water leaving a layer in a day. DC is a soil input parameter but in case it is not provided by the user, it is estimated from the hydraulic conductivity:  total drainage is limited by saturated conductivity: if :  water content is updated:  the excess water (EXCESS i ; cm) is redistributed in layers above (from the actual layer up to the 1 st one):  first the capacity of the layer above is filled and in the case of remaining excess moving one layer above  if at the end of the redistribution (in the top soil layer) is greater than 0, its value is added to pond water and/or runoff depending on the actual pond depth and maximal pond depth, which is an input soil parameter  infiltration into the next layer is the drainage from the layer above: If is less than or equal to :  water content update:  drainage is calculated only if soil moisture content is above field capacity, otherwise it is zero:  drainage is limited by the saturated conductivity:  state update moisture content with drainage:  infiltration into the next layer is the drainage of the above one: On rainless days the calculation is simpler: there is no infiltration at the top of the soil profile  Maximum of drainage (DRMX; cm) is greater than zero only if soil moisture content is above field capacity  HOLD is greater than zero only if soil moisture content is below field capacity  Drainage rate in the top layer:  Drainage rate in layers below:  Drainage is limited by saturated conductivity:  The SWC is updated from 9 th layer to 2 nd layer with saturation as the upper limit: if  Status update of SWC in the top (1 st ) layer:

Capillary water flow
The diffusion calculations follow the formulae of the 4M model (Fodor et al., 2014). The downward diffused water from the bottom active (9 th ) layer is a net loss for the soil system, while the upward diffusion flux from the passive (10 th ) layer is a net gain.
 water flow between the actual layer and the one below is estimated with the product of the gradient (GRAD) of the normalized water contents (THET) between two adjacent layers and the diffusivity at the border of the two layers (DBAR). Diffusivity is also assumed to be a function of the normalized water contents of the two adjacent layers.
weighted average water content of the two layers: where are the parameters of diffusion calculation, default values are 0.88, 35.4 and 100, respectively.
 gradient between the actual and the layer below:  water flow between the actual and the layer below:

Soil evaporation
The method calculates the cumulated soil evaporation (soilEvap) in two phases (Ritchie, 1972 In the second evaporation phase the cumulative evaporation is a quadratic function of the time (number of days) since the last rainfall event (dsr).
First evaporation phase: if : if : Second evaporation phase: If : if : if : If : limitation: potential evaporation: Soil evaporation cannot be larger than the available soil water content (soilwAVAIL) in the top layer which is the difference between the actual and the hygroscopic water content. The hygroscopic water content is a soil input parameter.

Transpiration
Transpiration demand ( where: is the proportion of the root length in the given soil layer, is the rooting distribution plant input parameter (detailed description in User's Guide).

FC-rising effect of groundwater
The FC-rising effect of groundwater ( ) for the layers above the water table (WT) is calculated based on the ratio of actual (GWdist) and critical distance ( ) from WT, but only after the layers below have already been charged up to their modified FC values.
is equal to zero for layers above the capillary fringe.
The modified FC due to the close groundwater table is calculated with a linear function: The water flux coming from the groundwater on a given day and in a given layer (dischargeGW; kgH 2 O m -2 day -1 ) is a net gain for the soil profile above the water table and is the function of drainage coefficient (DC) of the given layer.

Nitrogen mineralization function
The nitrogen mineralization fluxes of the litter pools are the function of the corresponding potential litter loss rate ( in kg C m -2 day -1 ), C:N ratio (kg C kg N -1 ) and respiration fraction. The nitrogen mineralization fluxes of the SOM pools are the function of carbon loss rate, which is the function of potential rate constant (k_base: reciprocal of residence time in day -1 ), and environmental integrated response function (Fr).
where: is the potential carbon loss of the source pool (A) in the actual layer, is the respiration factor of the source pool, and are the C:N ratios of the given source and target compartments in the actual layer, respectively is the actual rate constant scalar of the source pool in the i th layer, is the potential rate constant scalar of the source compartment (soil input parameter, see detailed in User's Guide), is the response function of the i th layer. A index denotes the source pools (L1, L2, L4, S1, S2, S3), B index denotes the target pools (S1, S2, S3, S4), AB index denotes the transformations between different pools: L1-S1, L2-S2, L4-S3, S1-S2, S2-S3 and S3-S4 (see Figure 4 in Section 4.2).
The decomposition of the passive SOM pool (S4) is special: only mineralization is assumed, and the respiration factor is supposed to be equal to 1, namely all released carbon is respired, no target pool is defined.
The response function ( ) of the transformation processes is the product of depth, water and temperature dependent factors ( , ,and , respectively) calculated layer by layer.
The factor of the i th layer is the exponential function of ratio the layer's midpoint ) and e-folding depth of decomposition rate's depth scalar (ED in m; see detailed in the User's Guide): The factor is zero below HW, is equal to 1 between FC and a critical SWC ( is defined by a soil input parameter (see User Guide).

;
The , is the bell function of the soil temperature in the given layer ( ):

Nitrification
The first step of nitrification calculation is to determine the net ammonia mineralization flux in the actual layer ( ). Net ammonia mineralization is the sum of mineralization fluxes detailed in Section 4.2) and the ammonia mineralization of the passive humus pool ( ).
Nitrification flux is calculated layer by layer based on the following Equation: where: is the net mineralization proportion of nitrification parameter soil input parameter, , , are the response functions of soil temperature, soil water content and pH in the actual layer, respectively. , is the same as , is also a trapezoidal function (such as ) depending on soil input parameter (detailed in the User Guide). is an empirical, exponential function of pH (Parton et al., 2001). If the value of the response function is zero (below a critical value of the given soil properties), the nitrification process is totally limited. If the value of the response function is 1 (above a critical value of the given soil properties), the process is not limited. In case the actual soil water content exceeds the optimal soil water content, the response function can decrease accounting for the saturation related stress (anoxic soil).

N 2 O-emission and N-emission
During both nitrification and denitrification N 2 O-emission occurs which (added to the N 2 Oflux originated from grazing processes) contributes to the total N 2 O-emission of the simulated ecosystem.
In Biome-BGCMuSo v6.2 a fixed part (set by the coefficient of N 2 O emission of nitrification soil input parameter; User Guide) of nitrification flux is lost as N 2 O and not converted to NO 3 .
During denitrification, nitrate transforms to N 2 and N 2 O gas depending on the environmental conditions: NO 3 availability, total soil respiration (proxy for microbial activity), soil water content and pH. The following equation is used to calculate the ratio of N 2 and N 2 O that is generated during denitrification: where is the ratio of N 2 and N 2 O gas (both expressed as N equivalent mass), represents the effect of NO 3 availability, is a total soil respiration related modifying factor (proxy for microbial activity), represents the effect of soil water content, and is included to represent the effect of soil pH (this latter was not included in Parton et al., 1996) in the given layer. Importance of pH is highlighted in Wagena et al. (2017).
It is important to note that the equations presented in Wagena et al. (2017) are incorrect. The exact formulae of the equations are much better described in Parton et al. (1996)  The equation for (Parton et al., 1996): where NO3 is the nitrate content of the soil expressed in μgN g -1 (=ppm).
Effect of soil respiration on the ratio in the given layer: where is total soil respiration (respiration of litter and soil pools) expressed in kg C ha -1 day -1 .
Effect of soil water content status (Parton et al., 1996) in the given layer: Effect of soil pH (Wagena et al., 2017) in the given layer: The ratio is used to calculate denitrification related N 2 O flux of the given layer from the amount of denitrification (DN) using this formula: Denitrification related N 2 /N 2 O ratio multiplier ( ) as a soil input parameter is used to take into account the effect of the soil type (Del Grosso et al., 2000).