Modeling the Effects of Spatial Variability of Irrigation Parameters on Border Irrigation Performance at a Field Scale

The interaction between surface and subsurface water flows plays an important role in surface irrigation systems. This interaction can effectively be simulated by the physical-based models, which have been developed on the basis of the numerical solutions to the Saint-Venant and Richards’ equations. Meanwhile, the spatial variability of field physical properties (such as soil properties, surface micro-topography, and unit discharge) affects the interaction between surface and subsurface water flows and decreases the accuracy of simulating surface irrigation events at large scales. In this study, a new numerical methodology is developed based on the physical-based model of surface irrigation and the Monte Carlo simulation method to improve the modeling accuracy of surface irrigation performance at a field scale. In the proposed numerical methodology, soil properties, unit discharge, surface micro-topography, roughness, border length, and the cutoff time for the unit discharge are used as the stochastic parameters of the physical-based model, while field slope is assumed as the constant value because of the same field tillage and management conditions at a field scale. Monte Carlo simulation is used to obtain the stochastic parameter sample combinations of the physical-based model to represent the spatial variability of field physical properties. The updated stochastic simulation model of surface micro-topography, which is developed to model the spatial distribution of surface elevation differences (SED), is used to obtain the surface micro-topography samples at a field scale. Compared with the distributed-parameter modelling methodology and the field experimental data, the proposed numerical methodology presents the better simulation performance.


Introduction
Representing surface irrigation processes and evaluating surface irrigation performance (e.g., distribution uniformity, irrigation efficiency) at large spatial scales is a major challenge in farmland hydrology. Difficulties arise mainly from the spatial variability of field physical properties which create satisfactorily capture the field physical characteristics and display a better distribution uniformity of model parameters than the distributed-parameter modeling. However, the stochastic parameter modeling method has been rarely used in border irrigation simulation of the field scale.
The purpose of this paper is to develop a stochastic parameter model for evaluating the impact of the spatial variability of field physical properties on border irrigation performance at a field scale. The proposed stochastic parameter model is developed on the basis of a hybrid coupled model of border irrigation [23] and Monte Carlo simulation framework. The proposed stochastic parameter model is applied and tested in the Mawan irrigation district in Guangrao Country, Shandong Province, China, and the Yehe irrigation district in Shijiazhuang City, Hebei Province, China.

Site Description
The first field experiment was conducted in the spring of 2012 at the Mawan irrigation district (MW) in Dongying City, Shandong Province, China and the second field experiment was conducted in the fall of 2012 at the Yehe irrigation district (YH) located in Shijiazhuang City, Hebei Province, China, as shown in Figure 1. The experimental field areas are 76.2 ha 2 in MW and 43.4 ha 2 in YH, respectively. MW is located in the coastal area of the lower reaches of the Yellow River, which has a typical warm temperate zone continental monsoon climate, with an average annual rainfall of 537.4 mm. The soils are alluvial sediments with primary textures of silt loam and sandy loam. YH is located in the North China Plain, which has a typical temperate zone monsoon climate, with the average annual rainfall of 444.1 mm. The experimental soil is classified as sandy loam in the root zone. The purpose of this paper is to develop a stochastic parameter model for evaluating the impact of the spatial variability of field physical properties on border irrigation performance at a field scale. The proposed stochastic parameter model is developed on the basis of a hybrid coupled model of border irrigation [23] and Monte Carlo simulation framework. The proposed stochastic parameter model is applied and tested in the Mawan irrigation district in Guangrao Country, Shandong Province, China, and the Yehe irrigation district in Shijiazhuang City, Hebei Province, China.

Site Description
The first field experiment was conducted in the spring of 2012 at the Mawan irrigation district (MW) in Dongying City, Shandong Province, China and the second field experiment was conducted in the fall of 2012 at the Yehe irrigation district (YH) located in Shijiazhuang City, Hebei Province, China, as shown in Figure 1. The experimental field areas are 76.2 ha 2 in MW and 43.4 ha 2 in YH, respectively. MW is located in the coastal area of the lower reaches of the Yellow River, which has a typical warm temperate zone continental monsoon climate, with an average annual rainfall of 537.4 mm. The soils are alluvial sediments with primary textures of silt loam and sandy loam. YH is located in the North China Plain, which has a typical temperate zone monsoon climate, with the average annual rainfall of 444.1 mm. The experimental soil is classified as sandy loam in the root zone.

Data Collection
The border fields with 272 plots in MW and 154 plots in YH were used to irrigate the winter wheat by the Yellow River water and Gangnan reservoir water, respectively. Before the beginning of border irrigation, soil samples were collected at four depths (0 cm to 20 cm, 20 cm to 40 cm, 40 cm to 70 cm, 70 cm to 100 cm) from the middle of border field with 100 plots selected randomly to analyze soil bulk density, soil particle size distribution, and initial soil water content. After the end of border irrigation with a period of 48 h, soil samples were collected at four depths (0 cm to 20 cm, 20 cm to 40 cm, 40 cm to 70 cm, 70 cm to 100 cm) from the top, middle, and bottom of the border fields, with 100 plots selected randomly to analyze soil water content.
The relative surface elevation values of border field with 100 plots selected randomly were observed using a water level gauge, and the space between two adjacent measured points was 10 m for all experiments. Three lines were set at the left, middle, and right sides along the border width ( Figure 2). The border length and width were measured using an engineering measure ruler. The observational tests for the surface flow advance processes were carried out at the beginning of surface irrigation ( Figure 2). The surface flow depth was measured using water depth measuring devices [24], which were placed at every observation point before the irrigation was initiated. The water depth was used to estimate Manning's roughness coefficient n [25]. The recorded data were then transferred to a computer at the end of the experiments. The border fields with 272 plots in MW and 154 plots in YH were used to irrigate the winter wheat by the Yellow River water and Gangnan reservoir water, respectively. Before the beginning of border irrigation, soil samples were collected at four depths (0 cm to 20 cm, 20 cm to 40 cm, 40 cm to 70 cm, 70 cm to 100 cm) from the middle of border field with 100 plots selected randomly to analyze soil bulk density, soil particle size distribution, and initial soil water content. After the end of border irrigation with a period of 48 h, soil samples were collected at four depths (0 cm to 20 cm, 20 cm to 40 cm, 40 cm to 70 cm, 70 cm to 100 cm) from the top, middle, and bottom of the border fields, with 100 plots selected randomly to analyze soil water content.
The relative surface elevation values of border field with 100 plots selected randomly were observed using a water level gauge, and the space between two adjacent measured points was 10 m for all experiments. Three lines were set at the left, middle, and right sides along the border width ( Figure  2). The border length and width were measured using an engineering measure ruler. The observational tests for the surface flow advance processes were carried out at the beginning of surface irrigation ( Figure 2). The surface flow depth was measured using water depth measuring devices [24], which were placed at every observation point before the irrigation was initiated. The water depth was used to estimate Manning's roughness coefficient n [25]. The recorded data were then transferred to a computer at the end of the experiments.

Statistical Analysis of Field Physical Properties at a Field Scale
Border irrigation performance at a field scale is mainly affected by the spatial variability of field physical properties and field management measures. According to the same field tillage and management conditions, the frequency distribution histograms of border length in MW and YH are shown in Figure 3. It can be seen that border length in MW and YH presents the normal distribution characteristic and can be fitted by the normal distribution function. The mean values and standard deviations for the normal distribution of border length in MW and YH are 120, 33.6 m and 200, 30.2 m, respectively. Similarly, two parameters of the fitting distribution of border width in MW and YH are 9, 1.98 m and 6, 1.44 m, respectively.

Statistical Analysis of Field Physical Properties at a Field Scale
Border irrigation performance at a field scale is mainly affected by the spatial variability of field physical properties and field management measures. According to the same field tillage and management conditions, the frequency distribution histograms of border length in MW and YH are shown in Figure 3. It can be seen that border length in MW and YH presents the normal distribution characteristic and can be fitted by the normal distribution function. The mean values and standard deviations for the normal distribution of border length in MW and YH are 120, 33.6 m and 200, 30.2 m, respectively. Similarly, two parameters of the fitting distribution of border width in MW and YH are 9, 1.98 m and 6, 1.44 m, respectively. The border fields with 272 plots in MW and 154 plots in YH were used to irrigate the winter wheat by the Yellow River water and Gangnan reservoir water, respectively. Before the beginning of border irrigation, soil samples were collected at four depths (0 cm to 20 cm, 20 cm to 40 cm, 40 cm to 70 cm, 70 cm to 100 cm) from the middle of border field with 100 plots selected randomly to analyze soil bulk density, soil particle size distribution, and initial soil water content. After the end of border irrigation with a period of 48 h, soil samples were collected at four depths (0 cm to 20 cm, 20 cm to 40 cm, 40 cm to 70 cm, 70 cm to 100 cm) from the top, middle, and bottom of the border fields, with 100 plots selected randomly to analyze soil water content.
The relative surface elevation values of border field with 100 plots selected randomly were observed using a water level gauge, and the space between two adjacent measured points was 10 m for all experiments. Three lines were set at the left, middle, and right sides along the border width ( Figure  2). The border length and width were measured using an engineering measure ruler. The observational tests for the surface flow advance processes were carried out at the beginning of surface irrigation ( Figure 2). The surface flow depth was measured using water depth measuring devices [24], which were placed at every observation point before the irrigation was initiated. The water depth was used to estimate Manning's roughness coefficient n [25]. The recorded data were then transferred to a computer at the end of the experiments. Border irrigation performance at a field scale is mainly affected by the spatial variability of field physical properties and field management measures. According to the same field tillage and management conditions, the frequency distribution histograms of border length in MW and YH are shown in Figure 3. It can be seen that border length in MW and YH presents the normal distribution characteristic and can be fitted by the normal distribution function. The mean values and standard deviations for the normal distribution of border length in MW and YH are 120, 33.6 m and 200, 30.2 m, respectively. Similarly, two parameters of the fitting distribution of border width in MW and YH are 9, 1.98 m and 6, 1.44 m, respectively.  The unit discharge is calculated on the basis of border width and inflow observed at the field experiments. The unit discharge presents the obvious statistical characteristics of normal distribution by the Kolomogorov-Smirnov test and the corresponding parameters and normal test results are shown in Table 1. The same field tillage and management methods are conducted at a field scale. Therefore, the standard deviations of surface elevation differences from the fields of 100 plots are assumed as the same in this study. The surface elevation differences (SED) values are checked for normality by the Kolmogorov-Smirnov test and a normal distribution with a significance level of 0.05 is obtained. The mean values and coefficients of the standard deviation for the two experiments are 4.06, 0.12 cm and 3.5, 0.21 cm, respectively, calculated through the method of Bai et al. [26]. The spatial variability of the SED values is analyzed through the experimental semivariogram function, which can be fitted by the spherical semivariogram model [26]. The parameters of the spherical semivariogram model are calculated from the two sets of field-observed data ( Table 2). These parameters are used for the stochastic modeling of surface micro-topography. In addition, the field slopes in MW and YH are about 0.6‰ and 2.4‰, respectively. The probability distribution characteristics of soil sandy loam, silt loam, clay loam, bulk density and organic carbon content can be obtained on the basis of the field-observed data. Each variable of soil properties in MW and YH is checked for normality by the Kolmogorov-Smirnov test. A log transformation is conducted if the statistical distribution of the raw data is highly skewed. The primary statistical parameters, such as the mean, maximum, minimum, and standard deviation values, are calculated. The statistical results of the soil properties, each variable of which represents the normal or log-normal distribution at a significance level of 0.05, are shown in Table 3. The unit discharge is calculated on the basis of border width and inflow observed at the field experiments. The unit discharge presents the obvious statistical characteristics of normal distribution by the Kolomogorov-Smirnov test and the corresponding parameters and normal test results are shown in Table 1. Table 1. Statistical characteristics and normal test results of unit discharge of border field at a field scale. The same field tillage and management methods are conducted at a field scale. Therefore, the standard deviations of surface elevation differences from the fields of 100 plots are assumed as the same in this study. The surface elevation differences (SED) values are checked for normality by the Kolmogorov-Smirnov test and a normal distribution with a significance level of 0.05 is obtained.

Experiment Site Minimum (L/(s·m)) Maximum (L/(s·m)) Mean (L/(s·m)) S.D. K-S Test
The mean values and coefficients of the standard deviation for the two experiments are 4.06, 0.12 cm and 3.5, 0.21 cm, respectively, calculated through the method of Bai et al. [26]. The spatial variability of the SED values is analyzed through the experimental semivariogram function, which can be fitted by the spherical semivariogram model [26]. The parameters of the spherical semivariogram model are calculated from the two sets of field-observed data ( Table 2). These parameters are used for the stochastic modeling of surface micro-topography. In addition, the field slopes in MW and YH are about 0.6 and 2.4 , respectively. Table 2. The parameter values of the pherical semivariogram model for surface relative elevation at a field scale. The probability distribution characteristics of soil sandy loam, silt loam, clay loam, bulk density and organic carbon content can be obtained on the basis of the field-observed data. Each variable of soil properties in MW and YH is checked for normality by the Kolmogorov-Smirnov test. A log transformation is conducted if the statistical distribution of the raw data is highly skewed. The primary statistical parameters, such as the mean, maximum, minimum, and standard deviation values, are calculated. The statistical results of the soil properties, each variable of which represents the normal or log-normal distribution at a significance level of 0.05, are shown in Table 3. According to the results of Agricultural Drought Research Center in the United States, the surface roughness coefficient can be determined by the crop growth situations and soil properties. The surface roughness coefficient values are checked for normality by the Kolmogorov-Smirnov test and a normal distribution with a significance level of 0.05 is obtained. The mean values and coefficients of the standard deviation for the two experiments are shown in Table 4.

Numerical Modeling Methodology
Stochastic parameter model for surface irrigation is usually developed on the basis of deterministic irrigation model and stochastic parameter modeling method. In this section, we mainly introduce the physical-based surface irrigation model proposed by Dong et al. [23], which is a surface irrigation system model for the border irrigation simulation and implemented by coupling a one-dimensional surface flow model that uses the complete hydrodynamic form of the Saint-Venant equation with a one-dimensional subsurface flow model, and Monte Carlo simulation method, which is used to obtain the stochastic parameter values of the surface irrigation model.

Coupled Modeling of Surface-Subsurface Water Flow for Border Irrigation
Numerical simulation can be used to quantify surface irrigation application efficiency and surface irrigation distribution uniformity from irrigating fields and to design surface irrigation scheme based on controlling variables such as soils, border length, unit discharge, and surface micro-topography. In the proposed paper, the equations that govern the flow of water for surface irrigation modeling are presented below and the processes that describe water movement in the surface flow and subsurface flow domains are discussed.
Border irrigation water flow can be characterized by the complete hydrodynamic form of the one-dimensional Saint-Venant equations, which can accurately describe the non-steady water flow of surface irrigation. The continuity and momentum equations in the complete hydrodynamic model are taken as the governing equations: ∂U ∂t where t is the time (s), x is the spatial coordinate (m), U is the vector of the conserved dependent variables, F is the physical flux, and S is the source vector that comprises bottom elevation vector S 1 , bed roughness vector S 2 and infiltration vector S 3 , which are expressed as: where h s is the water depth (m), q denotes the unit discharge along coordinate direction x (m 3 /s/m), g is the acceleration due to gravity (m/s 2 ), z is the bottom elevation (m), n denotes the Manning's roughness coefficient (m 1/6 ), and u is the vertically averaged flow velocity (m/s), I represents the infiltration rate (m/s).
Infiltration and soil water dynamics are simulated using the one-dimensional Richards equation, which describes vertical one-dimensional and non-hysteretic flow in a variably saturated porous medium where h is the pressure head (m), K(h) is the vertical hydraulic conductivity (m/min), C(h) = dθ/dh is the hydraulic capacity (−), z denotes the vertical coordinate(m), and t is the time (min). The curve of soil water retention characteristics is described using the van Genuchten (1980) equations: and where θ s is the saturated soil water content (−), θ r is the residual soil water content (−), α (1/m), n (−) and m (−) denote the soil water retention curve parameters, m = 1 − 1/n, and K S is the saturated hydraulic conductivity (m/min), and The hydraulic capacity is derived from van Genuchten's model:

Numerical Solution
The various components of the mass and momentum conservation equations are solved by a hybrid numerical method to improve simulation accuracy and reduce computational effort. The details on the numerical solution of the complete hydrodynamic form of the one-dimensional Saint-Venant equations were previously discussed by Zhang et al. [20] and Dong et al. [23]; hence, only a brief description of it is provided.
At the time discretization used for the governing equations, the Crank-Nicolson implicit time scheme is used for ∂U/∂t, ∂F/∂x, and S 2 at an unknown time (w + 1) and a known time w. Then, S 1 and S 3 are discretized using the explicit time scheme at a known time w. At the space discretization of the governing equations, the finite difference scheme is used to discretize the space derivative of the physical flux linear approximation generated by the implicit time scheme. The advection upstream splitting method is used to discretize the space derivative of the advection physical flux, S 2 and S 3 are then calculated according to the corresponding space-node values. After spatio-temporal discretization processes of the governing equation of the one-dimensional complete hydrodynamic model, the discretization equations of mass conservation and momentum conservation are distinguished shown in Dong et al. [23]. The final format of surface flow equations can be expressed by: where B is the coefficient matrix consisting of the coefficient of ∆U of cell (j − 1), j, (j + 1), respectively; at a known time w; B MA is the coefficient matrix consisting of the coefficient of h s of cell (j − 1), j, (j + 1), respectively, at a known time w; B MO is the coefficient matrix consisting of the coefficient of u of cell (j − 1), j, (j + 1), respectively, at a known time w; X is the matrix of ∆U of each cell at a known time w; X MA is the matrix of h s of each cell; at a known time w; X MO is the matrix of u of each cell at a known time w; D is the coefficient matrix of numerical flux vector of interface of cell j, bed roughness vector S 2 , infiltration vector S 3, and the averaged form of bottom elevation vector S 1 of cell j at a known w; D MA is the coefficient matrix about h s and u at a known w; and D MO is the coefficient matrix of numerical fluxes of interface of cell j, bed roughness, infiltration, and bottom elevation at a known w.
A finite volume method with fourth-order accuracy is used for the solution of the one-dimensional Richards equation in space, and the implicit second-order accurate finite difference scheme is used for the solution in time [23]. Then, the Picard iteration scheme [27] is applied to cope with the nonlinearity of the spatio-temporal discretization equation. Finally, the linearized algebraic equation obtained by spatio-temporal discretization processes is numerically solved. The discretized subsurface flow equation is as follows: where h k+1,m i is the approximate value of the h of the ith cell at the mth discrete time level of the (k + 1)th time step, and R P is the residual associated with the Picard iteration. If R k+1,m P ≤ 10 −5 and the iteration stops, the soil pressure heads of every spatial point at the (k + 1)th time step are calculated. The calculation then proceeds to the next time step.
For the surface flow system, the initial water depth h s,initial and water flow velocity u initial in each cell are zero before surface irrigation modeling is initiated. However, the numerical procedure requires all depths to be greater than zero to avoid undefined terms in Equation (1). Therefore, a small and positive value (0.0001 m) is initially assigned to the depth of flow at all cells in the flow domain. At the upstream end, the unit inflow rate is q = q 0 . If irrigation stops, the water flow velocity at the field head is zero. The field end is closed. The initial time step for the subsurface flow system is given by users and the size of time step changes in accordance with the one used for the surface flow model. Time-variable surface pressure head h s is set as the upper boundary condition; the bottom boundary conditions can be assumed as free drainage. The uniform pressure head throughout the soil profile is set as the initial condition in the computational domain.
This model can describe the impact of the interactions among soil, surface microtopography, field slope, roughness, border length and unit discharge on surface irrigation performance on a second time-step basis and exhibits the satisfactory simulation accuracy, low total volume balance error and high computational efficiency. Surface irrigation processes are described using the individual modules, which can be custom-configured to simulate soil water redistribution processes after surface irrigation. Further details on this model can be found in Dong et al. [23].

Monte Carlo Modeling
The proposed numerical methodology of modeling surface irrigation is also developed on the basis of stochastic parameter modeling method of controlling variables. In this section, Monte Carlo simulation method is used as the stochastic parameter simulation methodology.

Probability Distribution Functions of Model Parameters
Stochastic parameters, such as border length, unit discharge, soil hydraulic properties, surface micro-topography, field slope, and roughness, are input to the hybrid coupled model proposed by Dong et al. [23] as the initial conditions. Soil hydraulic properties are often nonlinearly related to soil sand content, clay content and organic carbon content. Therefore, soil hydraulic properties can be indirectly obtained from soil sand content, clay content, and organic carbon content using Pedo-Transfer Functions (PTFs). In order to assess border irrigation performance of field scale, the normal or logarithmic normal distribution functions are used to describe the distribution characteristics of controlling variables because of their spatial uncertainty. Border length, border width, unit discharge, soil sand content, soil clay content, and soil organic carbon content are of approximately normal distribution at a field scale. The probability distribution function of normal distribution can be expressed as follows: where x is the controlling variable, u is the mean value of the controlling variable, and σ is the standard deviation of the controlling variable. Surface micro-topography is one of the most important factors affecting the irrigation performance of a border irrigation system and it significantly affects the advance and recession processes of an irrigation event [28]. Micro-topography varies spatially along both the longitudinal and transversal directions of a border, thus the spatial distribution of surface elevation differences (SED) also affects the advance and recession processes and system performance. Therefore, the spatial variability of surface micro-topography includes the degree of unevenness of SED expressed by the standard deviation S d and the spatial distribution of SED between the actual and the target design elevations. Under the same field tillage and management conditions, all of the S d values of SED at a field scale are assumed as the same [26]. For a same S d the spatial distribution of SED may vary, which makes it difficult to describe the impacts of the spatial variability of SED on surface irrigation performance. Therefore, surface micro-topography at a field scale can be represented by a set of elevation distributions of SED for the same S d . Based on the field geometry (length and width) and the statistical characteristics of observed SED, the surface elevation nodes can be expressed by the following distribution: where z mt is SED (cm), u mt is the mean value of SED (cm), and σ mt is the standard deviation of SED (cm).

Monte Carlo Simulation
Monte Carlo method is developed on the basis of probability and statistics theory and can be used to analyze most of the uncertainty issues. It is assumed that the probability distribution function of the random variable is known and a set of input parameter samples are obtained using the random sampling method. Thus, numerical models run based on the input parameter samples and output the corresponding simulation results. Finally, the performance indicators are calculated on the basis of the output data and their probability distribution characteristics are fitted, respectively.
In the proposed paper, it is assumed that all of the border slope, surface roughness and the standard deviation S d of surface micro-topography at a field scale are the same under the same filed tillage and management conditions and the spatial variability of soil properties, border length, unit discharge, and the spatial distribution of SED are the main factors that affect surface irrigation performance of field scale and the spatial distribution characteristics. The Monte Carlo simulation process of border irrigation at a field scale shown in Figure 4 is discussed as follows: process of border irrigation at a field scale shown in Figure 4 is discussed as follows: (1) The probability and statistical analysis of soil sand content, soil clay content, soil organic carbon, border length, unit discharge and SED are conducted on the basis of the field test data of field scale. The roughness coefficients of 100 plots selected randomly are calculated based on the surface elevation data and field slope values of 100 plots selected randomly are obtained based on surface water depth and crop growth situation. (2) The probability distribution functions of soil sand content, clay content, silt content, and organic carbon content are obtained on the basis of the field test data of field scale and a set of random samples are obtained through their probability distribution functions by means of Monte Carlo simulation. If the parameter sample combinations of soil properties do not satisfy the given condition of Equation (12), the sample combination is replenished. Then, soil bulk density, field water holding capacity and wilt water content are obtained when the random samples of soil sand content, clay content and organic carbon content satisfying Equation (12)   (1) The probability and statistical analysis of soil sand content, soil clay content, soil organic carbon, border length, unit discharge and SED are conducted on the basis of the field test data of field scale. The roughness coefficients of 100 plots selected randomly are calculated based on the surface elevation data and field slope values of 100 plots selected randomly are obtained based on surface water depth and crop growth situation.
(2) The probability distribution functions of soil sand content, clay content, silt content, and organic carbon content are obtained on the basis of the field test data of field scale and a set of random samples are obtained through their probability distribution functions by means of Monte Carlo simulation. If the parameter sample combinations of soil properties do not satisfy the given condition of Equation (12), the sample combination is replenished. Then, soil bulk density, field water holding capacity and wilt water content are obtained when the random samples of soil sand content, clay content and organic carbon content satisfying Equation (12) are input to Soil Plant Atmosphere Water (SPAW) software [29]. Finally, the random samples of soil sand content, soil clay content, soil organic carbon, soil bulk density, field water holding capacity and wilt water content are input to Rosetta software and the corresponding samples of soil hydraulic properties are obtained.
S i + C + S + SOC = 100% (12) where S i is soil silt content, %, C is soil clay content, %, S is soil sand content, %, and SOC is soil organic matter content, %.
(3) A set of border length and unit discharge samples are obtained through their probability distribution functions by means of Monte Carlo simulation, respectively.
(4) The number of surface elevation node in a single border is determined by field length, field width and the length of surface cell which is divided along the basin length direction. Surface elevation distribution data are generated by the stochastic simulation model of basin micro-topography, which is developed to simulate two-dimensional surface micro-topography on the basis of Monte Carlo method and Kriging interpolation method [26]. To simulate the one-dimensional border micro-topography, Equation (13) is used to update the above-mentioned surface elevation distribution data for border irrigation simulation: where z y,j is the average value of SED at the grid j along the basin length direction, representing the new SED value (cm), z x,k is the SED value at the grid k along the basin width direction (cm), N x is the grid number along the basin width direction (−), and N y is the grid number along the basin length direction (−). (5) Steps (2) to (4) are repeated to obtain a set of random sample combinations, including border length, unit discharge, soil hydraulic properties, and surface micro-topography. These sample combinations and the constant parameters such as the mean values of field slope and roughness coefficient are input to the hybrid coupled model. Then, the hybrid coupled model runs repeatedly based on a set of stochastic parameter combinations and the border irrigation uniformity coefficient (CU) and border irrigation application efficiency (E a ) of each irrigation event are calculated based on the simulation values of soil water content along the border length direction. In the proposed paper, the maximum number of randomly selected samples is set to 272, which is the field-observed number.
(6) If the mean values of CU and E a of the iteration (n + 1) and mean values of CU and E a of the iteration n cannot satisfy the convergence criterion of Equations (14) and (15) (shown in Section 2.3), the iterative process continues. Otherwise, the iteration is terminated.
(7) The probability and statistical analysis method is used to analyze the changing trend of each evaluation indicator at a field scale, assess the impacts of spatial variability of field physical properties on border irrigation performance of field scale and optimize surface irrigation management plan at a field scale.

Evaluation Indicators of Surface Irrigation Performance
The application efficiency (E a ) and the uniformity coefficient (CU), as defined by Walker (1987), are selected as surface irrigation performance indicators. The application efficiency E a (%) is defined as: and the uniformity coefficient CU (%) is given by: where θ a and θ b are the average soil volumetric water content in the root zone before and after irrigation event, cm 3 /cm 3 , RD is the depth of root zone, cm, Z avg is the average infiltrated depth of water applied to the field (mm) and determined from the cut-off time and the inflow rate, Z i is the infiltrated depth of water in the ith field unit along the field length direction. N is the number of field units along the border length direction.

Model Performance Criteria of Field Scale
A numerical experiment is developed to assess the number of sets of generated stochastic parameter samples aiming at representing the spatial variability of field physical properties to be analyzed through simulation for assessing the impacts of the spatial variability of field physical properties on border irrigation performance at a field scale. Based on the probability distribution function of each variable of field physical properties, stochastic parameter samples are selected and input to the hybrid coupled model proposed by Dong et al. [23]. In these simulations with the hybrid coupled model, the same roughness, filed slope, and initial soil water content are used. When stochastic parameter samples are generated using the above-mentioned numerical methodology, more than one set of stochastic parameter sample can be generated. Different sets of stochastic parameter samples generated with the same condition will produce different values for the border irrigation performance at a field scale. Thus, it is necessary to determine how many stochastic parameter sample sets need to be generated for a given condition to appropriately analyze the impacts of the spatial variability of field physical properties on border irrigation performance at a field scale.
The number of stochastic parameter sample generations can be determined by analyzing the change trend of mean values of E a and CU resulting from the simulation data of the different irrigation event through a set of stochastic parameter sample combinations. M sets of irrigation performance indicators can be obtained by modeling the different irrigation event with M sets of stochastic parameter samples for a given condition. The number m (m < M) of stochastic parameter sample generations required to characterize the spatial variability of field physical properties may be determined by analyzing the changes of border irrigation performance with the number of stochastic parameter samples generations. The irrigation performance simulation experiment is conducted through the hybrid coupled model. Statistical analysis is carried out by using Equations (16) and (17) (m < 272), which can reflect the effect of the stochastic parameter distribution on the border irrigation performance. If the MARE k < 0.5% and RARE k < 5% of six consecutive data points are used as the discriminant indicators that the mean value and standard deviation reach the steady state, the mean values of border irrigation performance indicators and the corresponding standard deviations can be determined. MARE k and RARE k can be calculated as follows: where MARE k is the relative error between the (k + 1)th and the kth mean values of border irrigation performance indicators, RARE k is the relative error between the (k + 1)th and the kth standard deviations of border irrigation performance indicators, S k d represents the kth standard deviations of border irrigation performance indicators, and S k+1 d represents the (k + 1)th standard deviations of border irrigation performance indicators.

Data Analysis
The spatial variability of soil properties, border length, unit discharge and surface micro-topography results in the uncertainty of model input parameters. In order to represent the effects of the spatial variability of field physical properties on border irrigation performance at a field scale, the probability distribution and the classical statistical method are used in the proposed paper to verify the proposed stochastic model.

Results and Discussion
The proposed numerical methodology of assessing surface irrigation performance at a field scale is verified in this section.

Total Water Volume of Surface Irrigation
The total water volume of surface irrigation is defined as the sum of surface irrigation water volume applied to the experimental field at a field scale, calculated as, where V s is the total water volume of surface irrigation (m 3 ), V i is the irrigation water volume of the ith basin (m 3 ), A is the total area of the experimental field (m 2 ), L is the basin length (m), W is the basin width (m), and n is the number of basins in the experimental field (−).
The relative error of the total water volumes of surface irrigation between the simulation values and the observed values (RE V ) is calculated as, where V s is the simulation value of the total water volume of surface irrigation (m 3 ), V o is the observed value of the total water volume of the surface irrigation (m 3 ).
To assess the simulation performance of the proposed numerical methodology (PNM), the distributed numerical methodology (DNM) is employed. In distribution simulation, the experiment areas in MW and YH are divided into 16 and 12 sub-areas, respectively, according to soil properties, border length and unit discharge. The model parameters from each sub-area are input to the hybrid coupled model proposed by Dong et al. [23] and then the results of the sub-area level simulation are aggregated at a field scale to evaluate border irrigation performance in MW and YH, respectively. In the numerical modeling of the border irrigation, the considered surface spatial step along the border length is 2.5 m, and the considered subsurface spatial step along the vertical direction is 0.01 m. The time step of the surface and subsurface is 1 s. The initial soil water contents in MW and YH are 0.21 and 0.20 cm 3 cm −3 , respectively. Table 5 presents the relative error of the total water volumes of the border irrigation between the model-predicted values and the field-observed values. The RE V of the PNM is lower than that of the DNM. This result indicates that the PNM for scaling up the border irrigation performance at a field scale exhibits the more satisfactory simulation accuracy.

Surface Flow Advance
The surface flow advance time changes with the border length increasing and the maximum value appears at the bottom of the border. A power function curve such as described by Equation (19) is fitted successfully to the surface flow advance trajectories of the field-observed and the model-predicted, respectively. (20) where T is the time of the surface flow advance (min), X is the length of the surface flow advance (m), and a and b are fitting parameters. Parameter a represents the change range of the whole surface flow advance trajectories and parameter b reflects the change in the speed of the surface flow advance rate. The change trend of the surface flow advance trajectories within the surface irrigation processes can be expressed through parameters a and b of Equation (20). The parameter values obtained from each experimental field are checked respectively for normality through the Kolmogorov-Smirnov test statistics. They yielded a normal distribution with a significance level of 0.05. The relative errors of the average a and b values from the model-predicted and field-observed data were respectively calculated as where RE p is the relative error of the average p value from the model-predicted data and the field-observed data, p represents the fitting parameters a and b, and P s and P o are the average values of the fitting parameters from the field-observed data and the model-predicted data, respectively (−). Table 6 presents the relative error of the average a and b values calculated from the model-predicted values and field-observed values, showing that the PNM has the better simulation performance compared with the FNM and SNM.

The Irrigation Time of Surface Flow Advance
The irrigation time of surface flow advance (FT) is defined as the time of the surface flow arriving at the bottom of the border. The average absolute differences between the field-observed and the model-predicted values of the PNM and DNM are about 1.03 min, 1.46 min and 1.49 min, 2.1 min for the two experimental fields, respectively. The relative error of the average FT values of the fieldobserved and model-predicted data (RE t ) is calculated as, where t s is the average FT value under the model-predicted conditions (min) and t o is the average FT value under the field-observed conditions (min). The statistical characteristics of the FT derived from the field-observed data and model-predicted data of the PNF for the two experimental fields are shown in Figure 5. The statistical characteristics of the FT from the field-observed data and model-predicted data of the PNM present the similar results, which indicate that the PNM performs satisfactorily. The relative errors of the average FT values from the field-observed and model-predicted data of the PNM and DNM are shown in Table 7, which illustrates that the PNM has the better simulation accuracy compared with the DNM. model-predicted data of the PNF for the two experimental fields are shown in Figure 5. The statistical characteristics of the FT from the field-observed data and model-predicted data of the PNM present the similar results, which indicate that the PNM performs satisfactorily. The relative errors of the average FT values from the field-observed and model-predicted data of the PNM and DNM are shown in Table 7, which illustrates that the PNM has the better simulation accuracy compared with the DNM.

Soil Water Content
In this case, soil water content is used to assess the modeling accuracy of the proposed numerical methodology. The average soil water content of each sampling site is calculated based on those of three depths (0 cm to 20 cm, 20 cm to 40 cm, and 40 cm to 70 cm), using the weight average method. The spatial distribution of the average soil water content in MW and YH after the ending of the surface irrigation is intuitively presented in Figure 6. The spatial distribution data of the average soil water content were statistically analyzed to obtain the cumulative probability distribution of the

Soil Water Content
In this case, soil water content is used to assess the modeling accuracy of the proposed numerical methodology. The average soil water content of each sampling site is calculated based on those of three depths (0 cm to 20 cm, 20 cm to 40 cm, and 40 cm to 70 cm), using the weight average method. The spatial distribution of the average soil water content in MW and YH after the ending of the surface irrigation is intuitively presented in Figure 6. The spatial distribution data of the average soil water content were statistically analyzed to obtain the cumulative probability distribution of the average soil water content, as shown in Figure 6. The cumulative probability of the average soil water content of the model-predicted data with PNM, which is calculated using the soil water content data along the vertical direction of soil profile, is also shown in Figure 7. The agreement between the cumulative probability curves of the field-observed and the model-predicted data with the PNM is satisfactory. average soil water content, as shown in Figure 6. The cumulative probability of the average soil water content of the model-predicted data with PNM, which is calculated using the soil water content data along the vertical direction of soil profile, is also shown in Figure 7. The agreement between the cumulative probability curves of the field-observed and the model-predicted data with the PNM is satisfactory.

Conclusions
To assess the impact of the spatial variability of field physical properties on the surface irrigation performance at a field scale, a numerical methodology is proposed on the basis of the hybrid coupled surface-subsurface flow model of surface irrigation and Monte Carlo simulation method. The hybrid coupled surface-subsurface flow model is developed from the numerical solution of the one-dimensional Saint-Venant equations using an improved hybrid numerical method and from the one-dimensional Richards equation using a finite volume method with fourth-order spatial accuracy. Monte Carlo simulation method is used to obtain the stochastic samples of field physical properties, and the stochastic simulation model of surface micro-topography is used to obtain a set of surface micro-topography samples for the same standard deviation of surface elevation differences. The proposed numerical methodology is tested in terms of the total water volumes of the surface irrigation, surface flow advance trajectories, irrigation time of each border, and soil water contents from the model-predicted and the field-observed data. Compared with the distributed-parameter modeling methodology, the

Conclusions
To assess the impact of the spatial variability of field physical properties on the surface irrigation performance at a field scale, a numerical methodology is proposed on the basis of the hybrid coupled surface-subsurface flow model of surface irrigation and Monte Carlo simulation method. The hybrid coupled surface-subsurface flow model is developed from the numerical solution of the one-dimensional Saint-Venant equations using an improved hybrid numerical method and from the one-dimensional Richards equation using a finite volume method with fourth-order spatial accuracy. Monte Carlo simulation method is used to obtain the stochastic samples of field physical properties, and the stochastic simulation model of surface micro-topography is used to obtain a set of surface micro-topography samples for the same standard deviation of surface elevation differences.
The proposed numerical methodology is tested in terms of the total water volumes of the surface irrigation, surface flow advance trajectories, irrigation time of each border, and soil water contents from the model-predicted and the field-observed data. Compared with the distributed-parameter modeling methodology, the proposed numerical methodology shows the better simulation performance in total water volume of surface irrigation, surface flow advance trajectories and soil water content redistribution. In conclusion, the proposed numerical methodology presents the satisfactory simulation performance and can be used as a tool for analyzing surface irrigation performance and similar hydrologic patterns at large scales.