Modeling the Spatial and Seasonal Variations of Groundwater Head in an Urbanized Area under Low Impact Development

: Increasing impervious land cover has great impacts on groundwater regimes in urbanized areas. Low impact development (LID) is generally regarded as a sustainable solution for groundwater conservation. However, the effects of LID on the spatial-temporal distribution of groundwater are not yet fully understood. In this case study, a coupled Storm Water Management Model (SWMM) and Finite Element Subsurface FLOW system (FEFLOW) model was used to simulate surface and groundwater ﬂow in an urbanized area in Shenzhen, China. After veriﬁcation, the model was used to analyze the spatial-seasonal variations of groundwater head and hydrological processes under different LID scenarios. The results indicate that if the runoff from 7.5% and 15% of impervious area is treated by LID facilities, the annual surface runoff decreases by 5% and 9%, respectively, and the spatial average groundwater head relative to sea level pressure increases by 0.9 m and 1.7 m in the study area, respectively. The rise in groundwater head generally decreases from the recharge zones to the discharge zones surrounded by the streams and coastal waters. However, the groundwater head change is determined not only by the location in the catchment, but also by the hydraulic conductivity of underlying aquifer and LID inﬁltration intensity. Moreover, LID signiﬁcantly enhances groundwater recharge and aquifer storage in the wet seasons; in turn it increases aquifer release and groundwater discharge in the dry seasons. However, LID has the potential to increase the risk of groundwater ﬂooding during wet seasons in areas with poor aquifer drainage capacity and shallow groundwater depth. The ﬁndings from this study provide the basis for further assessing the beneﬁt and risk of LID inﬁltration for groundwater supplementation in the urbanized areas.


Introduction
Impervious land cover, for example in the form of roads, roofs, and parking lots, significantly decreases rainwater infiltration, restricts rainwater recharge to groundwater, and results in a groundwater decrease in urbanized areas [1,2]. Increasing groundwater deficit in urban aquifers causes serious risks to the ecological environment of wetlands and streams fed by groundwater [3]. Low impact development (LID) has been recommended as an innovative solution for storm water management [3,4]. LID facilities can store, infiltrate, and retain the storm water by applying decentralized and small-scale ecological systems such as bioretention basins, permeable pavements, vegetated swales, and green roofs [5]. LID also has capacity to enhance groundwater recharge and increase base flow during dry seasons [6,7]. Recharging groundwater through LID has been promoted as an alternative to centralize rainwater infiltration facilities in some urban areas with groundwater depletion issues [8,9].
The effects of LID on groundwater systems have been evaluated in many previous studies based on experimental observation or model simulation. Most of the studies focused on single LID facilities and investigated infiltration processes and groundwater mounding beneath the facility [10,11]. For example, based on groundwater monitoring beneath an experimental storm water infiltration basin, Machusick et al. [12] reported that the extent of groundwater mounding was localized to the infiltration basin and over a limited distance. Newcomer et al. [6] carried out an in situ observation of the infiltration process beneath an 11 m 2 best management practice (BMP) infiltration trench and then applied a Hydrus-2D model to predict the change of groundwater recharge under climate variability and change.
Due to the lack of observation data for large areas with extensive LID infiltration, hydrological models have been widely used to investigate the impacts of LID on groundwater at the urban catchment scale. Most catchment-scale model research shows that LID infiltration has the potential to increase groundwater level. Thomas et al. [9] developed a regional multivariate regression model to predict the changes of groundwater level after implementing recharge best management practices (BMPs) in the Back Bay region of Boston, and the results indicated that recharge BMPs have a limited effect on groundwater elevations. Endreny and Collins [13] used a MODFLOW model to analyze the rise magnitude of groundwater mounding in an 8-ha drainage area in New York considering different arrangements of bioretention basins and different hydraulic conductivity. They found that groundwater elevation could increase 1.1 m after 30 years of recharge. Based on a distributed hydrological model, Locatelli et al. [14] found that land use change associated with widespread storm water infiltration increased the local infiltration and caused an overall increase in groundwater in a catchment in Perth, Western Australia. Increasing groundwater storage by LID infiltration helps the restoration of base flow in urban rivers [15]. Kidmose et al. [16] applied a fully coupled MIKE URBAN and MIKE SHE model to explore the effects of forced infiltration on urban water balance in Silkborg City, Denmark. The results showed that the groundwater discharge to runoff network could increase by up to 45% due to the growth of groundwater recharge, and the water level could rise up to 2-3 m in the area with poor soil permeability. Göbel et al. [8] used Hydurs-2D, GwNeu, and SPRING to assess the effects of decentralized stormwater infiltration through swales and infiltration trenches on the water budget and groundwater surface in an urban area, Germany, predicting an increase of base flow by 40% with local groundwater mounding up to 2.9 m.
Although many previous studies investigated the effects of LID on the changes of groundwater head at a catchment or an urban scale, few of them explored spatial-temporal variations of groundwater head and hydrological processes under LID [13,17]. However, the infiltration enhancement and the recovery of groundwater level due to LID are determined by many factors such as the spatial arrangement of LID and spatial differences of hydrogeological characteristics. Moreover, the increase of groundwater recharge under LID has great variations in different seasons. In addition, rising groundwater levels due to LID possibly increase the risk of groundwater flooding and the risk shows spatial and temporal differences in urban areas [3,14,16]. Therefore, it is necessary to understand the spatial and temporal changes of groundwater under LID at the catchment scale.
This paper focuses on analyzing the spatial and seasonal variations of groundwater head and hydrological processes in an urbanized area, Shenzhen, China, where LID practices are designed to reduce urban runoff. A coupled Storm Water Management Model (SWMM) and Finite Element Subsurface FLOW system (FEFLOW) model was used to simulate surface and subsurface water flow in the area. Based on the model simulation, the objectives of this study are to investigate: (1) the changes in hydrological balance in the area, (2) the spatial variations of groundwater head and its influencing factors, and (3) the seasonal variations of groundwater systems under designed LID scenarios. This study provides an improved understanding on effectiveness of LID infiltration for groundwater supplement in the urbanized areas.

Study Area
The study area is located in the southwestern coast of Shenzhen, China ( Figure 1), with an area of 73.1 km 2 . It is adjacent to the Pearl River Estuary to the east and the Shenzhen Bay (also called Deep Bay) to the west. The time-average sea level in the Pearl River Estuary and Shenzhen Bay is about 0 m based on the Huanghai Vertical Datum. The Shuangjie River flows along the northwestern boundary of the area, with a length of 3.2 km and average water levels ranging from 7 m in the north to 0 m in the south. The Dasha River flows along the eastern boundary of the study area, with a length of 4.5 km and average water levels ranging 0-6 m. The area has a subtropical monsoon climate with very wet summer and dry winter. The mean annual rainfall is 1944 mm, and up to 85% of annual precipitation occurs during wet seasons from April to September.

Study Area
The study area is located in the southwestern coast of Shenzhen, China ( Figure 1), with an area of 73.1 km 2 . It is adjacent to the Pearl River Estuary to the east and the Shenzhen Bay (also called Deep Bay) to the west. The time-average sea level in the Pearl River Estuary and Shenzhen Bay is about 0 m based on the Huanghai Vertical Datum. The Shuangjie River flows along the northwestern boundary of the area, with a length of 3.2 km and average water levels ranging from 7 m in the north to 0 m in the south. The Dasha River flows along the eastern boundary of the study area, with a length of 4.5 km and average water levels ranging 0-6 m. The area has a subtropical monsoon climate with very wet summer and dry winter. The mean annual rainfall is 1944 mm, and up to 85% of annual precipitation occurs during wet seasons from April to September. The fluvial-marine plain spreads around the catchment with the elevation ranging from 0 to 30 m for most of areas. According to the aquifer lithology, the study area can be divided into four hydrogeological zones (Figure 2a). A fractured aquifer comprising heavily weathered biotite granite exists in Zones 1 and 2. The groundwater table depth in Zone 2 is larger than that in Zone 1 due to the high hill around the Zone 2. The porous Quaternary aquifer is mainly distributed around the river valley in Zone 3. The area has experienced four reclamation activities between 1983 and 2005, so a wide area of fill material is spread around the coastline and forms a new porous aquifer in Zone 4 [18]. The fluvial-marine plain spreads around the catchment with the elevation ranging from 0 to 30 m for most of areas. According to the aquifer lithology, the study area can be divided into four hydrogeological zones (Figure 2a). A fractured aquifer comprising heavily weathered biotite granite exists in Zones 1 and 2. The groundwater table depth in Zone 2 is larger than that in Zone 1 due to the high hill around the Zone 2. The porous Quaternary aquifer is mainly distributed around the river valley in Zone 3. The area has experienced four reclamation activities between 1983 and 2005, so a wide area of fill material is spread around the coastline and forms a new porous aquifer in Zone 4 [18]. There are four soil or rock formations from the ground to the bedrock (Figure 2b). Layer 1 is composed of manmade fill materials, muddy clay, and silts. The thickness of this layer is about 2-7 m, increasing from north to south. Layer 1 is regarded as an aquitard due to the poor permeability. Layer 2 is sand and gravel alluvium interbedded clay layers with an averaged thickness of 10 m. This layer is the main aquifer in the area and is partly confined because of the overlying aquitard. The hydraulic conductivity of this layer ranges from 2.5 to 6.5 × 10 −4 m/s. Layer 3 is made of less permeable silty clay, and Layer 4 is formed by decomposed granite with a hydraulic conductivity of 0.8 to 1.2 × 10 −4 m/s [19]. Groundwater in the area is mainly recharged by rainfall, flows generally from north to south, and eventually discharges to the surrounding sea and rivers. The average annual discharge of groundwater in the Dasha River catchment is about 4-5 million m 3 /year according to the local hydrological survey report [20]. The depth to groundwater table is highly related to the topography in the study area. The groundwater table depth in biotite granite is over 10 m around the mountains and is less than 3 m in the areas of Quaternary sediment and fill material near the rivers and bays.
The study area has undergone rapid urbanization in the last 30 years. Seventy-two percent of the land has been urbanized, with main land uses for residence and industry. Rapid urbanization has introduced a greater impervious land surface, which has substantially changed the hydrological processes and caused more surface runoff and less groundwater recharge in this area. In order to restore the natural hydrological cycle and improve the urban ecological environment, the local government is planning to promote LID, and a number of LID facilities such as bioretention basins, green roofs, and permeable pavements will be implemented in the study area by 2020.

Urban Storm Water Model
The US EPA Storm Water Management Model (SWMM version 5.1) [21] is used to simulate the rainfall-runoff process in the study area. SWMM has been widely used for modelling the quantity and quality of surface runoff in urban areas [22,23]. The SWMM model also provides a LID control module to evaluate the effects of LID on hydrological processes [24][25][26].
The study area is simplified to 29 subcatchments and 31 junctions for simulation according to the topography and drainage networks ( Figure 3). As shown by the data of local land use and cover, the percentage of impervious covers of these subcatchments range from 9.5% to 86%. The Dasha River and the Shuangjie River are set as open channels, and the stormwater sewers are set as conduits in the model. The model simulates the rainfall-runoff process, which includes infiltration, depression There are four soil or rock formations from the ground to the bedrock (Figure 2b). Layer 1 is composed of manmade fill materials, muddy clay, and silts. The thickness of this layer is about 2-7 m, increasing from north to south. Layer 1 is regarded as an aquitard due to the poor permeability. Layer 2 is sand and gravel alluvium interbedded clay layers with an averaged thickness of 10 m. This layer is the main aquifer in the area and is partly confined because of the overlying aquitard. The hydraulic conductivity of this layer ranges from 2.5 to 6.5 × 10 −4 m/s. Layer 3 is made of less permeable silty clay, and Layer 4 is formed by decomposed granite with a hydraulic conductivity of 0.8 to 1.2 × 10 −4 m/s [19]. Groundwater in the area is mainly recharged by rainfall, flows generally from north to south, and eventually discharges to the surrounding sea and rivers. The average annual discharge of groundwater in the Dasha River catchment is about 4-5 million m 3 /year according to the local hydrological survey report [20]. The depth to groundwater table is highly related to the topography in the study area. The groundwater table depth in biotite granite is over 10 m around the mountains and is less than 3 m in the areas of Quaternary sediment and fill material near the rivers and bays.
The study area has undergone rapid urbanization in the last 30 years. Seventy-two percent of the land has been urbanized, with main land uses for residence and industry. Rapid urbanization has introduced a greater impervious land surface, which has substantially changed the hydrological processes and caused more surface runoff and less groundwater recharge in this area. In order to restore the natural hydrological cycle and improve the urban ecological environment, the local government is planning to promote LID, and a number of LID facilities such as bioretention basins, green roofs, and permeable pavements will be implemented in the study area by 2020.

Urban Storm Water Model
The US EPA Storm Water Management Model (SWMM version 5.1) [21] is used to simulate the rainfall-runoff process in the study area. SWMM has been widely used for modelling the quantity and quality of surface runoff in urban areas [22,23]. The SWMM model also provides a LID control module to evaluate the effects of LID on hydrological processes [24][25][26].
The study area is simplified to 29 subcatchments and 31 junctions for simulation according to the topography and drainage networks ( Figure 3). As shown by the data of local land use and cover, the percentage of impervious covers of these subcatchments range from 9.5% to 86%. The Dasha River and the Shuangjie River are set as open channels, and the stormwater sewers are set as conduits in the model. The model simulates the rainfall-runoff process, which includes infiltration, depression storage, evaporation and surface runoff. The surface infiltration is estimated using the Green-Ampt equation and the infiltration rates range from 8 × 10 −6 m/s to 1 × 10 −5 m/s in the study area. The evaporation losses are estimated by the monthly-average potential evaporation rates which range from 2.7 mm/day to 3.9 mm/day during the dry months, and from 4.8 mm/day to 6.3 mm/day during the wet months, respectively. The other model parameters are set by referring to the existing studies and reports in this area [27,28]. The model can be used to evaluate the effects of the LID designs on the surface runoff and the infiltration. storage, evaporation and surface runoff. The surface infiltration is estimated using the Green-Ampt equation and the infiltration rates range from 8 × 10 −6 m/s to 1 × 10 −5 m/s in the study area. The evaporation losses are estimated by the monthly-average potential evaporation rates which range from 2.7 mm/day to 3.9 mm/day during the dry months, and from 4.8 mm/day to 6.3 mm/day during the wet months, respectively. The other model parameters are set by referring to the existing studies and reports in this area [27,28]. The model can be used to evaluate the effects of the LID designs on the surface runoff and the infiltration.

Distributed Groundwater Model
The Finite Element Subsurface FLOW system (FEFLOW version 6.2) [29] is used to simulate the groundwater head in the area. FEFLOW is a popular distributed groundwater model for the simulation of subsurface flow. According to the lithology data from 158 boreholes, the soil and rock formations in the study area are simplified to four layers and 95,516 elements ( Figure 4). The average area of these elements is about 3000 m 2 , and the elements around the monitoring wells or the model boundary were subdivided into a number of smaller elements, which have an area ranging from 10-120 m 2 . The second layer is the main porous aquifer with a composition of sand and gravel. Since the northern part of study area is the uplifted granite bedrock, the northern boundary (C-D) is set as a no-flow boundary. The Shuangjie River and the Dasha River insect with the porous aquifer and exchange with groundwater. Hence, the western boundary (A-D) and the eastern boundary (B-C) of model are both set to the head-dependent boundary according to the river stage. The exchange flux between the groundwater and river is calculated as the product of the transfer rate and water head difference in the FEFLOW [29]. The transfer rate is defined as the hydraulic conductivity of the river bed divided by its thickness. The inflow and outflow transfer rate at the river boundary is taken as 0.05 1/day and 0.08 1/day, respectively [30]. Groundwater discharges into Pearl River Estuary and Shenzhen Bay in the southern part of study area, so the southern boundary of the model (A-B) is set to the Dirichlet head boundary with a mean sea level of 0 m. The bottom of model is defined as noflow boundary due to the underlying impervious fresh bedrock. The hydraulic conductivity, specific yield, and specific storage of each model layer (Table 1) are acquired from related reports and research [18,19,30].
The parameter "inflow on top" in the FEFLOW model is applied to describe the net groundwater recharge from the top layer and is usually estimated by a specific percentage of receiving precipitation [19,30]. However, this method is unable to explicitly consider the effects of LID on the groundwater recharge. Assuming that the effects of periodic change of soil water content on the

Distributed Groundwater Model
The Finite Element Subsurface FLOW system (FEFLOW version 6.2) [29] is used to simulate the groundwater head in the area. FEFLOW is a popular distributed groundwater model for the simulation of subsurface flow. According to the lithology data from 158 boreholes, the soil and rock formations in the study area are simplified to four layers and 95,516 elements ( Figure 4). The average area of these elements is about 3000 m 2 , and the elements around the monitoring wells or the model boundary were subdivided into a number of smaller elements, which have an area ranging from 10-120 m 2 . The second layer is the main porous aquifer with a composition of sand and gravel. Since the northern part of study area is the uplifted granite bedrock, the northern boundary (C-D) is set as a no-flow boundary. The Shuangjie River and the Dasha River insect with the porous aquifer and exchange with groundwater. Hence, the western boundary (A-D) and the eastern boundary (B-C) of model are both set to the head-dependent boundary according to the river stage. The exchange flux between the groundwater and river is calculated as the product of the transfer rate and water head difference in the FEFLOW [29]. The transfer rate is defined as the hydraulic conductivity of the river bed divided by its thickness. The inflow and outflow transfer rate at the river boundary is taken as 0.05 1/day and 0.08 1/day, respectively [30]. Groundwater discharges into Pearl River Estuary and Shenzhen Bay in the southern part of study area, so the southern boundary of the model (A-B) is set to the Dirichlet head boundary with a mean sea level of 0 m. The bottom of model is defined as no-flow boundary due to the underlying impervious fresh bedrock. The hydraulic conductivity, specific yield, and specific storage of each model layer (Table 1) are acquired from related reports and research [18,19,30].
The net groundwater recharge and the ET from the soil, plants and groundwater are equal to fRI × Infiltration and (1 − fRI) × Infiltration, respectively. The amount of rainfall that eventually becomes recharge is mainly influenced by the lithology property and groundwater table depth [31]. In order to consider the spatial heterogeneity of hydrogeological conditions, the study area was divided into several geological zones with different aquifer lithology and groundwater table depth. The parameters in different zones are calibrated by using the measured groundwater head.

Model Coupling
A coupled SWMM and FEFLOW model is used to evaluate the effects of LID designs on groundwater. Firstly, the infiltration of each subcatchment with LID designs is calculated by the SWMM model. Secondly, the infiltration converted to net groundwater recharge can be calculated according to the parameter of fRI, and then fed into the FEFLOW model. Since the geological zones intersect with different subcatchments, the infiltration of each subcatchment is mapped to different geological zones according to the percentage of area. Then, the net groundwater recharge can be assigned to each FEFLOW element. In addition, the time step of the SWMM model simulation is set to minutes or hours to capture the dynamic variations of the surface runoff and infiltration process, while the FEFLOW model is run at monthly step. Thus, the infiltration is aggregated on a monthly basis to address temporal difference.

Verification Methods
There is no specific monitoring of groundwater in the study area. However, some drilling data are available from the local geological survey enterprises ( Figure 1). There are 15 accessible drills in the study area, and the groundwater head data of all these wells in September 2008 are available. Moreover, the groundwater head data of three wells were recorded once per month during the period  The parameter "inflow on top" in the FEFLOW model is applied to describe the net groundwater recharge from the top layer and is usually estimated by a specific percentage of receiving precipitation [19,30]. However, this method is unable to explicitly consider the effects of LID on the groundwater recharge. Assuming that the effects of periodic change of soil water content on the groundwater system can be ignored, the infiltration is equal to the net groundwater recharge plus the water loss due to evapotranspiration (ET) from soil, plants and shallow groundwater. In this study, the ratio of net groundwater recharge to infiltration (f RI ) has been used as the internal parameter in the coupled SWMM-FEFLOW model to connect the surface and groundwater modules. The net groundwater recharge and the ET from the soil, plants and groundwater are equal to f RI × Infiltration and (1 − f RI ) × Infiltration, respectively. The amount of rainfall that eventually becomes recharge is Water 2018, 10, 803 7 of 17 mainly influenced by the lithology property and groundwater table depth [31]. In order to consider the spatial heterogeneity of hydrogeological conditions, the study area was divided into several geological zones with different aquifer lithology and groundwater table depth. The parameters in different zones are calibrated by using the measured groundwater head.

Model Coupling
A coupled SWMM and FEFLOW model is used to evaluate the effects of LID designs on groundwater. Firstly, the infiltration of each subcatchment with LID designs is calculated by the SWMM model. Secondly, the infiltration converted to net groundwater recharge can be calculated according to the parameter of f RI , and then fed into the FEFLOW model. Since the geological zones intersect with different subcatchments, the infiltration of each subcatchment is mapped to different geological zones according to the percentage of area. Then, the net groundwater recharge can be assigned to each FEFLOW element. In addition, the time step of the SWMM model simulation is set to minutes or hours to capture the dynamic variations of the surface runoff and infiltration process, while the FEFLOW model is run at monthly step. Thus, the infiltration is aggregated on a monthly basis to address temporal difference.

Verification Methods
There is no specific monitoring of groundwater in the study area. However, some drilling data are available from the local geological survey enterprises ( Figure 1). There are 15 accessible drills in the study area, and the groundwater head data of all these wells in September 2008 are available. Moreover, the groundwater head data of three wells were recorded once per month during the period from September 2008 to December 2009. The rainfall data were acquired from the local meteorological department (Meteorological Bureau of Shenzhen Municipality).
Firstly, the FEFLOW model was performed in a steady state in order to set the initial condition for transient simulations. Secondly, the coupled model was verified against groundwater head data measured in 2008-2009. The SWMM model was run at a time step of 1 h and used to calculate the infiltration for each subcatchment. Then, the infiltration was converted to monthly net groundwater recharge fed into the FEFLOW model. Taking the groundwater head obtained in the first simulation as the initial condition, the FEFLOW model was run at a monthly time-step to simulate the spatial-seasonal variations of groundwater head. The model parameters (e.g., soil saturated hydraulic conductivity, the f RI of geological zones, hydraulic conductivity, specific yield, and specific storage) were tuned by the trial-and-error method to obtain a reasonable estimation of the measured data.

Verification Results
The calibrated model parameters are listed in the Table 1. The simulated average annual surface runoff is about 58% of precipitation and is consistent with the value estimated by other research in this area [32,33]. Hu and Jiao [30] estimated the submarine groundwater discharge of 1.41 × 10 4 m 3 /day in the western part of Shekou peninsula that covers a similar area with the study area. The simulated groundwater seaward discharge in this study is 1.48 × 10 4 m 3 /day and shows good agreement with local research. Figure 5 shows that the simulated groundwater head of the 15 monitoring wells is basically consistent with the measured data of September 2008. The root mean square error (RMSE) of groundwater head is 0.55 m, and the square correlation coefficient (R 2 ) is 0.87. The accuracy of the model is relatively poor at the lithological boundary of biotite granite and Quaternary deposit (e.g., RMSE of 0.5-1 m) due to the complex geology and the uncertainty of hydraulic conductivity.
The results indicate that the model can reasonably reflect the spatial distribution of groundwater head in the study area.

LID Scenarios
In terms of urban planning of Shenzhen, LID is promoted to reduce the runoff from 5-15% of impervious surface in the study area by 2020. However, there is no information about the potential locations of specific LID facilities in the existing LID planning of the study area [32]. In this study, the LID implementation intensity of a subcatchment is defined as the ratio of the impervious area from which the runoff is treated by LID facilities to the total area of the subcatchment. In addition, the runoff from the impervious area is assumed to be equally divided into three parts and respectively treated by bioretention basins, permeable pavements, and green roofs. Three LID scenarios were proposed according to the local planning and above assumptions:


Base case: The land use remains unchanged and no LID practices will be implemented in the area.  Scenario 1: The runoff from 7.5% of impervious area in each subcatchment will be treated respectively by bioretention basins (2.5%), permeable pavements (2.5%), and green roofs (2.5%). The bioretention is designed to have the capacity to treat the runoff from the impervious area that is 10 times its own size according to the Technical Guidance for Sponge City Construction-  The simulated data of Well A can fit well with the variations of monthly observed data, whereas the simulated groundwater head of Well B has a smaller fluctuation range than the observations. The possible reason is that Well B is located in an industrial agglomeration area and the groundwater head may be disturbed by human activities such as pumping activities or underground drainage. The trends of the simulated groundwater head of the two wells are generally in agreement with the observations.

LID Scenarios
In terms of urban planning of Shenzhen, LID is promoted to reduce the runoff from 5-15% of impervious surface in the study area by 2020. However, there is no information about the potential locations of specific LID facilities in the existing LID planning of the study area [32]. In this study, the LID implementation intensity of a subcatchment is defined as the ratio of the impervious area from which the runoff is treated by LID facilities to the total area of the subcatchment. In addition, the runoff from the impervious area is assumed to be equally divided into three parts and respectively treated by bioretention basins, permeable pavements, and green roofs. Three LID scenarios were proposed according to the local planning and above assumptions:


Base case: The land use remains unchanged and no LID practices will be implemented in the area.  Scenario 1: The runoff from 7.5% of impervious area in each subcatchment will be treated respectively by bioretention basins (2.5%), permeable pavements (2.5%), and green roofs (2.5%). The bioretention is designed to have the capacity to treat the runoff from the impervious area that is 10 times its own size according to the Technical Guidance for Sponge City Construction-

LID Scenarios
In terms of urban planning of Shenzhen, LID is promoted to reduce the runoff from 5-15% of impervious surface in the study area by 2020. However, there is no information about the potential locations of specific LID facilities in the existing LID planning of the study area [32]. In this study, the LID implementation intensity of a subcatchment is defined as the ratio of the impervious area from which the runoff is treated by LID facilities to the total area of the subcatchment. In addition, the runoff from the impervious area is assumed to be equally divided into three parts and respectively treated by bioretention basins, permeable pavements, and green roofs. Three LID scenarios were proposed according to the local planning and above assumptions:

•
Base case: The land use remains unchanged and no LID practices will be implemented in the area. • Scenario 1: The runoff from 7.5% of impervious area in each subcatchment will be treated respectively by bioretention basins (2.5%), permeable pavements (2.5%), and green roofs (2.5%). The bioretention is designed to have the capacity to treat the runoff from the impervious area that Water 2018, 10, 803 9 of 17 is 10 times its own size according to the Technical Guidance for Sponge City Construction-Low Impact Development Rainwater System Construction [34]. For example, Figure 7 shows the LID implementation intensity of each subcatchment under Scenario 1. • Scenario 2: The runoff from 15% of impervious area in each subcatchment will be treated respectively by bioretention basins (5%), permeable pavements (5%), and green roofs (5%).
Water 2018, 10, x FOR PEER REVIEW 9 of 17 Low Impact Development Rainwater System Construction [34]. For example, Figure 7 shows the LID implementation intensity of each subcatchment under Scenario 1.  Scenario 2: The runoff from 15% of impervious area in each subcatchment will be treated respectively by bioretention basins (5%), permeable pavements (5%), and green roofs (5%). The settings of different LID facilities in the SWMM model refer to the values reported by Qin et al. [27] and are illustrated in Table 2. The rainfall data of a normal year (2013) with an approximately 50% probability of exceedance was used to conduct the long-term continuous simulation. The corresponding annual rainfall amount is 1997 mm, and 84% of rainfall occurs during the wet seasons (Meteorological Bureau of Shenzhen Municipality). The response of the coastal aquifer system to external interference is a slow process and takes a decade for the system to reach a new equilibrium [30]. In order to investigate the long-term hydrological effects of LID on groundwater, the coupled model was run until the difference of groundwater head between two consecutive years was less than 0.1% under each scenario. In the base case, the model was run continually until the groundwater system reached the equilibrium. In Scenarios 1 and 2, the groundwater head in the base case was set as the initial condition, and then the model was continued to run under the LID scenarios until the groundwater system reached a new equilibrium. The spatial and seasonal variation of groundwater head and the relevant hydrological processes (e.g., surface runoff, infiltration and groundwater discharge) obtained from the simulation were used to investigate the effects of LID infiltration on groundwater in this area. The settings of different LID facilities in the SWMM model refer to the values reported by Qin et al. [27] and are illustrated in Table 2. The rainfall data of a normal year (2013) with an approximately 50% probability of exceedance was used to conduct the long-term continuous simulation. The corresponding annual rainfall amount is 1997 mm, and 84% of rainfall occurs during the wet seasons (Meteorological Bureau of Shenzhen Municipality). The response of the coastal aquifer system to external interference is a slow process and takes a decade for the system to reach a new equilibrium [30]. In order to investigate the long-term hydrological effects of LID on groundwater, the coupled model was run until the difference of groundwater head between two consecutive years was less than 0.1% under each scenario. In the base case, the model was run continually until the groundwater system reached the equilibrium. In Scenarios 1 and 2, the groundwater head in the base case was set as the initial condition, and then the model was continued to run under the LID scenarios until the groundwater system reached a new equilibrium. The spatial and seasonal variation of groundwater head and the relevant hydrological processes (e.g., surface runoff, infiltration and groundwater discharge) obtained from the simulation were used to investigate the effects of LID infiltration on groundwater in this area.

The Effect of LID on Hydrological Balance
When precipitation falls onto ground, some rainfall retains by surface depressions and evaporates back to atmosphere while some infiltrates through land surface. The rest becomes surface runoff. As for infiltration, one part is stored in the aquifer, one part becomes groundwater discharge, and the other part returns to the atmosphere by evapotranspiration (ET) from plants, soil water, and shallow groundwater. According to the annual average statistics, precipitation is equal to surface runoff plus groundwater discharge plus ET. However, storm water management practices such as LID may change the hydrological balance and disturb the original equilibrium status [8,15]. The simulation results show that the groundwater system requires around 10 years to reach a new equilibrium under LID scenarios. Table 3 shows the hydrological balance in the area when the groundwater system reaches the equilibrium under the three scenarios. The annual precipitation of the selected year is 1997 mm. Under base case, the surface runoff, groundwater discharge, and evapotranspiration values are 1207, 74, and 716 mm, respectively. The infiltration is 421 mm and the spatial average groundwater head is 3.0 m. Since the groundwater storage remains stable on an annual basis in the equilibrium state, the annual net groundwater recharge is equal to groundwater discharge. Compared with base case, the surface runoff decreases by 5% and 9% for Scenarios 1 and 2, respectively. As the infiltration increases, the groundwater discharge increases by 46% and 85% for Scenarios 1 and 2, respectively. Furthermore, the spatial average groundwater head relative to sea level pressure increases by 0.9 m and 1.7 m for Scenarios 1 and 2, respectively. In addition, the evapotranspiration also increases by 4-7% under LID scenarios.

The Effect of LID on Spatial Variation of Groundwater Processes
The 29 subcatchments can be divided into recharge zones (R1-R12) and discharge zones (D1-D17) according to the connectivity between the subcatchments and the receiving water. Under base case, the annual average groundwater head gradually decreases from 10 m at the upstream boundary to 0 m at the coastline (Figure 8). Compared with base case, the rise in groundwater head ranges from 0.2 m to 2.2 m, and 0.3 m to 3.9 m, respectively, under Scenarios 1 and 2. The groundwater head changes show similar spatial variability between Scenarios 1 and 2, which decreases gradually from recharge zones to discharge zones along the flow direction. However, there are several exceptions. For example, in recharge zones, R3 has less of an increase in groundwater head than nearby R4; in discharge zones, D9 has a greater increase in groundwater head than D10, and D14 has less of an increase as compared to D13.
Water 2018, 10, x FOR PEER REVIEW 11 of 17 by 5% and 9% for Scenarios 1 and 2, respectively. As the infiltration increases, the groundwater discharge increases by 46% and 85% for Scenarios 1 and 2, respectively. Furthermore, the spatial average groundwater head relative to sea level pressure increases by 0.9 m and 1.7 m for Scenarios 1 and 2, respectively. In addition, the evapotranspiration also increases by 4-7% under LID scenarios.

The Effect of LID on Spatial Variation of Groundwater Processes
The 29 subcatchments can be divided into recharge zones (R1-R12) and discharge zones (D1-D17) according to the connectivity between the subcatchments and the receiving water. Under base case, the annual average groundwater head gradually decreases from 10 m at the upstream boundary to 0 m at the coastline (Figure 8). Compared with base case, the rise in groundwater head ranges from 0.2 m to 2.2 m, and 0.3 m to 3.9 m, respectively, under Scenarios 1 and 2. The groundwater head changes show similar spatial variability between Scenarios 1 and 2, which decreases gradually from recharge zones to discharge zones along the flow direction. However, there are several exceptions. For example, in recharge zones, R3 has less of an increase in groundwater head than nearby R4; in discharge zones, D9 has a greater increase in groundwater head than D10, and D14 has less of an increase as compared to D13. The groundwater head difference relative to base case is mainly determined by the difference between groundwater recharge and discharge. In order to investigate the factors that influence the spatial variation of groundwater head changes, the cumulative increase in net groundwater recharge and discharge in each subcatchment during the 10-year period was calculated for Scenario 1 ( Figure  9). The results indicate that D8 and D15 have more of an increase in net groundwater recharge, while D9, D14 and D13 have less of an increase in net recharge (Figure 9a). Figure 10a shows that the increase in net groundwater recharge has positive correlation with the proportion of land for LID implementation with a correlative coefficient of 0.99 (sig = 0.01). The results implicate that the proportion of land for LID implementation is a key factor that influences the net groundwater recharge. The cumulative increase in groundwater discharge is usually affected by many factors such as increased aquifer recharge, the overlying soil property, and the distance from the receiving water. Figure 9b shows that D14, D15, and D10 have a greater increase in groundwater discharge while D9 and R2 have less of an increase in discharge. The increase in groundwater discharge has a positive The groundwater head difference relative to base case is mainly determined by the difference between groundwater recharge and discharge. In order to investigate the factors that influence the spatial variation of groundwater head changes, the cumulative increase in net groundwater recharge and discharge in each subcatchment during the 10-year period was calculated for Scenario 1 (Figure 9). The results indicate that D8 and D15 have more of an increase in net groundwater recharge, while D9, D14 and D13 have less of an increase in net recharge (Figure 9a). Figure 10a shows that the increase in net groundwater recharge has positive correlation with the proportion of land for LID implementation with a correlative coefficient of 0.99 (sig = 0.01). The results implicate that the proportion of land for LID implementation is a key factor that influences the net groundwater recharge. The cumulative increase in groundwater discharge is usually affected by many factors such as increased aquifer recharge, the overlying soil property, and the distance from the receiving water. Figure 9b shows that D14, D15, and D10 have a greater increase in groundwater discharge while D9 and R2 have less of an increase in discharge. The increase in groundwater discharge has a positive correlation with the hydraulic conductivity of the main underlying aquifer, with a correlative coefficient of 0.83 (sig = 0.01) (Figure 10b). The results imply that the hydraulic conductivity of underlying aquifer is the main factor that influences the groundwater discharge.
Water 2018, 10, x FOR PEER REVIEW 12 of 17 correlation with the hydraulic conductivity of the main underlying aquifer, with a correlative coefficient of 0.83 (sig = 0.01) (Figure 10b). The results imply that the hydraulic conductivity of underlying aquifer is the main factor that influences the groundwater discharge.  The above analysis indicates that the LID scenarios cause relatively less of an increase in groundwater head in discharge zones than in recharge zones. This is because the groundwater increased by infiltration due to LID in the discharge zone may directly discharge into the nearby rivers or offshore areas (Figure 9). However, the groundwater head change in each subcatchment is not only determined by the location in the catchment, but also the degree of LID implementation and the aquifer drainage capacity. For example, compared with the nearby R2 and R4 in recharge zones, R3 has a lower percentage of land for LID, larger aquifer hydraulic conductivity, and thus less of an increase in groundwater head. In discharge zones, compared with the nearby D8 and D10, D9 has less gain of net groundwater recharge; however, D9 has much less increase in groundwater discharge due to the poor hydraulic conductivity of underlying granite. Therefore, the rise in groundwater head of D9 is greater than in the nearby subcatchments. As for D14, it has significantly higher hydraulic conductivity in underlying sediments than D13 and D15, leading to less of an increase in groundwater head. Therefore, it can be concluded that the spatial variation of groundwater head changes is mainly controlled by drainage capacity of underlying aquifer, and also affected by many  correlation with the hydraulic conductivity of the main underlying aquifer, with a correlative coefficient of 0.83 (sig = 0.01) (Figure 10b). The results imply that the hydraulic conductivity of underlying aquifer is the main factor that influences the groundwater discharge.  The above analysis indicates that the LID scenarios cause relatively less of an increase in groundwater head in discharge zones than in recharge zones. This is because the groundwater increased by infiltration due to LID in the discharge zone may directly discharge into the nearby rivers or offshore areas ( Figure 9). However, the groundwater head change in each subcatchment is not only determined by the location in the catchment, but also the degree of LID implementation and the aquifer drainage capacity. For example, compared with the nearby R2 and R4 in recharge zones, R3 has a lower percentage of land for LID, larger aquifer hydraulic conductivity, and thus less of an increase in groundwater head. In discharge zones, compared with the nearby D8 and D10, D9 has less gain of net groundwater recharge; however, D9 has much less increase in groundwater discharge due to the poor hydraulic conductivity of underlying granite. Therefore, the rise in groundwater head of D9 is greater than in the nearby subcatchments. As for D14, it has significantly higher hydraulic conductivity in underlying sediments than D13 and D15, leading to less of an increase in groundwater head. Therefore, it can be concluded that the spatial variation of groundwater head changes is mainly controlled by drainage capacity of underlying aquifer, and also affected by many The above analysis indicates that the LID scenarios cause relatively less of an increase in groundwater head in discharge zones than in recharge zones. This is because the groundwater increased by infiltration due to LID in the discharge zone may directly discharge into the nearby rivers or offshore areas ( Figure 9). However, the groundwater head change in each subcatchment is not only determined by the location in the catchment, but also the degree of LID implementation and the aquifer drainage capacity. For example, compared with the nearby R2 and R4 in recharge zones, R3 has a lower percentage of land for LID, larger aquifer hydraulic conductivity, and thus less of an increase in groundwater head. In discharge zones, compared with the nearby D8 and D10, D9 has less gain of net groundwater recharge; however, D9 has much less increase in groundwater discharge due to the poor hydraulic conductivity of underlying granite. Therefore, the rise in groundwater head of D9 is greater than in the nearby subcatchments. As for D14, it has significantly higher hydraulic conductivity in underlying sediments than D13 and D15, leading to less of an increase in groundwater head. Therefore, it can be concluded that the spatial variation of groundwater head changes is mainly controlled by drainage capacity of underlying aquifer, and also affected by many complicated factors such as the degree of implementation of LID and the connectivity between groundwater and receiving water.

The Effect of LID on Seasonal Variation of Groundwater Processes
The seasonal variations of groundwater system in the whole area are illustrated in Figure 11 when it reaches equilibrium under the above three scenarios. In base case, taking Well A as an example, the groundwater head increases during the wet seasons (April to September) and decreases during the dry seasons (October to March) (Figure 11a). This is because the aquifer has a positive change in storage during the wet seasons and a negative change during the dry seasons (Figure 11b). Furthermore, the seasonal variation of groundwater head is determined by the state of aquifer water exchange. For example, the water absorption rate of aquifer in April is higher than that in other months, leading to the fastest groundwater head rise in the year, while the groundwater head has no significant decrease in December because the release rate of aquifer is close to zero. Figure 11b also indicates that aquifer can regulate hydrological processes like a sponge. It absorbs and stores the abundant infiltration of stormwater during wet seasons and releases the stored water in the dry months to maintain groundwater discharge. The state of aquifer water exchange is mainly affected by the net groundwater recharge. Generally, the aquifer absorbs more water when there is adequate recharge and releases more water in dry months with little rainfall (Figure 11c). In addition, the rate of aquifer water exchange is also affected by the antecedent dry period. For example, although the net groundwater recharge rates in April and August are nearly the same (0.8 million m 3 /month), the water absorption rate of aquifer in April is 0.11 million m 3 /month more than that in August. This can be explained by the greater storage capacity of aquifer in April due to the continuous water release in former dry seasons. Due to the "buffer" function of the aquifer, the seasonal change in groundwater discharge appears a relatively small fluctuation relative to the rainfall variations (Figure 11d), ranging from 0.32 million m 3 /month during the dry months to 0.55 million m 3 /month during the wet seasons. complicated factors such as the degree of implementation of LID and the connectivity between groundwater and receiving water.

The Effect of LID on Seasonal Variation of Groundwater Processes
The seasonal variations of groundwater system in the whole area are illustrated in Figure 11 when it reaches equilibrium under the above three scenarios. In base case, taking Well A as an example, the groundwater head increases during the wet seasons (April to September) and decreases during the dry seasons (October to March) (Figure 11a). This is because the aquifer has a positive change in storage during the wet seasons and a negative change during the dry seasons (Figure 11b). Furthermore, the seasonal variation of groundwater head is determined by the state of aquifer water exchange. For example, the water absorption rate of aquifer in April is higher than that in other months, leading to the fastest groundwater head rise in the year, while the groundwater head has no significant decrease in December because the release rate of aquifer is close to zero. Figure 11b also indicates that aquifer can regulate hydrological processes like a sponge. It absorbs and stores the abundant infiltration of stormwater during wet seasons and releases the stored water in the dry months to maintain groundwater discharge. The state of aquifer water exchange is mainly affected by the net groundwater recharge. Generally, the aquifer absorbs more water when there is adequate recharge and releases more water in dry months with little rainfall (Figure 11c). In addition, the rate of aquifer water exchange is also affected by the antecedent dry period. For example, although the net groundwater recharge rates in April and August are nearly the same (0.8 million m 3 /month), the water absorption rate of aquifer in April is 0.11 million m 3 /month more than that in August. This can be explained by the greater storage capacity of aquifer in April due to the continuous water release in former dry seasons. Due to the "buffer" function of the aquifer, the seasonal change in groundwater discharge appears a relatively small fluctuation relative to the rainfall variations ( Figure  11d), ranging from 0.32 million m 3 /month during the dry months to 0.55 million m 3 /month during the wet seasons. The seasonal variations of groundwater hydrological processes under LID scenarios are similar to the base case. However, comparing with base case, the groundwater head in each month under LID scenarios shows a significant increase due to the 10-year accumulative recharge (Figure 11a). The aquifer stores 0.56-1.02 million m 3 more water during the wet seasons under LID scenarios. The seasonal variations of groundwater hydrological processes under LID scenarios are similar to the base case. However, comparing with base case, the groundwater head in each month under LID scenarios shows a significant increase due to the 10-year accumulative recharge (Figure 11a). The aquifer stores 0.56-1.02 million m 3 more water during the wet seasons under LID scenarios. Correspondingly, the aquifer releases more water during the dry seasons under Scenarios 1 and 2 ( Figure 11b). Compared with the base case, the net groundwater recharge under Scenarios 1 and 2 during wet seasons increases by 49% and 89%, respectively, while no significant growth of groundwater recharge occurs in dry months (Figure 11c). However, the groundwater discharge still has a great growth in dry seasons due to the gain in aquifer water release under LID scenarios (Figure 11d). For the base case, the groundwater discharge during dry months is 0.39 million m 3 /month, and the discharge increases by 0.17 million m 3 /month and 0.31 million m 3 /month for Scenarios 1 and 2, respectively. The river discharge from groundwater increases from 1.2 million m 3 /year in base case to 2.5 million m 3 /year in Scenario 2. The increasing groundwater discharge to receiving water during dry seasons will help increase river base flow and improve the ecological status of rivers and bays.

The Impacts of Increasing Groundwater
LID infiltration brings a lot of benefits to urban ecological environment. Increasing groundwater head in urban coastal aquifers can increase the potential of local water supply. It may mitigate the inland movement of saltwater interface and reduce the risks of seawater intrusion [6]. However, it is also important to balance the relationship between natural groundwater restoration and urban underground infrastructure safety [13]. Increasing groundwater head in aquifers due to LID infiltration may trigger groundwater flooding. The flooding may induce the emergence of groundwater above low-lying surface terrain or subsurface infrastructure and cause the instability of infrastructure foundation, damage of building assets and the threats to drinking water safety [3,35,36]. Understanding the seasonal and spatial variations of groundwater head under LID scenarios can help to assess the risk of groundwater flooding. The recharge zones overlying biotite granite (such as R2, R11, R12) have low groundwater flood risk. This is because the increase in groundwater head (2-3 m) under LID scenarios is much lower than the groundwater table depth (>10 m) in these zones. However, subcatchments on top of Quaternary sediments in the recharge zones (such as R3, R4, R9) are located in a dense residential and commercial area and have an average groundwater table depth of about 2-3 m. Thus, an over 1-m increase in groundwater head will have negative impacts on the building base. In the discharge zones, the groundwater head only has slight increase due to the large drainage capacity of underlying aquifer, leading to low risk of groundwater flood under LID scenarios. Moreover, the groundwater head increases significantly after LID infiltration in the wet months, so there will be a higher risk of groundwater flooding. Although the groundwater flooding is investigated qualitatively due to model and data limitations, the findings from this analysis can help address the stormwater managers' concerns and provide the basis for further assessing the impacts of LID infiltration in urban areas.

Model Limitations
It should be noted that the study is limited to the model assumption. In this study, it is assumed that the LID implementation intensity is proportional to existing impervious area in each subcatchment. In fact, LID implementation is affected or limited by many factors, e.g., land use type, soil properties, and groundwater table depth. The effects of these factors on the spatial differentiation of LID implementation should be considered in LID scenario design in further study [32]. It is also assumed that the infiltration is not affected by the fluctuation in groundwater level beneath LID facilities. However, the infiltration capacity of LID facilities possibly decreases if LID is applied in areas with rising groundwater mounding [11]. Therefore, the model is likely to overestimate the groundwater recharge in areas with shallow groundwater under LID scenarios. Moreover, the value of f RI is likely to slightly increase with the increase of soil moisture in the unsaturated layer, so this model possibly underestimates net groundwater recharge during the wet seasons and overestimates net groundwater recharge during the dry seasons. In order to better simulate the interaction between the surface and subsurface hydrological processes, it is necessary to develop a two-way coupling model. The outputs (e.g., groundwater head, soil moisture) of the groundwater model at each time step need to be transferred into the surface water model, and to be considered as the basis to calculate the outputs (e.g., infiltration from LID and non-LID areas) of the surface water model in the next time step. In addition, it is necessary to collect more temporal and spatial variation data on the groundwater head for model calibration. Most importantly, the groundwater monitoring should be conducted in the catchment with widespread LID facilities in order to investigate the hydrological impacts of LID at a catchment or urban scale. As for spatial distribution, the monitoring points should cover the groundwater recharge and discharge zones as well as the different hydrogeological conditions and land use types. Furthermore, more than 10 years of monitoring data are required because of the slow response of groundwater system in the study area.

Conclusions
In this study, a coupled Storm Water Management Model (SWMM) and Finite Element Subsurface FLOW system (FEFLOW) model was used to evaluate the effects of LID scenarios on the spatial and seasonal variations of groundwater head and hydrological processes in an urbanized area in Shenzhen, China. The results obtained are summarized below: LID changes the hydrological balance in urbanized areas significantly. In this study area, the annual surface runoff is reduced by 5% and 9%, respectively, when the runoff from 7.5% and 15% of impervious area is treated by LID facilities. As a result of increased infiltration, the spatial average groundwater head relative to sea level pressure in the equilibrium state is predicted to rise by 0.9 m and 1.7 m, respectively.
In this area, the increase in groundwater head decreases generally from the recharge zones to the discharge zones near the receiving streams and coastlines under LID scenarios. The groundwater head changes are determined not only by the location in the catchment, but also by the hydraulic conductivity of underlying aquifers and LID infiltration intensity. LID significantly enhances groundwater recharge and aquifer storage in the wet seasons. In turn, it increases aquifer release and groundwater discharge in the dry seasons.
LID infiltration brings many benefits to the urban ecological environment. Increasing groundwater head in urban aquifers can increase the river base flow and improve the potential of local water supply. However, LID infiltration also has the potential to increase the risk of inland groundwater flooding in areas with poor aquifer drainage capacity and shallow groundwater depth, especially during the annual wet seasons. Therefore, further catchment-scale modelling and experiments are required to better evaluate the trade-off of LID infiltration.