Three-Dimensional Hierarchical Hydrogeological Static Modeling for Groundwater Resource Assessment: A Case Study in the Eastern Henan Plain, China

Abstract

Groundwater is closely related to hydrogeological structure and hydro-lithology, which mainly refers to the spatial distributions and properties of the environment where groundwater occurs. To analyze the constraints of hydrogeological structure and hydro-lithology on regional groundwater resources in the Eastern Henan Plain, China, we reconstructed the three-dimensional (3D) hierarchical models at two scales, hydrogeological structural models and hydro-lithological models, using hydrogeological cross-sections. First, the hydrogeological structural models of four aquifer groups, corresponding to four formations of the Quaternary in the study area, were reconstructed. Second, the hierarchical hydro-lithological model was built using SIS and IK estimation under the constraint of each aquifer group model space, respectively. Compared to global model, the variograms of hierarchical model captured more spatial characteristics of lithology in each aquifer group. The IK hierarchical model presents more continuities, clear boundaries, and realistic geometric shapes of the three lithologies, excluding the banding characteristics of the IK global model. The hierarchical SIS models reproduced the lithology distribution of each aquifer group and captured small changes in the lithology, with the smallest absolute percentage errors (APEs). Third, coupling the SIS hierarchical models and the groundwater levels, the groundwater resource in the study area was estimated to have a total volume of 1.2339 × 10⁴ m³. The shallow groundwater in the study area is mainly concentrated in Hebi City and the Puyang basin of the Yellow River, and deep groundwater is mainly concentrated in the northern Anyang City and Hebi City. Finally, the possible quantities of shallow and deep groundwater recharges were estimated for future groundwater management decision in the study area. The hierarchical hydrogeological model, groundwater resource assessment, and possible groundwater recharge estimation can also provide a basis for groundwater vulnerability, groundwater extraction, and land subsidence assessment.

1. Introduction

In arid or semi-arid areas, groundwater is the most important and even the only source of drinking water [1,2]. In recent years, due to the arid climate, reduced precipitation, increased extraction, and increasing intense human activities, groundwater extraction has also led to a series of environmental problems, such as ground subsidence, depletion of water resources, and water pollution [3,4].

Groundwater is closely related to the hydrogeological structure, which mainly refers to the spatial distributions and properties of the environment where groundwater occurs. The hydrogeological structure includes the spatial characteristics and hydro-lithological properties of various aquifers and aquicludes, which affect the characteristics of groundwater levels. To distinguish the hydrogeological structure, the key is to describe the spatial morphological characteristics and hydro-lithological properties of the aquifers. However, the definitions of aquifer and aquiclude depend on specific conditions, which are defined using lithological hydraulic conductivity. Generally, lithological properties can be used to express the spatial distribution and attribute of aquifers, because lithology is often more abundant and easily available in actual investigations. In addition, the visualization of the spatial structure of aquifers is also very important, and an integrated 3D hydrogeological and hydro-lithological model can be especially useful.

Three-dimensional (3D) geological modeling is a useful tool that can provide more realistic descriptions for subsurface geological phenomena and structures [5,6], and it is widely used in the development and use of urban underground space [7], geological disasters [8], metallogenic prediction [9,10,11], hydrocarbon exploration [12,13,14], hydrology and water resources assessment [15], and groundwater management. The 3D hydrogeological model enables the evaluation of groundwater flows under future scenarios, which is a basic support tool for decisions concerning the management of groundwater resources [16]. Using lithology data to construct a 3D hydro-lithological model of aquifers is easy to achieve. In addition, groundwater resource management needs the models to understand how the groundwater system works, to estimate available groundwater resources, and to simulate the impacts of groundwater exploitation [17].

The hydrogeological structure model is used to express the hydrogeological structure and spatial distribution of an aquifer, which is the research basis for groundwater and its environment [18]. Katsuaki et al. (1996) constructed a model of shapes of water-bearing strata and applied it to estimate the total volume of the groundwater in a sedimentary basin [19]. Whether for groundwater resource protection or urban development, the demand for hydrogeological models is constantly increasing. To meet the long-term water supply/demand in Turku, Artimo et al. (2003) built a 3D aquifer solid model, which represents the geometry, topology, and hydro-stratigraphy to support urban water resource planning [20]. Nasrin Nury et al. (2010) constructed a 3D geological model of the subsurface aquifer system of the Barwon Downs Graben in Victoria, Australia, using the borehole database, hydrogeological data, geological information, and surface topography. [21]. Three-dimensional hydrogeological modeling can also distinguish hydro-facies in which the lithology has the same hydraulic conductivity [22,23,24]. In this way, the aquifers and aquicludes can be distinguished by modeling from lithological data. Nuria et al. (2019) constructed a 3D geological model based on sequential indicator simulation (SIS) to identify six hydro-facies in Doñana National Park [17]. Medina Ortega et al. (2019) constructed a 3D model of four hydro-facies based on SIS to characterize the heterogeneous hydrogeological system of the Mexico City aquifer using borehole lithological records [25]. In particular, the sediments have strong spatial variability between different Quaternary formations, and there are strong heterogeneities within each aquifer group. However, previous works usually deemed the aquifer group as unitary and homogenous sediments, which would cause the deviations of groundwater resource and recharge estimation. In summary, the lithology data of boreholes or drillings can be directly used for reconstructing hydrogeological structural models. However, the boreholes are usually sparsely distributed within a region and are limited in quantity. A small number of boreholes or drillings as conditional data can only be used to construct an uncertain 3D hydrogeological model. Using a hydrogeological cross-section as the modeling data, connected by the boreholes with the regional geological knowledge, will highly improve the quality of the 3D hydrogeological model with more constraints. To facilitate research, hydrogeological models often use lithological data to express the spatial distributions and attributes of aquifers and aquicludes.

In 3D modeling, the Kriging estimation of geostatistical methods and conditional simulation are widely used to express the spatial variation characteristics of lithology and to construct hydrogeological models in a geological field. Geostatistical methods, for example, simple Kriging, ordinary Kriging, and indicator Kriging, have been widely used in hydrogeology to characterize spatial variations of the hydrogeological parameters on the basis of investigations and other relevant measurements [26]. The geostatistical characterization of a complex aquifer system and hydrogeological structural modeling at the regional scale have achieved a great deal [27,28]. Currently, the indicator Kriging (IK) algorithm is the dominant method for estimating categorical variables. Thus, it enables the production of simulations that can account for local variations of lithology distribution. The IK algorithm works well for non-parametric and categorical variables, such as the stratigraphical and lithological type, and it has been widely used in soil mapping [29], sedimentology [30], and hydrogeology [31]. In addition, the IK algorithm is commonly used to estimate the probability of distribution or single location uncertainty because it can obtain the conditional cumulative distribution function (CCDF). However, since IK is designed based on the kriging estimator, it has been criticized for its smoothing effect [32,33]. According to Goovaerts (1997), the smoothing effect is minimal when the locations of the observed data are closely spaced, and the smoothing effect increases as the distance among the observed data increases [29]. Unlike local estimation of interpolation, stochastic simulation considers the overall characteristics of the results and the statistical spatial correlation of the simulated values first, and the accuracy of the local estimates last. Therefore, stochastic realizations overcome the smoothing effect and are more realistic in representing spatial heterogeneity [34]. Sequential indicator simulation (SIS) is a conditional stochastic simulation algorithm that is straightforward in constructing a CCDF [35]. IK is designed to minimize local criteria, whereas the objective of SIS is to reproduce a global histogram and variogram. SIS is increasingly preferred over Kriging in the cases where spatial variation measured in the field must be preserved [32]. However, in practice, some studies show that the mean prediction error tends to be larger for simulated values than for kriging estimates [29,34,36]. Therefore, the selection of estimation or simulation to predict lithology properties may involve trade-offs in terms of errors of the results not only in prediction accuracy but also in the reproduction of spatial variability.

This study attempted to directly reconstruct a 3D hydrogeological model using hydrogeological cross-sections, and also compared IK and SIS to reconstruct a three-dimensional hydro-lithological model considering the spatial heterogeneity of these lithologies. In addition, the Quaternary sediments were divided into four aquifer groups from top to bottom, namely the Holocene Series (Q₄), Upper Pleistocene (Q₃), Middle Pleistocene (Q₂), and Lower Pleistocene (Q₁). This study also took the boundaries of the four aquifer groups as constraints to hierarchically construct the model of lithologies, and the details of the hydro-lithologies were better expressed in this hierarchical hydrogeological model. Coupling shallow and deep groundwater levels with the hierarchical SIS model, the total volumes of shallow and deep groundwater resources were estimated. Meanwhile, the possible quantitative recharges of shallow and deep groundwater were assessed on the 3D hydro-lithological and groundwater level models.

2. Study Area and Dataset

2.1. Physical Geographical Background

The Eastern Henan Plain is located at 112°31′~116°40′ E and 32°00′~36°12′ N, Henan Province, in mid-eastern China, with an area of about 8.53 × 10⁴ km². The Eastern Henan Plain is surrounded by mountains in the northwest, west, and south, and by vast plains in the east. The overall terrain is higher in the west and lower in the east, with altitudes of less than 100 m [37]. This region is a warm temperate zone with a semi-humid and semi-arid monsoon climate and with four distinct seasons. The annual precipitation is 600–1200 mm, decreasing from south to north. Precipitation is mainly concentrated from May to August, which accounts for about 60% in the annual precipitation. The annual potential evaporation in the area is generally 900–1400 mm, increasing sequentially from the west to the east, and the annual mean temperature is 14.2 °C [38]. The study area is located in the Eastern Henan Plain to the north of the Yellow River, which is also part of the Northern China Plain, as shown in Figure 1a. There are two main rivers flowing through the study area, namely the Wei River and the Yellow River. The Yellow River is the second largest river in China, and its mean annual flow at the Huayuankou Station is 4.7 × 10¹⁰ m³. The annual mean flow of the Weihe River at the Xinxiang Station is 1.7 × 10⁹ m³. The study area belongs to the alluvial plain of the Yellow River; therefore, its formation and development are closely related to the changes and flooding of the Yellow River. There are mostly saline–alkali depressions along the banks of the Yellow River, with low production levels. The study area and the layout lines of seven hydrogeological cross-sections are shown in Figure 1b.

2.2. Hydrogeological Setting

According to the intensity of the new tectonic movement, the geomorphology, genetic type, and lithofacies characteristics of the Quaternary, the study area is divided into the piedmont plain and the plain strata zones. The piedmont stratigraphic zone is distributed in the piedmont slope plain to the west of the boundary, and its material came from the western mountainous areas. It consists of mainly alluvial proluvial deposits and slope alluvial deposits, whose particles are thicker, thickness is thinner, and color is obviously different. The plain stratigraphic zone is distributed in the alluvial plain to the east of the boundary. The upper sediment mainly came from the Yellow River, and the lower part is fluvio-lacustrine deposits whose grains are relatively fine, thickness is large, and color is not distinctly different.

On the basis of the quartet method (classifying, grading, sorting, and partitioning), the Quaternary stratigraphy of the eastern Henan Plain is divided into the Wuzhi Formation (Q₁), Kaifeng Formation (Q₂), Taikang Formation (Q₃), and Puyang Formation (Q₄), which correspond to the Holocene Series, Upper Pleistocene, Middle Pleistocene, and Lower Pleistocene, respectively. Therein, the Q₄ and Q₃ contain shallow groundwater, while the Q₂ and Q₁ contain deep groundwater. The Q₄ and Q₃ aquifer groups are the early exploitation sources of groundwater in the eastern Henan Plain. With the increasing demand for groundwater, the aquifer groups Q₂ and Q₁ have become the main sources of groundwater extraction [1].

The study area contains two subsystems of groundwater, which are the Weihe River alluvial plain (I₂) and the Yellow River alluvial plain (II₂), as shown in Figure 2. Subsystem I₂ includes tributaries, piedmont alluvial–pluvial fan, and its front depression, where groundwater exists in porous and phreatic water. The middle and upper parts of the alluvial fan have good water abundance, and the lower and front areas have poor water abundance. In Subsystem II₂, groundwater exists in pore phreatic water and micro confined water. The water-bearing strata in the fan-shaped terrain and mainstream belt are mainly fine sand, medium sand, and siltstone, with abundant water.

2.3. Dataset

Seven hydrogeological cross-sections connected by drillings were obtained from the databases of the Institute of Hydrogeology and Environmental Geology, China Academy of Geological Sciences (CAGS) (2005) [40], as shown in Figure 3. The division of the aquifer groups of the hydrogeological cross-sections is mainly based on the four formations of the Quaternary, i.e., Q₄, Q₃, Q₂, and Q₁, which correspond to the Holocene Series, Upper Pleistocene, Middle Pleistocene, and Lower Pleistocene, respectively. The horizontal scale of the hydrogeological cross-sections is 1:250,000, and the vertical scale is 1:3000. The lithological categories consist of silt, clay, and aquifer sand.

We extracted the boundary lines of the aquifer groups and a total of 120,210 sampled lithological points from the seven hydrogeological cross-sections in the study area, of which there were 12,581, 36,046, 34,645, and 36,938 lithological points from Q₄ to Q₁, respectively. The horizontal distance between the sampling points along the cross-section layout lines was 250 m, and the vertical sampling distance was 3 m. Then, the aquifer group boundaries (Figure 4a) and the lithological sampling points (Figure 4b) were converted into 3D coordinates according to the match mapping between the cross-section layout lines and the cross-sections.

3. Methodology

3.1. Flowchart

We constructed a 3D hierarchical hydrogeological model of the aquifer and lithology from the hydrogeological cross-sections in the Eastern Henan Plain using the GOCAD^® software (www.pdgm.com (accessed on 1 March 2017), Paradigm, Houston, TX, USA). The key processes of this study consisted of the following steps. First, in data preprocessing, the lithologic sampling points and aquifer group boundaries were extracted from the seven hydrogeological cross-sections. Second, five interface models were constructed using the aquifer group boundaries, and the stratigraphical model was constructed. Third, variogram analysis was carried out on the lithological samplings, and then the global hydro-lithological model and four hierarchical hydro-lithological models constrained by the corresponding aquifer group models were constructed and compared using the IK and SIS methods, respectively. Finally, according to the characteristics of aquifers in the study area, the groundwater resources of the confined and unconfined aquifers were estimated using the hierarchical hydrogeological model. The workflow of this study is shown in Figure 5.

3.2. Indicator Kriging (IK)

The aim of IK is to estimate the CCDF belonging to any category $z_k$ under the condition data n [41], as follows:

$$I^{*}\left(u;z_k\right)=E^{*}\left\{I\left(u;z_k\right)\mid \left(n\right)\right\}={\mathrm{Prob}}^{*}\left(z\left(u\right)\in z_k\mid \left(n\right)\right)$$

(1)

where $I^{*}\left(u;z_k\right)$ is the estimated indicator variable of category $k$ at the location of u. $z\left(u\right)$ is the regionalized variable belonging to any category $z_k$. The IK algorithm assumes that the marginal probabilities$E^{*}\left\{I\left(u;z_k\right)\right\}$are constant and known for all $u_{\alpha }$and $z_k$ [29].

For each category k, the indicator variable is defined as follows:

$$I\left(u,z_k\right)=\{\begin{array}{cc}1& \mathrm{if} z\left(u\right)\le z_k\\ 0& \mathrm{otherwise} \end{array}$$

(2)

The probability $I^{*}\left\{\left(u;z_k\right)\right\}$ for$z\left(u\right)$ belonging to a category is estimated by simple kriging, as follows:

$$I^{*}\left(u;z_k\right)-E\left\{I\left(u;z_k\right)\right\}={\displaystyle \sum }_{\alpha =1}^N{\lambda }_{\alpha }\left(I\left(u_{\alpha };z_k\right)-E\left\{I\left(u_{\alpha };z_k\right)\right\}\right)$$

(3)

where weights ${\lambda }_{\alpha }$can be calculated from the simple kriging system, as follows:

$${\displaystyle \sum }_{k=1}^{\mathrm{n}\left(u_0\right)}{\lambda }_k\left(u_0;z_k\right)C_I\left({\mathrm{u}}_k-{\mathrm{u}}_j;z_k\right)=C_I\left({\mathrm{u}}_j-u_0;z_k\right), \mathrm{for} j=1 \mathrm{to} \mathrm{n}\left(u_0\right)$$

(4)

where $\mathrm{n}\left(u_0\right)$ weights are associated with the neighboring data, $z_k$ denotes the corresponding indicators,$C_I\left({\mathrm{u}}_k-{\mathrm{u}}_j\right)$ is the covariance between indicators ${\mathrm{u}}_k$ and ${\mathrm{u}}_j$, and $C_I\left({\mathrm{u}}_j-u_0\right)$is the covariance between the sample point ${\mathrm{u}}_j$ and the point being estimated $u_0$, with

$$C_I\left(h;z_k\right)=E\left(z_k\right)\left[1-E\left(z_k\right)\right]-r_I\left(h;z_k\right)$$

(5)

where the used indicator variograms $r_I\left(h;z_k\right)$is as follows:

$$r_I\left(h;z_k\right)=\frac{1}{2N\left(h\right)}{\displaystyle \sum }_{i=1}^{N\left(h\right)}{\left[I\left(u;z_k\right)-I\left(u+h;z_k\right)\right]}^2$$

(6)

where $h$ is the lag, and $N\left(h\right)$ is the number of data pairs separated by$h$.

Only neighboring data close to the estimated location $u_0$can actually be used; therefore, $\mathrm{n}\left(u_0\right)$<< N, and N is the total number of condition data [29]. The estimated category can be decided by the final probability distribution of the IK results.

3.3. Sequential Indicator Simulation (SIS)

Usually, the SIS algorithm relies on the IK to infer the probability density function (PDF) of the categorical variable z(u). It simulates the non-parametric distribution by combining the indicator formalism with the sequential paradigm [42]. Through stochastic simulation, a series of alternative and equally probable realizations of the distribution of an indicator variable z(u) are produced. For example, if z(u) belongs to a category k is simulated at a spatial site u, and the PDF to be estimated becomes

$$Prob\left\{I\left(u\right)=1\mid \left(n\right)\right\}=E\left\{I\left(u\right)\mid \left(n\right)\right\}$$

(7)

Using simple Kriging to estimate the probability of variables$z_k$on location u yields the following:

$${\mathit{Prob}}^{*}\left\{I\left(u;z_k\right)=1\mid \left(n\right)\right\}=p_k+{\displaystyle \sum }_{\alpha =1}^n{\lambda }_a\left[I\left(u;z_k\right)-p_k\right]$$

(8)

where $p_k=E\left\{I\left(u;z_k\right)\right\}\in \left[0,1\right]$is the marginal frequency of category$z_k$, and${\lambda }_a$ is weight of the simple Kriging system.

The detailed steps of SIS are as follows [43]:

(i)

Define a path for visiting all locations to be simulated.

(ii)

For each location u along the path, first, retrieve the neighboring categorical conditioning data $z\left(u_{\alpha }\right),\alpha =1,\dots ,N$; second, estimate the indicator random variable $I\left(u;z_k\right)$ for each of the k categories by solving a kriging system; third, define an estimate of the discrete conditional probability density function (CPDF) of the categorical variable z(u); then, sample a realization from CPDF and assign it as a datum at location u. The previous simulated results can be used as conditioning data for the later visited location. Finally, loop the above processes until all locations are visited.

(iii)

Repeat the previous two steps to generate other realizations.

4. Results

4.1. Aquifer Group Model

The study area was divided into four aquifer groups according to the formations of the Quaternary, and the boundaries of the four aquifer groups were extracted from the hydrogeological cross-sections. A 3D surface model of the five interfaces, i.e., the ground surface, Q₄ bottom surface, Q₃ bottom surface, Q₂ bottom surface, and Q₁ bottom surface, was constructed using the extracted boundaries, as shown in Figure 6.

Based on the interface model of the ground surface and four aquifer groups, four 3D stratigraphical models were constructed for the corresponding aquifer groups to serve as a framework for the hydro-lithological IK estimation and SIS simulation, as shown in Figure 7. The resolutions of the stratigraphical models in the X, Y, and Z directions were determined by the densities of the lithological sampling points, which were 250 m, 250 m, and 3 m, respectively. The grid numbers of the aquifer groups Q₄ to Q₁ were 1,249,875, 3,499,902, 4,916,116, and 5,916,004, respectively.

4.2. Variograms

Most geostatistical algorithms use to model properties take spatial autocorrelations into account, and this spatial autocorrelation is often analyzed and described by a variogram. In this study, the hierarchical hydrogeological modeling was carried out at two scales, and the lithological variograms at the two scales of the global and aquifer groups were plotted, respectively. A variogram is a cross-plot showing how the lithology changes with separation distance, in other words, how data points become uncorrelated as the distance between them increases. Martinius et al. [44] argued that the spatial structure of sedimentary facies should be considered comprehensively to determine the ellipse size used to search the simulated grid points as conditional data around an estimated grid node before stochastic simulation. Usually, plotting the variograms on the four azimuths related to the depositional direction can focus on simulating the lithological deposition direction and its vertical direction, eliminating the error caused using a single variogram curve to express the degree of variation in all directions [45].

After transforming the sampled lithology points into indicator data, variograms were calculated for the vertical and horizontal directions based on the dataset. GOCAD software was used for the plotting of variograms for both IK and SIS, and the scatter plot of variogram$h$distance was generated (Figure 8). The X-axis represents the separation distance between data points, and the Y-axis represents the variance of data pairs determined by the distance. Each point in the variogram represents the average variance computed for all the pairs of points with the same spacing. The theoretical variogram functions often used to fit the experimental values in geostatistics include the spherical, exponential, Gaussian, and linear-with-sill models in GOCAD [46]. In this study, the variograms in the horizontal and vertical directions were plotted for all the global and aquifer group sampling points, respectively, and the spherical function fitting was performed on the discrete points of the horizontal and vertical variograms with the best-fitting value. The purpose of variogram analysis is mainly to obtain the major horizontal, minor horizontal, and vertical ranges. In the lithofacies modeling, the range value is a variation of the extension scale of the lithology, which has certain effects on the prediction of the lithology scale and the division of the lithofacies’ guiding significance [47]. The variogram function shows that the hierarchical experimental variation values of the three types of lithologies fit better with the theoretical model than the global ones. The determination of these parameters in the theoretical variogram model of each aquifer group is reasonable because the fitted theoretical model reproduces the variation of the experimental variogram.

When fitting the variogram, the most important thing is to determine the parameters of the theoretical model, such as the range, nugget, still, and contribution. In this study, all variograms were fitted using a spherical model and a series of parameters of the fitted variograms of each lithology, for example, nugget C₀, contribution C, sill C₀ + C, and the nugget-to-sill ratio (C₀/(C₀ + C)) were calculated in the major horizontal, minor horizontal, and vertical directions. The range is the maximum correlation distance of the sampling points in the space and directly determines the lithology distribution distance in the hydro-lithological model. The nugget-to-sill ratio designates the degree of spatial heterogeneity arising from random components to that of the total spatial heterogeneity [48]. A nugget-to-sill ratio value close to 0 indicates that the variable has strong spatial autocorrelation. Conversely, a nugget-to-sill ratio value close to 1 indicates the spatial heterogeneity is dominated by randomness or the nugget effect [35]. The parameters of the global variograms of silt, clay, and sand are shown in

Table 1. Parameters of global indicator variograms.

Lithology	Nugget (C0)	Sill (C0 + C)	Contribution (C)	Major Range (m)	Minor Range (m)	Vertical Range (m)	C0/ (C0 + C)
Silt	0.05	0.09	0.04	13,499.2	42,682.2	301.21	0.56
Clay	0.02	0.25	0.23	15,173.6	1,858,690	289.53	0.08
Sand	0.15	0.22	0.07	11,253.2	1,874,810	317.20	0.68

. We can see that clay in the horizontal direction has a strong spatial auto-relation (C₀/(C₀ + C) = 0.08), reflecting regular sedimentary deposition. However, the nugget-to-sill ratio of the silt (0.56) and sand (0.68) in the horizontal direction exhibits some degree of nugget effect, indicating random heterogeneity and reflecting complicated formation conditions at the scale of sampling. The major horizontal, minor horizontal, and vertical ranges of the three lithologies are also shown in Table 1. The correlation range observed depends on the direction, explained by the lateral extent of the lithologies, which is usually greater than their thickness; hence, the vertical range is much smaller than the horizontal range [49]. The clay presents the longest horizontal integral scale, and the silt and sand resulted in shorter horizontal integral scales than the clay.

As shown in

Table 2. Parameters of indicator variograms for each aquifer group model.

Lithology	Aquifer Group	Nugget (C0)	Sill (C0 + C)	Contribution (C)	Major Range (m)	Minor Range (m)	Vertical Range (m)	C0/ (C0 + C)
Silt	Q4	0.15	0.22	0.07	10,510.60	55,095.40	17.15	0.68
	Q3	0.05	0.09	0.04	4568.09	52,030.40	51.73	0.56
	Q2	0.05	0.11	0.06	27,221.40	11,843.50	113.09	0.45
	Q1	/	/	/	/	/	/	/
Clay	Q4	0.10	0.22	0.12	21,971.60	41,831.60	45.73	0.45
	Q3	0.15	0.23	0.08	14,588.20	57,652.30	94.11	0.65
	Q2	0.15	0.24	0.09	14,150.40	13,203.30	377.32	0.63
	Q1	0.13	0.15	0.02	36,818.80	4631.67	44.30	0.46
Sand	Q4	0.15	0.22	0.07	11,058.20	41,399.60	18.39	0.87
	Q3	0.15	0.22	0.07	17,764.80	23,082.30	49.67	0.87
	Q2	0.15	0.21	0.06	43,156.30	24,813.30	78.02	0.71
	Q1	0.13	0.16	0.03	28,522.60	40,387.10	118.79	0.81

, the variograms of each aquifer group were fitted separately, and the nugget-to-sill ratio value of each aquifer group was calculated according to the parameters of the theoretical variograms. Compared to the global model, the nugget-to-sill ratio values of the silt and sand were similar in the hierarchical models, showing some degree of nugget effect in the horizontal direction. Conversely, the nugget-to-sill ratio values of the clay were 0.45 and 0.46 in the Q₄ and Q₁ aquifer groups, respectively, which indicates moderate spatial autocorrelation. The nugget-to-sill ratio value of the clay was higher in the Q₃ and Q₂ aquifer groups, which exhibits some degree of nugget effect in the horizontal direction, indicating random heterogeneity. The ranges of lithologies in each aquifer group, and the silt, clay, and sand had the biggest range in the Q_2, Q₁, and Q₂ aquifer groups, respectively. Although the proportion of sand is relatively large in the Q₃ aquifer group, its distribution is relatively concentrated so that the final fitting major range of sand is relatively small.

4.3. IK Models

An IK estimation was conducted on global grids 250 m long, 250 m wide, and 3 m high (10,248,975 voxels) using all lithological sampling points. To compare it with the hierarchical IK lithological model of the four aquifer groups, the global IK model was assigned to the Q₄–Q₁ aquifer groups, as shown in Figure 9. As expected, sand, corresponding to lacustrine sediments, was mostly horizontally distributed, and appeared almost continuously in the upper portion of each aquifer group. The clay occupied a relatively large volume in the model, which is also consistent with the original statistics of hydrogeological cross-sections. We will further discuss the difference between the proportion of each lithology in the IK model and the sampled data using absolute percentage errors (APEs) in Section 5.1. The silt was mainly distributed in aquifer group Q₄ of the model in the range of roughly 10–30 m. The silt and the clay formed a fine aquiclude or aquitard, and the interbedded sand aquifer occurred in the horizontal direction. However, the global IK model showed that the sand formed a series of regular bands in the Q₄ formation, as shown in Figure 9a. The global model analysis showed that the IK overestimated the spatial continuity and proportions of the sand and clay. The boundaries of the sand and clay were not clear and very messy in the Q₃, Q₂, and Q₁ aquifer groups, and sand occurred in huge blocks in the southeastern part of the study area, where few sampling points were located (Figure 9d). Moreover, the silt sampling points did not appear in the Q₁ aquifer group in the study area; however, there was silt in the Q₁ aquifer group of the IK global model, which proves that the global IK model cannot accurately reconstruct the distribution of the lithology in the Q₄, Q₃, Q₂, and Q₁ aquifer groups without aquifer group constraints.

Therefore, the lithological sampling points distributed in the Q₄, Q₃, Q₂, and Q₁ aquifer groups were extracted from the samplings of the study area. Sequentially, variogram analysis was performed on the sampling points of each aquifer group, and the theoretical variograms were fitted. Then, with the constraint of the corresponding aquifer group models (Figure 7), the IK estimator was used to hierarchically construct the 3D lithological models of the four aquifer groups, as shown in Figure 10. Compared to the IK global model, the 3D hierarchical IK model ensured the extension range of the lithology sampling points in each aquifer group. The spatial distribution of the Quaternary lithologies was expressed well in these constructed models. In the Q₄ aquifer group model (Figure 10a), the silt presented some continuity at the southeastern and western parts of the study area. The clay showed more continuity at the northwestern and southeastern portions of the Q₃ and was interbedded with the sand, which may have interrupted its continuity. In the Q₂ and Q₁ aquifer groups, the clay aggregated in a large area and presented more continuity. The sand bands were excluded from the four aquifer groups, and the boundaries of the three lithologies were very clear.

To analyze the distribution and shape of sand in a 3D space, sand was extracted from the four IK aquifer groups, as shown in Figure 11. All the interpolated sands had a good and realistic spatial continuity. The hierarchical IK lithological model estimated the least volume of sands and can reproduce a realistic shape and size of most of the sands. The horizontal distribution of sand was reproduced by IK hierarchical lithological model. In the Q₄ aquifer group, the sand was mainly distributed in the Weihui, Changyuan and southern Puyang. The sand has a large extent in the Q₃ aquifer groups. Most of the study area is occupied by sand except for the Qixian and Huaxian areas. In the Q₂ aquifer group, sand is lacking in the Weihui, Huaxian, Neihuang, northern Changyuan, and southwestern Puyang. In the Q₁ aquifer group, the sand was concentrated in the central part of the study area, i.e., the Weihui, Changyuan, Huaxian, and southern Puyang.

4.4. SIS Models

A SIS was also conducted on the same grid, and Figure 12 shows the global SIS model of lithologies in the study area was assigned to the Q₄–Q₁ aquifer groups. Obviously, the general pattern was captured in the simulated map, where the sand is distributed horizontally. However, the sand of the global SIS model was present to a large extent in the four aquifer groups, which is not in good agreement with the observed sand in the original hydrogeological cross-sections. Moreover, there were a series of band shapes of sand in the Qixian and Puyang areas. The global SIS model in the Q₄ aquifer group shows that all the silt was underestimated in the spatial continuity and proportions in the southeastern portion, which is near Puyang City, Changyuan, and Huaxian areas. The silt was also simulated in the Q₄ aquifer group, which means that the range of the lithology simulation will be misestimated without the constraints of the aquifer groups. In the global SIS model (Figure 12), the clay was interbedded with the other two lithologies, which interrupted its continuity. The clay was underestimated, especially in the Q₁ aquifer group. However, the lithologies in the global SIS model could not clearly show their 3D geometrical shapes.

A SIS hierarchical hydrogeological model was constructed using the constraints of the boundaries of the four aquifer groups of Q₄, Q₃, Q₂, and Q₁, as shown in Figure 13. The realizations obtained with the separate simulations of the four aquifer groups were by far more realistic (Figure 13). The hierarchical SIS lithological model takes into account the differences between the four aquifer groups that are evident in the lithological samples of the study area and worked better over individual aquifer groups than over the global study area. Compared to the global SIS model, the silt showed more continuity at the southern portion of the Q₁ aquifer group and had clear boundaries in the hierarchical SIS model (Figure 13a). In the Q₃ aquifer group (Figure 13b), sand presented a realistic shape and more continuity, which is well connected among the hierarchical models. The proportions of sands in the Q₂ and Q₁ aquifer groups were distinctly reduced, but the boundaries were not obvious. The banding feature did not appear in the Q₄ aquifer group model. The clay was almost all distributed in Q₂ and Q₁, and its continuity was very strong, and it was interbedded with the other two lithologies. This paper aims to highlight the possibilities of developing groundwater numerical models (i.e., noncontinuous lithofacies without smooth effect limits) from geostatistical estimations. Hence, we calculated the expected groundwater volume with a hierarchical SIS realization.

4.5. Validation

The correlation of the hierarchical IK lithological model of the four aquifer groups with boreholes lithological logging is shown in Figure 14. The lithology predicted by the hierarchical IK model has a roughly similar trend between boreholes. The rates of coincidence of the predicted lithology with the boreholes C2, G6, C8, and G8 are 70%, 66%, 73%, and 63%, respectively. In summary, the overall validation results of the model were relatively good.

5. Discussion

5.1. Absolute Percentage Error (APE) of Models

Since the proportion of each lithological category is an important evaluation indicator, herein, the absolute percentage error (APE) was calculated according to Equation (9), which is the difference between the sample proportion ${\mathrm{A}}_{\mathrm{t}}$ and the predicted model proportion ${\mathrm{A}}_{\mathrm{i}}$ of a certain category divided by ${\mathrm{A}}_{\mathrm{t}}$ [35].

$$\mathrm{APE}=\left|\frac{{\mathrm{A}}_{\mathrm{t}}-{\mathrm{A}}_{\mathrm{i}}}{{\mathrm{A}}_{\mathrm{t}}}\right|$$

(9)

To evaluate the APEs of the global and hierarchical models, the proportions of the original lithological sampling points were first calculated, as shown in

Table 3. The proportions of the lithological sampling points in the global space and each aquifer. group.

Statistic	Samples	Silt	Clay	Sand
Proportion (%)	Global model	10.30	56.40	33.30
	Q4	33.60	35.00	31.40
	Q3	10.90	36.90	52.20
	Q2	12.20	57.30	30.60
	Q1	/	81.90	18.10

.

The APE of the global IK model and the global SIS model were calculated, as shown in

Table 4. Proportions and absolute percentage errors (APEs) of global IK and SIS models.

Statistic	Model	Silt	Clay	Sand
	Samples	10.30	56.40	33.30
Proportion (%)	IK	5.60	60.8	33.70
	SIS	9.20	48.10	42.70
APE (%)	IK	36.15	7.8	1.20
	SIS	10.70	14.72	28.23

. The APEs of clay and sand in the IK model were smaller than those in the SIS model; however, the APE of silt was larger, which is related to the proportions of original samplings and the smoothing effect of the IK estimator for the local optimum. Overall, the lithologic reconstruction of the IK global model, except for silt, was more accurate than the SIS global model.

The APE of the global IK model and hierarchical IK model of the four aquifer groups are compared in

Table 5. The proportions and absolute percentage errors (APEs) of the global IK model and hierarchical IK model of the four aquifer groups.

Statistic	Model	Aquifer Group	Silt	Clay	Sand
Proportion (%)	Global model	Q4	20.7	46.7	32.6
		Q3	6.28	40.0	53.7
		Q2	5.61	66.4	28.0
		Q1	0.53	79.1	20.4
	Hierarchical model	Q4	29.3	46.7	24.00
		Q3	6.70	33.50	59.80
		Q2	8.40	65.50	26.10
		Q1	/	89.90	10.10
APE (%)	Global model	Q4	38.4	33.4	3.8
		Q3	42.4	8.4	2.9
		Q2	54.0	15.9	8.5
		Q1	/	3.4	12.7
	Hierarchical model	Q4	12.78	33.13	23.57
		Q3	38.53	9.21	14.56
		Q2	31.15	14.31	14.71
		Q1	/	9.77	44.20

. Compared to the global IK model, the silt performed better in the hierarchical IK model, especially in the Q₄ aquifer group. However, except for the clay in the Q₄ and Q₂ aquifer groups, the clay and sand have smaller AEPs in the IK global model, which is closer to the corresponding proportions of the original lithological sampling points.

The proportions and APEs of the global SIS model and hierarchical SIS model of the four aquifer groups are shown in

Table 6. Proportions and absolute percentage errors (APEs) of global SIS model and hierarchical SIS model of four aquifer groups.

Statistic	Model	Aquifer Group	Silt	Clay	Sand
Proportion (%)	Global model	Q4	17.9	39.4	42.7
		Q3	10.7	37.0	52.3
		Q2	8.6	52.0	39.4
		Q1	5.4	58.6	36.0
	Hierarchical model	Q4	35.00	31.90	33.10
		Q3	8.20	37.00	54.80
		Q2	10.70	59.00	30.30
		Q1	/	82.80	17.20
APE (%)	Global model	Q4	46.7	12.6	36.0
		Q3	1.8	0.3	0.2
		Q2	29.5	9.2	28.8
		Q1	/	28.4	98.9
	Hierarchical model	Q4	4.17	8.86	5.41
		Q3	24.77	0.27	4.98
		Q2	12.30	2.97	0.98
		Q1	/	1.10	4.97

. The SIS global model performed better than the hierarchical SIS model in the Q₃ formation. However, the APEs of the aquifer groups were smaller than the global model in the Q₄, Q₂, and Q₁ strata, which proves that the hierarchical SIS model of the four aquifer groups was more accurate than the global SIS model in lithology reproduction. In brief, this indicates that the hierarchical SIS lithological model proposed in this work is able to reproduce the basic statistics of the informed data.

The lithological proportions of each aquifer group in the IK and SIS hierarchical models were compared to the samplings (Figure 15), which showed that the lithological proportions of the hierarchical SIS model were closer to the original samplings and more accurate. This result is consistent with the calculated APEs of the hierarchical SIS model of the four aquifer groups.

5.2. Groundwater Resource Assessment

Comparing the global IK model, the hierarchical IK model of the four aquifer groups, the global SIS model, and the hierarchical SIS model, the hierarchical SIS is more realistic because the APE value was the smallest, and the proportion of each lithology was closest to the original samplings; therefore, the hierarchical SIS models of Q₄–Q₁ were used in a subsurface groundwater resource assessment. The groundwater level surfaces of the shallow groundwater system (mainly including Q₄ aquifer group) and deep groundwater system (main including Q₃, Q₂, and Q₁ aquifer groups) were coupled with the hierarchical SIS hydrogeological model to quantitively estimate the reserve and possible recharge of groundwater, respectively. According to the hierarchical SIS model, the quantities of possible shallow and deep groundwater recharges were estimated as the corresponding difference between the theoretical groundwater storage of sand and the reserve of groundwater, respectively, as shown in Figure 16.

The shallow and deep groundwater level lines were extracted from the cross-sections, and interpolated into 3D surfaces, as shown in Figure 17.

The main aquifer, sand, is distributed into four aquifer groups. Therein, Q₄ is an unconfined aquifer group, about 60 m thick; Q₃ is a shallow confined aquifer group, 60 m thick; Q₂, more than 90 m thick, is a confined aquifer group; Q₁ is 50–60 m thick and is the deep confined aquifer group [1].

The groundwater volume of the unconfined aquifer group Q₄ was estimated according to the volume of sand voxels multiplied by the effective porosity (or specific yield) S [50]. According to the definition of effective porosity, the S value of aquifer group Q₄ is 0.3.

$${\mathrm{V}}_{\mathrm{w}}={\mathrm{V}}_{\mathrm{uc}}\times \mathrm{S},$$

(10)

where ${\mathrm{V}}_{\mathrm{w}}$ is the volume of groundwater (milliliter), ${\mathrm{V}}_{\mathrm{uc}}$ is the volume of unconfined aquifer voxels (m³), and$\mathrm{S}$ is the specific yield.

The groundwater volume (elastic storage) of the confined aquifer groups, i.e., Q₃, Q₂, and Q₁, were estimated according to Equation (10) [21]. The ${\mathrm{S}}_{\mathrm{c}}$ is the storage co-efficient related to the thickness and the porosity of the aquifer. Because the thicknesses and porosities of the three confined aquifer groups in the study area are different [1,51], the storage co-efficient of aquifer groups Q₃, Q₂, and Q₁ are set as 6.915 × 10⁻², 1.034 × 10⁻², and 1.013 × 10⁻², respectively.

$${\mathrm{V}}_{\mathrm{w}}={\mathrm{V}}_{\mathrm{c}}\times {\mathrm{S}}_{\mathrm{c}},$$

(11)

where ${\mathrm{V}}_{\mathrm{w}}$ is the volume of groundwater (milliliter), ${\mathrm{V}}_{\mathrm{c}}$ is the volume of confined aquifer voxels (m³), and ${\mathrm{S}}_{\mathrm{c}}$ is the storage co-efficient.

The total reserves of groundwater resources in shallow and deep aquifer groups were estimated and mapped, as shown in Figure 18. The shallow groundwater in the study area is mainly concentrated in Hebi City and the Puyang basin of the Yellow River. During the Holocene, these two areas were alluvially deposited by the Weihe River and Yellow River to form thicker sand aquifer. The deep groundwater is mainly concentrated in northern Anyang City and Hebi City. The strata experienced severe deformation activities in the Middle Pleistocene, and a huge depression was formed in the Q₂ aquifer group in the northern Anyang City. In the late Pleistocene, the topography of the Q₃ aquifer became flat and received abundant sediments to form thicker sands.

The quantities of possible shallow and deep groundwater recharges were calculated, as shown in Figure 19. The main possible recharge area of shallow groundwater is distributed in the Changyuan-Puyang area, which is along the Yellow River. The sediments carried by it precipitated and accumulated here to form a series of sands which provides a possible recharged reservoir of shallow groundwater. The main possible recharge areas of deep groundwater are distributed in northern Anyang City and southern Hebi City. In the Late Pleistocene, the strata received a lot of sediments to form thicker sand aquifers, which made it possible to host a large amount of groundwater.

Table 7. The volume of each lithology and groundwater reserve and recharge in the study area.

Code	Lithology	Hydraulic Property	Volume $\left({\mathbf{m}}^3\right)$	Reserve $\left({\mathbf{m}}^3\right)$	Recharge $\left({\mathbf{m}}^3\right)$
1	Silt	Aquitard	7.903 × 1010
2	Clay	Aquitard	5.178 × 1011
3	Sand	Aquifer group Q4	2.731 × 1010	5.187 × 103	3.008 × 103
		Aquifer groups Q3–Q1	2.277 × 1011	7.152 × 103	1.645 × 103

gives the volume of each lithology and reserve and recharge of groundwater in the 3D SIS hierarchical hydrogeological model. We estimated the aquifer volume of Quaternary as 2.5501 × 10¹¹ m³ for 5070.45 ${\mathrm{km}}^2$ using the lithological models of four aquifer groups. The aquifer thickness and porosity were fully considered when evaluating the groundwater volume. The groundwater volume is totally estimated to be 1.2339 × 10⁴ m³ (including 5.187 × 10³ m³ shallow groundwater and 7.152 × 10³ m³ deep groundwater) in the study area, and the possible groundwater recharge is calculated as 4.653 × 10³ m³ (including 3.008 × 10³ m³ shallow groundwater recharge and 1.645 × 10³ m³ deep groundwater recharge).

6. Conclusions

We present a methodology to integrate accessible stratigraphical and lithological information of hydrogeological cross-sections for a 3D hierarchical hydrogeological model in the Eastern Henan Plain, China. A 3D hierarchical hydrogeological model was constructed at two scales, i.e., a stratigraphical model and a lithological model, which analyzed the hydrogeological characteristics from the spatial structure and the lithological statistics, respectively. From the perspective of spatial structure, the spatial distribution of three lithologies was analyzed in the global model and hierarchical model. Comparing the global IK model with the hierarchical IK model, the sand of the global model presents band-like distributions in the study area. However, the hierarchical model clearly shows the sand’s geometric shape and spatial distribution characteristics, which proved that the hierarchical model can accurately control the extension range of lithology to reflect the lithological characteristics in each aquifer group. The sand also appears with striped distributions in the global SIS model; however, these incorrect characteristics do not appear in the hierarchical SIS model. From the lithological statistics, the APEs of the 3D global and hierarchical models by the IK and SIS methods were calculated, respectively. Although the APEs of the global model are smaller, this is related to the smoothing effect of the IK estimator due to the original proportion of lithological samplings. The APEs in the hierarchical SIS model of the three aquifer groups were smaller than the global model. Therefore, the proportion of each lithology in the hierarchical SIS model was closer to the original samplings. Moreover, comparing the hierarchical IK model and the hierarchical SIS model, the lithology proportion of the SIS model was closer to the original samplings, which also proves that the hierarchical SIS model is more accurate.

The 3D hydrogeological modeling of the study area offers a new perspective for sustainable groundwater management. These aquifer groups are important sources of fresh groundwater and have complex stratigraphical and lithological settings. Coupling shallow and deep groundwater level surfaces with the hierarchical SIS model, the amount of shallow and deep groundwater reserves in the study area was assessed to be 1.2339 × 10⁴ m³, taking into account the influence of each aquifer groups’ type, thickness, and porosity on groundwater resources, and the volumes of the four aquifer groups were calculated in the hierarchical SIS model. The shallow groundwater in the study area is mainly concentrated in Hebi City and the Puyang basin of the Yellow River, and deep groundwater is mainly concentrated in northern Anyang City and Hebi City. The volume of possible shallow and deep groundwater recharge in the study area were also calculated as 4.653 × 10³ m³. The main possible recharge area of shallow groundwater is distributed in the Changyuan-Puyang area, and the main possible recharge areas of deep groundwater are distributed in northern Anyang City and southern Hebi City. The hierarchical hydrogeological model can be further used for groundwater flow modeling, vulnerability assessment, and environmental analysis, which would provide information to complement the development and management policy for sustainable water extraction from the aquifers. In the future, the static hydrogeological model reconstructed needs to integrate the hydraulic properties, boundary conditions, and groundwater calibration to dynamically simulate the groundwater flow and to analyze the groundwater sensitivity in the study area.