Simulating advance distance in border irrigation systems based on the improved method of characteristics

Improving the simulation accuracy of the advance distance based on the method of characteristics is essential to develop numerical solutions for simulating surface irrigation. Instead of volume balance in the traditional method of characteristics (T-MC), the position of critical flow is determined to simulate the advance distance in the improved method of characteristics (I-MC), which is used in border irrigation systems with rapid variation in inflow discharge in the current research. Specifically, the zones of both subcritical and supercritical flow were firstly distinguished to determine the position of the critical flow point, which was also the upstream boundary of the wetting front region, and then the advance distance was calculated by applying the diffusive wave equation in the wetting front region. The results showed that the I-MC accurately simulated the advance distance with high determination coefficients (0.984-0.998) and low errors (root mean square error of 0.35-1.56 min and coefficient of residual mass of 0.01-0.06), which performed much better than the T-MC. The I-MC provided a suitable and simple numerical simulation tool to improve the establishment of numerica l surface irrigation models.


Introduction 
Surface irrigation is the most common method for irrigating farmland across the world. Although surface irrigation is easy to operate, it is complex to design and optimize, especially under varied-discharge conditions, such as the increased-discharge irrigation [1] and the real-time control irrigation [2] . Numerical analysis methods are critical tools for evaluating and simulating the surface flow of irrigation systems. The mathematical models of surface irrigation mainly include the complete hydrodynamic model, the zero-inertia model and the kinematic-wave model [3,4] . The complete hydrodynamic model, known as the Saint-Venant equations, is the most comprehensive and accurate one for simulating surface irrigation [3] .
Saint-Venant equations, the hyperbolic partial differential equations, can be derived through numerical solutions, such as the finite difference method, the method of characteristics, the finite element method and the finite volume method. Although the finite difference method is relatively intuitive [5] , it leads to numerical vibrations in the flow front region [6] . The finite element method and the finite volume method, which divide the computational domain into irregular units, are advantageous in solving two-dimensional flow calculations with complex boundary conditions and wide water surface [7,8] , while the procedure of computations is complex [9,10] .
The method of characteristics has numerous advantages, including rigorous mathematical analysis, clear physical concept, simplicity and high precision, and it is especially suitable for hyperbolic partial differential equations [11,12] . Moreover, the method of characteristics is more accurate in calculating flow transients than other methods [13] . The method of characteristics has been widely used for various engineering purposes, such as calculating the time-dependent head and velocity of fluids in a complex pipeline network [14] , determining the turbine and compressor boundary conditions [15] and simulating hydraulic transients and flow propagation in hydraulic systems [12] .
However, the method of characteristics is rarely used to solve Saint-Venant equations in surface irrigation systems due to the following two reasons: (1) The method of characteristics is usually explicit schemes with limited time, which results in long operation time for simulating long border or furrow irrigation; (2) The flow depth rapidly approaches 0 near the wetting front of surface flow and its characteristics are strongly curved. When using the method of characteristics to calculate the advance distance, it is difficult to converge the numerical analysis and a special numericalprocess is required [16] .
The implicit method of characteristics has removed the limitations of time step [17] . Deriving the advance distance based on the shape-factor method [18] makes it feasible to apply the method of characteristics for numerical solutions of surface irrigation models. This method used a simple shape factor and assumed a normal depth at key points along the surface (such as the upstream end of the field). The amount of irrigation water then was approximated by the flow and infiltration depths at the key points of the field, and the advance distance was derived through the volume balance approach. However, the shape factor is difficult to obtain directly and is usually recorded as a fixed empirical value (generally 0.7-0.8).
The accuracy of the shape-factor method is controversial due to the high variability of the shape factor [19] . Many researchers have attempted to modify this method to improve its accuracy [20][21][22] . However, when the inflow discharge changes during the advance phase, for instance, when the inflow is cut off prior to the water flow reaching the border endpoint, which is an important irrigation strategy in border irrigation systems [23] , this method will no longer be viable due to its indeterminable shape factor. Another approach, named as a depth-cumulation method based on the method of characteristics (D-MC), implemented by Liu and Hui [24] , calculated the surface water by accumulating flow depth at each node. D-MC almost gets rid of the limitation of the shape factor and derives the advance distance when inflow discharge changes. However, this method leads to numerical vibrations due to inaccurate calculations of the surface water by the depth of each node. To improve its accuracy, the space step needs to be very small, and the calculation will be more extensive.
The purposes of this study are: (1) to present an improved method of characteristics (I-MC) in calculating the advance distance in the varied-discharge border irrigation system; (2) to evaluate the proposed method by comparing with field data, results simulated from WinSRFR and the traditional method D-MC.

Basic equations and boundary conditions
Border irrigation is characterized by a shallow flow over a gentle slope bed. Since the width of the border is small relative to the length, the flow conforms to one-dimensional long-wave equations. The mass conservation equation where, x is the distance along the border, m; t is the time from the beginning of the inflow, s; h is the flow depth, m; V is the mean velocity in the cross-section, m/s; i is the infiltration rate, m/s. The common approximation to the momentum equation is where, g is the acceleration due to gravity, m/s 2 ; s 0 is the bed slope.
In the border irrigation system, since the border width is much larger than the water depth, the hydraulic radius is approximated as the flow depth and s f is approximately calculated as where, Q is volume flow rate, m 3 /s; A is the cross-sectional area of flow, m 2 ; R is the hydraulic radius, m; n is the Manning roughness coefficient.
Kostiakov equation is the one of most frequently used infiltration equation for the numerical simulation of surface irrigation [25] and is written as: where, I is the cumulative infiltrated depth, m; k is the infiltration coefficient, m/s α ; α is the infiltration index.
Considering the wide application of closed-ended border irrigation in China, boundary conditions of the closed-ended border irrigation system are described. The upstream and downstream boundary conditions are written as Equations (5) and (6), respectively.
where, q i is the initial discharge per unit width, m 2 /s; q c is the discharge per unit width after changing, m 2 /s; t c is the time of discharge changing, s; t s is the cut-off time, s; t z is the time of the water depth at the upper end of the border falls to 0, s; x t is the distance of the trailing edge, m; x f is the distance of wetting front, also named as advance distance, m; L is the border length, m; t d is the time of the flow reaches the end of the border, s.

Computational schemes of I-MC
In border irrigation systems, the slope is generally small and the irrigated surfaces are relatively rough, resulting in subcritical ). However, near the wetting front, the flow depth rapidly approaches 0 while the velocity negligibly changes, which creates supercritical flows ( V gh  ) [24] . There will be a critical point, named critical flow point ( V gh  ), during the transition from subcritical flow to supercritical flow. The area from the critical flow point to the distance of the wetting front is defined as the wetting front region. In this region, the flow depth rapidly drops to 0, and it is difficult to converge numerically using the method of characteristics.
Since the system of equations (Equations (1) and (2)) are hyperbolic, the numerical solution of surface flow except for the wetting front region follows the method of characteristics. At every point in the x-t plane, characteristic curves are drawn backward from each node until they intersect with the previous timeline, as shown by point P in Figure 1. Except for the boundary points, each point has two characteristic curves, one positive (PM in Figure 1) and one negative (PN in Figure 1).
Note: The horizontal lines indicate time, the vertical lines indicate distance, the curves indicate characteristic curves, and the letters indicate calculation nodes.

Figure 1 Schematic characteristic curves of nodes
The Kostiakov equation (Equation (4)) is substituted into Equation (1) and the basic equations are reduced to equations that involve the total differentials as Equation (7a) for positive characteristic curve and Equation (7b) for negative characteristic curve [26] . Figure 1 shows an example of calculating the period of t n based on the period of t n−1 with a rectangular time-space computation grid system, where h and V are known at the points J, K, and L at the time t n , while the unknown points M and N can be solved by linear interpolation formulas, then h and V at point P and the time t n+1 can be calculated by Equation (7). The calculation of the upstream boundary point is similar to the calculation for the interior points, except that there is only one characteristic curve, as shown by point C in Figure 1. In this case, Equation (7) needs to be solved simultaneously with the upstream boundary condition Equation (5).
In the process of node calculation, if a node is calculated that corresponds to flow depth h<V 2 /g, then the point is written as x E and the previous point is considered to be the downstream boundary point x D . A finer grid spacing Δx' is then taken from x D to x E . These refined nodes are sequentially solved by Equation (7) from x D to x E until the special node (a node is calculated to correspond to flow depth h≤V 2 /g and its previous node is calculated to correspond to flow depth h>V 2 /g) is found. This special point was considered as a critical flow point (x C ) with Δx' precision. In the wetting front region, since the flow depth quickly drops to 0, the absolute values of ∂h/∂x and s f are far greater than other items [27] . Therefore, the momentum equation is reduced to a diffusive wave equation [28] .
According to Chanson and Liem and Köngeter, the flow velocity in the wetting front region is not rapidly changed [28,29] . Therefore, the mean velocity in the cross-section in the wetting front region is considered homogeneous and written as V instead of V in Equation (8). Then the relationship between h and x in the wetting front region can be obtained by the integral of Equation (8) and written as: where, f is the integration constant. At the point of critical flow, f can be obtained with the following equation: where, x C is the distance of the point of critical flow, m; and h C is the critical depth, m.
Equation (10) is substituted into Equation (9) and flow depth h can be written as: Since the flow depth of the flow front is 0, V can be obtained by solving for Equation (11).
In the wetting front region, Equation (8) is simplified by the momentum equation. 4 4 3 3 where, h f is the flow depth of the wetting front, m; h C is the critical depth, m; V C is the mean velocity in the cross-section at point x C , m/s; V f is the mean velocity in the cross-section at point x f , m/s; and the other variables are as previously defined.
Since the flow depth of the wetting front is 0 and the mean velocity in the cross-section at point x f and x C are all equal to V , the length of the wetting front region is written as: Since x C , V and h C have been derived in the previous section, the advance distance x f can be calculated with Equation (14). It should be noted that there is no initial inflow discharge at t=0. Therefore, the initial conditions are 0 for the cross-section, and h and V are also equal to 0 for every x. The initial inflow discharge is constant. To start the calculation, the surface shape factor r Y and the subsurface shape factor r Z are introduced [20] as follows: Z 0 where, W S is the total amount of surface water per unit width, m 2 ; W I is the total amount of infiltration water per unit width, m 2 ; h 0 is the flow depth at the upper end of the border, m; z 0 is the infiltration depth at the upper end of the field, m.
According to the volume balance approach, W S plus W I equals the amount of irrigation water. Then the advance distance of the first-time interval Δt is written as: The shape-factor method is only used in the first calculation step, which means that r Y and r Z are not used in the subsequent calculation steps. Therefore, the values of the shape factors have little effect on the simulation results. The r Y and r Z values are used directly as constants. In this study, these values are r Y = 0.8 and r Z = 0.7. In addition, before the change in inflow discharge (t≤t c ), the advance distance, x f , can still be derived by the simultaneous solutions of Equations (5), (7) and (17), and this is the shape-factor method based on the method of characteristics.

WinSRFR and D-MC
To evaluate the accuracy of I-MC, the calculation results were compared with the predicted results by WinSRFR4.1 simulation model [30] and D-MC [24] . When the discharge changes rapidly, the surface flow profile will change accordingly, so W S is hard to be calculated directly by the shape-factor method.
D-MC approximates the surface water from 0 to x c by the accumulation of the water depth over the space step, as shown by the dotted line in Figure 2. The surface water from x c to x f is calculated by r Y and h c . Then the advance distance is obtained by the volume balance: Vol. 14 No. 3 159 WinSRFR is an integrated software for analyzing surface irrigation systems and developed by the USDA-Agricultural Research Service. WinSRFR4.1 simulation model uses finite volume method to discrete control equations, then solves equations by double-sweep technique [31] . At present, WinSRFR is the most commonly used software to simulate surface irrigation systems.

Field experiments
The ability of the proposed method to simulate border irrigation was verified with eleven sets of field test data, and each test was distinct in several technical elements, including border length, Manning roughness coefficient and rapid discharge variations (cut-off and change of discharge). The experimental field is located at the Nanpi Ecological Agricultural Experiment Station of Chinese Academy of Sciences, Hebei Province, China (116°40ʹE, 38°06ʹN). This area has a continental monsoon climate. Winter wheat is the main crop to be cultivated, and closed-ended border irrigation is the main irrigation method. At this site, the groundwater depth is below 3 m, the surface soil is silt loam (67.02% silt, 25.19% sand, and 7.79% clay, on average), and the dry bulk density is 1.48 g/cm 3 . As rainfall is insufficient and uneven (mainly concentrated in summer) in this area, it is urgent and essential to increase the border irrigation water use efficiency.
The experimental data from Salahou1 to Salahou5 were obtained from the literature [22] . Supplementary experiments with different border lengths (B1 and B2) and changing discharge during irrigation (B3, B4, B5 and B6) were conducted at the same experimental field. Consistent with local irrigation forms, all test borders were closed-ended, and the cut-off and change of discharge were based on distance [23] . The soil infiltration tests were carried out during the irrigation process, and the empirical coefficients k and α of the Kostiakov equation were calculated by the infiltration time and cumulative infiltration depth of a series of observation points [32] . The Manning roughness coefficient n was determined using the trial and error approach. To be specific, enter different n in WinSRFR simulation model from 0.01 to 0.30 with steps of 0.01 to predict advance and recession curves. And the Manning roughness coefficient of the border was obtained when the best fit between the predicted and the observed curves [23] .
The time-averaged inflow rate and the distance-averaged slope were measured to reduce the impact of any errors. 14 Note: L is the border length; D is the border width; s0 is the average border slope; qi is the initial discharge per unit width; qi is the initial discharge per unit width; xc is the distance of wetting front when discharge changed; qc is the discharge per unit width after changing; CR is the cut-off distance ratio; k is the infiltration coefficient; α is the infiltration index; n is the Manning roughness coefficient.

Assessment statistics
Model simulations were carried out for estimating advance trajectories and flow depths. Several statistics were computed to evaluate the accuracy of I-MC, as follows: coefficient of determination, R 2 [21,33] , root mean square error, RMSE [34] , coefficient of residual mass, CRM [33] . R 2 , RMSE and CRM were used to compare measured and predicted advance trajectory, while R 2 and RMSE were also used to compare measured and predicted flow depths. R 2 varies from 0 to 1, and the closer it gets to 1, the more valuable simulation results are. RMSE has a minimum value of 0, and the closer it gets to 0, the more accurate the simulation results are. CRM can be negative or positive, and it indicates the underestimation or overestimation amounts of predictions.
where, P i is the calculated value; O i is the observation value; O ave is the average of the observation values; n is the number of evaluated points.

Results and discussion
In order to solve the Saint-Venant equations with the method of characteristic, D-MC was developed to calculate the advance distance [24] . For accuracy, this method requires small time-space steps. The simulated advance trajectories generated by D-MC (Δx = 1 m, Δt = 1 s) and I-MC (Δx = 2.5 m, Δt = 3 s) are shown in Figure 3. The RMSE values of the advance trajectories of D-MC and I-MC were 1.643 min and 0.544 min, respectively. There was a prominent numerical vibration in the advance distance with D-MC. The accuracy of I-MC was significantly higher than D-MC. In addition, the time-space step size had a large influence on the efficiency of the calculation. With the current time-space step, the computational times for D-MC and I-MC were 7.98 min and 0.75 min, respectively. The absolute error of simulated advance time with D-MC gradually increased, while this trend did not occur with I-MC ( Figure 4). D-MC is essentially a method of volume balance, which has inherent inaccuracies [35] . The inaccuracies are further exacerbated when approximating the amount of irrigation water by the accumulating flow depth of each node. As shown in Figure 2, when calculating the surface water from 0 to x c by the accumulation of the water depth over the space step, the calculated value is slightly larger than the actual value. Therefore, the value of advance distance calculated by Equation (18) is slightly smaller than the actual value. This error is systematic and accumulates as the advance distance increases. To reduce the numerical vibration and inaccuracies, the time-space steps need to be refined, which in turn will lead to a significant increase in computational time. D-MC becomes even less practical when used with a finer grid. Therefore, I-MC is more accurate, convenient and practical than D-MC. For further verification, I-MC was used to predict the irrigation advance trajectories for all borders and the results were compared with the field measured data and WinSRFR simulated results. The predicted advance trajectories were in good agreement with the field data ( Figure 5). The average absolute errors and average relative errors were 0.28 min and 4.6% for B1, 0.94 min and 4.3% for B3, respectively. And the WinSRFR software also simulated advance trajectories very similar to I-MC.
More detailed statistical data were shown in Table 2   In addition, the accuracy of I-MC was verified by comparing the measured flow depth with the simulated flow depth. It should be noted that due to the flow fluctuates near the upper boundary, the flow depth observation position was not 0 m but 5 m. The flow depth simulated by WinSRFR is very limited and does not include this observation position, so there was no comparison with WinSRFR. Flow depth was measured by Odyssey Capacitance Water Level Logger (resolution is 1 mm). Figure 6 shows that I-MC simulated flow depth with reasonable accuracy. The R 2 and RMSE values for B5 were 0.987 and 4.032 mm and those for B6 were 0.957 and 6.206 mm, respectively.  Zhang et al. [9] presented a hybrid numerical method, which used the finite-difference method, the finite element method and the finite volume method to discretize the different vector terms of Saint-Venant equations. They verified the method with four border irrigation tests, and the average relative errors of flow advance times were 4.0%, 3.2%, 6.3% and 6.1%, respectively. Sayari et al. [31] developed a zero-inertia finite element method to simulate the flow in surface irrigation. And the RMSE values of advance times were 0.38 min, 0.54 min and 2.36 min for 50 m, 72 m and 72 m furrow, respectively. These two methods were similar in accuracy to I-MC. Soroush et al. [36] simplified the Saint-Venant equations and solved them using the finite difference method. They compared their method with WinSRFR on the simulated advance times, and the results showed that for 71% of the total number of irrigation tests, WinSRFR was more accurate. Therefore, it was concluded that the accuracy of I-MC is satisfactory. However, the method of characteristic is more convenient and intuitive than the finite element method and the finite volume method (WinSRFR also used this method) [12,13] . Furthermore, the method of characteristics is also advantageous for dealing with convection-dominated diffusion problems [37] . And future surface irrigation simulation models will need to simulate surface water flow and solute transport [38,39] . Therefore, I-MC has the potential to be integrated into solute transport models and better serve the surface irrigation simulations.

Conclusions
A new method, based on the method of characteristics, was developed to calculate the advance distance in border irrigation systems with rapid discharge variations. The inapplicable area of the method of characteristics in simulating surface flow was identified, and the diffusive wave equation was used in this area to derive the advance distance. This method is highly accurate and simplifies the numerical simulation of border irrigation. In all examined cases, the proposed method provided predictions that were in good agreement with the field data. Additionally, this method effectively avoids issues with numerical vibration. Compared with other methods, including WinSRFR, the hybrid numerical method, the zero-inertia finite element method and the simplified Saint-Venant equations finite difference method, the numerical simulation accuracy of the proposed method is equivalent or better.
Considering that the method of characteristics is convenient and advantageous for solving Saint-Venant equations and dealing with convection-dominated diffusion problems, the proposed method provides a suitable numerical approach to the numerical simulation of border irrigation while laying the foundation for the application of the method of characteristics in future surface irrigation simulations to model surface water flow and solute transport.