A semianalytical solution of the modified two‐dimensional diffusive root growth model

Accurate estimation of the temporal and spatial root water uptake patterns in root zone is needed for an improved understanding of water and chemical transport dynamics in vadose zone. Rooting system is of great importance to describe the plant root water uptake directly. In this study, the diffusive root growth model was coupled with the Environmental Policy Integrated Climate (EPIC) crop growth model through changing the boundary condition, and a new semianalytical solution of the modified root diffusive model was derived considering the dynamic root length density distribution simulated by the EPIC crop growth model. To test the modified root growth model, the field‐measured data of maize (Zea mays L.) and tomato (Solanum lycopersicum L.) root length density distribution were used for parameter optimization. A MATLAB program was developed by coupling the modified root diffusive model with the genetic algorithm to facilitate the parameters optimization. Results showed that the simulated root length density distribution at different times was in a good agreement with the observed values. The RMSE and bias values ranged from 0.22 to 0.25 cm cm–3 and from −3.0 to 24.5%, respectively. The modified diffusive root growth model can therefore be used to simulate the two‐dimensional root growth during the crop growing period.

management strategy (Wöhling & Schmitz, 2007). It is well known that the interaction between soil water dynamics and crop growth is very complicated; however, both soil water process and crop growth process can be quantified with mathematical models. The interaction between the two processes can be described with the term "root water uptake." Therefore, accurate estimation of the temporal and spatial root water uptake patterns is needed for an improved understanding of the magnitude of soil water transport to atmosphere and groundwater (Vrugt et al., 2001).
Actually, root water uptake in field condition is a threedimensional (3D) problem. However, many crops (e.g., maize [Zea mays L.], tomato [Solanum lycopersicum L.], and pepper [Capsicum annuum L.]) are planted in rows. The distance between two plants is quite small, whereas the distance between two rows is much larger than the former one. Plant roots develop and connect each other along the row direction, whereas they develop but do not connect along the direction perpendicular to the row. Thus, for simplicity, soil water flow and root water uptake are treated as two-dimensional (2D) (Li et al., 2015;Skaggs et al., 2004;Wang et al., 2014;Wöhling & Schmitz, 2007). Many models and software packages (i.e., SWMS_2D [Šimůnek, Vogel, & van Genuchten, 1994], CHAIN_2D , Nitrogen-2D [Lu et al., 2004], and HYDRUS-2D [Šimůnek et al., 1999]) have been developed to quantify the 2D soil water flow. However, the root water uptake term in most models does not account for the root dynamics. For example, the HYDRUS model (Šimůnek et al., 2006) used the root water uptake term of Vrugt et al. (2001) without considering the plant root developing process. Some other models have included the root distribution pattern and root development in the root water uptake terms (Roose & Flowler, 2004;Wang et al., 2014). For example, Wang et al. (2014) developed a model to simulate the root growth and root water uptake, and in turn to simulate soil water flow and crop yield on the basis of coupling CHAIN_2D with crop growth model EPIC (Environmental Policy Integrated Climate). In this developed model, the root depth growth model was included in the Vrugt root water uptake term (Vrugt et al., 2001) with simplified assumption that the horizontal root length is a proportion of rooting depth. Furthermore, Wöhling and Schmitz (2007) developed a physically based coupled model for simulating one-dimensional (1D) surface-2D subsurface flow and plant water uptake. Within this model, an empirical exponential function with uniform root density distribution was applied to describe the root water uptake. However, these simplifications cannot quite well describe the root growth in the horizontal direction. Therefore, a multidimensional root growth model for describing the dynamic root distribution in vertical and horizontal directions is of great importance to accurately model soil water movement in vadose zone.

Core Ideas
• A novel diffusive root growth model was developed by coupling the EPIC crop growth model. • A semianalytical solution of the modified diffusive root growth model was derived. • The modified root growth model was tested using field experimental data of maize and tomato.
The root system architecture model and diffusive-type root growth model are the two main types of models capable of describing the 2D/3D dynamic root distribution in vertical and horizontal directions. The architecture model can explicitly describe and visualize the root system and analyze plantenvironment interactions, especially in 3D cases (Dunbabin et al., 2002;Dupuy et al., 2010;Pagès et al., 2004). However, the root system architecture model is complex with incorporating the space and considering the biophysical interactions between the root and its environment. Consequently, it requires a large number of input parameters that are difficult to obtain (Dupuy et al., 2010). Furthermore, the root architecture model requires sophisticated algorithms as coupling with soil process models (Dupuy & Vignes, 2012), whereas the diffusive-type root growth model is a simple model and provides good spatial description of root growth with fewer parameters (Dupuy et al., 2010). Page and Gerwitz (1974) first developed the diffusion root growth model in analogy to the equations for describing diffusion or heat flow. They assumed that the root proliferation was driven by the gradient of root density, that the root growth rate was a constant, and that the downward proliferation was independent of soil state. With the assumptions of variable diffusivity and root growth rate, Hayhoe (1981) developed a more general diffusion model. Acock and Pachevsky (1996) proposed a generic 2D convective-dispersive root growth model. They added the model with a convection term including the convection-like downward propagation caused by geotropism. In addition, they also assumed that the roots were classified as two types: young roots and mature roots. Elsewhere, de Willigen et al. (2002) considered the anisotropic coefficients in the diffusion model, which can better characterize the root development in different directions. Although it is difficult to attribute explicit biological meaning to the parameters, the diffusive root growth model has been proven to be powerful in describing the root growth and proliferation, as well as water distribution and uptake (Acock & Pachevsky, 1996;Hayhoe, 1981;Heinen et al., 2003;Reddy & Pachepsky, 2001). As a plant develops, its biomass accumulates and is partitioned into the root system and shoot system; all these biological processes can be described by a crop growth model. However, plant biomass accumulation was not counted in the diffusive root growth model. Thus, it is necessary to integrate the diffusive root growth model and crop growth model to provide better physical insight into the root proliferation process and obtain more accurate estimation of root distribution.
Several have tried to integrate the diffusive root growth model with the crop growth model (Dathe et al., 2014;Mollier et al., 2008;Pronk et al., 2005). For example, Dathe et al. (2014) described potato (Solanum tuberosum L.) root growth using a diffusive root growth model that was coupled with a specific crop growth model (SPUDSIM). However, the root depth growth dynamics were not considered in the above coupled crop growth models. At present, there are many simple and universal crop growth models including the root depth growth process, such as CropSyst (Stöckle et al., 2003), WOFOST (Spitters et al., 1989), the Daisy crop production model (Abrahamsen & Hansen, 2000), and the EPIC crop growth model (Williams et al., 1989). Among these models, the EPIC crop growth model with less demanding data input, would be preferred to be coupled with the diffusive root growth model.
In general, the diffusive model can be solved by using either a numerical method or an analytical method. With the assumptions of variable diffusivity and root growth rate, the numerical solution was obtained using the finite element method (Acock & Pachevsky, 1996;Hayhoe, 1981;Heinen et al., 2003;Pronk et al., 2005). The analytical solution can also be derived in more limiting situations, such as those with a constant diffusivity or a constant root growth rate (de Page & Gerwitz, 1974). As the diffusive model was coupled with the crop growth model, it will be rather difficult to obtain its analytical solution similar to that of de . A semianalytical solution is an option for the coupled model. Therefore, the objectives of this study were (a) to develop a modified diffusive root growth model integrating with the EPIC crop growth model; (b) to derive a semianalytical solution of the modified diffusive root growth model with considering the constrained rooting depth; and (c) to test the modified root growth model using field experimental data of maize and tomato.

Diffusive root growth model
Root growth can be considered as a diffusion process driven by the gradient of root density, and the diffusive root growth model was developed in analogy to the equations for describing diffusion or heat flow (Page & Gerwitz, 1974). With the assumptions of neglecting the geotropic velocity, the impacts of water stress and overwetting on root growth and soil heterogeneity, and so on, de  developed a 2D diffusive root growth model. In this model, the root decay was treated as a first-order sink term. The root biomass accumulation at soil surface was treated as a flux boundary condition at a rate of Q. The diffusivity coefficient (D) was used to reflect the effect of soil variables on root growth. D was assumed as a constant. To describe the root growth in a row planting system, the root zone was conceptualized as a 2D rectangular geometry. In the horizontal direction (X direction) perpendicular to the row, the domain is finite. In the vertical direction (Z direction), the region is considered to be infinite . Then, the governing diffusive equation for root length density in 2D Cartesian coordinates can be described as The boundary conditions are The initial condition is where C is the root length density (cm cm −3 ), X and Z are the spatial coordinates in the horizontal and vertical directions, respectively, T is time (day), λ ′ is the decay rate constant (d −1 ), X 1 and L are the half distance between two crops and two rows (cm), respectively, D x and D z are the diffusivity coefficients (cm 2 d −1 ) in the X and Z directions, respectively, and Q is the fine root growth rate located at the surface in root system (cm cm −2 d −1 ).

Modification of boundary conditions
In this study, to better predict the root distribution, the diffusive root model was coupled with the EPIC crop growth model (Williams et al., 1989). The EPIC model needs less data input and uses a unified approach to simulate the growth for >80 types of crops (Williams et al., 2006). The EPIC model includes crop phonological development based on daily accumulated heat unit, a harvest index for partitioning grain yield, the Monteith's approach for potential biomass accumulation, and water and temperature stress adjustments. With the coupled model, the rooting depth estimated by both the diffusive model and the EPIC crop growth model should be the same at any given time. Therefore, the bottom and lateral boundary conditions of the diffusive root growth model should be changed as follows: where RD is the root depth, which was considered as a constant value on a certain day; it can be derived from the EPIC crop growth model (Williams et al., 1989): where RD max is the maximum root depth, RD i is the root depth of ith day, ΔRD is the daily change of root depth, and HUI i is the heat unit index and is given as (Williams et al., 1989) where PHU j is the potential heat unit required for maturity of crop j. PHU may be set as a default value or calculated by the model from planting to harvest. HU k is the value of heat unit, which is specified by (Williams et al., 1989) where T mx,k and T mm,k are maximum and minimum air temperatures for the kth day (˚C), and T b,j is the crop-specific base temperature of crop j (˚C). In this study, the parameter T b,j for maize and tomato was set as default values referring to the manual of EPIC model (Williams et al., 2006).

Semianalytical solution of modified diffusive root growth model
For simplicity, we define the following dimensionless terms to transform Equations 1-3 and Equations 6-9 into the dimensionless ones: where ref is the reference root length density, which is 1 cm cm −3 in this study. Substituting Equation 14 into Equations 1-3 and Equations 6-9) leads to Applying the Laplace transform to Equation 15, one obtains where the Laplace transform of c is denoted by u, and the Laplace parameter is s. The boundary conditions retain their original forms except for Equation 16, which becomes in which φ( ) is the Laplace transform of ( ). Then, the finite cosine transform with respect to x is taken, and this transformation is defined as (Churchill, 1972, p For the boundary condition of Equation 18, the general solution of Equation 25 should bē Therefore, the suitable solution of Equation 25 should bē ) , for = 1, 2, 3, … (33) As the root growth diffusive model was coupled with a crop growth model, then ( ) can be considered as a constant value in a certain time (i.e., 1 d). Eventually, A 0 ,ū 0 , A n , andū n become (Equations 30-33): Furthermore, the solution for u is given by (Churchill, 1972, p. 355): Then dimensionless root length density c can be obtained as follows where L −1 denotes the inverse Laplace transform. In Equation 39, 200 summations were used to obtain the analytical solution. However, there are no analytic expressions for the inverse Laplace transformation of u 0 and u n . Therefore, the root length density in the real-time domain was calculated by the numerical inverse Laplace transform algorithm such as Stehfest (1970), which is presented as (41) where n must be a positive even integer, and P is a Laplace function. The V i is coefficient matrix and calculated by using the Fortran program of Zhan and Zlotnik (2002).

Field description and measurements
The modified diffusive root growth model was tested by using the observed maize and tomato root length density distribution data obtained from field experiments (Gao et al., 2010;Zotarelli et al., 2009 (Gao et al., 2010). Maize was sown on 16 April and harvested on 18 August in 2007 and 20 August in 2008. Root samples were collected from the soil profile by pressing sharp-edged iron boxes vertically into the soil surface. The soil cores were collected from the top layer at 10-cm intervals along the north-south axis and the east-west axis. Root length was measured with the modified Newman-line-intersect method (Tennant, 1975). Only a part of the roots was measured directly, and the rest of the roots were measured by using the regression of root length and root weight. In this study, the treatment of sole maize was selected to test the modified diffusive root growth model (Gao et al., 2010 Note. D x , diffusivity coefficient in the x direction; D z , diffusivity coefficient in the z direction; λ′, decay rate constant; Q, ; PHU, potential heat unit.

Parameter optimization and model test
Genetic algorithm (GA) was applied to get the optimized parameters of the modified diffusive root growth model based on the observed root length density of the soil core samples. For a relatively fast convergence to the global optimum, we used a relatively high crossover fracture value of 0.85 (Vrugt et al., 2001). In addition, the population size and generations were 1,000 and 30, respectively. The fitness function for GA was described as follows (Vrugt et al., 2001): where OF(s) is the objective function, n is the number of measurements; b′(x, z) and b(x, z; s) are the measured and predicted root length density at the point (x, z) according to the semianalytical solution of the modified model, respectively; and s is the parameter vector representing the fitting parameters. Smaller values of the objective function OF(s) imply that better fitting is obtained. The allowable ranges of the parameters included in s are shown in Table 1 according to Heinen et al. (2003). The GA optimization was carried out using MATLAB (version 7.12.0 [R2011a], the Math Works). We developed a MATLAB program by coupling the modified root diffusive model with GA for facilitating the computation. A GA is a powerful tool for parameter identification, when the number of fitted parameters is large (Vrugt et al., 2001). Although GA is an effective way to determine the global minimum region, it is not the most efficient way to find the exact optimum location (Vrugt et al., 2001). In order to obtain the optimization parameters, three steps were conducted in this study. Firstly, GA was used to determine the global minimum region. Then, 100 series of the parameters of the modified diffusive root growth model were tested to search the exact optimal location. Finally, the fitted parameters were determined as the RMSE being the smallest in the F I G U R E 1 Simulated maize root length density vs. observed values 100 tested series. To calibrate the diffusive root growth model, the observed maize and tomato root length densities were used to optimize the parameters. The observed root length densities were compared with the simulated ones. Then, the fitted value, the medium value, and the 95% confidence interval for each parameter were analyzed. The RMSE was used to evaluate the model performance for predicting root length density according to the criteria proposed by Willmott (1982). It represents the discrepancy between observations and predictions. The closer RMSE is to 0, the more accurate the model is. The RMSE was described as follows: where O i and P i are the observed and simulated values respectively, and n is the number of the pair values. In addition, bias was used to identify the model over-or underpredicted root length density and calculated as It becomes positive when the model over predicts and negative when the model under predicts observed values.
Statistical analysis of optimized parameter series was performed by the descriptive statistical of SPSS 16.0 software (SPSS).

RESULTS AND DISCUSSIONS
To calibrate the modified root growth model, the observed root length density values were compared with the simulated ones based on the fitted parameters. As shown in Figure 1, F I G U R E 2 Simulated tomato root length density vs. observed values a good agreement can be found between the simulated and measured maize root length densities. However, large discrepancy appeared when the root length density was high. In addition, statistical test results showed that the bias was 24.5%, which was lower than the result of Dathe et al. (2014), in which the bias ranged from −50.7 to 226.3% for potato root density at harvest. In addition, the bias between the predicted tomato root length densities and the observed ones was −2.8% ( Figure 2). It indicated that the prediction produced a good agreement with the observation results. The fitted parameters were shown in Tables 2 and 3 for the maize and tomato root growth model. To show the stability of the optimization, the predicted medium value and 95% confidence interval for each parameter of maize in 2007 and 2008 are presented in Table 2. As shown in Table 2, the decay rate constant λ′ has a relatively higher uncertainty as compared with the other parameters. It may be because decay rate is a sensitive parameter. If the observed decay rate value was considered as an input value, the predicting accuracy of root growth and distribution can be improved. However, the lowest CV value was observed for PHU. The fitted values of parameters D x , D z , λ′, Q, and PHU were in the ranges of 0.4320-0.6073 cm 2 d −1 , 3.2086-3.6974 cm 2 d −1 , 0.0001 d −1 , 0.1682-0.225 cm cm −3 d −1 , and 3,447-3,460°C, respectively. Similar results were obtained by de , who estimated D z = 5.6 cm 2 d −1 when solving the analytical solution of diffusive root growth model for maize. Heinen et al. (2003) got a numerical solution and the diffusive coefficient in horizontal and vertical directions were 2.234 and 4.552 cm 2 d −1 for maize-row, respectively. For tomato, the predicted medium value and 95% confidence interval as well as the fitted values for each parameter are shown in Table 3. The parameter λ′ had the greatest CV value. The fitted values of parameters D x , D z , λ′, Q, and PHU T A B L E 2 The fitted parameters, the 95% confidence intervals and the standard deviation of the parameters of the diffusive root growth model for maize in optimization within and without considering potential heat unit (PHU) Note. D x , diffusivity coefficient in the x direction; D z , diffusivity coefficient in the z direction; λ′, decay rate constant; Q, fine root growth rate located at the surface in root system.

T A B L E 3
The fitted parameters, the 95% confidence intervals and the standard deviation of the parameters of the diffusive root growth model for tomato in optimization within and without considering potential heat unit (PHU) When the diffusive root growth model was coupled with the crop growth model, the parameter PHU can be calibrated by the leaf area index, crop biomass, and yield. To reduce the number of optimized parameters and increase the optimization accuracy and efficiency, the parameter PHU was consid-ered as a constant value calibrated by the crop growth model in this study. Based on the above optimized parameters, the values of PHU were set as the medium of the optimization values (i.e., 3,244°C in 2007, 3,183°C in 2008 for maize, and 1,020°C for tomato).

Parameters
The predicted medium value and 95% confidence interval and the fitted values for each parameter of maize in  Table 2. The fitted values of parameters D x , D z , λ′, and Q were in the ranges of 0.4320-0.6073 cm 2 d −1 , 3.2086-3.6974 cm 2 d −1 , 0.0001 d −1 , and 0.1682-0.225 cm cm −2 d −1 , respectively. There was little discrepancy of the fitted parameter values with or without con-sideration of PHU in optimization. It indicated that the three steps' optimization scheme is a powerful tool for parameter calibration. Similar results were obtained for tomato (Table 3).
The fitted values of tomato root growth parameters D x , D z , λ′, and Q were 1.0448 cm 2 d −1 , 5.6928 cm 2 d −1 , 0.0001 d −1 , and 0.4526 cm cm −2 d −1 , respectively. With these fitted parameters, the bias values between simulated and observed for maize and tomato root length densities were 17.8 and −3.0%, respectively.
The simulated root length density distributions with fitted parameters at different time along the growing period were compared with the observed ones. As shown in Figure 3, a good agreement was obtained between the simulated and the observed root length density distributions on 5, 17, and 28 June 2007. However, the greatest discrepancy was found on 25 May 2007. The simulated root length density was higher than the measured values, especially in the topsoil of 0-10 cm. A similar result was obtained for maize in 2008 ( Figure 4) and tomato ( Figure 5). As shown in Figure 5, the greatest difference between the predicted and the measured root length densities was found in the topsoil (0-20 cm) on 24 DAT. Thereby, it can be concluded that the predicted root length density of the upper root zone was greater than the observed values during the early growing period. This may be attributed to the constant value of Q (Heinen et al., 2003). In this study, Q was considered as an optimized parameter, similar to Pronk et al. (2002) and Heinen et al. (2003). However, as a linkage between the diffusive root growth model and the crop growth model, Q can be transformed from the rate of root biomass input (g cm −2 d −1 ) Q M , which was estimated by the crop growth model: where R is the mean root radius (cm), ρ is the bulk density of root (g cm −3 ), d is the dry matter content of root (g dry g −1 fresh), L surface is the surface length occupied by a single plant (cm), and A is the sum of area where root input occurs (cm 2 ). Mollier et al. (2008) integrated the diffusive root growth model with the crop growth model, and the input rate of the root length density Q was calculated by the crop growth model predicted root biomass input rate Q M . Furthermore, Heinen et al. (2003) indicated that the root biomass input Q M changed during the growth stage. It is a very small value during the early growth stage and then increases rapidly after middle growth stage. Therefore, if the root biomass input variation during the growing period was considered in the modified diffusive root growth model, the predicted root length distribution would fit well with the observed values.
The RMSE values for maize and tomato root length density during 100 runs are presented in Figure 6a. The RMSE values ranged from 0.22 to 0.39 cm cm −3 , and the highest F I G U R E 6 The RMSE of the root length density for maize and tomato in optimization (a) with and (b) without considering potential heat unit (PHU) value was obtained for tomato root length density. This may be because root samples were smaller for tomato optimization. The number of samples for maize and tomato were 84 and 36, respectively. Without considering PHU, the RMSE values for maize in 2007 and 2008 were in the range of 0.22-0.27 cm cm −3 and 0.22-0.28 cm cm −3 , and that for tomato was 0.25-0.39 cm cm −3 , respectively (Figure 6b). Similar results for tomato were obtained by Heinen et al. (2003). They indicated that the RMSE values for numerical simulation of tomato root length density were in a range of 0.21-4.49 cm cm −3 , whereas they were about 0.10 cm cm −3 for maize. The bias of maize roots ranged from 17.7 to 24.5% in this study. The greater RMSE and bias for maize predicted by the semianalytical solution of the modified root growth model may be attributed to the two reasons. One is the measurement error. Dathe et al. (2014) indicated that several factors influenced experimentally obtained root distribution. When roots are washed from the surrounding soil, medium fine roots will get lost, yielding a lower total root mass than that actually present in the soil, especially in deeper soils. The other reason may be the limitation of the modified diffusive root growth model. In this study, the modified diffusive root growth model was developed with several assumptions, without considering the regulated or controlled factors, including soil temperature, soil water content, soil resistance to penetration, soil structure, soil chemical transport and advection condition, and soil microbiological conditions. In addition, Dathe et al. (2014) simulated the effect of water stress on potato root growth in diffusive process by using the root biomass input rate. A similar approach was adopted by Mollier et al. (2008). Nevertheless, although the semianalytical solution of the root diffusive model was easily used to describe the rooting pattern and distribution, it cannot handle complex problems such as the situa-tion that root input occurs inside the rooting medium (Heinen et al., 2003). If abovementioned factors were all integrated into the proposed root distribution function, the numerical model developed by coupling the crop growth to soil water and solute transport dynamics could be improved. Then, the modified diffusive root growth model is expected to perform better in simulating root length density distribution under a 2D domain.

CONCLUSIONS
To establish the process-based simulation tools for modeling water flow in the surface-soil-crop system, rooting system is of great importance to accurately estimate the temporal and spatial root water uptake patterns in root zone. In this study, integration of the diffusive root growth model and EPIC crop growth model was carried out by changing the boundary conditions. In the modified model, the diffusive root growth was constrained by the root depth simulated with the EPIC crop growth model. Then, the semianalytical solution of the modified diffusive root growth model was derived. Moreover, the modified root growth model was tested by using field experiment data of root length density for maize and tomato. A MATLAB program was developed to integrate the modified root diffusive model and the GA for facilitating the parameters optimization. With the optimized parameters, the modified root growth model performed well in predicting root length density, with RMSE and bias values ranging from 0.22 to 0.25 cm cm −3 and from −3.0 to 24.5%, respectively. It can be concluded that the modified root growth diffusive model can be used to simulate 2D root growth conditions during the crop growing period.