Introduction

Shenzhen, in Guangdong province, is an important megacity in southern China which has experienced a huge and rapid development. Its area extends along the southern Chinese Sea coastline between the Zhujiang River (Pearl River) Estuary and Daya Bay and shares its southern boundary with Hong Kong (Fig. 1). The rapid development has caused a large reduction in agricultural lands and green areas (Li et al. 2010; Hao et al. 2011). Urbanization increased from 9 to 40% of the city area between 1970 and 2008 (Chen et al. 2014; Hong and Guo 2017). Urban expansion continues and has reached 42% of the study area, as shown by analysis of satellite images (Fig. 1). Construction of reservoirs for drinking and industrial water uses (Wang et al. 2004) has greatly altered the morphology of the river valleys (Qin et al. 2013). The water supply (2000 × 106 m3/year) is mostly imported (~75%), as the remaining contribution (only ~25%) comes from reservoirs and water-wells distributed inside the study area (WRBSM 2016). The groundwater resources of Shenzhen are underutilized and the amount of pumped water is only 5.9 × 106 m3/year (WRBSM 2016). Extensive pollution affects surface waters and the aquifers, with high costs for wastewater treatment systems (Li et al. 2016). Notwithstanding the link between economic development and the environment, the hydrogeological setting of Shenzhen has not been thoroughly investigated.

Fig. 1
figure 1

Satellite image of Shenzhen and Zhujiang River estuary (from CNES-2017). The urban areas were defined by manual analysis of satellite imagery at 1:5.000 scale

Full size image

The earliest available compilation of hydrogeological information was GBGP (Geological Bureau of Guangdong Province, China, Unpublished report on regional hydrogeological survey, scale 1:200.000, 1980) which compiled many borehole logs from throughout Shenzhen along with recovery tests and hydrogeochemical data. A decade later, Yang et al. (1991) described the geology of the Shenzhen area, attributing the current structural setting to strike-slip tectonics with a main element oriented NE–SW. Campbell and Sewell (1997, 1998) reconstructed the different plutonic stages and their relationship with the strike-slip tectonics. Xu et al. (2015) investigated the structural setting from a geophysical point of view. The most comprehensive hydrogeological description of Shenzhen is Yang et al. (2007).

Groundwater is strongly influenced by the relationship between the aquifers and anthropic infrastructure. In urban areas, groundwater sources and pathways are numerous and complex. Buildings, roads and other impervious surfaces impede the natural aquifer recharge, increasing surface runoff, while on the other hand, leakage from water mains and sewers provide unconventional groundwater recharge (Lerner 2002). Many studies discuss the impacts of urbanization on the hydrogeological setting (Held et al. 2006; Passerello et al. 2012; Kruse et al. 2013; Di Salvo et al. 2014). Due to the economic importance of the study area, groundwater is of great interest to local water and environmental management agencies. The key science question is how a megacity such as Shenzhen should manage and protect its groundwater resources as a strategic resource and environmental asset. Because of the lack of literature and its geological and urban complexity, this study aims to be a first step to investigate this important key scientific question through:

  • Constructing the three-dimensional (3D) distribution of the aquifers and defining a conceptual circulation model

  • Evaluating the recharge and permeability of the hydrogeological units

  • Identifying and quantifying the groundwater resources of Shenzhen

Hydrogeological setting and conceptual model

Shenzhen is located on the continental margin of the Cathysian block (Zhang et al. 2013) and experienced the creation of a basin and range structure, affected by the complex geodynamic setting of the Guangdong region and, more generally, of Southeast Asia (Li and Li 2007). In Shenzhen, the stratigraphic basement consists of Lower Paleozoic or older clastic rocks with low-grade metamorphism (Fig. 2). Overlying the basement, the sedimentary succession consists of Upper Paleozoic and Mesozoic clastic sequences characterized by a continental to shallow marine-water environment. The stratigraphic sequence varies from siltstone to sandstone with layers of gravel and breccia (Fig. 2). Carbonate sedimentation in the Lower Carboniferous is represented by a sequence of dolostone and limestone, preserved in the Pingshan and Longgang area (Fig. 2b). The thickness of the succession is large and difficult to characterize in its entirety, due to the heterogeneous sedimentation with several hiatus and transgressive events. Yang et al. (2007) estimate a maximum thickness around 15 km for the whole succession. Late Jurassic–Early Cretaceous granitic intrusions related to the subduction of the paleo-Pacific slab (Fig. 2a,b) are widespread and interrupt the sedimentary succession (Li and Li 2007; Xia et al. 2011). This regional magmatic event is known as the Yanashanian event (Zhang et al. 2013). As a result of the intrusive activity, the sedimentary succession is intersected by contact metamorphism, related to two phases of intrusion that outcrop in the Bao’An and Nanshan areas and along the Mirs Bay, Dapeng peninsula (Fig. 2b). An acid magmatism resulted in a volcanic edifice of rhyolitic lava flows and tuffs located at Wutong Mt. (Fig. 2b). Cretaceous tuffs associated with loose ash are scattered along the Dapeng peninsula as well as in the Pingshan and Longgang basins (Fig. 2). Younger reworked tuff in a fluvial-lacustrine environment also outcrops in the Pingshan basin (Fig. 2b). Yang et al. (2007) date these sequences from the Cretaceous to the Eocene. Quaternary deposits are thickest in the main and secondary river valleys and along the coast. The thickness of these deposits varies as a function of the inherited morphology; between 30 and 60 m in the coastal sector and main river valleys, and 15–40 m in the secondary river valleys (Lu and Sun 1991).

Fig. 2
figure 2

a Geological and hydrogeological setting of the western Guangdong Province, b the Shenzhen city with sections, and c a review of aquifer characteristics. Key to the legend: (1) Quaternary clastic unit from coastal plain and alluvial valleys; (2) Tertiary reworked tuff unit; (3) Upper Mesozoic acid lava flow unit; (4) Jurassic granite unit, more fractured (a) or more compact (b); (5) Upper Paleozoic-Mesozoic sedimentary unit; (6) Carboniferous carbonate unit; (7) Paleozoic sedimentary unit; (8) Lower Paleozoic (?) basement unit; (9) normal, contractional and/or strike slip fault; (10) spring, number indicates the elevation m a.s.l. (a), hot spring numbers indicate the elevation and the temperature (b); (11) main river; (12) main reservoir, number indicates the average water level; (13) main domestic water-well or piezometer, number indicates the elevation m a.s.l.; (4) hydrogeological section; (15) limit of the study area of Fig. 2b (represented in Fig. 2a)

Full size image

The distribution of the Paleozoic and Mesozoic rocks is controlled by the different tectonic stages that also influence the trend of the river valleys. The Lower Paleozoic tectonics is characterized by contractional faults and thrust planes, NE–SW oriented and NW verging, created during the phases of the Caledonian orogenesis (Pang et al. 2016). During the early Mesozoic, the orogen experienced a distension associated with high rate of subsidence, due to roll back of the paleo-Pacific plate. The study area was intersected by strike-slip tectonics (Allen et al. 1997) with the Meixian-Shenzhen fault as important tectonic line (Sun et al. 2007). The latter divides the Bao’An area, characterized by an Upper Paleozoic hiatus, from the Paleozoic sequence of the Longgang-Pingshan area. These faults were active up to the Quaternary and in some cases up to the Middle part of the Upper Pleistocene (Lu and Sun 1991).

For the conceptual model the great thickness of sedimentary units, associated with the strike-slip tectonics, is schematized as a block structure, whereby every block is delimited by high-angle faults or granite intrusions. Cross-sections in Fig. 2b are representative of the regional geological setting. The literature, especially the borehole stratigraphy from GBGP(Geological Bureau of Guangdong Province, China, Unpublished report on regional hydrogeological survey, scale 1:200.000, 1980) and the hydrogeological map from Yang et al. (2007), supports a division into eight main hydrogeological units:

  • The Quaternary clastic unit (legend 1 in Fig. 2) consists of a multi-layered sequence with the hydraulic conductivity strongly related to the predominant size of the clasts. Deposits are heterogeneous; layers consist of coarse to silty sands, with intercalations of sandy gravels along the river valleys or silty clays in coastal areas as the result of marine transgressions and/or temporary lacustrine events. In addition, layers of detritus and talus, derived from the weathered bedrock, have been observed. This hydrogeological unit has a thickness varying between 30 and 60 m along the main valleys and in the coastal sectors.

  • The Tertiary reworked tuff unit consists of loose ash and tuffs (legend 2 in Fig. 2) which belong to the upper portion of the volcanic sequence. This hydrogeological unit also includes reworked deposits with intercalations of coarse gravel observed in outcrop in the Pingshan area.

  • The Upper Mesozoic lava unit (legend 3 in Fig. 2) constitutes Wuton Mt. in the southern part of Shenzhen. An interconnected fracture system characterizes the lava flows. Individual fracture orientation is directly correlated to the movement of the flow. Intercalations of paleo-soils or weathered deposits locally reduce the hydraulic conductivity.

  • The Jurassic granite unit (legend 4 in Fig. 2) is an aquitard. Facture density is very low and RQD (rock quality designation) index records excellent values (90–100%), beyond 30 m from the ground surface. Towards the land surface, a blanket of more fractured and weathered deposits is present. This band typically extends to 30 m below the land surface based on the analysis of the stratigraphy of many water-wells and boreholes. Usually the shallowest 15-m-thick portion of this unit consists of a very permeable coarse quartz sand. This shallowest part represents the main aquifer, where the groundwater circulation is more developed. In the eastern sector, the granite aquitard (legend 4a in Fig. 2) is more fractured (Yang et al. 2007) but with a thinner weathered band, gravitationally mobilized because of the steep morphology. From a hydrogeological point of view, this granite (legend 4a in Fig. 2) is grouped with the Upper Mesozoic lava unit.

  • The Upper Paleozoic–Mesozoic unit (legend 5 in Fig. 2) consists of a sequence of quartz sandstone, siltstone, sandy siltstone and shale. Deposits are widely intersected by tectonic sediments. Sediments are lithified and the primary porosity is very low even in the coarser deposits. Fractures are abundant but the apertures were strongly reduced by dynamic and regional metamorphic processes as proven by the presence of andalusite, chlorite and other metamorphic minerals (Yang et al. 2007).

  • The Carboniferous carbonate unit (legend 6 in Fig. 2) is strongly karstified and it is present at the bottom of the main Quaternary depressions in the Pingshan and Longgang areas. Numerous sinkholes occur where this unit is present. Sporadic outcrops also suggest metamorphic facies of marble, dolomitic marble and crystalline limestone.

  • The Paleozoic unit (legend 7 in Fig. 2) consists of metamorphic sandstone, phyllite and slate with serpentine. Porosity is very low due to the huge lithostatic load that these rocks experienced. Due to the different tectonic cycles, fracture systems are interconnected and more developed than in the Upper Paleozoic–Mesozoic unit. The Basement unit (legend 8 in Fig. 2) is metamorphosed Lower Paleozoic rocks of quarzitic sandstone and muddy shale. Porosity is very low and fracture apertures were strongly reduced by metamorphic processes. From a hydrogeological point of view, this unit is an aquitard and is grouped with the Upper Paleozoic–Mesozoic unit (legend 5 in Fig. 2).

The stratigraphy is summarized as loose Quaternary clastic aquifers superimposed on old Paleozoic and Mesozoic bedrock aquitards, covered by a layer of alluvium and colluvium. To illustrate the resulting hydrogeological conceptual model, three sections have been drawn. Figure 3 represents the shallow groundwater flow paths. Eluvium and colluvium, combined with Quaternary deposits, constitute the main aquifer, recharged in the rural areas by rainfall. Additional recharge occurs across the urban areas, in accordance with Lerner (Lerner 2002); Fig. 3a,b). Groundwater flows toward and discharges into rivers, reservoirs or directly to the ocean as suggested by Wang et al. (2017). The hydrogeology of the Dapeng peninsula (Fig. 3c) also has a thin aquifer overlying a bedrock aquitard, with recharge from rainfall and urban sources.

Fig. 3
figure 3

Hydrogeological conceptual model of the Shenzhen area, described via three hydrogeological sections through a the Zhuijian River Estuary-Guanglang R., b Pingshan-Longgang and c across the Dapeng peninsula. Section locations are represented in the mini-map (d). Sections are not to scale to emphasise the hydrogeological circulation through the shallow eluvium/colluvium soils, the Quaternary clastic sequence and the tuff layers. Recharge is from rainfall and urban sources and groundwater directly discharges towards the reservoirs, rivers or the sea. Granite and the sedimentary successions are aquitards; thus, hydraulic catchments and hydrogeological basins coincide

Full size image

The Lower Carboniferous carbonate unit is an exception to this conceptual model; it is the most permeable unit of the study area but it is only found in the bottom of the Pingshang and Longgang River valleys, buried by the Quaternary aquifers, as shown in Fig. 3b. As a result of granite intrusions there is no continuity of the karst unit from one hydraulic basin to the next. Despite the high hydraulic conductivity and the ongoing karstification, evidenced by sinkholes (Yang et al. 2007), there is little groundwater movement in this hydrogeological unit except for local flows caused by groundwater pumping.

In this conceptual model, groundwater and surface watersheds coincide and thus groundwater circulation can be studied at the catchment scale. Bedrock is essentially an aquitard but the hydraulic conductivity can be locally altered by the numerous faults. Across the fault plane, the cataclastic bands increase the hydraulic conductivity, functioning as a conduit. On the contrary, the fault zones with a mylonite core have significantly lower permeability (Caine et al. 1996). The effects of faults on groundwater flow are small and localized, as observed by Ran et al. 2013 for normal displacements and Cilona et al. 2015 for strike-slip systems. In the Mazhou River basin (Fig. 2), a single thermal spring illustrates upward heat migration from the bedrock to the shallow zone along a fault, affecting the local groundwater path. The thermal occurrence is isolated and favored by the cataclastic band; however, there is no impact of the faults on the bedrock circulation at the regional scale.

Methodology

Construction of the numerical model

A steady-state groundwater model of Shenzhen was developed using the finite difference code MODFLOW (Harbaugh 2005) and the graphical user interface GMS. The subsurface was discretized with five model layers (Fig. 4a), set as convertible (GMS 2017), with each layer containing 110,000 cells, with dimensions of 500 m in the east–west direction and 100 m in the north–south direction. The Upstream Weighting (UPW) flow package and the Newton solver (NWT) were used (Niswonger et al. 2011). The top of the first layer is the land surface, derived from a digital elevation model (DEM) (BIGEMAP database), with a 2.0 m resolution, while the bottom of the first layer is specified as 30 m below land surface in areas where Quaternary, granite, volcanic and sedimentary hydrogeological units are present (Fig. 4b). The bottom of the second layer was set at 60 m below land surface. This layer represents the deeper Quaternary basin of the Longgang area, the coastal area of Bao’An, Nanshan and Futian districts as well as the tuff banks in Dapeng. The bases of the third, fourth and fifth model layers were set at elevations of −150, −300 and − 500 m a.s.l., respectively. The latter layers were used to characterize the distribution of the karst aquifer at depth. The bottom of the model domain is −500 m a.s.l., a depth at which groundwater movement is negligible at a regional scale. Geological structures such as folds and monoclines are represented, due to the primary importance of the karst aquifer. Hydraulic conductivity and recharge rates were specified for each cell based on the hydrogeological unit mapped at the location of the cell.

Fig. 4
figure 4

a Three-dimensional grid of Shenzhen model, displaying the elevation (m a.s.l.); b schematic discretization of the grid covering five layers; c boundary conditions of the Shenzhen area (map and section). Key to the legend: (1) constant head (IBOUND = −1) for the sea level (h = 0 m a.s.l.); (2) river (RIV package); (3) no flow boundary (IBOUND = 0); (4) recharge from rainfall and urban areas via Unsaturated Zone Flow package (UZF); and (5) General Head Boundary Condition package (GHB), in which the conductance has been calculated assuming an average hydraulic conductivity, in accordance with the conceptual model

Full size image

Constant head boundary conditions (IBOUND = −1, Harbaugh 2005) were specified along the coast with a hydraulic head equal to zero. Rivers have been represented via the river boundary package (RIV 1, Harbaugh et al. 2000). River stage has been derived by DEM. Along the river beds, values of Carc (conductance for unit length, GMS 2017) vary between 0.02 and 0.0001 m2/s, as a function of the width of the rivers and assuming the hydraulic conductivity of the river sediments equal to 0.0001 m/s. The unsaturated zone flow package (UZF) was used to represent the upper model surface (Niswonger et al. 2006). The package can represent thick portions of unsaturated soils, as it is typical for rough and steep morphologies. Recharge rate is the input variable of this package and depends on the hydrogeological units that outcrop. The Brooks and Corey (1966) model has been assumed constant to describe the groundwater flow in the unsaturated soils. Epsilon (EPS) has been set equal to 3.5 and the water content of the unsaturated zone (THTS) to 0.33 (GMS 2017).

The northern boundary of the model consists of many sub-catchments that perfectly follow the watershed lines, as for the Guanlang and Pingshan river valleys. Because the groundwater does not flow beyond the watershed boundary, no-flow boundary conditions (IBOUND = 0, Harbaugh 2005) were specified (Fig. 4c). Nevertheless, the Longgang River catchment extends beyond the boundary of the studied area. Here, the limit of the studied area cuts across a granite mountain section with narrow gorges. Because of the particular hydrogeological setting across the boundary, continuous aquifers are not present. In accordance with the constructed hydrogeological conceptual model, no-flow boundary conditions (IBOUND = 0, Harbaugh 2005) have been extended along the north-west boundary (Fig. 4c).

However, the model boundary cuts across the Longgang River Valley and the groundwater level varies with a horizontal gradient. The general head boundary (GHB) package (Harbaugh et al. 2000) has been set to better define the groundwater setting of the area (Fig. 4c). The conductance term of the GHB package has been calculated, analyzing the elevation of the river and its distance from the boundary cells.

The boundary between Hong Kong and the study area follows the Shenzhen River and the Shatoujian stream (Fig. 2). The downstream portion of the Shenzhen River has been included in the River package. The upper stream portion of the Shenzhen River and the Shoutoujian Stream are less important and characterized by a seasonal discharge. This area is characterized by volcanic lavas with a shallow groundwater flow directed towards these two rivers. Hydraulic exchanges between Hong Kong and Shenzhen do not occur and the no-flow boundary condition (IBOUND = 0, Harbaugh 2005) is appropriate.

Model calibration

Homogeneous hydraulic conductivity and recharge were assumed for each hydrogeological unit. Analysis of recovery tests of water wells (Geological Bureau of Guangdong Province, China, Unpublished report on regional hydrogeological survey, scale 1:200.000, 1980) showed a wide range of hydraulic conductivities. For each hydrogeological unit, a starting value was set at the average of the recovery test results. Recharge in undeveloped areas was based in the first step of calibration on average recharge rates from Yang et al. (2007), with a zero value for urban recharge. Recharge and hydraulic conductivity for each block were calibrated using a trial-and-error procedure.

The manual calibration investigates how changes in model parameters vary the fit between simulated values and targets. The calibration targets were 158 water-wells and piezometers. The data were collected via a classic water level meter in March 2017, at the end of the dry season, and in July 2017, during the wet season. The two datasets typically differ by 3 m in the rural areas and less than 1 m in the flat urban areas. Calibration started by varying the hydraulic conductivity of the Quaternary aquifers that outcrop in the coastal plain and in the alluvial aquifers. Next, the conductivity of weathered soils and then the bedrock units were calibrated. Initial results were satisfying from a qualitative point of view but simulated hydraulic heads differed from the observed by up to 30 m, indicating that recharge needed to be included as a model calibration parameter.

For the rural areas, recharge values were constrained to the literature range of Yang et al. (2007). In urban areas, simulated hydraulic heads were lower than the observed. Therefore unconventional water recharge was considered (i.e. leakage from water mains, sewer pipes, septic tanks and storm drains (Lerner 1990, 2002). For urban areas, a best fit was found for a recharge of 250 mm/year. The final calibrated hydraulic conductivity and recharge datasets are shown in Table 1. All the values are in accordance with the conceptual model and the literature ranges. The calibration plot shows an excellent distribution of the simulated versus the observed hydraulic heads (Fig. 5a), with a residual value ranging ±6 m (Fig. 5b).

Table 1 Hydraulic conductivity and recharge datasets used for the calibration of the numerical model
Full size table
Fig. 5
figure 5

Calibration plots showing a simulated vs. observed hydraulic head and b the residual vs. observed hydraulic head

Full size image

However, calibration of a ground model does not usually yield a unique solution. A reduction in hydraulic conductivity will raise the simulated heads and can be counterbalanced by a reduction in recharge. Therefore, a sensitivity analysis was performed; the entire hydraulic conductivity and recharge datasets were varied by ±10%. An increase in the error was observed in both cases, without varying the groundwater flow direction. This informal sensitivity analysis gives some confidence that the numerical model is satisfactory.

Results and discussion

Numerical results

The model calibration results are excellent considering the multi-catchment scale and a conceptual model with homogeneous parameter values for each hydrogeological unit. In the Zhuijiang River Estuary sector, the flat morphology, the nearness to the coastline and the high hydraulic conductivity of the Quaternary deposits depresses the hydraulic heads to less than 15 m a.s.l. (Fig. 6). The hydraulic gradient steepens to 5% in the piedmont transition between the coastal plain and the less permeable sedimentary and intrusive rocks. Simulated hydraulic heads match the elevation of some springs in the Bao’An area (Fig. 2b). The highest hydraulic heads are observed under Yangtai Mt. (maximum head of 150 m a.s.l.) and in the Yangtian Mts. area (maximum head 192 m a.s.l.).

Fig. 6
figure 6

Simulated hydraulic head distribution of Shenzhen, the result from the top layer. Key to the legend: (1) groundwater flow direction, (2) main river considered in the simulation (RIV package), (3) Applied General Head Boundary (GHB package); (4) single/clustered measured water levels

Full size image

From this area of maximum hydraulic heads, the simulated groundwater fluxes are toward the NNW into the Maozhou River basin and NNE to the Guanlang River (Fig. 6). Shenzhen Bay receives groundwater from the granite intrusions of the Yangtai Mt. and Futian area, via fluxes directed from North to South. In the Shenzhen River basin (Fig. 6), the flatter morphology and the clastic morphology contribute to flatter gradients, with a variation of 50 m in 15 km (average hydraulic gradient of 0.3%). In the southern part of the catchment, groundwater directly discharges to the Mirs Bay. The hydraulic head variation is related to the topography, with a gradient varying from 3 to 5%. Fluxes are directly perpendicular to the coastline and are strongly influenced by the fluvial morphology. The groundwater of the Dapeng peninsula discharges into the Daya Bay in a similar way. The irregular coastline favors the discharge of the groundwater directly to the sea, via hydraulic gradients from 0.3 to 1% and hydraulic heads typically between 0 and 30 m a.s.l. Finally, the Pingshang and Longgang river basins have similar settings where groundwater discharges directly to the river valleys. Hydraulic gradients vary strongly as a function of the predominant lithology along the flow path and range from 0.1% in the alluvial aquifer to 4% in the granite. Carboniferous karstified layers, located at the bottom of the Pingshang and Longgang valleys, are fully saturated, with a groundwater flow that reflects the shallow groundwater paths.

As expected from the hydrogeological conceptual model, the simulated groundwater basins coincide with the surface watersheds. The estimated groundwater balance is summarized in Table 2, where ‘average groundwater discharge’ is the annual flux directed into the rivers, reservoirs or the ocean, according to the black arrows of Fig. 6. The model simulations demonstrate that urban water recharge is an important part of Shenzhen’s water balance, consistent with Lerner (1990, 2002).

Table 2 Estimated groundwater balance from the Shenzhen numerical model
Full size table

Limits and uncertainties of the numerical analysis

The conceptual model incorporates the most plausible geological and hydrogeological settings, based on the stratigraphic and the structural analysis from borehole data. All the available information has been used for the model construction and an uncertainty analysis of alternative conceptual models (Reefsgaard et al. 2006) cannot be performed. The parameter values of the calibrated model are consistent with the literature ranges indicated by GBGP (Geological Bureau of Guangdong Province, China, Unpublished report on regional hydrogeological survey, scale 1:200.000, 1980) and Yang et al. (2007), and are in accordance with the aquifer and aquitard definition of the hydrogeological conceptual model, and have hydrosense (Hunt and Zheng 2012); however, many combinations of parameter values can give similar results in groundwater modelling. A sensitivity analysis was performed that showed that small variations in hydraulic conductivity and recharge (±10%) slightly increase the error in head calibration but do not change the groundwater flow paths, whereas, on the contrary, the hydraulic head distribution is not sensitive to the conductance specified in the river package. The study has the objective of identifying and quantifying the groundwater resources at a regional scale from the current constrained dataset. Many limitations are inherent from the regional character of this study and a formal uncertainty analysis of the numerical results (Refsgaard et al. 2007) is not possible. Parameter values are suitable to describe the regional flow paths; downscaling to analyze a single catchment, though, would require a new calibration process. Shenzhen is one of the pilot areas for the Chinese “Sponge City Initiative” (Shenzhen Committee 2016), involving the introduction of LID (low impact development) elements to store rainfall in the subsoil and ease its reuse (Chang et al. 2018). In order to assess the interaction between urban drainage and groundwater, the model would require finer discretization and the inclusion of a sewer network as a boundary condition.

Despite the high number of well-distributed hydraulic head measurements, there is a lack of measured hydraulic heads in the areas with the highest heads near Yangtai Mt. and the Yantian Mts. In these areas, the simulated hydraulic head suggests a thickness of unsaturated granites from 400 to 800 m. Supplementary hydraulic heads in these areas would help the calibration process by providing more data at locations away from the boundary conditions (i.e. main rivers and coastline); nevertheless, since the granite is an aquitard, changing the hydraulic head maxima would not influence the flux directions.

The diffuse urbanization and the construction of numerous reservoirs make unreliable the calibration of the model via discharge measures. Indeed, the estimation of the base flow cannot be performed in the areas located downstream the reservoirs, as here, the flux is controlled by the water authorities and can also vary in a small lapse of time. Obviously, reliable baseflow data also cannot be inferred in the urban areas. Data can be collected just in the rural portion, where small sub-catchments close to the watersheds have not been urbanized. Discharge measures would not be sufficient to homogeneously cover the study area and provide a reliable calibration for the whole model.

Discretization also creates inaccuracies, and more cells would better define morphology as well as the geometries of the aquifers. The top layer has a thickness of 30 m and was introduced to represent the Quaternary clastic complex and the eluvium/colluvium overlying the bedrock. In reality, the thickness is not homogeneous and has an inverse relationship with the slope angle; indeed, eluvium/colluvium is usually gravitationally mobilized on slopes steeper than 30°, leaving the bedrock exposed. A better definition of the eluvium/colluvium thickness would improve the model, defining the small perched aquifers that feed the ephemeral mountain streams. The role of these aquifers has been omitted from this regional-scale study.

Conclusions

A hydrogeological conceptual model of the Shenzhen area (Figs. 2 and 3) and numerical analysis identified and quantified the Shenzhen groundwater resources (Fig. 6). Despite a lack of measurements of hydraulic conductivity and recharge, the model fits well to the observed head data. The Quaternary coastal plain, alluvial valleys and eluvium/colluvium blankets represent the main aquifer. Tropical humid conditions have created the eluvium-colluvium blanket, which facilitates infiltration. The dense vegetative cover, in the nonurban areas, increases the porosity of the soil which slows the runoff and also facilitates infiltration. Bedrock is a regional aquitard with relatively unimportant circulation. According to the conceptual model watersheds and groundwater basins coincide. The numerical simulation based on this conceptual model constrained the hydraulic conductivity and recharge dataset.

Over 40% of the land area of Shenzhen is densely urbanized, and in this area, infiltration from precipitation is small but leakage from water mains, sewers and irrigation excess provides a significant source of groundwater recharge. Estimated urban recharge from model calibration is approximately 250 mm/year, a value which is significant, and consistent with Lerner (2002) who found that in urban areas, recharge is equal to or higher than in natural areas. In natural areas, recharge varies as a function of the outcropping hydrogeological units.

The conceptual model and associated numerical simulation provide a comprehensive framework for understanding groundwater conditions in Shenzhen and supporting future management of the city. The shallow aquifer is important in supporting streams which are environmentally significant. A sustainable development and management of the water resources is essential for the Shenzhen area and this study represents an important first step for pollution remediation of the aquifers, and the water-table restoration. An adaptation of the model could be useful to analyze the “sponge city” concept. The same approach of a block-structured hydrogeological model and calibrated dataset can be reapplied to describe groundwater flow in other areas of Guangdong and nearby Fujian province which have similar geological and climate features, and high population density.