A Comprehensive Quasi-3 D Model for Regional-Scale Unsaturated-1 Saturated Water Flow 2

Abstract. For computationally efficient modeling of unsaturated–saturated
flow in regional scales, the quasi-three-dimensional (3-D) scheme that
considers one-dimensional (1-D) soil water flow and 3-D groundwater flow
is an alternative method. However, it is still practically challenging for
regional-scale problems due to the highly nonlinear and intensive input data
needed for soil water modeling and the reliability of the coupling scheme.
This study developed a new quasi-3-D model coupled to the UBMOD 1-D soil water balance
model with the MODFLOW 3-D hydrodynamic model. A new implementation
method of the iterative scheme was developed in which the vertical net
recharge and unsaturated zone depth were used as the exchange information. A
modeling framework was developed to organize the coupling scheme of the soil
water model and the groundwater model and to handle the pre- and
post-processing information. The strength and weakness of the coupled model
were evaluated by using two published studies. The comparison results show
that the coupled model is satisfactory in terms of computational accuracy
and mass balance error. The influences of spatial and temporal discretization
as well as the stress period on the model accuracy were discussed.
Additionally, the coupled model was used to evaluate groundwater recharge in
a real-world study. The measured groundwater table and soil water content
were used to calibrate the model parameters, and the groundwater recharge
data from a 2-year tracer experiment were used to evaluate the recharge
estimation. The field application further shows the practicability of the
model. The developed model and the modeling framework provide a convenient
and flexible tool for evaluating unsaturated–saturated flow systems at the
regional scale.


2 Abstract: For computationally efficient modeling of unsaturated-saturated flow in 22 regional scales, the Quasi three-dimensional (3-D) scheme that considering one-23 dimensional (1-D) soil water flow and 3-D groundwater flow is an alternative method. 24 However, it is still practically challenging for regional-scale problems due to the high 25 non-linear and intensive input data needed for soil water modeling and the reliability of conductivity Ks, saturated water content θs, field capacity θf, and residual water content 208 θr) is used to describe the hydraulic diffusivity D(θ). The heterogeneity of soils is also 209 taken into account by adding a correction item in the right side, which makes the model 210 applicable to heterogeneous situations. With the help of the diffusive term, the UBMOD 211 can consider upward soil water movement, which is ignored by most water balance 212 models. The details of D(θ) are shown in the Appendix. 213 The original UBMOD is a soil water balance model, which cannot consider 214 groundwater

The Brief Introduction of MODFLOW and Two Peripheral Tools 217
MODFLOW is a computer program that numerically solves the 3-D groundwater 218 flow equation for a porous medium using a block-centered finite-difference method 219 (Harbaugh, 2005). The governing equation solved by MODFLOW is, 220 ArcGIS for Python (Toms, 2015), which provides a useful and productive way to 234 perform geographic data analysis, data conversion, data management, and map 235 automation with Python. 236 237   12   The procedures of the modeling framework are composed of three major parts,  238 including the pre-processing, the coupled model, and the post-processing. The 239 preparation of geographic input information of the model shown in Fig. 1(a)  partitioned into a number of sub-areas in the horizontal direction mainly according to 257 the spatially distributed inputs (each sub-area is considered to be homogeneous in 258 13 horizontal). There are l sub-areas, j layers for a specific soil column shown in Fig. 1(b). 259

The Process of Geographic Input Information
Soil water flow of each sub-area is simulated by using one 1-D soil column. The 260 recharge at the bottom boundary calculated by UBMOD is treated as the upper 261 boundary condition of MODFLOW. The whole saturated zone is discretized into a grid 262 with cells, and there are m rows and n columns cells of the saturated zone as shown in 263 Since the independent variable of UBMOD is the soil water content and the 276 independent variable of MODFLOW is the hydraulic head, this study uses the vertical 277 net recharge and the unsaturated zone depth to couple the unsaturated zone and 278 saturated zone. The domain shown in Fig. 1(b) is used as an example to illustrate the 279 14 spatial and temporal coupling methods in the study. The vertical net recharge is 280 represented by vector R with mn elements, and the unsaturated zone depth by vector 281 Du with l elements, as illustrated in Fig. 1(b). Scalar R is used to denote the specific net 282 recharge of a soil column to the corresponding saturated sub-area, and scalar du denotes 283 the depth of the soil column. Figure 2(a) shows the spatial coupling method of a soil 284 column connected with groundwater system. The water table locates in the j-th layer. 285 The net recharge R from soil zone is calculated by UBMOD as follows, 286 These four terms are corresponded to the four major components in UBMOD, as 292 described in Sect. 2.1. Specifically, the infiltration water is allocated first according to 293 Eq. (1) if there is precipitation or irrigation. When there is residual infiltration water 294 across the water table in the j-th layer, the amount of residual infiltration is denoted as 295 qI. Then the advective flow qA across the water table driven by gravitational potential 296 is calculated by Eq. (2). The direction of these two terms is downward. The qS term is 297 caused by evapotranspiration. When the critical depth of evapotranspiration is 298 shallower than the groundwater table depth, the groundwater can be consumed by 299 evapotranspiration and it causes an upward qS term. A virtual layer is needed when 300 15 calculating the diffusive movement driven by matric potential across the water table 301 based on Eq. (4). As shown in Fig. 2 (a), the virtual layer will be added under water 302 table, numbered as layer j+1. The thickness, Mj+1 [L], of the layer is set as, 303 The UBMOD can give acceptable results when Δtu is shorter than 10 d for assumed 315 cases and 1 d for a real-world case (Mao et al., 2018). Δts is set as the technical report 316 described by Harbaugh (2005) and can be changed during the calculation. 317 The implementation of iterative coupling scheme is shown in Fig. 2(c), which 318 shows the calculation period from t to t+ΔT. At the time t, the saturated hydraulic head 319 is known, marked as H t (mn dimension). When the model runs from t to t+ΔT, firstly, 320 the initial saturated hydraulic head H t+ΔT at t+ΔT is set to be equal to H t , and then the 321 16 average unsaturated depth from t to t+ΔT is calculated according to H t+ΔT , marked as 322 for one soil column is 323 calculated as follows, 324 where εH is a user-specified tolerance [L]. If the criterion is met, the iteration stops, and 344 H t+ΔT, p is the convergent results at time t+ΔT, and the model proceeds to the next stress 345 period. Otherwise, the iteration continues to p+1 and H t+ΔT, p will be used to calculate 346 the average unsaturated depth shown in Eq. (8). The above procedures will be repeated 347 until the convergence criterion of Eq. (11) is met. 348

Model Evaluation 349
In this section, two test cases were designed to evaluate the model accuracy and   The groundwater table depth calculated by scenarios with different temporal 458 discretizations (S1-S3) are compared with those from HYDRUS-1D in Fig. 7(a) scenarios with different spatial discretizations (S1, S4 and S5) are compared with those 464 23 from HYDRUS-1D in Fig. 7(b). The ARE values of the three scenarios are smaller than 465 25%. The maximum RMSE value is 0.215 m. The IA and R 2 values are larger than 0.95. 466 The water table depth calculated by scenarios with different stress periods (S1, S6, S7 467 and S8) are compared with those from HYDRUS-1D in Fig. 7(c). It should be noted 468 that the model collapsed at the time of 227 d when the stress period is 15 d (S8). The 469 statistical index values for S1, S6 and S7 are shown in Table 2 Fig. 4(b). This is 480 caused by a significant lateral flow in the unsaturated zone during the early period due 481 to the relatively low initial soil water content condition. Therefore

Study Site and Input Data 511
The coupled model was used to calculate the regional-scale groundwater recharge 512 in a real-world case, where the shallow groundwater has significant impact on the soil 513 water movement. Figure 8 which are farm land, villages and bared soil, as shown in Fig. 8(b). The crop types in 527 26 the farmland were not considered for determining the sub-areas. The surface digital 528 elevation model (DEM) is shown in Fig. 8(c). 529 The were divided into 6 layers for numerical simulation. Ten groundwater monitoring wells 548 27 were set in this district, and the groundwater tables were observed every 6 days. Well 549 1, well 2, well 3, well 5 and well 6 are located in farm land areas, well 4 and well 8 in 550 villages, and well 7, well 9 and well 10 in bared soil areas. Additionally, there are 5 soil 551 water content monitoring points in the farm land and 2 points in the bared soil area, as 552 shown in Fig. 8(a). Soil water contents within 1 m depth were observed 1-3 times every 553 month from May to October in 2004. 554 Five GIS files are prepared as the shapefile files of the study domain, the land 555 usage types, the boundary conditions, and raster files of the surface DEM and initial 556 hydraulic head. There were 150 rows and 50 columns used in the MODFLOW model. 557 The spatial discretization of UBMOD was set to be 0.1 m. The stress period ΔT was set 558 as 5 d, and the MODFLOW time step Δts and UBMOD time step Δtu were set as 1 d. 559

Model Calibration Results 560
There are two soil types in the first layer as loamy sand and loam. The unsaturated 561 hydraulic parameters of the two soils are listed in Table 4. The hydraulic conductivity 562 of the top aquifer in MODFLOW was set as the same as the unsaturated layer, and the 563 hydraulic conductivity of the bottom sand aquifer was set as 3.5 m/d during the 564 calibration, and the specific yields of the top and bottom were set as 0.08 and 0.1, 565 respectively. Figure 10 shows the comparison of the simulated and observed water table 566 depth for the whole area and locations of different monitoring wells. The statistical 567 index values are listed in Table 5. It can be found that the ARE, RMSE, IA and R 2 values 568 are 9.9%, 0.203 m, 0.869 and 0.71 for the regional average water table depth. Larger 569 28 deviations of simulated water table depth can be found for the locations of monitoring 570 wells, with RMSE values ranging from 0.25 m-0.39 m. Figure 11 further shows the 571 spatial distribution of the simulated water table depth at different output times. The 572 increasing trend is obviously found in Fig. 11(a) to Fig. 11(c) in the crop growing season, 573 during which the groundwater was consumed by crop transpiration and strong soil 574 evaporation. When the intensive autumn irrigation happened after the 160 th day, the 575 water table depth in the farm land decreased rapidly, as shown in Fig. 11(d) The computational cost of the real-world application is 120 s, which is efficient 604 considering the scale of the problem. 605

Regional Groundwater Recharge 606
In the tracer experiment, bromide (Br) was used as the tracer for calculating 607 groundwater recharge. The tracer was injected at 1 m depth at two locations shown in 608