Modeling Root Growth , Crop Growth and N Uptake of Winter Wheat Based on SWMS _ 2 D : Model and Validation

Simulations for root growth, crop growth, and N uptake in agro-hydrological models are of significant concern to researchers. SWMS_2D is one of the most widely used physical hydrologically related models. This model solves equations that govern soil-water movement by the finite element method, and has a public access source code. Incorporating key agricultural components into the SWMS_2D model is of practical importance, especially for modeling some critical cereal crops such as winter wheat. We added root growth, crop growth, and N uptake modules into SWMS_2D. The root growth model had two sub-models, one for root penetration and the other for root length distribution. The crop growth model used was adapted from EU-ROTATE_N, linked to the N uptake model. Soil-water limitation, nitrogen limitation, and temperature effects were all considered in dry-weight modeling. Field experiments for winter wheat in Bouwing, the Netherlands, in 1983-1984 were selected for validation. Good agreements were achieved between simulations and measurements, including soil water content at different depths, normalized root length distribution, dry weight and nitrogen uptake. This indicated that the proposed new modules used in the SWMS_2D model are robust and reliable. In the future, more rigorous validation should be carried out, ideally under 2D situations, and attention should be paid to improve some modules, including the module simulating soil N mineralization.


INTRODUCTION
'SWMS' is a computer program for simulating water and solute movement in two-or three-dimensional, variably saturated media and is one of the most widely used mechanistic models.The model can simulate, in either two or three dimensions, soil water movement and pollutant transport under transient flow conditions in heterogeneous media using a finite element method (FEM) (Šimůnek et al., 1994, 1995;Palla et al., 2009).There have been many applications and hydrology-related case studies using this model.This model was applied to simulate the variably saturated flow within a green roof system, in order to characterize the hydrological behavior of the system (Palla et al., 2009).It has also been found to be useful inexamining and evaluating other models, such as the WAVE model (El-Sadek, 2009).The SWMS_3D version of the model was applied by Dages et al. (2008) to evaluate a three-dimensional Richards' equation-based modeling approach for simulating groundwater recharge.As the SWMS code is in the public domain and can be used freely, researchers have been adapting, extending, and developing the SWMS model by adding new components, introducing new boundary conditions and others features.
Crop models have been constructed to include transpiration, water and nitrogen limitations, and partitioning of biomass for assimilation between shoots and roots.Soil water content, N uptake, and root and crop growth are the main components for simulations using crop models and on-site measurements, which are of particular concern to both researchers and farmers.Corre-Hellouet al. (2009) used the STICS intercropping model to simulate crop growth and N accumulation in pea-barley intercropping.They found simulation results were in agreement with the observations from the experiment dataset.Improvements in crop growth, soil water, and groundwater modules in SWAT were implemented by Luo et al. (2008).These modifications improved simulations of crop evapotranspiration and biomass, as well as soil water dynamics, under dry soil profile conditions.An evaluation of the accuracy of the Opus model was made by Wegehenkel and Mirschel (2006) regarding the simulation of crop growth, soil water, and N balance dynamics under the site conditions of the northeast German lowlands.
Widely used nutrient response models that cover a range of crops are the EPIC models (Balkovic et al., 2013) and DSSAT models (Dokoohaki et al., 2016).These two models have been used to study the effects of climate and management on growth and yield (Zhang et al., 2009).However, these models are generally complex and, therefore, the model parameters are not easily obtained.In addition, the 1D nature of the models also means that they cannot be used for wide row crops or ridged crops satisfactorily.The major difference between this study and previous research is that the model proposed here attempts to deal with 2D situations, which are common in the real world.Further, the input parameters of the proposed model are easily available.
The study incorporated crop growth, root growth, and N uptake modules into SWMS_2D.This newly developed model can be used in a wide area of research.SWMS_2D was chosen as the base model for the following reasons: first, the source code of SWMS_2D is free and open to the public, while many other similar models are not.Second, the study attempted to simulate for 2D agro-hydrological cases.And third, SWMS_2D performs particularly well in modeling of soil water movement.There are indeed some other crop N models dealing with similar processes (soil water movement, crop growth, root growth, N uptake, and others), but they are not process based.The main objective of this study was to develop such a model for agro-hydrological simulations for various crops.In farmland simulations, there are wide row crops and ridged crops, which ideally should be considered as 2D cases.The model developed is process based, and thus complex.It involves many of the processes (including soil water movement, root and crop growth, and N uptake) that govern water and N dynamics in the soil-crop system.
Since winter wheat is one of the most important cereal crops globally, research into simulation and monitoring of this crop have been widely conducted.For example, Rev Bras Cienc Solo 2017;41:e0150064 Bechini et al. (2006) carried out experiments in northern Italy from 1986 to 2001 to quantify the dynamics of aboveground biomass, plant N concentration, and N uptake and used the measurements for testing the CropSyst model.Other similar studies can be found in a large body of literature (Groot and Verberne, 1991;Li et al., 2009;Zhang et al., 2009;Biernath et al., 2011).
The objectives of this paper were to develop the SWMS_2D model by adding the root growth model, crop growth model, and nitrogen uptake model; and to validate the newlydeveloped model by comparison between simulation and real measurement.

Governing equation for soil water movement and root water uptake
The governing equation for soil water movement and root water uptake in the 2-D dimension within the soil profile is expressed in terms of soil water content, θ, the pressure head, h, and root water uptake, S (Šimůnek et al., 1994): where θ (m 3 m -3 ) is the volumetric water content, h (m) is the pressure head, S (d -1 ) is the rootwater uptake rate, x i (I = 1.2) (m) are the spatial coordinates, t (d) is time, K ij A are components of a dimensionless anisotropy tensor K A , K (m d -1 ) is the unsaturated hydraulic conductivity function, h 0 (x,z) (m) is the initial soil pressure head in the profile, E (m d -1 ) is the maximum potential rate of infiltration or evaporation under the current atmospheric conditions, n i are the components of the outward unit vector normal to boundary, and ψ (m) is prescribed functions of x, z, and t.Equations 3 and 4 are the upper and lower boundary conditions, respectively.
Equation 1 is a non-linear differential equation and the finite element method (FEM) based on the source code of SWMS_2D (Šimůnek et al., 1994) is used.Details of SWMS_2D can be seen in Šimůnek et al. (1994).
The soil hydraulic functions are defined according to van Genuchten (1980) and Mualem (1976 Eq. 5 where Θ is the relative saturation, θ s and θ r (m 3 m -3 ) are the saturated and residual soil water contents, respectively, α (m -1 ) and n are the shape parameters of the retention and conductivity functions, m=1-1/n, and K s (m d -1 )is the saturated hydraulic conductivity.
The rootwater uptake rate S is defined according to Feddes et al. (1978) and Zuo et al.(2004): Eq. 7 Rev Bras Cienc Solo 2017;41:e0150064 Eq. 9 where S max (d -1 ) is the maximal specific water extraction rate under optimal soil water conditions, γ(h) is a dimensionless reduction function simulating water stress, L r (m) is root penetration depth, T pot (m d -1 ) is daily potential transpiration, calculated by the Penman-Monteith method (Tegos et al., 2015), L nrd (z r ) is the root length density distribution function, L d (z r ) (m m -3 ) is root length density at z r depth, z r is normalized z, ranging from 0 to 1 (z r = z/L r ), h 1 (m) and h 4 (m)are the soil matric potentials at the anaerobiosis point and wilting point, respectively, h 2 (m) is the value of the pressure head, below which roots start to extract water at the maximum possible rate, and h 3 (m) is the value of the limiting pressure head, below which the roots cannot extract water at the maximum rate.

Root growth model
The root growth dynamic model is divided into 2 sub-modules.One is for root penetration calculation, based on the study of Pedersen et al. (2010).The other is used to model root length density distribution, based on the study of Wu et al. (1999).Root penetration calculation in the 2-D dimension within the soil profile is according to equations 10-12: Eq. 10 Eq.11 where DD (°C) is cumulative day-degrees, DD lag (°C) is the lag phase for initial root growth, R z and R x (m) are vertical and horizontal root penetration depths, respectively, R z-min and R x-min (m) are initial root lengths, R z-max and R x-max (m) are maximum root lengths, k rz and k rx (m d -1 °C-1 ) are vertical and horizontal root growth rates, respectively, T air (°C) is daily average temperature, and T min (°C) and T max (°C) are the minimum and maximum temperatures for root growth.
Rev Bras Cienc Solo 2017;41:e0150064 The newlycoupled root penetration model, which is based on DD (°C), has been shown in field trials to be applicable to both monocot and dicot crops (Smit and Groenwold, 2005;Kirkegaard and Lilley, 2007;Pedersen et al., 2010).For a 1D case, root penetration can be simplified to equations 10 and 12.
Normalized relative root density distributions at different growth stages are quite similar, so Wu et al. ( 1999) introduced a third-order polynomial equation to describe L nrd of winter wheat as follows: where R i (i = 0, 1, 2, 3) are the polynomial coefficients and z r is normalized z, ranging from 0 to 1 (z r = z/L r ).

Crop growth and N uptake model
For crop growth: where W (Mg ha -1 ) is the dry weight per unit area of plant dry matter excluding fibrous roots, K 1 (Mg ha -1 )is a growth constant which is equal to 1 Mg ha -1 for all crops, K 2 is a growth coefficient characterizing crop growth when there is ample mineral N, ample available water, and a comfortable temperature for crop growth, and G N , G T , and G W are growth coefficients (ranging from 0 to 1) dependent on crop % N, DD, and water supply, calculated as described below (Rahn et al., 2007).Equation 14is the crop growth model developed by Greenwood (2001), which is adopted in the EU-ROTATE_N model (Rahn et al., 2007).
G N , G T , and G W are calculated in the following equation, respectively: Eq. 15

T pot
Eq. 16 Eq. 17 where %N is the actual %N in the dry matter of the whole plant (excluding fibrous roots), %N crit is the critical %N, T act and T pot (m d -1 ) are the daily actual and potential transpiration rate, and T max ' and T min ' (°C) are the maximum and minimum temperature for crop growth.
For equation 17, we revised G T to 0.3 instead of 0.0 when T air is less than T min '.
Eq. 18 %N crit = a.(1+b.e -0.26W ) Eq. 19 where N g is %N of the growth related tissue and N s is %N of the storage tissue, a and b are two parameters for critical N calculation, and ∆N (kg ha -1 ) is the one-day N increment of the crop.Equation 19 is used for calculating N uptake between t and t+∆t, and ∆t is daily in this model.
For N transformation in soil, a simple approach was adopted, based on the study of Greenwood (2001) and Rahn et al. (2007).This is partly because the simple method adopted in our study was based on experiments on agricultural soils, and this approach has successfully been used for many such situations.

Evaluation criteria
The performance of the model was evaluated using the following statistical indices (Nash and Sutcliffe, 1970;Willmot, 1981): the Nash-Sutcliffe coefficient of model efficiency (E NS ), the coefficient of determination (R 2 ), the index of agreement (D), the root of the mean squared errors (RMSE), and the mean error (ME).
Eq. 21 Eq. 22 Eq. 23 Eq.25 where y i and y i ' are the measured value and predicted value, y i-a and y i-a ' are the average value of the measured value y i and predicted value y i ' , and N is the number of measurements.

Validation case
During the 1982-1983 and 1983-1984 growing seasons, field experiments with winter wheat were conducted in Bouwing, the Netherlands.These experiments were to obtain data on crop growth and development, nitrogen uptake, soil nitrogen dynamics, and soil water dynamics under different nitrogen treatments at three sites (Bouwing, Eest, and PAGV) under natural weather conditions (Groot and Verberne, 1991;Yang et al., 2009).The data from the Bouwing experiment, undertaken during 1983 to 1984, was chosen for the current study.A summary of the experiment is given in table 1.The atmospheric data (rainfall, and minimum and maximum temperature) for the Bouwing experiment during 1983 to 1984 are presented in figure 1.There was no irrigation and only one fertilizer treatment in this case.Other details of the experiments can be seen in Groot and Verberne (1991).

Model parameterization
The model parameterization includes data for weather, soil, crop and others factors.The study used the measured values of soil water content and mineral N down the soil profile and the crop dry weight on Feb 14, 1984, as the initial condition (Table 2).Weather conditions, including rainfall and daily maximum and minimum temperature, were obtained from the closest meteorological station at Wageningen, which is 7 km away from the experimental farm at Bouwing (Groot and Verberne, 1991).The fitted values of   3.These values were obtained by applying RETC software (van Genuchten et al., 1991) to estimate soil water retention and conductivity curves for the 0.00-0.40 and 0.40-1,00 m layers on the experimental farm.
Values of R 0 , R 1 , R 2 , and R 3 are obtained from the study of Wu et al. (1999).From that study, it can be concluded that normalized relative root density distributions at different growth stages are quite similar when equation 13 is used.R z-min is set as rooting depth at sowing or planting, and R z-max is rooting depth at harvest.R z-min is 0.10 m according to Pedersen et al. (2010), and R z-max is 1,00 m from the experiment data (Groot and Verberne, 1991).The other four parameters, DD lag , K rz , T min , and T max , were obtained and evaluated from Pedersen et al. (2010), as there are too few parameter studies for winter wheat in this area.
As for the crop growth model, equation 14 was used differently from the study of Greenwood and Draycott (1989).Here we divided crop growth into two growth phases.In phase 1, K 2 = 0.3 d -1 and in phase 2, K 2 = 0.55 d -1 .K 1 was set at 1 Mg ha -1 during the whole growing season and initial dry weight was 0.033 Mg ha -1 .N s and N g were 1 and 5.5 % respectively for N uptake modeling, while a and b were 1.35 and 3.0 for calculating critical %N.T max ' and T min ' were set to 20 °C and 4 °C in equation 17 according to Rahn et al. (2007).If mean temperature is greater than T max ', GT is 1.0, indicating that the crop grows very well with no temperature restriction, while if average temperature is less than T min ', the crop grows slowly and G T = 0.3 is set.
Mineralization rate is influenced by many factors, but as there were not many mineralization conditions to consider in this study, the mineralization rate was setat an average rate of 0.6 kg ha -1 d -1 .
The soil domain was calculated to a depth of 1,10 m, and the top and lower boundary condition of water considered as variable flux and free drainage, respectively.The spatial discretization of the 0.012 m depth of the upper and lower boundaries is 0.04 m, while the other layer within the soil profile is from 0.01 to 0.04 m.

Soil water contents in different layers
The measured and simulated soil water contents according to time in various layers in this study are in figure 2. Generally speaking, good agreement was achieved between measurement and simulation.In the top 0.00-0.20 m layer, both the simulated curve and the measurement fluctuated from day to day, and changed frequently.This indicated that the upper layer was influenced more acutely than the other five layers by boundary conditions (infiltration and evaporation).For the 0.20-0.30and 0.30-0.40m layers, both measured and simulated soil water content decreased until a peak occurred in the 140 th to 180 th day period.This was because of heavy and intense rainfall events during this period.The simulated curve well fits the measurement for the other three layers from 0.40-1,00 m.Both the measured and simulated soil water content changed steadily and decreased slowly day to day, even when heavy and intense rainfall occurred.One possible reason was that water uptake of roots in these layers increased as root penetration and root length density increased.It should be noted that the layers from 0.00-0.40m are more easily influenced by upper boundary conditions than those from 0.40-1,00 m.
Evaluation of model performance was carried out for soil water content (Figure 3 and Table 4).Although the coefficient of determination R 2 and E NS are only 0.49 and 0.016 respectively, the fitted curve ŷ = 0.005 + 0.9753 x is very close to y = x.The values of D, RMSE, and ME are 0.81, 0.034, and 0.00013 m 3 m -3 respectively, indicating that the simulated results are acceptable.The overall comparison between measurement and simulation of soil water content are not very satisfactory.This is mainly because two different layers were simulated in this study.

Depth of root penetration and normalized relative root density distribution
A comparison between measurement and simulation for normalized root length density distribution is shown in figures 4 and 5. z r is normalized z, ranging from 0 to 1 (z r = z/R z ), while L nrd means normalized root length distribution, which can be transformed by equation 9. Simulation results of L nrd were obtained by coding the root dynamic model (Equations 10, 12, and 13) into SWMS_2D.On April 25, 1984, good agreement was achieved between simulation and measurement.Although good agreements were still obtained between simulation and measurement down through the 0.20-1,00 m soil the simulated values varied from the measured values in the 0.00-0.20 m layer during the latter two phases.Underestimated root length density distribution led to a lower water root uptake value for simulation; thus, soil water content in the 0.00-0.20 m layer was overestimated in simulation during the latter period (Figure 2).This can also be seen in figure 5, and there are two points located far from the 1:1 line.As a result, more attention should be     root length density distribution.Roots penetrated to 0.50, 0.90, and 1,00 m depths for three phases, while 0.527, 0.865, and 1,00 m were simulated in this study.

Dry weight growth simulation
Comparison between measurement and simulation of dry weight over time for this study can be seen in figure 6. Equations 14-20 were used for simulation and two strategies different from the study of Greenwood and Draycott (1989) were used.One was that K 2 was 0.3 t -1 before the 144 th day and 0.55 t -1 afterwards, indicating that the growing season was divided into two phases and each phase had a different potential dry weight growing rate.The second was that growth coefficients include G N , G T , and G W, which were recorded in the manual of the EU-ROTATE_N model (Rahn et al., 2007) and differed from the study of Greenwood and Draycott (1989).
Different gradients were found for simulated dry weight growth rate, which was induced by a different value for K 2 in two growing periods.From the measurement, we found the same phenomenon, and the modeling curve fit well.At harvest, dry weight was 15,162 Mg ha -1 for the field experiment and 15,300 Mg ha -1 for simulation.However, dry weight is overestimated before the 144 th day, and underestimated afterwards.High values for the coefficient of determination of 0.97 and an E NS of 0.96 were obtained.This, together with the gradient of the fitted line close to 1, indicates that the model performed well in dry weight simulation (Figure 7).
The variation in the G W , G T , and G N growing coefficients with time, for this study, can be seen in figure 8. Generally, G W changed together with variationin rainfall.It declined from 1.0 to 0.36 from the 44 th day to the 85 th day because there was insufficient rainfall, and soil water content in the soil domain decreased (Figure 2).However, G W was above 0.8 from the 88 th day to the 104 th day and from the 148 th day to the 170 th day period, as heavy and intense rainfall events occurred in both these periods.After about the 170 th day, G W declined with time, although there were still rainfall events during this period.The reason for these results was that crop potential transpiration increased as crop dry weight grew, thus the net water in the soil domain decreased every day.G T increased with rising average temperature as of day 44, although fluctuations were found throughout the whole growing season.In the beginning, G T was as low as the minimum value 0.3, and increased to as high as 1.0 at intervals in the end.
G N ranged from 0.65 to 1.0, with an average value as high as 0.82, indicating that growth rate was not greatly influenced by N. As crop dry weight increased, %N and %N crit of the crop decreased at different rates.From the 44 th day to the 151 st day, %N decreased more rapidly, which led G N to decrease from 1.0 to 0.65.After the 151 st day, G N increased with time, and reached 0.89 at harvest.The inflexion of G N for time was on the 151 st day, very near the 144 th day, when K 2 changed from 0.3 to 0.55.
The comparison between measurement and simulation of cumulative N uptake by the plant over time is in figure 9. Before the 124 th day, the simulated cumulative N uptake by plants was overestimated, compared with the experiment data, while after the 124 th day, it was underestimated.This was mainly caused by the simulation results of dry weight (Figure 6).Simulation results were divided into two phases by the 144 th day, and the different N uptake rate, seen from figure 9, is due to a different dry matter growing rate.Better agreements were achieved in the latter period of the growing season, compared with simulation for the beginning.At harvest, 194.36 kg ha -1 was simulated for cumulative N uptake by the plant, very close to the 194.70 kg ha -1 of the experiment data.The overall comparison between measurement and simulation of cumulative N uptake (Figure 10 and Table 4) also shows very high R 2 and E NS , and the best fit line is very close to the 1:1 line.
Indices of how good the fit for this model is can be seen in table 4. Apart from soil water content, the other simulated values are very acceptable, with a very high E ns and R 2 .For soil water content, the E ns and R 2 values are relative low, partly due to the fact that the soil water content measured was not markedly variable.This phenomenon, that the simulated soil water content is associated with low R 2 and E ns values, can be seen elsewhere (Kröbel et al., 2010;Nendel et al., 2012).
Compared with other crop nutrient models, such as STICS (Corre-Hellou et al., 2009)  Rev Bras Cienc Solo 2017;41:e0150064 using SWMS_2D as the core is 2D, which makes it more widely applicable.In addition, as the model uses relatively easily available data, it has the potential for use across a much wider field.While the model requires improvement in some areas, such as the simulation of soil N transformation, future studies will include further validation of the model, as well as strengthening of those weak modules.

CONCLUSIONS
The newly developed SWMS_2D is reliable for agro-hydrological simulation, especially for root growth, dry weight growth, and nitrogen uptake by the plant.Important points of the study are considerations of soil water limitation, nitrogen limitation, and temperature influence for dry weight growth modeling, which are now calculated by simple but reliable methods.
The study is the first step in validating the proposed 2D model, but more rigorous validation should be carried out, ideally under 2D situations.Attention should also be paid to improving some modules, including the module simulating soil N mineralization.

Figure 1 .
Figure 1.Atmospheric data (rainfall, the minimum and maximum temperature) for Bouwing experiment during 1983 to 1984.

Figure 2 .
Figure 2. Comparison between the measured and simulated soil water contents in different layers (The dashed line stands for rainfall).

Figure 3 .
Figure 3. Overall comparison of soil water content between measurement and simulation.

Figure 4 .
Figure 4. Comparison of normalized root length distribution between measurement and simulation.

Figure 5 .Figure 6 .
Figure 5. Overall comparison of normalized root length density distribution between measurement and simulation.

Figure 7 .
Figure 7. Overall comparison of dry weight between measurement and simulation.

Figure 8 .
Figure 8. Variety of growing coefficient G W , G T and G N with time.

Figure 9 .Figure 10 .
Figure 9.Comparison of cumulative nitrogen uptake by plant with time between measurement and simulation.

Table 1 .
Summary of the Bouwing experiment

Table 2 .
Initial values of soil water content (θ) and mineral nitrogen

Table 3 .
Fitted van Genuchten parameter values for the Bouwing experiments using the RETC software(van Genuchten et al., 1991, Yang et al., 2009) θ s and θ r are the saturated and residual soil water contents, respectively; α and n are the shape parameters of the retention and conductivity functions, K s is the saturated hydraulic conductivity.

Table 4 .
Statistical index values of performance of the model