资源环境科技发展态势分析平台

Entropy-Guided Multivariate Groundwater Network Design for Multi-Source Data Assimilation Under Observational Uncertainty

https://doi.org/10.1029/2025GL117466
2025-09-27
查看原文
Geophysical Research Letters . Volume 52 , issue 19
Mingxu Cao, Zhenxue Dai, Junjun Chen

Abstract

Observational uncertainty poses major challenges to groundwater model calibration. As the primary source of information for multi-source data assimilation, monitoring network design is critical for accurately characterizing subsurface dynamics. Under limited measurement accuracy or cost constraints, monitoring networks must remain robust to observational errors. This study develops a multivariate network design framework that quantifies the uncertainty of multicomponent responses using joint entropy and employs deep learning to accelerate computations. Case study results show that the framework reliably estimates non-Gaussian permeability fields even under high-noise observations. The calibrated reactive transport model demonstrates strong capability in reproducing historical data and predicting system responses. This work advances the understanding of multi-source data fusion and supports the development of groundwater monitoring networks under observational uncertainty. Moreover, the proposed approach can be extended to the design of geophysical survey lines that integrate geophysical data.

Plain Language Summary

Groundwater monitoring is essential for understanding subsurface flow and contaminant transport. To obtain more reliable estimates of hydrogeological parameters, researchers often integrate diverse field observations with numerical models. However, observational uncertainty critically limits the reliability of groundwater modeling. We developed a monitoring network design framework grounded in the joint uncertainty of multivariate responses, employing an entropy-based data discretization strategy to guide the optimization algorithm in selecting monitoring stations. This framework enables robust estimation of non-Gaussian permeability fields under high-noise observations at optimal networks. Overall, the proposed approach provides a new perspective on groundwater monitoring and data utilization in the presence of observational uncertainty.

Key Points

  • A multivariate groundwater network design framework is developed by combining information entropy theories and deep learning methods

  • The developed framework can identify non-Gaussian permeability fields by integrating multi-source observations with high noise

  • The calibrated model honors both historical and predictive responses of water flow and multicomponent transport

1 Introduction

Accurate hydrological modeling is essential for predicting the spatiotemporal evolution of water flow and contaminant transport in aquifer systems (Sheng et al., 2024). Nevertheless, limited knowledge of system dynamics and the parameter uncertainty restrict the accurate simulation of real-world hydrological processes (M. Cao et al., 2024; Olorunsaye & Heiss, 2024). Data assimilation (DA) frameworks, which integrate observational data to update model states (Zhan et al., 2022), hydraulic parameters (C. Cao et al., 2024), as well as initial (Meyer et al., 2019) and boundary conditions (Chen et al., 2021), provide a powerful means of reducing uncertainty.

Within the Bayesian probabilistic framework, DA algorithms explicitly characterize discrepancies between model predictions and observations, enabling optimal integration and posterior uncertainty quantification. The deterministic inversion methods can produce a unique solution (Chen & Dai, 2024), while stochastic DA provides a comprehensive representation of uncertainty in both model states and parameters (Parrish et al., 2012). DA is therefore particularly well suited to complex hydrological systems characterized by substantial parameter uncertainty and limited data availability (Ghorbanidehno et al., 2020). The advancement of generalized DA methods for nonlinear, high-dimensional, and non-Gaussian problems critically depends on the quality of observational information and the performance of the inversion algorithms (J. Zhang et al., 2024). On the algorithmic front, substantial progress has been made, evolving from basic sampling approaches, such as Markov Chain Monte Carlo (Vrugt et al., 2008), to more efficient ensemble-based methods, including the ensemble Kalman filter (Xue & Zhang, 2014) and the ensemble smoother (Bao et al., 2020). These methods have been further enhanced through local update strategies to better address nonlinearity, such as the Iterative Local Update Smoother (Zhang et al., 2018), and through adaptive techniques for parameter estimation and uncertainty reduction (Ju et al., 2018).

Despite advances in DA methods, their utilization remains constrained by the observation data value (Zheng et al., 2018). Multi-source data fusion is a vital strategy for enhancing the diversity of response information (Gettelman et al., 2022). However, the acquisition of observational data is inherently subject to multiple sources of uncertainty, including instrument-related errors, temporal and spatial variability in measurements, environmental noise, and human-induced biases during data collection or preprocessing. Moreover, different data types often exhibit varying accuracies and sensitivities to system dynamics, leading to complex error structures during data fusion (da Silveira Barcellos & de Souza, 2022).

Optimizing monitoring networks is crucial for enhancing data utilization, particularly under observational uncertainty (Meggiorin et al., 2024). This requires defining an optimization criterion to quantify the information value of each station. Common sensitivity-based design criteria (Hsu & Yeh, 1989; Knopman et al., 1991; Sciortino et al., 2002) optimize monitoring stations iteratively by maximizing sensitivity and minimizing correlation. However, these criteria, which rely on the first two-order moments of parameter distributions, are not well suited for complex systems with multimodal parameter distributions. Entropy-based metrics provide a quantitative framework for evaluating the uncertainty associated with monitoring variables and the information transfer among stations (Shannon, 1948). Commonly used measures include marginal entropy, joint entropy, conditional entropy, transinformation, and total correlation. Keum et al. (2017) reviewed the application of entropy theory in designing monitoring networks for precipitation, streamflow, water level, soil moisture, groundwater, and water quality. Building on this foundation, multivariate network designs based on conditional entropy have been proposed by Keum and Coulibaly (2017), further demonstrating the potential of information-theoretic approaches to support integrated hydrological observations.

Incorporating optimized networks with hydrological models and DA frameworks is essential for validating their effectiveness and enhancing credibility (Ma et al., 2025). Careful consideration must also be given to data discretization strategies (Bosserelle & Hughes, 2024), objective functions (Foroozand & Weijs, 2020), temporal dynamics (Wang et al., 2018), and boundary constraints (Leach et al., 2016). Because entropy-based information estimation replaces probabilities with frequencies, the analysis requires data sets that accurately represent system uncertainty. Constructing such data sets in complex hydrological systems poses major challenges, including the “curse of dimensionality” and high computational costs. Recent advances in deep learning provide alternatives by reducing computational demands and enhancing practical applicability (Zhi et al., 2024). For instance, Generative Adversarial Networks (Ling & Jafarpour, 2024) and Denoising Diffusion Probabilistic Models (X. Zhang et al., 2024) enable the parameterization of high-dimensional model parameters into low-dimensional feature vectors, while deep learning-based surrogate models can efficiently replace CPU-intensive forward models (Zhan et al., 2025).

Our previous work (Cao, Dai, Chen, et al., 2025; Chen et al., 2022) has incorporated these recent advances in information entropy and deep learning into the monitoring network optimization framework, and reliable inversion results were achieved under ideal observations. However, as discussed earlier, measurement data often suffer from limited accuracy, or cost considerations necessitate higher tolerance to observational errors in practical applications. Therefore, the core objective of this study is to develop monitoring networks that can still ensure reliable estimation under high-noise observations. To this end, we developed a multivariate monitoring network design framework that integrates joint entropy with deep learning.

Accordingly, a multivariate, multi-objective optimization algorithm was established for optimal network design, in which the uncertainties of hydraulic heads, pH, and multicomponent concentrations were quantified through joint entropy. To improve computational efficiency, an integrated surrogate modeling framework was introduced, combining deep convolutional generative adversarial networks for reconstructing non-Gaussian permeability fields with deep convolutional residual networks for predicting multivariate system responses. A data discretization strategy grounded in information entropy was further employed to guide the optimization algorithm in monitoring station selection, thereby enhancing the robustness of multi-source DA under observational noise. The framework was validated with a coupled hydrogeochemical model, successfully identifying non-Gaussian permeability fields through assimilation of high-noise, multi-source observations. Overall, this study provides a novel theoretical foundation for integrating multi-source response data and designing multivariate groundwater networks with enhanced tolerance to high-noise observations, offering new insights for practical applications in complex groundwater modeling.

2 Methodology

This section introduces the methodology developed in this study, with the overall framework illustrated in Figure 1 and further details provided in Sections 2.1-2.3.

Details are in the caption following the image
Figure 1
Open in figure viewerPowerPoint

(a) Integrated surrogate modeling framework, combining a Deep Convolutional Generative Adversarial Network (DCGAN) for parameterization and a Deep Convolutional Residual Network (DCRN) as the surrogate model. (b) Schematic illustration of the multivariate groundwater network design method. The acronym MOO denotes multi-objective optimization. (c) Workflow diagram of network design strategies under high-noise observations.

2.1 Integrated Surrogate Modeling

This section introduces a framework of deep learning–based integrated surrogate modeling, as illustrated in Figure 1a. The framework couples a deep convolutional generative adversarial network (DCGAN) for the parameterization of non-Gaussian permeability fields with a deep convolutional residual network (DCRN) for the replacement of forward models.

DCGAN, proposed by Radford et al. (2015), consists of two adversarial networks: a generator G $G$ and a discriminator D $D$ . The generator G ( z ) $G(\mathrm{z})$ takes a noise vector z ∼ p z ( z ) $\mathrm{z}\mathit{\sim }{p}_{\mathrm{z}}(\mathrm{z})$ , sampled from a prior distribution, and reconstructs synthetic data resembling real data in the training set. Then, the discriminator D ( x ) $D(\mathrm{x})$ evaluates the generated data against real instances x ∼ p data ( x ) $\mathrm{x}\mathit{\sim }{p}_{\text{data}}(\mathrm{x})$ . The generator is trained to maximize the probability of the discriminator misclassifying its outputs as real. DCGAN improves training stability by employing strided convolutions in the discriminator and fractional-strided convolutions in the generator, making it particularly effective for image-like data. In this study, DCGAN is trained to capture complex spatial features and reconstruct non-Gaussian permeability fields from low-dimensional latent vectors, enabling efficient sampling. Furthermore, we employ DCRN to approximate the mapping of the forward model from permeability fields to system responses. The residual structure of the DCRN network mitigates the vanishing gradient issue in high-dimensional mappings. The training data set takes reconstructed non-Gaussian permeability fields as inputs and normalized multivariate response fields as outputs. Specifically, all response outputs from the numerical model are uniformly processed using min–max normalization (Chen et al., 2021), ensuring that different types of response data share the same scale. This facilitates surrogate model training and enables simultaneous multivariate prediction. The framework enables efficient sampling of non-Gaussian permeability fields and rapid prediction of multivariate dynamic responses, thereby substantially reducing the computational burden for multivariate network design and performance evaluation under uncertainty.

2.2 Multivariate Network Design

As shown in Figure 1b, let P = { X i j P ∣ i = 1 , 2 , … , M ; j = 1 , 2 , … , m } $\mathcal{P}=\left\{{X}_{ij}^{P}\mathit{\vert }i=1,2,\mathit{\text{\ldots }},M;j=1,2,\mathit{\text{\ldots }},m\right\}$ represent an initial hydrological monitoring network comprising M densely distributed preselected stations, and each recording m hydrological variables. X i P , i = 1 , 2 , … , M ${X}_{i}^{P},i=1,2,\mathit{\text{\ldots }},M$ , denote the multivariate random variables at station i. Unlike the marginal entropy in the univariate case, the uncertainty of X i P ${X}_{i}^{P}$ is quantified by their joint entropy, as follows:
H X i P = H X i 1 P , X i 2 P , … , X i m P = − ∑ r 1 = 1 n 1 ⋯ ∑ r m = 1 n m p r 1 , . . . , r m log 2 p r 1 , . . . , r m $H\,\left({X}_{i}^{P}\right)=H\,\left({X}_{i1}^{P},{X}_{i2}^{P},\mathit{\text{\ldots }},{X}_{im}^{P}\right)=-\sum\limits _{{r}_{1}=1}^{{n}_{1}}\mathit{\cdots }\sum\limits _{{r}_{m}=1}^{{n}_{m}}\,{p}_{{r}_{1},...,{r}_{m}}\,{\mathit{log}}_{2}\ {p}_{{r}_{1},...,{r}_{m}}$ (1)
where p r 1 , . . . , r m ${p}_{{r}_{1},...,{r}_{m}}$ denotes the joint probability of the m hydrological variables, estimated from empirical frequencies. To ensure computational feasibility, we adopted the merging algorithm (Alfonso et al., 2010), which reformulates multivariate entropy into lower-dimensional forms without information loss (e.g., H ( X , Y , Z ) = H ( 〈 X , Y 〉 , Z ) $H(X,Y,Z)=H(\langle X,Y\rangle ,Z)$ ).
The optimal monitoring networks are determined using a multi-objective optimization algorithm based on a multivariate MIMR criterion, which is formulated as:
max multivariate - MIMR = λ 1 H X 1 S , X 2 S , … , X K S + ∑ i U = 1 L U T X 1 S , X 2 S , … , X K S ; X i U − λ 2 C X 1 S , X 2 S , … , X K S $\underset{\text{multivariate}\mbox{-}\text{MIMR}}{\max }\,{=\lambda }_{1}H\,\left({X}_{1}^{S},{X}_{2}^{S},\mathit{\text{\ldots }},{X}_{K}^{S}\right)+\sum\limits _{{i}_{U}=1}^{{L}_{U}}T\,\left({X}_{1}^{S},{X}_{2}^{S},\mathit{\text{\ldots }},{X}_{K}^{S};{X}_{i}^{U}\right)-{\lambda }_{2}C\,\left({X}_{1}^{S},{X}_{2}^{S},\mathit{\text{\ldots }},{X}_{K}^{S}\right)$ (2)
where, K and L denote the numbers of selected and unselected monitoring stations, respectively. T X 1 S , X 2 S , … , X K S ; X i U = ∑ i 1 = 1 n 1 ⋯ ∑ i K = 1 n K ∑ i ′ = 1 n ′ p i 1 , . . . , i K ; p i ′ log 2 p i 1 , . . . , i K ; p i ′ p i 1 , . . . , i K p i ′ $T\,\left({X}_{1}^{S},{X}_{2}^{S},\mathit{\text{\ldots }},{X}_{K}^{S};{X}_{i}^{U}\right)=\sum\limits _{{i}_{1}=1}^{{n}_{1}}\mathit{\cdots }\sum\limits _{{i}_{K}=1}^{{n}_{K}}\,\sum\limits _{{i}^{\mathit{\prime }}=1}^{{n}^{\mathit{\prime }}}{p}_{{i}_{1},...,{i}_{K}};{p}_{i}^{\mathit{\prime }}\,{\mathit{log}}_{2}\frac{{p}_{{i}_{1},...,{i}_{K}};{p}_{i}^{\mathit{\prime }}}{{p}_{{i}_{1},...,{i}_{K}}{p}_{i}^{\mathit{\prime }}}$ defines the transinformation between the set of selected stations and each individual unselected station. C X 1 S , X 2 S , … , X K S = ∑ i S = 1 K S H X i S − H X 1 S , X 2 S , … , X K S $C\,\left({X}_{1}^{S},{X}_{2}^{S},\mathit{\text{\ldots }},{X}_{K}^{S}\right)=\sum\limits _{{i}_{S}=1}^{{K}_{S}}\,H\,\left({X}_{i}^{S}\right)-H\,\left({X}_{1}^{S},{X}_{2}^{S},\mathit{\text{\ldots }},{X}_{K}^{S}\right)$ quantifies redundancy among selected stations. The weighting parameters λ1 and λ2 balance the trade-off between effective information and redundant information, and are set to 0.8 and 0.2 following Li et al. (2012).

The multivariate groundwater monitoring network design procedure developed in this study is summarized as follows. First, specify the information ratio R s p ∈ ( 0 , 1 ) ${R}_{sp}\in (0,1)$ and initialize the selected network set S $\mathcal{S}$ . Then, compute the multivariate joint entropy by Equation 1 for each preselected station and add the one with the highest uncertainty to S $\mathcal{S}$ . The set S $\mathcal{S}$ is then expanded according to Equation 2. After each addition, the joint entropy H ( S ) $H(\mathcal{S})$ is updated, and the process terminates when H ( S ) ≥ R s p × H ( P ) $H(\mathcal{S})\mathit{\ge }{R}_{sp}\times H(\mathcal{P})$ . Subsequently, perform random sampling of the permeability field, repeat the above selection procedure, and incorporate additional stations into S $\mathcal{S}$ . The ensemble network is finalized when successive trials fail to add new stations.

2.3 Network Design Amid High-Noise Observations

2.3.1 Data Discretization Setups

The data discretization influences the probability distribution of model responses and thus determines the calculation of information entropy (Cellucci et al., 2005). In this study, multivariate random variables were normalized to the interval [0,1] and discretized according to the following expression:
x ′ = ⌊ x a ⌋ $\begin{array}{c}{\mathrm{x}}^{\prime }=\lfloor \frac{x}{a}\rfloor \end{array}$ (3)
where a $a$ denotes the bin width, finer discretization (i.e., a smaller a $a$ ) provides a closer approximation to the true data distribution, but limited sample sizes may introduce statistical fluctuations that compromise the accuracy of probability estimation.

With a sufficiently large sample size, applying finer discretization to continuous data enables more accurate quantification of data uncertainty. In contrast, coarser discretization is preferable for entropy estimation when samples are limited or data is noisy. In this study, three types of bin-width settings were defined: a $a$  = 0.1 for fine discretization, a $a$  = 0.25 for moderate discretization, and a $a$  = 0.5 for coarse discretization. These correspond to 10, 4, and 2 categorical divisions of the model responses, respectively, reflecting different levels of tolerance to data errors.

2.3.2 Multi-Source Data Assimilation

We employ the iterative local updating ensemble smoother (ILUES) developed by Zhang et al. (2018) to update parameters of interest by integration of multi-source observations, as illustrated in Figure 1c. Compared with the standard ensemble smoother, ILUES offers improved stability and convergence in nonlinear problems, making it more suitable for geological modeling and subsurface flow simulations under high uncertainty. The ILUES update equation is as follows:
p k n + 1 = p k n + K n d k ∗ − f G p k n $\begin{array}{c}{p}_{k}^{n+1}={p}_{k}^{n}+{K}^{n}\left({d}_{k}^{\ast }-f\left(G\left({p}_{k}^{n}\right)\right)\right)\end{array}$ (4)
where, the Kalman gain is computed as K n = C P M n C M M n + N iter · C D − 1 ${K}^{n}={C}_{PM}^{n}{\left({C}_{MM}^{n}+{N}_{\text{iter}}\mathit{\cdot }{C}_{D}\right)}^{-1}$ , with C P M n ${C}_{PM}^{n}$ and C M M n ${C}_{MM}^{n}$ denote the cross-covariance between parameters and model responses, and the auto-covariance of model responses, respectively. C D ${C}_{D}$ denotes the observation error covariance matrix. The perturbed observation is defined as d k ∗ = d + N iter · e k n ${d}_{k}^{\ast }=d+\sqrt{{N}_{\text{iter}}}\mathit{\cdot }{e}_{k}^{n}$ , with e k n ∼ 0 , C D ${e}_{k}^{n}\sim \left(0,{C}_{D}\right)$ representing the stochastic perturbation of the observation vector.

3 Numerical Example

The primary purpose of the modeling case is to demonstrate the feasibility of the developed multivariate groundwater monitoring network design framework under high observational noise. To balance computational efficiency, a two-dimensional multicomponent reactive transport model was employed to validate the framework, while maintaining enough complexity in the setup. It should be emphasized, however, that the proposed approach is generally applicable and can be readily adapted for three-dimensional models.

The model domain covers an 80 × 80  km 2 ${\text{km}}^{2}$ area discretized into a uniform 80 × 80 grid, and simulations are conducted with TOUGHREACT. Groundwater flow is governed by prescribed boundary conditions, including constant hydraulic heads on the left and right sides with a head difference of ΔH = 50 m, while the remaining boundaries are treated as no-flow. In addition to boundary conditions, the model incorporates the following assumptions. Groundwater flow is considered saturated and incompressible, following Darcy's law. The reactive transport processes include advection, dispersion, and chemical reactions, which were predefined based on previous studies, with reaction types and associated parameters treated as determined.

To represent non-Gaussian heterogeneity, the channelized field is generated from a training image (Figure 2a) using multiple-point statistics (MPS) simulations. The implementation follows the code provided by Zhang et al. (2020), while further details of the MPS approach can be found in Mariethoz et al. (2010). Figure 2b shows several random realizations of the channelized field. We followed the method proposed by Mo et al. (2020) to generate both the reference value and random realizations of the non-Gaussian permeability field. In TOUGREACT, the permeability of each grid was defined as k i = α i k ref ${k}_{i}={\alpha }_{i}{k}_{\text{ref}}$ , with k ref = 1.0 × 10 − 14 m 2 ${k}_{\text{ref}}=1.0\times {10}^{-14\,}{\mathrm{m}}^{2}$ in this study. The permeability modifier α i ${\alpha }_{i}$ of each facies follows a log-Gaussian distribution and is assigned distinct heterogeneous properties: ln α i 1 ∼ N 2 , 0.5 2 ${\mathit{ln}\left({\alpha }_{i}\right)}_{1}\sim \mathcal{N}\left(2,{0.5}^{2}\right)$ for low-permeability matrix, and ln α i 2 ∼ N 5 , 0.5 2 ${\mathit{ln}\left({\alpha }_{i}\right)}_{2}\sim \mathcal{N}\left(5,{0.5}^{2}\right)$ for high-permeability channels. Following the methodology outlined in Section 2.1, non-Gaussian permeability fields were reconstructed by a DCGAN model trained on a prior ensemble of 10,000 samples, implemented in PyTorch. Representative realizations are shown in Figure 2c. The performance of the generative model was assessed by comparing one-dimensional marginal probability density functions and two-dimensional scatter plots at randomly selected stations in the model domain with those derived from the training image (Figure S1 in Supporting Information S1). The reference permeability field used in this study is presented in Figure 2d, along with the uniformly distributed preselected monitoring stations across the model domain.

Details are in the caption following the image
Figure 2
Open in figure viewerPowerPoint

(a) Training image features channelized structures with a size of 250L × 250L. (b) The binary facies samples, sized 80L × 80L, capture the spatial distribution of high-permeability channels within a low-permeability matrix. (c) Non-Gaussian permeability fields generated using DCGAN. (d) The model domain and the spatial distribution of the permeability modifier α i ${\alpha }_{i}$ ( k i = α i × 1.0 × 10 − 14 m 2 ${k}_{i}={\alpha }_{i}\times 1.0\times {10}^{-14\,}{\mathrm{m}}^{2}$ ). Black dots indicate the preselected monitoring stations. (e) Multivariate responses of the reference field, including hydraulic heads, PH, and multicomponent concentrations (mol/L).

The reactive transport model is adapted from Dai and Samper (2006) and Chen et al. (2021), accounting for the main geochemical reactions in the aquifer system, including aqueous complexation, cation exchange, and calcite dissolution/precipitation. The reaction equations are listed in Table S1 in Supporting Information S1. The primary species in the model include H 2 O ${\mathrm{H}}_{2}\mathrm{O}$ , H + ${\mathrm{H}}^{+}$ , HC O 3 − $\text{HC}{\mathrm{O}}_{3}^{-}$ , S O 4 2 − $\mathrm{S}{\mathrm{O}}_{4}^{2-}$ , Ca 2 + ${\text{Ca}}^{2+}$ , Mg 2 + ${\text{Mg}}^{2+}$ , Na + ${\text{Na}}^{+}$ , K + ${\mathrm{K}}^{+}$ and XOH. Due to the manuscript length limitation, the associated model parameters, including cation exchange capacity, cation exchange coefficients, and the initial/boundary concentrations of the primary species, are provided in Table S2 in Supporting Information S1. TOUGHREACT calculates the concentrations of secondary species at each time step by solving coupled geochemical equilibrium/kinetic reactions and mass transport equations, where their concentrations are derived from the chemical speciation of primary species based on thermodynamics or reaction kinetics. Note that secondary species concentrations are intermediate variables only and not directly used in subsequent calculations. The multivariate dynamic responses included hydraulic heads, pH, and the concentrations of Ca2+, HCO3−, K+, Mg2+, and Na+, with the reference case simulations shown in Figure 2e.

The surrogate for the forward model was trained on a data set constructed from non-Gaussian permeability fields generated by DCGAN and multivariate responses simulated with TOUGHREACT. The performance of the surrogate model was evaluated by comparing its outputs with those of the numerical simulations. Details of the training parameters and error evaluation are provided in Text S1 in Supporting Information S1.

4 Results and Discussions

This section presents the results of the optimal monitoring networks and their application to multi-source DA under noisy observations for identifying non-Gaussian permeability fields. Section 4.1 investigates how data discretization influences the spatial distribution of monitoring stations. Section 4.2 assesses inversion results for non-Gaussian permeability fields across different noise levels, with emphasis on posterior estimates of historical and predictive responses under high-noise observations.

4.1 Optimal Monitoring Networks

Uniformly spaced preselected monitoring stations (361 in total, arranged in a 19 × 19 grid) were established across the modeling domain. Since information entropy is calculated by approximating probabilities through frequency counts, a sufficiently large sample data set of model responses is required to accurately quantify uncertainty. We evaluated the relationship between the information content at preselected stations and the data size (Figure S2 in Supporting Information S1). The results show that both the mean and variance of information entropy fluctuate significantly when the data size is small, but become stable when the sample size exceeds 10,000. Accordingly, a data size of N = 10 , 000 $N=10,000$ was adopted for subsequent analyses. For a uniform distribution, the theoretical upper limit of information entropy is log2 (N), while the actual entropy is determined by the specific data distribution and variable variances.

Figures S3–S5 in Supporting Information S1 present the stochastic simulation results of optimal single networks under the three data distributions. The results show that only 3–6 stations are needed to reach the 99% (the specific Rsp defined in this study) information coverage requirement. This is because the data information arises from the joint monitoring of multivariate variables, allowing each station to convey substantial information. Particularly, under finer discretization, each station transmits richer and more detailed information. Because single optimal networks derived from random simulations are inherently uncertain, our previous work (Cao, Dai, Chen, et al., 2025) introduced the concept of an ensemble monitoring network by aggregating all simulation outcomes. However, this approach is limited in that stations appearing only in a few simulations are often redundant and may impair parameter inversion. To address this, the present study proposes a practical termination criterion: the ensemble network is finalized once 6–8 consecutive trials fail to add any new stations. The range was determined empirically. Specifically, Figure S6 in Supporting Information S1 in the Supporting Information shows the number of newly added stations in each stochastic simulation, while Figure S7 in Supporting Information S1 illustrates the spatial distribution of the final optimal network and subsequently added stations. The results indicate that only a few stations were added in the later stages, and they were mainly located in non-critical areas. Accordingly, the final optimal networks comprise 36, 30, and 49 monitoring stations for the three data discretization schemes, respectively (Figure 3a). Spatially, the optimization based on a coarser discretization results in a larger number of monitoring stations with a more concentrated distribution. In contrast, the optimization based on a finer discretization yields fewer monitoring stations that are more widely dispersed.

Details are in the caption following the image
Figure 3
Open in figure viewerPowerPoint

(a) Optimal monitoring networks corresponding to three types of discretization bin widths: a = 0.1, 0.25, and 0.5. (b) Inversion results of non-Gaussian permeability fields using noisy observations with four noise levels (5%, 10%, 25%, and 50%) under the three optimal network scenarios. Black dots in the true field indicate the optimal monitoring stations. “Mean” denotes the ensemble mean of the posterior fields, while “Diff” shows the difference between the posterior mean and the reference field. (c) Historical data-matching results of Bayesian prior, posterior, and true responses at non-assimilated monitoring stations over the past 100,000 years (with 10,000-year time steps). Orange lines represent posterior responses under high-noise observations (50%) for the three optimal networks. The gray shaded area shows the prior uncertainty range, and the green line represents the true responses. (d) Comparison of predictive multivariate responses at 10,000 years under high-noise conditions (50%) across the three optimal network scenarios.

4.2 Inversion and Model Calibration Under Noisy Observations

In this case study, synthetic noise-free observations are first generated from the reference field, to which Gaussian noise ϵ ∼ N 0 , σ 2 ${\epsilon}\mathit{\sim }\mathcal{N}\left(0,{\sigma }^{2}\right)$ is added at four relative levels (5%, 10%, 25%, and 50%) to represent potential measurement uncertainties. In theory, a larger ensemble size for inversion provides a better representation of parameter uncertainty. However, computational efficiency must also be balanced against reliability. Based on sensitivity analysis, the ensemble size in this study was set to N e = 1000 ${N}_{e}=1000$ . In the case study, the inversion parameters essentially converged to the true values after the 5th iteration, therefore, the number of iterations was set to N i t e r = 10 ${N}_{\mathrm{i}\mathrm{t}\mathrm{e}\mathrm{r}}=10$ . The local smoother factor was set to α = 0.1 $\alpha =0.1$ , following the findings of Zhang et al. (2018).

Figure 3b presents the inversion results of the non-Gaussian permeability field under different levels of noisy observations, including the posterior ensemble means and their deviations from the true field for three optimal monitoring network scenarios. For the same data discretization scenario, increasing observational noise consistently reduces inversion reliability, whereas optimal monitoring networks designed with a larger discretization bin width enhance the tolerance of parameter inversion to noisy observations. At a fine discretization scenario (a = 0.1), the estimated accuracy is sensitive to higher noise, whereas it exhibits improved noise tolerance at a moderate discretization scenario (a = 0.25). It is noteworthy that under the high-noise scenario with a 50% observation error, the monitoring network designed with a discretization bin width of a = 0.5 yields the most reliable posterior estimates. This is because, with a sufficiently large sample size, finer discretization of the response data at a monitoring station can capture more detailed system information. Nevertheless, it should be emphasized that all optimal monitoring networks in this study provide the same total amount of information, which approaches the theoretical upper limit determined by the sample data set. Thus, monitoring networks derived from fine discretization consist of fewer stations and therefore require high-precision observations to ensure reliable inversion. In contrast, networks designed with coarse discretization involve a larger number of stations, which are often distributed in relatively concentrated areas. Such a setting enhances tolerance against observational noise, allowing reliable inversion results even when measurement accuracy is limited. This advantage provides a practical solution for real-world applications, where observational data are often subject to considerable uncertainty.

To further evaluate the posterior estimates under high observational noise (50%), we compared the calibrated model performance. Specifically, we reconstructed the long-term groundwater system evolution by fitting the spatial distributions of hydraulic heads, pH, and multicomponent concentrations over the past 100,000 years with 10,000-year time steps. A systematic comparison of the breakthrough curves for Bayesian prior, posterior, and reference responses was further carried out at monitoring stations that were not included in the DA. For clarity, Figure 3c illustrates the data-matching results at two randomly selected stations. The results indicate that, overall, the proposed inversion framework effectively reduces the uncertainty of posterior ensembles of model responses even under high observational noise. In contrast, cases with lower reliability in the posterior permeability fields exhibit larger historical fitting errors in the calibrated models. The results clearly indicate that when monitoring data accuracy is limited, optimizing the monitoring network with a larger discretization bin width can maximize the utility of available observations and thereby improve modeling accuracy. Calibrated models were further evaluated by comparing predicted hydraulic head, pH, and multi-component concentration fields over the next 10,000 years against the reference system responses. As shown in Figure 3d, the ensemble means of model responses accurately capture both high-permeability channels and low-permeability matrices across all scenarios. Examination of deviations from the reference fields further confirms the earlier conclusion that, under high-noise conditions, optimal monitoring networks designed with a larger discretization bin width can more effectively integrate observational information, thereby improving parameter inversion and model accuracy.

5 Summary and Conclusions

This study develops an entropy-based framework for multivariate groundwater monitoring network design. Under limited or low-precision observational data, the optimal network derived from this framework exhibits greater tolerance to observational errors. The framework is therefore advantageous for practical applications of network optimization and multi-source DA: with high-quality data, the number of monitoring stations can be reduced to lower costs without compromising inversion accuracy, whereas under restricted monitoring conditions, it still ensures reliable inversion estimation. Overall, the proposed approach thus offers valuable guidance for the strategies of multivariate data integration and the design of optimal monitoring networks amid high-noise observations.

Despite these contributions, this framework has some limitations for applications, including the reliance of the model on prior knowledge of the aquifers, the use of borehole data that may be sparse in practice, and the validation of the case study under the assumption of a deterministic process model. Nonetheless, the proposed approach and key findings are broadly applicable and can be extended to situations with uncertainty in geostatistical parameters or transport processes. Moreover, when borehole data are limited, geophysical data can be integrated within this framework to guide survey grid design.

Acknowledgments

This work was funded by the National Natural Science Foundation of China (NSFC: U2267217, 42402241, 42141011) and Shandong Key Water Conservancy Science and Technology Project (2024370203001957).

    Conflict of Interest

    The authors declare no conflicts of interest relevant to this study.

    Data Availability Statement

    Data are available at (Cao, Dai, & Chen, 2025).

      相关资源