Modified soil hydrological schemes for process-based ecosystem model Biome-BGC

Modified soil hydrological schemes were implemented in process-based ecosystem model Biome-BGC to improve water cycle simulation. The new schemes employed a 2-layer soil model and introduced new outflow components. Surface flow is input water exceeding soil infiltration capacity, and a scheme to predict it in the shorter duration than 1 day was proposed. Drainage from bottom soil layer and soil water redistribution between the layers driven by matric potential gradient and gravity were also modeled. The modified model was tested using the meteorological data and site parameters distributed with the original version. In comparison with the original, the modified model calculated 14% larger mean annual outflow. It also calculated a rapid outflow under dry antecedent soil condition, and reduced outflow peaks with slow drainage recession after soil wetted up, which were not simulated by the original model. The ecosystem carbon cycle responded to the change in soil water budget and net primary production (NPP) decreased by the modified model. Newly introduced parameters were determined automatically, so that users can use the same input data and parameter set for the original version. These simple schemes can be applied in other process-based ecosystem models.


INTRODUCTION
A number of process-based models for terrestrial ecosystems have been developed and applied in ecological, environmental, geophysical and socioeconomic researches aiming to investigate ecosystem energy and material cycles (Jung et al., 2007), landscape and succession (Seidl et al., 2011), biomass production and ecosystem services (Ooba et al., 2010), etc. The most significant advantage of those models is the capacity to simulate complex systems consisting of subsystems which are interdependent with each other. This feature makes it possible to predict the dynamic responses of ecosystems to both internal and external forcing, and the unforeseen future states of ecosystems especially under the influence of climate change (Melillo et al., 1993).
Biome-BGC (Running, 1993;Thornton, 1998) is one of the process-based ecosystem models focusing on carbon, nitrogen and water cycles at an ecosystem scale, and has been used in the ecological applications of forests (e.g., Ueyama et al., 2010) and grasslands (e.g., Milesi et al., 2005), as well as human-managed ecosystems after modifications (e.g., Hidy et al., 2012). One of the distinguished characteristics of the model is the relatively short calculation interval of 1 day for all processes, and it enables the explicit simulation of fast responding processes such as leaf phenology, water stress, and hydrological cycles including rainfall interception, evapotranspiration and runoff.
Although an advantage of Biome-BGC is its ability to simulate fast processes, its soil hydrology sub-model is too simple. The single layer bucket model considers only outflow from soil occurring when soil water content exceeds saturation or field capacity as runoff components, and does not include base flow. Some previous studies pointed out that part of the inaccuracy of the model was attributed to the soil sub-model and modified the hydrological schemes for their applications. For example, Pietsch et al. (2003) took the process of flooding and ground water flux into account for a forest in a floodplain. Hidy et al. (2012) developed a four-layer soil model calculating water percolation and diffusion across the layers to simulate water stress of crops. Gordon et al. (2004) assessed the performance of six ecosystem models in monthly runoff simulation by comparison with observation, and reported that runoff by Biome-BGC correlated well with the observation but showed relatively large negative bias (i.e., underestimation) among the models probably because base flow was not considered.
Aiming to improve the performance of water cycle simulation of Biome-BGC, this study modified the soil hydrological schemes. The basic development policy was to implement exact hydrological processes in soil including surface flow, base flow and water redistribution in soil with minimal modification and parameters added to the original model code to maintain the utility for users. A 2-layer soil model was selected rather than the multi-layer because of the consistency between temporal and physical dimensions, that is, a large water flux in shallow soil layers in a multi-layer model occasionally induces instability and inaccuracy in numerical computation with the calculation interval of 1 day. The modified model was tested using an example meteorological data and site parameter set distributed with the original model. Figure 1 illustrates the soil hydrological schemes in the Correspondence to: Takashi Machimura, Graduate School of Engineering, Osaka University, 2-1, Yamadaoka, Suita, Osaka 565-0871, Japan. E-mail: mach@see.eng.osaka-u.ac.jp original and modified Biome-BGC models. In the original model, water flows into the soil by throughfall, snow melt and canopy drainage, where canopy drainage is water temporally retained on wet canopy and comes in at the end of a day. Here, the term "throughfall" used by the original model developer means rain water not intercepted during a day consisting of throughfall and stem flow in commonly used terms. Inflow is first stored in the soil, and after subtracting evapotranspiration, water exceeding over-saturation or field capacity at the end of a day will flow out during the following day.

Model modification outline
In the modified model having two soil layers, water inflow, evapotranspiration and outflow over-saturation occur only from the top layer. Outflow over field capacity was removed. Instead of it, drainage from the bottom layer representing both over field capacity and base flows was introduced. Soil water moves across the top and bottom layer boundary both upward and downward according to the water potential gradient. The soil water flux F of drainage and inter-layer redistribution (kg m -2 day -1 ) in the vertical (z direction) is calculated based on Darcy's Law, where, θ is volumetric water content (VWC, m 3 m -3 ), K is hydraulic conductivity (m day -1 ), ρ is the density of water (kg m -3 ), ɡ is the acceleration of gravity (m s -2 ), ψ is soil matric potential (Pa) and z is vertical length (m). Another major modification is the introduction of surface flow that occurs when water inflow rate exceeds the infiltration capacity of soil, and this component is subtracted from inflow instantaneously before entering the soil. Soil water status affects various physical and biological processes in the model such as soil evaporation, stomatal conductance and soil organic matter decomposition. In order to maintain the consistency among related sub-models, the whole soil VWC and matric potential used in the original model was kept valid also in the modified model. Soil texture, and therefore the hydrological properties of soil were assumed to be homogeneous in whole soil. The depth of the top soil layer was determined by K s t c , where K s is saturated hydraulic conductivity (m day -1 ) and t c is calculation interval (= 1 day). This represents the maximum vertical reach of infiltrated water from the soil surface within the calculation interval. This approach may result K s t c exceeding the whole soil depth in the cases of very porous soil and/or very shallow soil layer. These extreme soil conditions are out of validity for the modified soil scheme and require different approaches such as simulating over shorter calculation intervals.

Hydraulic conductivity
Hydraulic conductivity is a fundamental soil parameter that controls water movement in soil, and the value depends both on soil texture and moisture. In the original model, saturated VWC θ s and matric potential are estimated by the percentages of sand, silt and clay following Cosby et al. (1984). This reference also proposed a generalized equation of saturated hydraulic conductivity K s as follows, where, sand % and clay% are the percentages of sand and clay by the US Department of Agriculture (USDA) classification procedure, respectively. Note that Equation (2) gives K s in inches h -1 and needs conversion into m day -1 . Following Clapp and Hornberger (1978) and Hidy et al. (2012), the hydraulic conductivity of unsaturated soil K (m day -1 ) was modeled as a function of VWC θ, where, θ s is saturated VWC (m 3 m -3 ) and b is the slope of logarithmic retention curve, and its dependency on soil texture proposed by Cosby et al. (1984) was also used in the original model.

Surface outflow
Water moves through the soil pore of saturated soil skin to enter into the soil layer, so that, the maximum infiltration water flux is restricted by saturated hydraulic conductivity. Surface outflow F out,surf (kg m -2 day -1 ) occurs when water inflow F in (kg m -2 day -1 ) exceeds the infiltration capacity of soil as, Because rainfall intensity of a storm is larger in the shorter observation periods, surface outflow may occur in short durations in a day even if daily inflow does not exceed the infiltration capacity. Figure 2 illustrates that a temporally concentrated storm generates surface outflow in shorter durations than the calculation interval (1 day = 24 h). The effect of duration on rainfall intensity is empirically modeled and Equation (5) is one of them (Sherman, 1931).
where, I(t) is rainfall intensity (mm day -1 ) in duration t (day), a is the scale parameter of storm magnitude (mm day -1 ) and n (0 < n < 1) is site specific shape parameter. Based on the formula, rainfall intensity in any duration can be calculated from daily rainfall I(1) (mm day -1 ), Among the components of water inflow F in , snow melt F in,sm and canopy drainage F in,cd are not changeable by duration. On the other hand, the temporal change of throughfall F in,tf is identical to that of rainfall I(t) after plant canopies wet up to saturation in large storms, and therefore it should respond to duration t in the same way of I(t). Assuming that F in,sm and F in,cd can infiltrate into soil prior to F in,tf , and that the intensity dependency of F in,tf on duration is similar to that of rainfall intensity I(t), surface outflow F out,surf (kg m -2 day -1 ) in any duration t is defined by the difference between F in,tf (kg m -2 day -1 ) and the infiltration capacity of soil reduced by F in,sm (kg m -2 day -1 ) and F in,cd (kg m -2 day -1 ), where, F in,tf (t) and F in,tf (1) are throughfall intensities (kg m -2 day -1 ) in the durations of t and 1 day, respectively, and t c is calculation interval (= 1 day).
Here, daily throughfall is scaled by infiltration capacity as, where, m is a scaling parameter. By assigning calculation interval t c = 1 day, the ratio of surface outflow in duration t to daily throughfall r is then defined as follows, r changes by t and reaches its maximum within the range of 0 < t ≤ 1 day when, The solution of Equation (10) F out,surf is then calculated by multiplying F in,tf (1) by maximum r and is used in the modified model instead of Equation (4). Equation (12) calculates a very small surface flow in tiny storms when Equation (11) gives very small t, for example shorter than 1 minute, which is however ridiculous and is outside the valid range of the rainfall intensity formula, Equation (5). In order to eliminate the surface flow in small storms, the use of Equation (12) should be constrained by the lower limit of t (e.g., t > 15 minutes). On the other hand, in large storms, optimum t may exceed 1 day by Equation (11), however it is also invalid within the calculation interval of the model. The upper limit of t ≤ 1 day gives the proper range of r ≤ n in Equation (12).

Soil water drainage and redistribution
Over-saturation outflow in the modified model is the same as in the original model but is calculated only in the top layer. The other soil outflow, drainage, is from the bottom layer and is driven by gravity. Assuming that matric potential gradient at the bottom of soil is zero, Equation (1)  where, θ bottom is VWC in the bottom layer (m 3 m -3 ). Field capacity is not used explicitly because it is not an exact indicator of drainage cessation, and instead, the rapid decrease of K by decreasing soil moisture can substantially control drainage cessation from dry soil. After calculating water input and output to/from soil in a day, the soil water redistribution between the top and bottom layers F rd (kg m -2 day -1 ) is calculated by the following equation based on Darcy's Law, Equation (1), where, θ and θ top are VWCs (m 3 m -3 ) of whole and top soil layers, respectively, and D is whole soil layer depth (m).

Implementation and simulation
The modified soil hydrological schemes were coded based on version 4.2 (final release) of Biome-BGC developed by Numerical Terradynamic Simulation Group (NTSG). In the implementation, the 4th order Runge-Kutta algorithm was used to improve the accuracy in drainage (Equation 13) and redistribution (Equation 14) water fluxes. The modified C codes were compiled using GCC 4.9.3.
The modified model was tested using example meteorological data at Missoula, MT, USA, a simulation and site parameter set, and default eco-physiological parameters for deciduous broadleaf forest, which were distributed with the original model. The test site was located in a cool temperate climate with latitude 46.8° N, and elevation 977 m a.s.l. Soil texture was 30% sand, 50% silt and 20% clay, and soil depth was 1 m. The rainfall intensity dependence on duration at the site was unknown, and therefore that the parameter n in Equation (5) was assumed to be 0.5, and the lower limit of t was set to be 0.01 day.
A simulation run of 71 years beginning in 1931 was performed following a spin up run by using the original and modified models, therefore, the results were for an ecosystem reaching its climax state.

Comparison of soil water budget
In the modified model run, the given site soil texture had a saturated hydraulic conductivity of K s = 0.272 m day -1 and therefore, a soil layer of 1 m deep was automatically divided into a top layer of 0.272 m and a bottom layer of 0.728 m. Table I compares the mean annual water budget produced by the original and modified models. Mean annual outflow was 234.6 mm year -1 by the original model, whereas the 267.7 mm year -1 by the modified model showed an increase of 14%. Consequently, evapotranspiration decreased in the modified model. Considering that Gordon et al. (2004) reported that runoff calculated by Biome-BGC was negatively biased from observation, outflow increase by the modified model may be beneficial in improving long-term runoff predictions.
Drainage accounted for 94% of mean annual outflow by the modified model, and the fractions of surface and oversaturation outflows were small. However, surface and oversaturation outflows contributed to the response improvement in individual large storms as shown later. A large difference was seen in the number of days with outflow; only 71 day year -1 by the original and 177 day year -1 by the modified, which was comparative to precipitation days of 181 day year -1 .
Daily outflow calculated by the original and modified models in the period from DOY 295 to 365 in 1952 is shown in Figure 3. Rainfall of 417.4 mm across five storms was observed in the period. This period was near the beginning of a rainy season in winter and the antecedent condition of soil was relatively dry. The largest outflow in the period appeared on DOY 300 by the modified model responding to the largest daily rainfall on DOY 298, where almost 100% of the outflow consisted of surface and over-saturation flows. In contrast, the original model did not calculate outflow during the first storm because soil VWC did not reach the field capacity. Thus, the modified model could respond to the large storm by considering rainfall intensity and saturation in the shallow top soil layer even when the VWC of the whole soil layer was low.
The original model calculated the high peaks and quick recession of outflow in the rest of the period indicating soil water reached the field capacity in the storm on DOY 315. In comparison, outflow peaks calculated by the modified model were relatively reduced and recessed more slowly after the storms. Drainage from the bottom soil layer accounted for 97% of total outflow, and over-saturation flow did not occur after DOY 315 because inflow in the top layer was quickly transferred to the bottom layer by the redistribution scheme, and was then drained out. This change in outflow components before and after DOY 315 calculated by the modified model was the result of hydraulic conductivity change by increased VWC, and a similar runoff component change depending on antecedent soil moisture conditions was observed in in situ studies (e.g., Hopper et al., 1990).
The mechanism of surface outflow in the modified model is identical to so called Hortonian surface flow, and is commonly recognized not to occur in forest soils. The modified model rarely calculated surface outflow in the forest because it calculated soil hydraulic conductivity based on Cosby et al. (1984) who determined the empirical Equation (2) from soils collected mainly from pastures and cultivated lands, and therefore its validity may not be guaranteed for very porous forest top soils. This mechanism of surface outflow should be implemented in the model for utility in a variety of non-forest or degraded forest ecosystems. Replacing Equation (2) for target ecosystem types may improve the simulation and is not difficult for the model users. Table II shows mean carbon contents in pools and net primary production (NPP) calculated by the original and modified Biome-BGC models. NPP, and consequently vegetation carbon storage (biomass carbon) decreased by 5% by the modified model, and they are consistent with evapotranspiration which decreased by 4% (see Table I). Note that longterm mean annual heterotrophic respiration consisting of litter and soil organic matter decomposition was also depleted by the modified model because it should be quasi-balanced with NPP in a climax forest. Whereas, soil carbon storage was slightly larger by the modified model. This suggested that the depletion of soil organic matter decomposition was more sensitive than the depletions of NPP and litter production in the dryer soil conditions calculated by the modified model in this case.

CONCLUSION
Modified soil hydrological schemes were implemented in process-based ecosystem model Biome-BGC to improve water cycle simulation. A single soil layer in the original model was divided into two layers, top and bottom, and soil water movement was modeled based on Darcy's Law. Surface outflow was defined as inflow exceeding the soil infiltration capacity, and a scheme to predict the effect of throughfall intensity change by duration on surface outflow was proposed. Drainage from bottom soil layer and soil water redistribution between the layers driven by matric potential gradient and gravity were also modeled.
The modified model was tested using the meteorological data and site parameters at Missoula, MT, USA. In comparison with the original model, the modified model calculated 14% larger mean annual outflow. It also calculated a rapid outflow under dry antecedent soil conditions, and reduced outflow peaks with slow drainage recession after soil wetted up, which were not simulated by the original model. The ecosystem carbon cycle responded to the change of soil water budget and NPP decreased by 5% by the modified model.
As demonstrated in this study, the modified model improved outflow calculation both in annual amount and short-term change from the original model, and therefore is potentially capable of improving the accuracy of ecosystem water cycle simulation. A model validation on soil moisture and runoff by comparison with observations was not conducted in this paper, and should be done for the actual cases in applications. The simple schemes could be implemented by minimal modification to codes from the original version. Almost all newly introduced parameters were determined  automatically from existing parameters and variables used in the original model, so that users can use the same input data and parameter set for the original version without modification. The simple soil hydrological schemes can be applied in the other process-based ecosystem models which currently have single layer soil schemes. This study did not focus on the full biogeochemical processes other than soil hydrology, and the further improvement of schemes is now considered to assess the soil water related processes. For example, soil evaporation, water uptake by plants, stomatal conductance control and soil organic matter decomposition should depend on the soil moisture and matric potential of the top soil layer more sensitively than on the whole soil layer, and nitrogen leaching by over-saturation flow and drainage may be different. The implementation of these processes is not difficult based on the modified Biome-BGC model presented in this study and will be done in coming reports.