Regional Groundwater Flow Assessment in a Prospective High-Level Radioactive Waste Repository of China

: The production of nuclear energy will result in high-level radioactive waste (HLRW), which brings potential environmental dangers. Selecting a proper disposal repository is a crucial step in the development of nuclear energy. This paper introduces ﬁrstly the hydrogeological conditions of the Beishan area in China. Next, a regional groundwater model is constructed using a multiphase ﬂow simulator to analyze the groundwater ﬂow pattern in the Beishan area. Model calibration shows that the simulated and observed hydraulic heads match well, and the simulated regional groundwater ﬂow pattern is similar to the surface ﬂow pattern from the channel network, indicating that the groundwater ﬂow is mainly dependent on the topography. In addition, the simulated groundwater storage over the period from 2003 to 2014 is similar to the trend derived from the Gravity Recovery and Climate Experiment satellite-derived results. Last, the established model is used to evaluate the inﬂuences of the extreme climate and regional faults on the groundwater ﬂow pattern. It shows that they do not have a signiﬁcant inﬂuence on the regional groundwater ﬂow patterns. This study will provide a preliminary reference for the regional groundwater ﬂow assessment in the site of the HLRW in China.


Introduction
Nuclear energy is considered as an additional dependable and clean energy source all over the world [1]. However, high-level radioactive waste (HLRW), which are the highly radioactive materials produced as a byproduct of the reactions that occur inside nuclear reactors, will probably bring potential dangers. Therefore, the proper disposal of nuclear waste is highly challenging. The safety of high-level radioactive waste disposal has gained much attention in the world. In America, Yucca Mountain was regarded as a potential repository site of high-level nuclear waste in 2002, with many experimental and numerical studies performed on groundwater flow. The initial assessment of the regional groundwater system at Yucca Mountains was carried out in 1997 [2], and then the comprehensive evaluation and simulation of the system was performed in 2010 by using MODFLOW-2000 [3]. In addition, a site-scale model for fluid flow in the unsaturated zone of Yucca Mountain, including fracture flow [4][5][6][7] and parallel computing [8], was developed by Transport Of Unsaturated Groundwater and Heat 2 (TOUGH2) to reproduce the overall system behavior. Finland was the first country to propose an ultimate disposal way for HLRW. In 2001, Finland selected a high-level radioactive waste disposal site based on field investigation and numerical simulation [9]. A kind of crystalline rock, which is tectonically and geochemically stable, was selected

Model boundaries and their representations
To obtain the drainage divide, the USGS/NASA Shuttle Radar Topography Mission (SRTM) data in 3 arc seconds is used to extract channel network and watershed information. The main procedures in Arc Hydro Tools in ArcGIS software include Digital Elevation Model (DEM) pretreatment, determination of the direction of flow, and extraction of flow accumulation, channel network, watersheds and watershed division. The colormap of the DEM and generated channel network are shown in Figure 2. The study area is surrounded by a natural surface divide in the north, the west and the partial east boundaries ( Figure 2). Considering human activities, such as groundwater pumping in Yumen City, the south boundary is set between Yumen City and Gobi. There are nearly no rivers in the study area. The top layer of this model is set as the atmosphere layer. The bottom layer is assumed to be a no-flow boundary. No long-term observation well data are available, and hydraulic heads at 33 open wells are collected from 2008 to 2009. The depth to groundwater in open wells almost varies from 5 m to 50 m, occasionally over hundreds of meters. Field investigation from 2012 to 2015 found that the inter-annual variation of water table in open wells is very small, usually less than 0.1 m, and yearly average of water table remain almost the same. According to the measured results [16], the depth to the groundwater level in BS09, BS12, BS14, BS05 are 59.5 m, 5.2 m, 36.6 m, 65 m, respectively. Vertical changes of soil moisture are not measured. Vertically, soil water usually changes from unsaturated state to fully-saturated state. As investigated by the Beijing Research Institute of Uranium Geology (2010) [16], the metamorphic rock and intrusion rock extend from the surface to the depth of over 600 m. The

Model Boundaries and Their Representations
To obtain the drainage divide, the USGS/NASA Shuttle Radar Topography Mission (SRTM) data in 3 arc seconds is used to extract channel network and watershed information. The main procedures in Arc Hydro Tools in ArcGIS software include Digital Elevation Model (DEM) pretreatment, determination of the direction of flow, and extraction of flow accumulation, channel network, watersheds and watershed division. The colormap of the DEM and generated channel network are shown in Figure 2. The study area is surrounded by a natural surface divide in the north, the west and the partial east boundaries ( Figure 2). Considering human activities, such as groundwater pumping in Yumen City, the south boundary is set between Yumen City and Gobi. There are nearly no rivers in the study area. The top layer of this model is set as the atmosphere layer. The bottom layer is assumed to be a no-flow boundary. No long-term observation well data are available, and hydraulic heads at 33 open wells are collected from 2008 to 2009. The depth to groundwater in open wells almost varies from 5 m to 50 m, occasionally over hundreds of meters. Field investigation from 2012 to 2015 found that the inter-annual variation of water table in open wells is very small, usually less than 0.1 m, and yearly average of water table remain almost the same. According to the measured results [16], the depth to the groundwater level in BS09, BS12, BS14, BS05 are 59.5 m, 5.2 m, 36.6 m, 65 m, respectively. Vertical changes of soil moisture are not measured. Vertically, soil water usually changes from unsaturated state to fully-saturated state.

Model boundaries and their representations
To obtain the drainage divide, the USGS/NASA Shuttle Radar Topography Mission (SRTM) data in 3 arc seconds is used to extract channel network and watershed information. The main procedures in Arc Hydro Tools in ArcGIS software include Digital Elevation Model (DEM) pretreatment, determination of the direction of flow, and extraction of flow accumulation, channel network, watersheds and watershed division. The colormap of the DEM and generated channel network are shown in Figure 2. The study area is surrounded by a natural surface divide in the north, the west and the partial east boundaries ( Figure 2). Considering human activities, such as groundwater pumping in Yumen City, the south boundary is set between Yumen City and Gobi. There are nearly no rivers in the study area. The top layer of this model is set as the atmosphere layer. The bottom layer is assumed to be a no-flow boundary. No long-term observation well data are available, and hydraulic heads at 33 open wells are collected from 2008 to 2009. The depth to groundwater in open wells almost varies from 5 m to 50 m, occasionally over hundreds of meters. Field investigation from 2012 to 2015 found that the inter-annual variation of water table in open wells is very small, usually less than 0.1 m, and yearly average of water table remain almost the same. According to the measured results [16], the depth to the groundwater level in BS09, BS12, BS14, BS05 are 59.5 m, 5.2 m, 36.6 m, 65 m, respectively. Vertical changes of soil moisture are not measured. Vertically, soil water usually changes from unsaturated state to fully-saturated state. As investigated by the Beijing Research Institute of Uranium Geology (2010) [16], the metamorphic rock and intrusion rock extend from the surface to the depth of over 600 m. The As investigated by the Beijing Research Institute of Uranium Geology (2010) [16], the metamorphic rock and intrusion rock extend from the surface to the depth of over 600 m. The permeability of original metamorphic and intrusion rock is low. However, because faults and fractures are occasionally distributed in the study area, hydraulic conductivity at the intersections of faults and fractures is relatively high. The metamorphic rock and intrusion rock can be simplified into three types according to logs of deep drillings [16], i.e., intensely weathered, middle weathered, and weakly weathered or even unweathered layers, where the thicknesses of the three layers are 50, 100, and 500 m, respectively. To better present the change of groundwater level and soil moisture in the unsaturated zone, the intensely weathered layer is refined to 15 model layers, where the thickness of each layer is 0.1, 1.4, 0.5, 1.0, 1.0, 2.0, 2.0, 2.0, 2.5, 2.5, 3.0, 3.0, 9.0, 10.0 and 10.0 m. The middle weathered layer is refined to five layers with the same thickness (20 m). The weakly weathered or unweathered layer is refined to three layers, and the thicknesses of the three layers are 100, 200 and 200 m ( Figure 3). In total, the model is divided into 23 model layers, so that the processes of the unsaturated-saturated interface can be identified by the transition from two-phase flow to single-phase flow in the numerical simulation.
Water 2017, 9,551 4 of 16 permeability of original metamorphic and intrusion rock is low. However, because faults and fractures are occasionally distributed in the study area, hydraulic conductivity at the intersections of faults and fractures is relatively high. The metamorphic rock and intrusion rock can be simplified into three types according to logs of deep drillings [16], i.e., intensely weathered, middle weathered, and weakly weathered or even unweathered layers, where the thicknesses of the three layers are 50, 100, and 500 m, respectively. To better present the change of groundwater level and soil moisture in the unsaturated zone, the intensely weathered layer is refined to 15 (Figure 3). In total, the model is divided into 23 model layers, so that the processes of the unsaturated-saturated interface can be identified by the transition from two-phase flow to singlephase flow in the numerical simulation.

Numerical Model
The parallel version of multiphase flow simulator is adopted in this study. TOUGH2 is a model based on the integral finite difference method, which could offer the advantage of being applicable to regular or irregular discretization in one, two and three dimensions [17]. The governing equations of TOUGH2 are established from mass and energy balance. Fluid flow is described with a multiphase extension of Darcy's law; in addition there is diffusive mass transport in all phases. Heat flow is governed by conduction and convection, also including sensible as well as latent heat effects. Time will be discretized fully implicitly as a first-order backward finite difference. The nonlinear equations in the residual form are solved by Newton/Raphson iteration [18]. The EOS3 module in TOUGH2 can simulate water and air movement in isothermal and nonisothermal conditions. The single and twophase conditions can be easily known from the range of a primary variable [17]. This module was employed to simulate the regional groundwater flow in the saturated-unsaturated zone of the study area, and it could describe fluid and heat flow in multiphase, multi-component systems [18,19]. The mass and energy equations for the two-phase flow of two components (water, air) in the TOUGH2/EOS3 module are summarized in Table 1 [20]. The unsaturated-saturated interface can be judged by the transition from the state of two-phase flow to that of single-phase flow.

Numerical Model
The parallel version of multiphase flow simulator is adopted in this study. TOUGH2 is a model based on the integral finite difference method, which could offer the advantage of being applicable to regular or irregular discretization in one, two and three dimensions [17]. The governing equations of TOUGH2 are established from mass and energy balance. Fluid flow is described with a multiphase extension of Darcy's law; in addition there is diffusive mass transport in all phases. Heat flow is governed by conduction and convection, also including sensible as well as latent heat effects. Time will be discretized fully implicitly as a first-order backward finite difference. The nonlinear equations in the residual form are solved by Newton/Raphson iteration [18]. The EOS3 module in TOUGH2 can simulate water and air movement in isothermal and nonisothermal conditions. The single and two-phase conditions can be easily known from the range of a primary variable [17]. This module was employed to simulate the regional groundwater flow in the saturated-unsaturated zone of the study area, and it could describe fluid and heat flow in multiphase, multi-component systems [18,19]. The mass and energy equations for the two-phase flow of two components (water, air) in the TOUGH2/EOS3 module are summarized in Table 1 [20]. The unsaturated-saturated interface can be judged by the transition from the state of two-phase flow to that of single-phase flow.

Description Equation
Conservation of mass and energy Mass accumulation where M K are the accumulation terms of the components and energy, kg m −3 for components and J m −3 for energy; k is the index for the components, =1 (water), 2 (air), and 3 (temperature); v n is the volume of the n th grid cell, m 3 ; F is the flux, g m −2 s −1 ; n is the normal vector on the surface element, dimensionless; Γ n is the area of closed surface, m 2 ; q k is the source/sink terms for mass or energy components, kg m −3 s −1 or J m −3 s −1 ; φ is the porosity, dimensionless; S β is the saturation of phase β , dimensionless; ρ β is the density of phase β, kg m −3 ; β = G for gas and = L for liquid; X k β is the mass fraction of component k in fluid phase β, dimensionless; u β is the Darcy velocity in phase β, m s −1 ; λ is the thermal conductivity, W K −1 m −1 ; u β is the viscosity, Pa s; h β is the specific enthalpy in phase β, J kg −1 ; C R is the heat conductivity, W K −1 m −1 ; T is temperature, • C or k; k γβ is the relative permeability to phase β; g is the gravitational acceleration constant, m s −2 ; and k is the absolute permeability, m 2 .

Elevation and Model Discretization
The topography is important for both regional and local groundwater flow patterns [21]. The elevation is derived from the SRTM data in 3 arc seconds, where Reuter recreated products by adopting void-filling interpolation methods to obtain higher-quality SRTM data [22]. The elevation changes from over 2300 m in the west to 900 m in the eastern part. By using IGMESH [23], which is a convenient irregular-grid-based pre-and post-processing tool for the TOUGH2 simulator, the study area is discretized into 12,390 gridblocks in the plane and 23 model layers in the vertical direction. There are 284,970 gridblocks and 1,120,107 connections in total ( Figure 4). The prospective site of the disposal repository is refined for better characterization of the groundwater flow.

Hydraulic Properties
Since 2005, the Beijing Research Institute of Uranium Geology has performed many field investigations via borehole drilling as well as hydrogeochemical and isotope sampling in an effort to identify a feasible site for high level radioactive waste disposal. The institute summarized hydraulic conductivity data from many pumping tests [16]. Hydraulic tests can be divided into three types: pumping tests, double ring infiltration tests and double packer hydraulic tests. There are 37 shallow wells with traditional pumping tests, and the locations of these tests are shown in Figure 1, where the range of the estimated hydraulic conductivity is from 0.0015 to 3.82 m/day. Double ring infiltration tests were conducted to estimate the hydraulic conductivity of the shallow soil in 2013 and 2014, and there are 113 test data points in total. The range of hydraulic conductivity estimated by the double

Hydraulic Properties
Since 2005, the Beijing Research Institute of Uranium Geology has performed many field investigations via borehole drilling as well as hydrogeochemical and isotope sampling in an effort to identify a feasible site for high level radioactive waste disposal. The institute summarized hydraulic Water 2017, 9, 551 6 of 15 conductivity data from many pumping tests [16]. Hydraulic tests can be divided into three types: pumping tests, double ring infiltration tests and double packer hydraulic tests. There are 37 shallow wells with traditional pumping tests, and the locations of these tests are shown in Figure 1, where the range of the estimated hydraulic conductivity is from 0.0015 to 3.82 m/day. Double ring infiltration tests were conducted to estimate the hydraulic conductivity of the shallow soil in 2013 and 2014, and there are 113 test data points in total. The range of hydraulic conductivity estimated by the double ring infiltration tests is 0.11-18.8 m/day, and the arithmetic mean value of hydraulic conductivity is approximately 4.27 m/day. Double packer hydraulic tests were performed at deep wells BS05 and BS06, where the depth of the drilling borehole is approximately 600 m below the surface. According to a large number of double packer hydraulic tests, the mean values of the hydraulic conductivities at BS05 and BS06 are 0.0088 and 0.00018 m/day, respectively. Considering all the test results, the change of permeability with the depth to the surface obtained from all test data [16] is summarized in Figure 5. The permeability is found to have an almost exponentially deceasing trend with the depth. The permeability obtained from these wells are in the range of 1 × 10 −11 to 1 × 10 −18 m 2 .

Hydraulic Properties
Since 2005, the Beijing Research Institute of Uranium Geology has performed many field investigations via borehole drilling as well as hydrogeochemical and isotope sampling in an effort to identify a feasible site for high level radioactive waste disposal. The institute summarized hydraulic conductivity data from many pumping tests [16]. Hydraulic tests can be divided into three types: pumping tests, double ring infiltration tests and double packer hydraulic tests. There are 37 shallow wells with traditional pumping tests, and the locations of these tests are shown in Figure 1, where the range of the estimated hydraulic conductivity is from 0.0015 to 3.82 m/day. Double ring infiltration tests were conducted to estimate the hydraulic conductivity of the shallow soil in 2013 and 2014, and there are 113 test data points in total. The range of hydraulic conductivity estimated by the double ring infiltration tests is 0.11-18.8 m/day, and the arithmetic mean value of hydraulic conductivity is approximately 4.27 m/day. Double packer hydraulic tests were performed at deep wells BS05 and BS06, where the depth of the drilling borehole is approximately 600 m below the surface. According to a large number of double packer hydraulic tests, the mean values of the hydraulic conductivities at BS05 and BS06 are 0.0088 and 0.00018 m/day, respectively. Considering all the test results, the change of permeability with the depth to the surface obtained from all test data [16] is summarized in Figure 5. The permeability is found to have an almost exponentially deceasing trend with the depth. The permeability obtained from these wells are in the range of 1 × 10 −11 to 1 × 10 −18 m 2 . Vertically, there are 23 model layers in total. Table 2 lists the vertical model layers and the zonation of hydraulic properties for each layer. The lithology zoning map at the 2nd layer and the Vertically, there are 23 model layers in total. Table 2 lists the vertical model layers and the zonation of hydraulic properties for each layer. The lithology zoning map at the 2nd layer and the 3rd layer ( Figure 6) are made from the geological zones in Figure 1. The rock grain density, formation heat conductivity and pore compressibility are assumed to be uniform, at 2650 kg·m −3 , 1 W·m −1 · • C −1 and 3.7 × 10 −10 Pa −1 , respectively. The heat capacity of the top layer is set as 8.0 × 10 99 J·kg −1 · • C −1 , and the values of the other layers are set as 1000 J·kg −1 · • C −1 . The relative permeability and capillary pressure were calculated by the van Genuchten-Mualem model [24]. In addition, related parameters and the pore compression coefficient of the model were obtained according to the empirical value of soil. Because the study area is almost composed of low-permeability media, such as metamorphic and intrusive rocks, parameters for estimating relative permeability in the van Genuchten-Mualem model remain the same for each type of rock. Among them, λ (a parameter in Van Genuchten-Mualem function) is equal to 0.492. Residual liquid saturation S lr and residual gas saturation S gr are both set as 0.01, and saturated liquid saturation S ls in the Van Genuchten-Mualem Model is 1.0. Initial capillary pressure P 0 and maximum capillary pressure P max are set as 1.0 × 10 3 Pa and 1.0 × 10 12 Pa also based on the empirical value, respectively. and intrusive rocks, parameters for estimating relative permeability in the van Genuchten-Mualem model remain the same for each type of rock. Among them,  (a parameter in Van Genuchten-Mualem function) is equal to 0.492. Residual liquid saturation Slr and residual gas saturation Sgr are both set as 0.01, and saturated liquid saturation Sls in the Van Genuchten-Mualem Model is 1.0. Initial capillary pressure P0 and maximum capillary pressure Pmax are set as 1.0 × 10 3 Pa and 1.0 × 10 12 Pa also based on the empirical value, respectively.  Table 2.

Regional Faults
Geological investigation [15] verified that fractures and regional faults are developed in the Beishan area ( Figure 1). Usually, there are granite and quartz fractures with filling materials. The field investigation in 2013 and 2014 and previous geological studies (scale is 1:500,000) found that the regional faults extend along the south-eastern, north-eastern, and east-western directions, and the inclination of those faults is mainly between 60° and 70°. In the model, the regional faults are assumed

Regional Faults
Geological investigation [15] verified that fractures and regional faults are developed in the Beishan area ( Figure 1). Usually, there are granite and quartz fractures with filling materials. The field investigation in 2013 and 2014 and previous geological studies (scale is 1:500,000) found that the regional faults extend along the south-eastern, north-eastern, and east-western directions, and the inclination of those faults is mainly between 60 • and 70 • . In the model, the regional faults are assumed to extend vertically to a depth of approximately 80 m. In addition, the widths of the faults are usually greater than 0.01 m on surface rocks, and then are reduced to 0.001-0.003 m in the depth of 650 m. However, the properties of faults are not fully recognized. In this model, three groups of regional faults (Figure 1) are modeled by the method of Equivalent Permeability Medium, and key parameters are listed in Table 2 (zone 10). They are considered to analyze their influences on the regional groundwater flow pattern.

Initial Conditions and Boundary Conditions
The initial pressure field is assumed to be under hydrostatic equilibrium by running the model over hundred thousands of years with yearly average rainfall infiltration (90 mm). The temperature field is determined from a surface temperature of 15 • C and a common geothermal gradient of 30 • C/km. This model is assumed to be isothermal. The hydraulic head can be converted from pressure according to Equation (1). The fluxes at the southern and partial eastern boundaries are difficult to ascertain. Thus, the model sets gridblocks at the boundary as constant pressure in TOUGH2, which is estimated by field investigations.
where P is the water pressure, Pa; P 0 is the standard atmospheric pressure, Pa; ρ is the water density, kg/m 3 ; g is the gravitational acceleration, m/s 2 ; and H is the hydraulic head, m.

Source and Sinks
The main groundwater recharge item is precipitation infiltration. Groundwater discharge mainly includes evapotranspiration and lateral outflow. Monthly precipitation data over the study area are obtained from data of the Global Precipitation Climatology Project (GPCP), which is a section of the WCRP (the World Climate Research program). The GPCP dataset integrated infrared and microwave satellite observations, and 6000 global conventional observation data points are proven to be reliable precipitation data sources [25][26][27]. The change of monthly rainfall from 2003 to 2014 is shown in Figure 7, and the yearly average rainfall is approximately 90 mm. Most rainfall is from June to September.
The methods for recharge estimation include in situ experiments of tritium migration (TM) and chloride mass balance (CMB) [28]. According to the tests, the estimated infiltration recharge rates are approximately 2.55 to 6.57 mm/year based on in situ experiments of TM, 0.05 to 1.57 mm/year based on the CMB method in unsaturated zones, and 0.287 mm/year based on the CMB method in saturated zones [14]. So the infiltration rates from the precipitation on the top boundary are set as 0.01 and 0.004 for Quaternary and the weakly permeable strata respectively. groundwater flow pattern.

Initial Conditions and Boundary Conditions
The initial pressure field is assumed to be under hydrostatic equilibrium by running the model over hundred thousands of years with yearly average rainfall infiltration (90 mm). The temperature field is determined from a surface temperature of 15 °C and a common geothermal gradient of 30 °C/km. This model is assumed to be isothermal. The hydraulic head can be converted from pressure according to Equation (1). The fluxes at the southern and partial eastern boundaries are difficult to ascertain. Thus, the model sets gridblocks at the boundary as constant pressure in TOUGH2, which is estimated by field investigations.
where P is the water pressure, Pa; P0 is the standard atmospheric pressure, Pa; ρ is the water density, kg/m 3 ; g is the gravitational acceleration, m/s 2 ; and H is the hydraulic head, m.

Source and Sinks
The main groundwater recharge item is precipitation infiltration. Groundwater discharge mainly includes evapotranspiration and lateral outflow. Monthly precipitation data over the study area are obtained from data of the Global Precipitation Climatology Project (GPCP), which is a section of the WCRP (the World Climate Research program). The GPCP dataset integrated infrared and microwave satellite observations, and 6000 global conventional observation data points are proven to be reliable precipitation data sources [25][26][27]. The change of monthly rainfall from 2003 to 2014 is shown in Figure 7, and the yearly average rainfall is approximately 90 mm. Most rainfall is from June to September.
The methods for recharge estimation include in situ experiments of tritium migration (TM) and chloride mass balance (CMB) [28]. According to the tests, the estimated infiltration recharge rates are approximately 2.55 to 6.57 mm/year based on in situ experiments of TM, 0.05 to 1.57 mm/year based on the CMB method in unsaturated zones, and 0.287 mm/year based on the CMB method in saturated zones [14]. So the infiltration rates from the precipitation on the top boundary are set as 0.01 and 0.004 for Quaternary and the weakly permeable strata respectively.

Model Calibration
The time period for model calibration is set from January of 2003 to December of 2014. In this study, calibration targets include observed hydraulic heads, and GRACE (Gravity Recovery and Climate Experiment)-derived groundwater storage.

Comparison of Simulated and Observed Hydraulic Heads
The simulated pressure from TOUGH2 software is converted to the hydraulic head by Equation (1). A comparison of the simulated and observed hydraulic heads at 33 open wells is shown inFigure 8. Here, the relative error is defined as the ratio of the absolute error of observed and simulated values to the observed values. Of the simulated data, 91% are in the range of the relative error in 2%. Overall, the fit between the simulated values and observed values is reasonable. The poor match between the observed and simulated results may be caused by the representations in the models of coarse resolution of hydrogeology zonation, averagely allocation of sources and sinks over a region, and heterogeneity of rock.

Comparison of Simulated and Observed Hydraulic Heads
The simulated pressure from TOUGH2 software is converted to the hydraulic head by Equation (1). A comparison of the simulated and observed hydraulic heads at 33 open wells is shown in Figure 8. Here, the relative error is defined as the ratio of the absolute error of observed and simulated values to the observed values. Of the simulated data, 91% are in the range of the relative error in 2%. Overall, the fit between the simulated values and observed values is reasonable. The poor match between the observed and simulated results may be caused by the representations in the models of coarse resolution of hydrogeology zonation, averagely allocation of sources and sinks over a region, and heterogeneity of rock. The change of soil moisture with depth at the end of 2014 is shown in Figure 9. From Figure 9, the depth of water table is basically in accordance with the measured depth to the groundwater level [16], which reveals that the simulated values fit the measured values well. In addition, some scholars performed tests on the soil permeability of the aeration zone in the study area [14], and they found that the soil permeability varied significantly at the spatial scale. As the graph shows, there is an obvious sudden change at the depth of 2 m for BS12 and approximately 10 m for the other boreholes. Because the shallow layer is composed of sand, clay and some intensely weathered granite and the permeability of these rocks are different, there is also some difference among the liquid saturation levels. For BS12, the depth of sand and clay is approximately 2 m; thus, a change occurs at that depth, and then, the liquid saturation reached saturation quickly. Regarding the other boreholes, the granite weathering degree is more serious than that below 10 m; thus, an abrupt change occurred at approximately 10 m below the surface. The change of soil moisture with depth at the end of 2014 is shown in Figure 9. From Figure 9, the depth of water table is basically in accordance with the measured depth to the groundwater level [16], which reveals that the simulated values fit the measured values well. In addition, some scholars performed tests on the soil permeability of the aeration zone in the study area [14], and they found that the soil permeability varied significantly at the spatial scale. As the graph shows, there is an obvious sudden change at the depth of 2 m for BS12 and approximately 10 m for the other boreholes. Because the shallow layer is composed of sand, clay and some intensely weathered granite and the permeability of these rocks are different, there is also some difference among the liquid saturation levels. For BS12, the depth of sand and clay is approximately 2 m; thus, a change occurs at that depth, and then, the liquid saturation reached saturation quickly. Regarding the other boreholes, the granite weathering degree is more serious than that below 10 m; thus, an abrupt change occurred at approximately 10 m below the surface.  Figure 10 shows the groundwater flow pattern in the Beishan area. The groundwater level in Figure 10 ranges approximately from 700 m to 2300 m (a.s.l.), with the lowest in the northeast and the highest in the mid-west. Then, the groundwater system in the study area can be divided into three local groundwater systems. Local groundwater in unit I and II flow towards the southeast and the northeast respectively, and then flow into the downstream of the Heihe River Basin. Groundwater in unit Ⅲ flows to the south boundary. Meanwhile, the prospective site for high-level radioactive waste disposal is located in the unit Ⅲ

Groundwater Flow Patterns
, so the hydrostratigraphic units should be further investigated for finalizing site selection. Comparing Figure 10 and Figure 2, the groundwater flow is found to have almost the same pattern as surface flow, suggesting that groundwater flow is mainly dependent on  Figure 10 shows the groundwater flow pattern in the Beishan area. The groundwater level in Figure 10 ranges approximately from 700 m to 2300 m (a.s.l.), with the lowest in the northeast and the highest in the mid-west. Then, the groundwater system in the study area can be divided into three Water 2017, 9, 551 10 of 15 local groundwater systems. Local groundwater in unit I and II flow towards the southeast and the northeast respectively, and then flow into the downstream of the Heihe River Basin. Groundwater in unit III flows to the south boundary. Meanwhile, the prospective site for high-level radioactive waste disposal is located in the unit III, so the hydrostratigraphic units should be further investigated for finalizing site selection. Comparing Figures 2 and 10, the groundwater flow is found to have almost the same pattern as surface flow, suggesting that groundwater flow is mainly dependent on the topography.  Figure 10 shows the groundwater flow pattern in the Beishan area. The groundwater level in Figure 10 ranges approximately from 700 m to 2300 m (a.s.l.), with the lowest in the northeast and the highest in the mid-west. Then, the groundwater system in the study area can be divided into three local groundwater systems. Local groundwater in unit I and II flow towards the southeast and the northeast respectively, and then flow into the downstream of the Heihe River Basin. Groundwater in unit Ⅲ flows to the south boundary. Meanwhile, the prospective site for high-level radioactive waste disposal is located in the unit Ⅲ

Groundwater Flow Patterns
, so the hydrostratigraphic units should be further investigated for finalizing site selection. Comparing Figure 10 and Figure 2, the groundwater flow is found to have almost the same pattern as surface flow, suggesting that groundwater flow is mainly dependent on the topography.

Change of Groundwater Storage
Limited point-based observation well data make it difficult to sufficiently calibrate the numerical model in the study area. In this situation, the simulated groundwater storage changes are compared with the GRACE-derived results. The Gravity Recovery and Climate Experiment (GRACE) satellite can monitor the time-dependent earth's gravitational field and then obtain the groundwater storage variation by inversion calculation [29,30]. The GRACE-derived terrestrial water storage (TWS) includes surface water, groundwater, soil water and snow water. Recent studies verified that TWS

Change of Groundwater Storage
Limited point-based observation well data make it difficult to sufficiently calibrate the numerical model in the study area. In this situation, the simulated groundwater storage changes are compared with the GRACE-derived results. The Gravity Recovery and Climate Experiment (GRACE) satellite can monitor the time-dependent earth's gravitational field and then obtain the groundwater storage variation by inversion calculation [29,30]. The GRACE-derived terrestrial water storage (TWS) includes surface water, groundwater, soil water and snow water. Recent studies verified that TWS can be directly estimated from GRACE data for certain large basins when appropriate smoothing is applied to reduce striping [31], and the TWS changes agree reasonably well, at least for long wavelengths, with those estimated from hydrology models and observation networks [29,[32][33][34]. In data-poor regions of the world, such as desert and some remote mountainous regions, GRACE may represent the only option for regional groundwater assessments [35]. Because the study area is in an arid region, and the infiltration rate is so small that there are nearly no rivers; thus, the changes of both surface water and snow water are set as zero [36]. The GRACE-derived data in the study area represent the sum of groundwater and soil water. The numerical model in the saturated-unsaturated zone of the study area can output the changes of groundwater and soil water. A comparison of the groundwater and soil water storage changes with the GRACE-derived values is shown in Figure 11, which shows that the simulated storage trend is accordance with the GRACE-derived TWS variations. can be directly estimated from GRACE data for certain large basins when appropriate smoothing is applied to reduce striping [31], and the TWS changes agree reasonably well, at least for long wavelengths, with those estimated from hydrology models and observation networks [29,[32][33][34]. In data-poor regions of the world, such as desert and some remote mountainous regions, GRACE may represent the only option for regional groundwater assessments [35]. Because the study area is in an arid region, and the infiltration rate is so small that there are nearly no rivers; thus, the changes of both surface water and snow water are set as zero [36]. The GRACE-derived data in the study area represent the sum of groundwater and soil water. The numerical model in the saturated-unsaturated zone of the study area can output the changes of groundwater and soil water. A comparison of the groundwater and soil water storage changes with the GRACE-derived values is shown in Figure 11, which shows that the simulated storage trend is accordance with the GRACE-derived TWS variations. Figure 11. Comparison between the simulated and GRACE-derived TWS variations.

Groundwater Budgets
Groundwater budgets play important roles in analyzing regional groundwater flow. Groundwater budgets can be used to evaluate the generation of the groundwater and guide the Figure 11. Comparison between the simulated and GRACE-derived TWS variations.

Groundwater Budgets
Groundwater budgets play important roles in analyzing regional groundwater flow. Groundwater budgets can be used to evaluate the generation of the groundwater and guide the groundwater management, especially in a large-scale region. In the Beishan area, the main groundwater recharge item is precipitation infiltration, and groundwater discharge mainly includes evapotranspiration and lateral outflow. The Beijing Research Institute of Uranium Geology conducted some field tests to obtain the precipitation infiltration recharge rate of unsaturated soil in the Beishan area.
According to the principle of water balance, the water resources budget can be expressed in Equation (2).
where P is the annual average precipitation, m; W is the study area, m 2 ; ∆Q is the total groundwater storage variation (including soil water and saturated groundwater), m 3 ; O is the outflow volume at the Dirichlet boundaries, m 3 ; and S denotes the surface runoff and evaporation, m 3 . In the numerical model, the yearly averaged precipitation is 90 mm. The yearly averaged outflow at the Dirichlet boundaries is approximately 1.03 × 10 6 m 3 . Changes in the soil water and saturated groundwater are approximately −3.22 × 10 5 m 3 /year. Thus, from Equation (2), surface runoff and evaporation account for approximately 99.99% of precipitation. Over the study area, approximately 0.01%-0.05% of the precipitation infiltrates into the groundwater.
For assessing the influence of the south boundary on groundwater flow in the Beishan area, different yearly annual precipitations (Table 3) are set to estimate the groundwater budgets. The results suggest the south boundary is not sensitive for the groundwater flow pattern in the Beishan area.

Influences of Climate Change on Groundwater Flow Pattern
The safety of the disposal repository in the Beishan area is the primary concern. As discussed earlier, the main recharge of groundwater is precipitation. As a result, the climate changes under extreme conditions, such as wet and dry conditions, should be considered to evaluate their influences on the regional groundwater flow [37]. A long-term calculation is thus required for the variation of climate; as a result, the prediction period is set as 10,000 years. Five precipitation scenarios are set with yearly value at 90, 40, 200, 300 and 450 mm for S1, S2, S3, S4 and S5 scenario, respectively. S1 is the base scenario. S2-S5 represent different annual mean precipitation scenarios. Figure 12 reveals the liquid saturation at the depth of 2 m in two boreholes. The result suggests that liquid saturation causes a change in a small range for each scenario, and higher precipitation will incur an increase of the liquid saturation level. The highest difference of liquid saturation between S5 and S1 is 7.0 × 10 −4 . The changes of liquid saturation with time at the two boreholes are very similar. From the pattern of the water table over the study area, the maximum difference of the water table relative to the base case is approximately 0.015 m for 10,000 years. The average groundwater velocity at the prospective site is approximately 4.66 × 10 −13 m/s, and it shows no obvious change under five climate scenarios, suggesting that climate changes will have a limited effect on the groundwater flow in the Beishan area.
the base scenario. S2-S5 represent different annual mean precipitation scenarios. Figure 12 reveals the liquid saturation at the depth of 2 m in two boreholes. The result suggests that liquid saturation causes a change in a small range for each scenario, and higher precipitation will incur an increase of the liquid saturation level. The highest difference of liquid saturation between S5 and S1 is 7.0 × 10 −4 . The changes of liquid saturation with time at the two boreholes are very similar. From the pattern of the water table over the study area, the maximum difference of the water table relative to the base case is approximately 0.015 m for 10,000 years. The average groundwater velocity at the prospective site is approximately 4.66 × 10 −13 m/s, and it shows no obvious change under five climate scenarios, suggesting that climate changes will have a limited effect on the groundwater flow in the Beishan area.

Influences of Faults on Regional Groundwater Flow
Because the properties of three faults are not fully understood, the method of scenario analysis is used to ascertain the influences of regional faults on the groundwater flow patterns. Four scenarios (shown in Table 4) are set by changing the permeability of different faults (S2-S5), and S1 is set as the base scenario. All models are run for 10,000 years. The changes of liquid saturation with time at the BS05 and BS14 wells under the four scenarios are nearly the same. To determine the influence on the regional groundwater level, the predicted contour maps of the water table difference in these four scenarios relative to the base scenario are shown in Figure 13. The result reveals that Fault 3 plays a significant role in the groundwater flow pattern, with the maximum variation of 30 m, and S5

Influences of Faults on Regional Groundwater Flow
Because the properties of three faults are not fully understood, the method of scenario analysis is used to ascertain the influences of regional faults on the groundwater flow patterns. Four scenarios (shown in Table 4) are set by changing the permeability of different faults (S2-S5), and S1 is set as the base scenario. All models are run for 10,000 years. The changes of liquid saturation with time at the BS05 and BS14 wells under the four scenarios are nearly the same. To determine the influence on the regional groundwater level, the predicted contour maps of the water table difference in these four scenarios relative to the base scenario are shown in Figure 13. The result reveals that Fault 3 plays a significant role in the groundwater flow pattern, with the maximum variation of 30 m, and S5 considering all three faults together is also similar with S4. However, the permeability values of Fault 1 and Fault 2 have a small influence on the water table, with the greatest variation of 6.5 m. According to the flow directions of groundwater system (Figure 10), faults 1 and 2 occupy the topographically high area between regions I and II such that groundwater drains away from them. Only fault 3 lies transverse to regional flow paths, and so the influence of fault 3 on the regional groundwater flow is significantly. Overall, the influences of the permeability of the three faults are mostly located at the region near the faults, and they will not change the pattern of regional groundwater flow. Table 4. Permeability distribution under four scenarios (m 2 ).

Conclusions
The proper disposal of high-level radioactive wastes is highly challenging. Numerical simulation is a wise method to evaluate the safety of a high-level radioactive waste repository [38]. The numerical method is also an effective approach to ascertain the groundwater flow patterns and the influences of climate change on the groundwater flow. In this study, the regional groundwater flow model in the Beishan area, a prospective disposal site for HLRW, was constructed based on the TOUGH2-MP/EOS3 module, which could describe fluid flow in multiphase, multi-component systems. The arid study area covers approximately 68,900 km 2 . This model was divided into 23 layers in the vertical direction to allow the unsaturated and saturated zones to be simulated accurately. The parameters of the model were identified by a series of investigations and statistical analysis. Initial pressure field was obtained by running the established model over 10,000 years after comparing a large period, and the differences of the pressure field between two periods were very small, so the initial pressure had been considered as a steady state.
The calibration targets included the hydraulic head from observation wells, and the GRACEderived groundwater storage changes. The good match between the simulated and observed results proved the reliability of this model. In addition, the model revealed that the regional groundwater flow pattern is quite similar to the surface flow from the channel network given by the DEM model, Figure 13. Contour map of the water table difference in the scenarios relative to the base case: (a) water table difference between S2 and S1; (b) water table difference between S3 and S1; (c) water table difference between S4 and S1; and (d) water table difference between S5 and S1.

Conclusions
The proper disposal of high-level radioactive wastes is highly challenging. Numerical simulation is a wise method to evaluate the safety of a high-level radioactive waste repository [38]. The numerical method is also an effective approach to ascertain the groundwater flow patterns and the influences of climate change on the groundwater flow. In this study, the regional groundwater flow model in the Beishan area, a prospective disposal site for HLRW, was constructed based on the TOUGH2-MP/EOS3 module, which could describe fluid flow in multiphase, multi-component systems. The arid study area covers approximately 68,900 km 2 . This model was divided into 23 layers in the vertical direction to allow the unsaturated and saturated zones to be simulated accurately. The parameters of the model were identified by a series of investigations and statistical analysis. Initial pressure field was obtained by running the established model over 10,000 years after comparing a large period, and the differences of the pressure field between two periods were very small, so the initial pressure had been considered as a steady state.
The calibration targets included the hydraulic head from observation wells, and the GRACE-derived groundwater storage changes. The good match between the simulated and observed results proved the reliability of this model. In addition, the model revealed that the regional groundwater flow pattern is quite similar to the surface flow from the channel network given by the DEM model, suggesting that groundwater is mainly dependent on the topography. The study demonstrated that only 0.01%-0.05% of the precipitation infiltrates into the groundwater, which is mainly caused by the wide distribution of low permeability media. Model predictions showed that extreme climate and the permeability of three groups of faults will not have a noticeable effect on the regional groundwater flow in the Beishan area within a period of 10,000 years.
In this paper, field investigation data were integrated to establish the numerical model to fully analyze and evaluate the regional groundwater flow pattern under extreme climate and uncertain fault properties, although the information on hydrogeological conditions is limited. This study is focused on analyzing the regional groundwater flow patterns, refining the groundwater flow systems and testing the influences of hypothesis that both climate change and the fault zone hydraulic properties. Results show three local groundwater systems are subdivided from the regional groundwater flow. To finalize the site selection, a series of chemical and physical experiments, and much data of observation well are highly required. It is suggested that experiments and refined field investigations could be carried out in the region of local flow system, such as the unit III. The study also found that local groundwater flow is sensitive to faults developed in the study area. A detailed investigation on fracture distribution, and development of a refined model, such as dual porosity or dual permeability model in the local area, are required to be carried out.