Comparative Assessment of Methods for Coupling Regional and Local Groundwater Flow Models: A Case Study in the Beijing Plain, China

: A coupled regional and local model is required when groundwater ﬂow and solute transport are to be simulated in local areas of interest with a ﬁner grid while regional aquifer boundary and major stresses should be retained with a coarser grid. The coupled model should also maintain interactions between the regional and local ﬂow systems. In the Beijing Plain (China), assessment of managed aquifer recharge (MAR), groundwater pollution caused by rivers, capture zone of well ﬁelds, and land subsidence at the cone of depression requires a coupled regional and local model. This study evaluates three methods for coupling regional and local ﬂow models for simulating MAR in the Chaobai River catchment in the Beijing Plain. These methods are the conventional grid reﬁnement (CGR) method, the local grid reﬁnement (LGR) method and the unstructured grid (USG) method. The assessment included the comparison of the complexity of the coupled model construction, the goodness of ﬁt of the computed and observed groundwater heads, the consistency of regional and local groundwater budgets, and the capture zone of a well ﬁled inﬂuenced by the MAR site. The results indicated that the CGR method based on MODFLOW-2005 is the easiest to implement the coupled model, capable of reproducing regional and local groundwater heads and budget, and already coupled with density and viscosity dependent model codes for transport simulation. However, the CGR method inherits shortcomings of ﬁnite difference grids to create multiple local models with inefﬁcient computing efforts. The USG method based on MODFLOW-USG has the advantage of creating multi-scale models and is ﬂexible to simulate rivers, wells, irregular boundaries, heterogeneities and the MAR site. However, it is more difﬁcult to construct the coupled models with the unstructured grids, therefore, a good graphic user interface is necessary for efﬁcient model construction. The LGR method based on MODFLOW-LGR can be used to create multiple local models in uniform aquifer systems. So far, little effort has been devoted to upgrade the LGR method for complex aquifer structures and develop the coupled transport models.


Introduction
Groundwater models are widely applied for groundwater management purposes. The finite-difference computer program MODFLOW is the most commonly applied code to analyse groundwater flow system, predict groundwater level changes, and calculate water balance components. The high capacity of the modern computers makes groundwater studies on large basin-scale much easier in recent decades. Regional groundwater models provide insights into the groundwater system on large scale and decision-makers can consider the bigger picture when choosing a good pathway towards the long-term groundwater sustainability of the region [1]. However, as local scale groundwater problems increase, a finer grid resolution is required to investigate problems like analyzing well fields [2], tracing groundwater contaminants [3], delineating the capture zones of pumping wells, etc. The complexity of regional groundwater models is normally insufficient to represent the local heterogeneity [4]. Using a large scale regional model to simulate the local scale groundwater problems can be inaccurate and not cost-beneficial [5] and it is common that the natural hydraulic boundary is located far from to the target local study area so that the extent of the model area becomes much larger than the area of interest [6]. In these situations, using unified global refined grids would be computationally expensive [7]. This issue can be solved using the emerging regional and local model coupling technique, which allows various grid size in one model domain. Finer grids in the target research area can be embedded into the regional model with coarser grids [8,9].
In general, the coupling of regional and local models can be achieved by three schemes: full coupling, one-way coupling and two-way coupling [5]. The conventional grid refinement method (hereafter referred to as the CGR method) in MODFLOW is the most stable and viable fully coupling scheme [10]. It integrates the regional and local model in one assembled matrix. The implementation of the CGR method is the most straightforward method and requires the least modification of the original structured regional model. However, the CGR method loses its advantage when there is more than one local model embedded with the regional model because it requires the refinement of entire rows and columns where the local model is located.
A typical one-way coupling scheme in MODFLOW is the telescopic mesh refinement (TMR) method [8]. The TMR method is usually numerical stable and computationally efficient, and compatible to operate some independently developed computer codes [11,12]. However, the drawback of this method is the lack of numerical rigor [9]. It only uses flux or heads information from the regional model as the local model boundary conditions. Without a feedback mechanism between regional and local grids, huge sub-model error can be introduced when local parameters or stresses are different from the regional scale.
The two-way coupling scheme links the regional grid and local refined grid and provides feedbacks between the two grids. The local grid refinement (LGR) method was developed based on the MODFLOW-2005 code as a two-way coupling scheme. MODFLOW-LGR allows to extract local model boundary condition from the regional model [13][14][15] and it can keep consistent groundwater heads at shared nodes but requires block-shaped refined area and normally requires a longer simulation time for convergence. There are also case studies combining the TMR and LGR method for multi-scale modelling [16]. However, LGR method is not suitable when the target refined area is not a centralized area, for example, when an aquifer storage recovery wells system is distributed over a large area [17].
In recent years, the unstructured grid method (hereafter referred to a USG) in MODFLOW-USG was applied by many studies. The USG method is a fully coupling method that generates a more flexible mesh [18]. The USG method is numerically rigorous and has been applied for different kinds of modelling purposes. Since the design of the grid is very flexible, the grid refinement is no longer restricted by the shape grid and refinement area. Multi-scale models can be constructed by the USG method either when the boundary of the local model area is irregular or the hydrological stresses change in the local area might have a significant influence on the larger region [19,20]. Finer grid resolution can be used to refine the area near river channels [21,22] or only on the model top layer [23] to investigate surface and groundwater interaction. Except for the ortho-grid, the Voronoi grid can be applied to simulate the area of interest with higher resolution [24][25][26]. The USG method can also be used to simulate axisymmetric problems by using coaxial cylindrical cells [27]. The USG method potentially results in asymmetry and irregular banding matrix in solving equations. It also requires more complex discretization files so that the model construction can be arduous.
Few studies have systematically compared these methods. Sbai [28] evaluated the improvement of the performance of using USG to represent the aquifer heterogeneity and compared it with the homogeneous coarser grid model. Lux [29] analyzed the difference between using the structured grid and the unstructured grid to simulate the multi-lateral wells. The unstructured grid could provide higher resolution keeping cell numbers relatively low by only refining around the wells. Some synthesis cases were analyzed to evaluate the accuracy and CPU time [7,30].
In this study, the regional groundwater flow model of the Beijing Plain was coupled with a local refined grid model in the target area for MAR simulation. Like all the other mega-cities, with rapid urbanization and population growth, the city of Beijing is facing many environmental problems, such as a continuous groundwater level decline [31], land subsidence [32] and groundwater contamination. To evaluate the impact of those problems, a coupled local groundwater model with higher resolution is required. We applied the CGR method, the LGR methods and the USG method to refine the area with a planned MAR project and a large groundwater well field. The refined local models will be used to design and simulate the MAR scheme and evaluate the effectiveness of groundwater storage recovery at a groundwater well field near the MAR site. Three refinement methods were compared in terms of the goodness of fit of the computed and observed groundwater heads, the consistency of regional and local groundwater budgets, and the capture zone of the well filed. The advantages and limitations of three methods were discussed. The results of this study are helpful for the modelers to select a suitable refinement method according to the complexity of the study area and processes to be simulated.

The Regional Groundwater Flow Model
Beijing city (39 • 28"~41 • 05" N, 115 • 25"~117 • 30" E) is located in the northwest corner of the North China Plain. The total area of the city is 16,800 km 2 with 62% of mountainous area and 38% of plain area. With the rapid development and the increasing population of Beijing city, the water shortage has become one of the most important environmental problems. To assist in managing the groundwater resources in the Beijing Plain, a regional groundwater flow model has been constructed, which aims to simulate the alluvial groundwater system [33]. In this study, the grid orientation of the regional model has been rotated 33 degrees clockwise toward the north-east to minimize the number of column and rows to be refined in the CGR method and at the same time maximize the refined area of the LGR method, which will be described in the next session. Figure 1 shows the physical structure, boundary conditions, and flow processes of the hydrogeological conceptual model and the regional numerical model grid. Four aquitards were included in the model to represent the leakage between the five major aquifers. The top aquifer is distributed in the whole plain area while the extents of the lower aquifers are truncated by the bedrock. Hydraulic properties for nine model layers were assigned to parameter zones. The mountain front boundary in the west and north were defined as lateral inflow boundary while the administrative boundary in the south and east was simulated with the general head boundary. The hydrological stresses include areal recharge, river and canal leakage, ET and abstractions.
Grid cell size of the regional model were designed as 1000 m by 1000 m. A steady-state groundwater flow was simulated and calibrated with average fluxes in a long period of wet years, using MODFLOW-2005 [34]. In this study, Groundwater Modeling System (GMS) version 10.4 [35] was used to create the coupled models. With the conceptual model approach in GMS10.4, boundary conditions, parameters and stresses were created and stored in coverages, which are independent from the numerical model grid. When a coupled regional and local model grid was designed, data in the conceptual model were transferred automatically to the numerical model. In this way, three different coupled models were created more efficiently.

Coupling Regional and Local Models
Different local grid refinement methods were applied on the basis of this regional model, to incorporate a coupled local groundwater flow model the planned MAR site in the Chaobai River catchment. The location of the MAR site is close to the No.8 Well Field, where groundwater has been over-exploited since the continuous drought in 1999. To evaluate the change in the local groundwater flow system in response to the MAR project implementation, CGR method, LGR method and USG method were applied to create a local refined model in the MAR area coupled in the regional steady-state flow model. The area of the local model is around 240 km 2 and defined as the area in the vicinity of the MAR site and along the direction of the Chaobai River.

Conventional Grid Refinement (CGR) Method
The CGR method refines the local model area using a variably spaced grid, which can be realized by further dividing the rows and columns into smaller cells in the local model area. In this study, each model row and column of the regional model grid was divided into 10 rows and columns, so that the local grid resolution became 100 m by 100 m in the local model area. The refined model grid by the CGR method is shown in Figure 2. The refined grids go through rows and columns of the entire mode area and cut through all model layers. After mapping the conceptual model to the CGR model grid, the model was executed by MODFLOW-2005.

Local Grid Refinement (LGR) Method
LGR method can be realized by the computer program MODFLOW-LGR share node [36] method developed by USGS, which refines the local model area both horizontally and vertically. The number of layers to be refined can also be specified by the user. However, the shape of the local model is restricted to a blocked-shaped area and needs to be unified through all the refined model layers. The horizontal refinement factor of the local model must be an odd integer so that the edges of the local model grid can be located inside the shared nodes of the regional grid. Cells in the regional model that overlapped with the local model are abolished. However, some data input (that depends on the model cell area directly or indirectly) needs to be modified for the cells at the interface of the regional and local grids to prevent the double-counted sources and sinks. In this study, due to the decreasing layer extent of the lower model layers, only the first three layers of the regional model were refined. Figure 3 shows the shape of the refined local model of the LGR method. The refinement factor was set as 9 so that the size of the local grid is 111.11 m. The regional model cells that overlapped with the local model area were inactivated only in the first two layers. The local grid's layer bottom replaces 1/2 of the thickness of the regional model layer. Thus, the overlapped cells in the third layer of the regional model were remained active to allow the communication between the regional and local grids. When transferring the conceptual model data to the refined local grids using GMS 10.4, the sources and sinks of the local model area in layer 3 were also transferred to the regional grid, which were removed manually to avoid double account of those stresses. The cells on the interface of the regional and local models were defined as specified head boundary for the local model and the head values were obtained from the computed heads from the regional model. The LGR model can be run separately or jointly with the regional model. The simultaneous solution scheme provides feedback between the regional and local models [7], keeping consistent groundwater heads on the adjoining interface of the two models.

Unstructured Grid Refinement (USG) Method
The USG method can be applied by using the MODFLOW-USG version 1 [18] to simulate the groundwater flow that supports various structured and unstructured grid types. As a fully coupling refinement method, the USG method is the most flexible grid refinement scheme among the three methods. It does not have restrictions on the shape of the refined area and is also very flexible for vertical refinement. It only requires faces of horizontal top and bottom cells. The discretization can be the combination of an arbitrary number of nested grids or grids with other irregular shapes. And the flexible vertical refining capability of USG makes it easier to delineate the discontinuous confining units in the model domain. In this case, the refinement factor of the USG method can only be an even integer because the finer and coarser grids are embedded by shared faces. Thus, to acquire a relatively similar refined level with other methods in this study, the refinement factor of the local model by the USG method was chosen as 8. The quadtree grids used to create the local model have a 125 m grid size (Figure 4). All the model layers have been refined in the area of interest. Due to the decreasing layer extent of the deeper layers, the numbers of refined grids are less in the lower layers. To evaluate the planned MAR system in the Chaobai River, it is necessary to compare the groundwater travel time computed by local models. MODPATH version 3.0 [37] and Mod-PATH3DU [38] were applied to the CGR, and USG models to calculate groundwater travel times and delineate the well field capture zones. In each refined local model, 100 particles were placed in each extraction wells at the No.8 well field near the MAR site for backward particle tracking. The pathlines report provided information on the groundwater capture zone and travel time from the recharge area to the well field.

Results
The regional groundwater flow model was constructed and calibrated with alternative model method [33]. The model created with the true-layer approach was found most suitable for the simulation of MAR scheme. This model was used in this study to develop the coupled regional and local models.

Comparison of Computed and Observed Heads
First, three coupled models with CGR, LGR and USG methods were compared with the original regional model [33]. Figure 5 shows the scatter plots of the computed groundwater heads versus the observed heads from three coupled models. The coefficient of determinations of the original regional model was 0.830 for the entire model domain and 0.693 for the observation wells in the local model area [33]. The coefficient of determinations of three coupled models are similar with the original model while the coefficients of determinations of three local models are higher than the original regional model, indicating that the calibration accuracy of the original regional model is maintained by three coupled models. Table 1 lists statistics of three coupled models which show that both the mean error (ME) and the root mean squared errors (RMSE) of the local models are much smaller than the regional models, which confirm the necessity to create a local refined model for MAR simulation with high accuracy. Although the LGR model gives the lowest errors of the local model, other models are also acceptable for MAR simulation.

Comparison of Computed Groundwater Budgets
The groundwater balance components of three coupled models are shown in Table 2. In general, all major inflow and outflow components are the same as the original regional flow model, no changes were found comparing with the original regional model. Only very small differences were found in the areal recharge due to a small change in the model area by grid refinement. Some comparatively large differences were found in flows through general head boundary (GHB) and ET values since they are head dependent. However, both flow components are relatively small and do not influence the total groundwater budget. The groundwater balance components of three local models are shown in Table 3. The computed flow components in the local model area (239.75 km 2 ) from the CGR and USG methods are very close. The LGR model computed relatively smaller water budget since the local model area from LGR refinement is smaller (212.19 km 2 ). During the LGR grid refinement, cells in the regional model that overlapped with the local model were inactivated. However, the local model area is smaller than the inactivated area because the local model is embedded with the regional grid by shared-nodes (centroid of the cell). Thus, some recharge and abstractions located at the interfaces of the regional and local models were not counted in the local model. Differences in other head-dependent flow components are the consequences of smaller recharge and abstractions.  Figure 6 shows the computed groundwater head contour maps in layer 1 and layer 3 of the local model area from three coupled models.

Comparison of Computed Contour Maps
Generally, the computed groundwater head contour lines are similar from three local models. The shape of the contour lines of the USG model is smoother than the other two models. All three models computed a cone of depressions at the local model area, which is caused by the No.8 well field. Small differences in the computed heads were found at the centre of the cone of the depression. In the first model layer, the computed minimum groundwater heads of the CGR, LGR and USG model are 25

Comparison of Computed Capture Zones
As shown in Figure 7, the capture zones of the No.8 Well Field were delineated separately for the shallow (layer 1) and deeper (layer 3) aquifers. The areas of the capture zone delineated by the CGR method and USG method in layer 1 are 63.9 km 2 and 53.2 km 2 , which are 2% and 18% smaller than the capture zone generated by the original regional model (65.1 km 2 ). The capture zone delineated by the CGR method and the USG method in layer 3 are 310 km 2 and 296 km 2 , which are 27% and 1.7% larger than the capture zone generated by the original regional model (291 km 2 ). However, the shape of the capture zone generated the CGR and USG models are quite different. Groundwater in wells located in the shallow aquifer mainly comes from river leakage and precipitation recharge. Groundwater abstracted from the deeper aquifer comes mainly from leakage of the shallow aquifer. It seems that the river leakage has a larger influence on the capture zone delineation in the shallow aquifer with the particle tracking model mod-PATH3DU linked to the USG model. Pathlines generated with mod-PATH3DU cannot cross the river while MODPATH linked to the CGR model treats the river as a weak source and generated pathlines crossed the river. For the deeper aquifer, pathlines generated with mod-PATH3DU ended up at the water table while pathlines generated with MODPATH reached inflow boundary as source water. MODPATH (version 3) did not work with the LGR model so that it was not possible to delineate the capture zone with the LGR model.

Discussion
Based on the results shown above, we summarize the pros and cons of the different methods with regards to the following aspects: the complexity of model construction, the limitations and the conditions for application, and the shortcomings of the current version of GMS interface (GMS10.4).
Among the three refinement methods, the CGR method is the easiest and most straightforward to implement and requires the least modification to an existing regional model. No extra data or model input files need to be prepared. Under the GMS environment, each row and column of the regional model grid in the local model area can be selected and divided into a number of sub-grids to create a coupled regional and local model grid. Afterwards, coverages in the conceptual model can be mapped and transferred to MODFLOW packages to establish a coupled regional and local flow model. With the same convergence criteria as the regional model, the CGR coupled model can reproduce the same results as the original regional model. Moreover, the CGR method is very flexible to modify the existing local model and to add a new local model. Like in this study, the cone of depression around the No. 8 well field has been expanded larger and deeper with further over-exploitation. Thus, it is necessary to enlarge the local model area to include the enlarged cone of the depression for a transient simulation. This modification can be easily implemented in the transient model construction. Furthermore, the coupled model was developed from an existing regional model in this study, elevations of all model layers were re-interpolated to have more accurate elevations for the refined grids. However, the refined grids with the CGR method extends over the entire model area and cut through all model layers, a coupled model with much larger number of model cells is created resulting in longer computation time. Although the extra simulation time for a steady-state flow model is negligible, the computation time might increase a lot with transient model simulation, especially when coupled with solute and heat transport simulation. Another shortcoming is that the CGR method cannot refine model layers in the local model area. Thus, the CGR method is more efficient when there is only one local model, the local model area is relatively small, and no vertical refinement is required. By far, the CGR method is still widely used for the coupled modelling of groundwater flow, solute and heat transports since all standard model codes such as MODFLOW, MODPATH, MT3DMS and SEAWAT can be run with the CGR girds without any modifications.
The construction of the LGR model needs more modifications than the CGR model. Similar to the CGR method, the local model was constructed based on the regional model in this study. Thus, the top and bottom elevations of local model layers were interpolated again to obtain more accurate model structure. The advantage of the LGR method is that a number of the local models with desired refinements can be created and embedded into the regional model so that the regional model grid is not changed. The coupled LGR model can be run jointly or sequentially from the regional model to the local model. Thus, computation efficiency can be achieved by only running the local model while the regional model is not influenced by the local model. However, the construction of the LGR model has a number of restrictions. Firstly, the LGR method requires a blocked-shape of the local model area with the same horizontal refinement for selected model layers. A problem occurs when the local model area covers the irregular boundary or the areal extents of the model layers are different. Like in this study, the local model area is located in the upstream of the Chaobai River where the model boundary is irregular and the extent of the model layer decreases with an increase of the depth as the basement uplifting towards the boundary. Under this circumstance, the extent of the local model area was not able to cover the entire influence area of the No. 8 well field. Secondly, the hydrological stresses at the interface of the local and regional grid may be double counted and need to be checked and revised manually by the user. If there are large number of sources and sinks located in the interface cells, a large discrepancy of the groundwater balance can be expected. Thirdly, MODFLOW-LGR [7] couples the regional and local model at the interface cells. The interface cells at the edge of the local model were simulated as specified head cells with the heads being computed by the regional model. In the regional model, the interface cells are treated as specified flow cells with the flows being computed by the local model. Therefore, a large number of iterations might be required to achieve consistence of computed heads and flows at the interface cells. For coupled models with complex hydrogeological conditions, it is more difficult to meet the convergence criteria in the MODFLOW-LGR. So far, only MODPATH version 5 is coupled with MODFLOW-LGR to perform particle tracking, which has not been accommodated to the GMS 10.4. No solute, salt and heat transport models have been developed for coupling with MODFLOW-LGR. Thus, the application of LGR method is still limited to the coupled regional and local groundwater flow simulations in less complex hydrogeological conditions.
Although the USG method supports a wide variety of structured and unstructured grid types, including nested grids and grids based on prismatic triangles, rectangles, and hexagons shapes, to simulate rivers and wells with flexible refined grids, the construction of the USG model requires a big effort. The unstructured grid needs to be set up separately from the regional model grid. For a large regional groundwater model, creating a desirable unstructured grid local model takes a lot of time and effort and once it is constructed, further modification is not possible. At the same time, pre-processing of the model input data also requires that the modeller has a deeper understanding of the MODFLOW-USG software. However, the USG model is the most flexible method that can accommodate a variety of demands. The shape and location of the local model area have fewer restrictions compared with the LGR method. Especially with the problems that require multi-scale model simulation simultaneously, the USG model does not have any restriction on the number of embedded local models and the interested model area can be refined into different levels of resolution. In this study, the characteristics of the USG model provide great convenience to the simulation of the MAR system in the Chaobai River catchment. A multi-scale model can be created which couples the regional scale, local scale and MAR site scale models and they interact with each other simultaneously. However, with the currently USG version 1, some flow packages are not supported yet. The ability to simulate the unsaturated zone flow in the USG model is relatively weak. Only the SRF2 package is currently supported, which can simulate the unsaturated flow beneath a hydraulically disconnected stream and aquifer. Furthermore, a particle racking code (mod-PATH3DU) and a solute transport model code (MODFLOW-USG Transport) have been developed for coupling MODFLOW-USG model. These coupled models enable the USG models to simulate flow and solute transport with coupled multiple scale models. However, a density and viscosity dependent model has not yet been developed for simulating saltwater and heat transport with the USG method. Thus, the USG method has good potential for a wide range of applications in simulation of flow and solute transport in multiple scales.
Since all of the abovementioned models require text input files and save model results in data files, a good graphical user interface (GUI) for the model construction and results analysis is required. GMS10.4 is one of better commercial modelling environments and was used in this study. The conceptual model approach in GMS is well suited to construct the coupled regional and local models. GMS supports all available models (MODFLOW, MODPATH, MT3DMS, RT3D, and PHT3D) with the CGR method. GMS also supports flow and transport models (MODFLOW-USG, mod-PATH3DU, and MODFLOW-USG Transport) with the USG method. In the current GMS 10.4 version, only MODFLOW-LGR is supported for coupled flow simulation with the LGR method. A latest particle tracking code MODPATH version 5 was not supported yet so that particle tracking cannot be performed in GMS10.4. Furthermore, some sources and sinks were double counted at interface cells with the LGR method and have to be removed manually.

Conclusions
Based on the regional groundwater flow model of the Beijing Plain, a coupled local groundwater flow model with finer grid resolution was constructed with three different refinement methods: the conventional grid refinement (CGR) method based on MODFLOW-2005, the local grid refinement (LGR) method based on MODFLOW-LGR, and the unstructured grid (USG) method based on MODFLOW-USG. All three methods were able to create a coupled local model for the MAR site in the Chaobai River. The computed regional groundwater contour maps and flow budgets from three methods are very close to the original regional flow model. However, the LGR model produced a slightly smaller flow budget from the local model since a smaller local model area was created with the refinement and coupling method. Large differences were found in capture zone delineation for the No. 8 well field. These differences may be caused by using different particle tracking codes. MODPATH in the CGR method can track pathlines across weak sources and sinks. mod-PATH3DU coupled with MODFLOW-USG in the USG method terminates pathlines at weak sources and sinks. It was not possible to delineate the capture zone with the LGR method since a latest MODPATH code was not linked in GMS 10.4.
The comparative assessment of the three methods found that the CGR method is capable of constructing coupled regional and local models for flow, solute, salt and heat transport simulations, but inherits the shortcomings of regular finite difference grids. The USG method is very flexible to construct multi-scale models with unstructured grids fitting to rivers, wells, irregular boundaries, and heterogeneities. So far, the USG method can be used to simulate flow and solute transport. A good graphic interface is very important for the USG model construction and result analysis due to the complexity to create unstructured grids and prepare model inputs. The LGR method has limited applications because the corresponding transport model code has yet to be developed and the interfacing problem is still to be solved. Thus, the selection of a method should consider the requirement of the local model, complexity of the hydrogeological conditions, flow and transport processes to be simulated, and the available graphic user interface.
Finally, it is recommended that a density and viscosity dependent model code should be developed for the USG method. The popular graphic user interfaces should be timely updated to include latest model codes.

Conflicts of Interest:
The authors declare no conflict of interest.