Modelling soil water dynamics and root water uptake for apple trees under water storage pit irrigation

Water storage pit irrigation is a new method suitable for apple trees. It comes with advantages such as water saving, water retention and drought resistance. A precise study of soil water movement and root water uptake is essential to analyse and show the advantages of the method. In this study, a mathematical model (WSPI-WR model) for 3D soil water movement and root water uptake under water storage pit irrigation was established based on soil water dynamics and soil moisture and root distributions. Moreover, this model also considers the soil evaporation, pit wall evaporation and water level variation in the pit. The finite element method was used to solve the model, and the law of mass conservation was used to analyse the water level variation. The model was validated by experimental data of the sap flow of apple trees and soil moisture in the orchard. Results showed that the WSPI-WR model is highly accurate in simulating the root water uptake and soil water distributions. The WSPI-WR model can be used to simulate root water uptake and soil water movement under water storage pit irrigation. The simulation showed that orchard soil water content and root water uptake rate centers on the storage pit with an ellipsoid distribution. The maximum distribution region of soil water and root water uptake rate was near the bottom of the pit. Distribution can reduce soil evaporation in the orchard and improve the soil water use efficiency in the middle-deep soil.


Introduction
Soil water movement and root water uptake are the key processes of the substance transport and energy transfer in farmland [1][2][3] .A quantitative study of soil water and root is the basis of formulating irrigation system and improving water use efficiency.Soil moisture and root interact restrict each other, root water uptake is the main method for soil moisture consumption in the root zone.On the other hand, soil moisture is the main factor that restricts root water uptake [4][5][6] .A reasonable distribution of soil moisture in the root zone can promote root water uptake and improve water use efficiency.
Water storage pit irrigation is a new method for simultaneously solving the problems of drought and soil erosion.Water storage pit irrigation with the advantages of water saving, water retention, drought resistance and soil and water conservation is suitable for apple trees [7][8][9][10] .In the field, several small cylindrical storage pits (pit diameter 25-30 cm and pit depth 40-60 cm) are dug along half of the radius of the apple tree canopy.During irrigation, water flows into the annular groove through the pipeline or field channel, then enters the pits and infiltrates into the rhizosphere soil along the pit wall.During water storage pit irrigation, the soil infiltration interface changes.The water content in middle-deep soil is high, whereas that in topsoil is low.Water content is different in surface, sprinkler and micro irrigations.The distribution of soil moisture under water storage pit irrigation can reduce ground evaporation between plants, boost the growth of roots in middle-deep soil, change the distribution of root system and promote the root water uptake [11][12][13][14] .
The mathematical model of soil water movement under root water uptake can be built by adding the root water uptake item to the right end of the basic equation of soil water movement, while considering boundary conditions such as irrigation, evaporation and rainfall.In general, the root water uptake model has two types, namely, micro and macro models [15] .The macro model is more common than the micro model, and the macro model can be subdivided into two categories [16] : the first category is built by distributing plant transpiration to different soil depths according to root density, and the second category is built on the basis of the spatial distribution of root water uptake intensity, which is obtained by the soil water dynamics.By contrast, the first category is easier and more popular than the second category.In recent years, many mathematical models for orchard soil water movement and root water uptake have been developed.Most of the models are 1D [17][18][19] or 2D [20][21][22] , only a few are 3D [23] .Numerical solutions, such as finite difference, finite element and finite volume methods, were adopted to solve the models [24][25][26][27][28] .Finite element method is widely used because of its favourable adaptability to irregular boundary conditions.Hydrus 2D/3D is developed based on the finite element method [29] .
This software performs well in simulating water-heat-solute transport in farmland [30,31] .However, Hydrus 2D/3D is a closed software with weak capability to handle complex boundary conditions.The water level in the pit, infiltration head and infiltration area that act on the pit wall infiltration interface during water storage pit irrigation declines with infiltration time.Moreover, the infiltration head has a wide range of variation.The soil water movement under water storage pit irrigation is 3D.Thus, factors such as root water uptake and complex boundary conditions must be considered in the simulation.
This study aims to establish the soil water movement and root water uptake model (WSPI-WR model) for apple trees under water storage pit irrigation, calibrate and validate the model by using the measured data in the orchard and show the spatial distribution characteristics of soil moisture and root water uptake through water storage pit irrigation method.

Study area
The field study was conducted at the Fruit Research Institute of Shanxi Academy of Agricultural Sciences (latitude 37°23′N, longitude 112°32′E) from May 23 rd to October 10 th , 2013.The region is located in the southwest of Taigu County, Shanxi Province.The region has a typical temperate continental climate with a long-term annual mean air temperature and precipitation of 9.8 °C and 462.9 mm, respectively.The frost period is from early October to mid-April of the following year, and the frost-free period spans 175 d.The area of the region is 450 ha, and the soil type is mainly loam.The apple trees (Malus pumila Mill) discussed in this study were planted with a row and plant spacing of 4 and 2 m, respectively.

Experimental design and test items 2.2.1 Experimental design
The field engineering of water storage pit irrigation included storage pits, annular groove, conveyance ditch (or a delivery pipe) and facilities to firm the pit wall (Figure 1).The storage pit (diameter 30 cm and pit depth 60 cm) was arranged under the crown and 60 cm away from the trunk.Four storage pits are under each apple tree.The bottom of the pit was waterproofed to reduce deep percolation, and metal mesh was used to firm the pit wall.The annular groove (depth 20 cm and width 30 cm) was used to connect all the storage pits around the tree.The delivery pipe is a fixed field pipe that connects the irrigation backbone system with the annular groove.Water flowed into the annular groove through the delivery pipe during irrigation and then entered the pits and infiltrated into the soil through the pit wall.The irrigation events were conducted on May 23 rd and September 11 st , 2013.The irrigation quota per tree was 300 L.

Measurement of soil evaporation
The orchard soil evaporation under water storage pit irrigation, included the pit wall and surface evaporations, was measured using micro-lysimeters [32] .Micro-lysimeters were arranged along the pit wall.Figure 3 depicts that each storage pit requires three micro-lysimeters with depths of 10, 30 and 50 cm.Figure 4 displays that micro-lysimeters for surface evaporation were arranged under the crown 50, 100 and 150 cm away from the trunk.The sap flow rate of the apple tree was measured by using a thermal diffusion probe sap flow meter (TDP, Dynamax Company, USA).Waterproof daub was applied at the junction of the probe and bark, and aluminium film and silver paper were wrapped around the trunk to eliminate the influence of the sun on the thermal probe.The end of the probe was connected to a data collector, and data were collected and recorded every 15 s and 30 min, respectively.

Measurement of root length density
Root length density was measured through the root drilling method (drill diameter 7 cm) in September 2013.The trunk was considered as the centre, and the measuring points were placed 30, 60, 90 and 120 cm away from the trunk (Figure 2).The soil samples were collected at a depth interval of 20 cm.The total sampling depth was 200 cm.Soil samples with roots were removed, immediately placed in sealed plastic bags and then transported back to the laboratory.The roots were tiled in organic glass dishes (20 cm × 25 cm) after washing and separation using clean water and tweezers, and no overlap existed among the roots.Subsequently, the roots were scanned using an Epson Perfection 4870 scanner.Finally, software WinRHIZO (Regent Instruments, Canada) was used to analyse the scanned results and obtain the root length density.

Measurement of leaf area index (LAI)
Leaf area index (LAI) was measured monthly from May 2013 to September 2013 using a LAI-2200 vegetation canopy analyser (LI-COR, USA).

Collection of meteorological data
Meteorological data, including precipitation, atmospheric temperature, soil temperature, wind speed, wind direction, pressure atmosphere, relative humidity and solar radiation, were collected from the Adcon_Ws weather station every 15 min.Figure 5 presents the variation of precipitation and reference crop evapotranspiration (ET 0 ) during the experiment.

Soil water movement and root water uptake model under water storage pit irrigation 2.3.1 Governing equation
Three-dimensional features are obvious in soil water movement under water storage pit irrigation.The computational region was determined according to the symmetrical placement of storage pits (Figure 6).The 3D equation of water movement that considers root water uptake was expressed as Equation ( 1). (

) [ ( ) ] [ ( ) ] [ ( ) ]
(1) where, h is the soil water pressure, cm; K(h) is the unsaturated hydraulic conductivity, cm/h; x, y and z are the space coordinates, cm; t is time, h; and S is the root water uptake rate, L/h.

Figure 6 Computational region of water storage pit irrigation
In order to run model, the hydraulic parameters θ s , θ r , K s , α and n were measured by specific experimental methods.First, a soil profile with length 2.0 m, width 1.0 m and depth 2.0 m was excavated at the experimental site.Then, according to the soil texture, color, tightness, etc., the soil profile was divided into 5 layers 0-40 cm, 40-70 cm, 70-120 cm 120-170 cm and 170-200 cm, and six soil samples were collected in each layer of soil using soil cutting ring.Then the soil samples were brought back to the laboratory and the soil water content corresponding to different pressures was measured using a 1500F1 pressure plate extractor (Soilmoisture Equipment Corp., Goleta, CA, US).According to the measured soil water content and pressure data, the equation ( 2) was fitted with MALAB2016 and the values of parameters θ s , θ r , α and n were obtained.The saturated hydraulic conductivity K s was measured by a TST-55 permeameter (Nanjing Soil Instrument Co., Ltd., China).Table 1 lists the soil hydraulic parameters.The root water uptake term in Equation ( 1) is presented in Equation (4) [23] .max ( , , , ) ( ) ( , , , )

Table 1 Soil hydraulic parameters of experimental plot soil
where, γ(h) is the coefficient of water stress (Equation ( 5)), and S max is the maximum root water uptake rate (L/h) (Equation ( 6)) [23] .

( ) 0 h h h h h h h h h h γ h h h h h h h
where, h 0 , h 1 , h 2 , h 3 are thresholds of soil water potential that affects root water uptake.For the apple tree, h 0 , h 1 , h 2 , h 3 are 0, −0.1, −10 and −150 m, respectively [18] .
where, T pot is the potential transpiration intensity (cm/d) (Equation ( 8)) [15] ; β(x, y, z) is the function of root distribution shape; X m , Y m , Z m are the maximum extension depths (cm) of root in x, y and z directions, respectively; A xy is the area of computational region on the surface, cm 2 ; and p x , p y , p z , x * , y * , z * are the fitting parameters of Equation ( 7) using the measured root length density (Table 2).
where, LAI is the leaf area index; a, b, k and A are constants with values of 0.2432, 0.448, 0.402 and 0.976, correspondingly [20] ; and ET 0 is the reference crop evapotranspiration calculated using the Penman-Monteith formula [33] .2.3.2Initial condition h(x, y, z, t) = h 0 (x, y, z) t = 0 (9) where, h 0 (x, y, z) is the suction head that corresponds to the initial water content, cm.

Boundary conditions 1) Surface boundary (AQJ, MEDL and PIRFN)
For surface boundary, the upper boundary is considered the second boundary when precipitation intensity is less than infiltration intensity; the upper boundary is transformed from the second boundary into the first boundary when precipitation intensity is higher than infiltration intensity, and excess rainfall becomes surface runoff and flows into the water storage pit (Equations ( 10) and ( 11)).

( ( ) ( )) h K h K h P P i z
where, P is the net precipitation intensity, cm/h; i is the infiltration intensity, cm/h; z u is the z coordinate value of the surface boundary and h s is the surface ponding depth, cm.The upper boundary is the second boundary during evaporation; the upper boundary is transformed into the first boundary when the ground becomes dry.
( , , ) ( , , ) where, e s is the surface evaporation intensity, cm/h; and h d is the suction head when the ground is dry, cm.
2) Lateral wall boundary of the annular groove (EFNM and

QJIP)
The lateral wall boundary of the annular groove is the second boundary during evaporation; the boundary is transformed into the first boundary when the ground becomes dry.
3) Pit wall boundary (IRFGOH) The pit wall boundary is the first boundary during irrigation.h(x, y, z) = h i (16) where, h i is the pressure head of each point on the pit wall, which has the value of the distance from the point to the water level in the pit.
The pit wall boundary is considered the second boundary when the infiltration in the pit is completed.( )( )

4) Pit bottom boundary (HOG)
The pit bottom boundary is impermeable.

5) Lower boundary (BCK)
The lower boundary is considered a free drainage boundary given the deep calculation depth and bury of groundwater.
h(x, y, z) = h 0 (19) 6) Boundary ABCD Boundary ABCD is considered a zero flux plane because of its symmetry.
7) Boundary ABKL Boundary ABKL is also considered a zero flux plane because of its symmetry.

8) Boundary CDLK
No horizontal water exchange is assumed for boundary CDLK because of the large calculation area.

Model solving
This study uses the Galerkin finite element method to solve the 3D Richards equation because the finite element method is easy to use under complex boundary conditions.The flow region was divided into a network of tetrahedral elements (Figure 7).Equation ( 1) is solved by using the finite element method with the Galerkin weighted residual method. where [ ] 36 where, where, i, j, k and m are the node numbers of the tetrahedral element e, V e is the volume of element e, and K is the average hydraulic conductivity over element e. L e is the area of flux boundaries in tetrahedral elements (cm 2 ), and q is the flux on the flux boundary of the tetrahedral element (cm 3 /min).Equation ( 28) is only significant to the flux boundary.The internal and Dirichlet boundary nodes are 0, and the boundary of the tetrahedron element is jkm, which is the area of the interface jkm.

Figure 7 Mesh generation results
An implicit backward difference scheme is used for the time term in Equation ( 23).
where, j 0 +1 denotes the current time level at which the solution is considered, j 0 refers to the previous time level and Δt j0 = t j0+1 -t j0 .The 'mass conservation' method [34] was adopted to handle the water content term ( t θ Δ Δ ); the first item in Equation ( 30) was divided into the following parts in the iterative process: The first term on the right-hand side can be expressed in terms of pressure head; thus, Equation ( 31) becomes the following one: (32) Equation ( 33) will be obtained by substituting Equation (32)  into Equation (30), as follows:

.5 Solution of water depth changing process in the water storage pit
The water level in the pit continuously declines during water storage pit irrigation, and the boundary conditions of the pit wall continuously change with time.Therefore, an accurate calculation of water level-changing process in the pit is the key step in the numerical simulation of soil water movement.According to the principle of mass conservation, the amount of water that reduces in the pit is equal to the amount of water that infiltrates into the soil by the pit wall during irrigation (Equation ( 34)).
where, Vol is the amount of irrigation water in a single pit, cm 3 ; Vol t is the amount of water in a single pit at time t (cm 3 ) (Equation (36)); and Q′ n is the boundary section flow controlled by the node n on the pit wall boundary, which is the product of the flow of the boundary node and the area of the control boundary section.N equations will be obtained by substituting the head value into Equation (33) when the real-time head value {h} of each node t is calculated.The node n on the boundary of the pit wall corresponds to the nth equation of N equations, in which Q n (the nth of {Q}) is unknown, and the remaining nodes are known.Q′ n can be calculated using Equation (35) as follows: where, N is the total number of nodes in the calculation region.Vol t = πr 0 H t (36) where, H t is the water depth in the pit at moment t, cm.Equation (37) will be obtained by substituting Equations ( 35) and (36) into Equation (34), as follows: The changing process of the water level in the pit can be obtained using Equation (37).The boundary point of the pit belongs to the first boundary when this point is below the water level; the boundary point belongs to the second boundary when this point is above the water level.
where, V S is the simulated value, V R is the measured value, and l is the number of measuring points.

Results and discussions
3.1 Comparison of measured and simulated root water uptake by the single apple tree Figure 8 illustrates the comparison of measured and simulated root water uptake.The measured root water uptake was expressed by the stem sap flow of a single apple tree obtained by a stem flow meter (TDP).The simulated root water uptake was obtained using the WSPI-WR model.The measured and simulated root water uptake demonstrated the same variation trend, and the simulated root water uptake agreed well with the measured water uptake.Figure 9 depicts the correlation analysis of measured and simulated root water uptake.A significant correlation was found between the measured and simulated root water uptake obtained using the WSPI-WR model.The determination coefficient was 0.93.The WSPI-WR model has RMSE, MRE, ARE and MAE values of 237 cm 3 /d, 9.88%, 5.42% and 205 cm 3 /d, respectively, thereby indicating the high accuracy in simulating root water uptake.

Comparisons of measured and simulated soil water content dynamics
Figure 10 presents the comparisons between measured and simulated soil water content.The measured and simulated soil water content had the same variation trend.
Figure 11 demonstrates the correlation analysis of measured and simulated soil water content.A significant correlation was found between the measured and simulated soil water content obtained by using the WSPI-WR model.The determination coefficient was 0.85.The RMSE, MRE, ARE and MAE values were 0.014 cm 3 /cm 3 , 14.26%, 3.85% and 0.014 cm 3 /cm 3 , correspondingly.
The WSPI-WR model performed well in simulating soil water content under the condition of root water uptake.

3D Distribution of Soil Moisture Content Simulated by the WSPI-WR Model
Figure 12 illustrates the 3D distribution of soil moisture content in 1 d, 3 d and 5 d after irrigation that simulated by the WSPI-WR model.The soil moisture contents in the topsoil (0-20 cm) were low, and soil moisture contents in the middle-deep soil (20-140 cm) were relatively high.The soil water centred ellipsoidal on the water storage pit, and the soil moisture content was getting lower when further distant from the pit.This condition is due to the irrigation water infiltrated into the soil through the pit wall and then moved around the pit wall under the influence of soil matrix and gravitational potentials.This kind of soil moisture distribution can effectively reduce soil evaporation and increase soil water use efficiency.Furthermore, the soil water in 1 day after irrigation was mainly distributed in the depth of 20-130 cm and then distributed in the range of 20-140 cm in 3 d after irrigation.The soil water distribution in 5 d after irrigation became uniform due to redistribution, root water uptake and evaporation.Moreover, no obvious wetting range was observed in the distribution.The results showed that the water movement model built in this study performed well in simulating the soil water dynamics under water storage pit irrigation.

Root water uptake distribution simulated by the WSPI-WR model under different pit depths
Figure 13 depicts the root water uptake distribution simulated by the WSPI-WR model 1 d, 3 d and 5 d after irrigation.The root water uptake under water storage pit irrigation mainly distributed in the range of 0-130 cm.This distribution centred on the storage pit and decreased along the pit wall, which was consistent with the distribution of soil moisture under water storage pit irrigation.The maximum region of root water uptake rate was near the bottom of the pit.The distribution of root water uptake rate 1 d, 3 d and 5 d after irrigation became increasingly uniform with time, and the root water uptake rate gradually decreased because the soil moisture becomes uniform and reduces with time under the influence of root water uptake and soil water redistribution.A mathematical model (WSPI-WR model) for a 3D soil water movement and root water uptake under water storage pit irrigation was established based on soil water dynamics, and soil evaporation, pit wall evaporation and water level variation in the pit had been considered.The model was solved through the finite element method.The experimental data of the sap flow of apple trees and soil moisture in the orchard were used for validating the model.The results showed that the WSPI-WR model performed well in simulating root water uptake and soil water movement.This model can be used to simulate the characteristics of root water uptake and soil water movement under water storage pit irrigation.
The simulation results of the WSPI-WR model showed that the orchard soil water content and root water uptake rate under water storage pit irrigation centred on the storage pit with an ellipsoid distribution.The water content of the surface soil was low, and the moisture content of the middle-deep soil was high.The root water uptake of apple trees was mainly distributed in the range of 0--130 cm, centred in the water storage pit and decreased around the pit wall.The spatial distribution of root water uptake rate was consistent with that of soil moisture, thereby reducing soil evaporation and improving the soil water use efficiency in the middle-deep soil.

Figure 1
Figure 1 Field engineering diagram of water storage pit irrigation 2.2.2 Measurement of soil moisture contentSoil moisture content was measured using the TRIME-PICO IPH measuring system every 5-10 d.An additional measurement should be performed after rain or irrigation.Figure2 illustrated

Figure 2
Figure 2 Layout of the measuring points of the soil moisture and root 2.2.3 Measurement of soil evaporationThe orchard soil evaporation under water storage pit irrigation, included the pit wall and surface evaporations, was measured using micro-lysimeters[32] .Micro-lysimeters were arranged along the pit wall.Figure3depicts that each storage pit requires three micro-lysimeters with depths of 10, 30 and 50 cm.Figure4displays that micro-lysimeters for surface evaporation were arranged under the crown 50, 100 and 150 cm away from the trunk.

Figure 3 Figure 4
Figure 3 Diagram of the pit wall evaporation test

Figure 5
Figure 5 Variation of precipitation and ET 0 during the experiment j m b c d bb cc dd bb cc dd bb cc dd b b c c d d b c d b b c c d d b b c c d d K A V b b c c d d b b c c d d b c d b b c c d d b b c c d d b b c c d d b

Figure 8 Figure 9 Figure 10 Figure 11 Figure 13
Figure 8 Comparison of measured and simulated root water uptake 3