Enhancing Coastal Aquifer Characterization and Contamination Inversion with Deep Learning

Abstract

Coastal aquifers are critical freshwater resources that face increasing threats from contamination and saltwater intrusion. Traditional approaches for characterizing these aquifers are challenged by complex dynamics, high-dimensional parameter spaces, and significant computational demands. This study presents an innovative method that combines an Auto-Regressive Convolutional Neural Network (AR-CNN) surrogate model with the Iterative Local Updating Ensemble Smoother (ILUES) for the joint inversion of contamination source parameters and hydraulic conductivity fields. The AR-CNN surrogate model, trained on synthetic data generated by the SEAWAT model, effectively approximates the complex input–output relationships of coastal aquifer systems, substantially reducing computational burden. The ILUES framework utilizes observational data to iteratively update model parameters. A case study involving a heterogeneous coastal aquifer with multipoint pollution sources demonstrates the efficacy of the proposed method. The results indicate that AR-CNN-ILUES successfully estimates pollution source strengths and characterizes the hydraulic conductivity field, although some limitations are observed in areas with sparse monitoring points and complex geological structures. Compared to the traditional SEAWAT-ILUES framework, the AR-CNN-ILUES approach reduces the total inversion time from approximately 70.4 h to 16.2 h, improving computational efficiency by about 77%. These findings highlight the potential of the AR-CNN-ILUES framework as a promising tool for efficient and accurate characterization of coastal aquifers. By enhancing computational efficiency without significantly compromising accuracy, this method offers a viable solution for the sustainable management and protection of coastal groundwater resources.

1. Introduction

Coastal aquifers play a crucial role as essential freshwater sources. However, they are extremely susceptible to contamination resulting from both natural and human-induced factors. The complex dynamics of saltwater intrusion, coupled with the heterogeneous properties of aquifers, present substantial challenges for precise characterization and efficient pollution tracing [1,2,3]. Coastal areas are particularly vulnerable due to interactions between freshwater and saltwater systems, influenced by factors such as tides, seasonal water demand, and climate change. These interactions can lead to saltwater intrusion, significantly impacting freshwater availability and quality—a critical issue for communities relying on these resources for drinking water, agriculture, and industry [3,4].

Effective management of coastal aquifers necessitates detailed spatial information regarding hydraulic conductivity and contamination source characteristics. Conventionally, numerical models like SEAWAT are utilized for this objective, linking groundwater flow with solute transport processes [5,6]. While powerful, these models can be computationally intensive, especially when dealing with inverse problems that require multiple iterations for parameter estimation. Such computational demands often hinder real-time monitoring and management, particularly in large-scale or highly dynamic systems [7,8].

Inverse problems in groundwater modeling involve estimating aquifer properties, such as hydraulic conductivity fields, and identifying contamination sources by utilizing observational data [9]. Data assimilation techniques are frequently employed to iteratively adjust model parameters, minimizing the differences between observed and simulated states [10,11]. However, repeatedly running complex numerical models like SEAWAT to address these inversion problems creates a significant computational bottleneck, particularly in coastal aquifers where density variations and complex boundary conditions further add to the complexity [12,13,14].

Recent advancements in machine learning, particularly deep learning, have provided researchers with the means to tackle these computational challenges in subsurface modeling [15,16,17]. Deep learning methods have the ability to capture intricate, nonlinear relationships within the data, making them highly suitable for either replacing or supplementing traditional numerical simulations. Convolutional neural networks (CNNs), a category of deep learning models commonly utilized in image recognition, have demonstrated potential in approximating the behavior of complex physical systems. By leveraging their capacity to learn spatial features from training data, CNNs can serve as surrogate models that approximate the output of numerical simulations at a significantly reduced computational cost.

Several studies have successfully employed CNNs and other machine learning techniques to develop surrogate models for groundwater applications [18,19,20,21,22,23]. For example, Jardani et al. (2022) utilized CNNs with an encoder–decoder structure to develop a surrogate model for hydraulic tomography applications, achieving comparable accuracy to traditional models while significantly reducing computational times [20]. Zhou et al. (2022) demonstrated the use of CNN-based surrogate models for the simultaneous identification of hydraulic and transport properties, offering a viable alternative to conventional methods that require extensive computational resources [19]. Zheng et al. (2023) integrated Generative Adversarial Network (GAN) and CNN techniques within an ILUES framework for the efficient estimation of contaminant source and hydraulic conductivity, showing that the approach could effectively reduce computational costs while maintaining accuracy [18]. Similarly, Mo et al. (2020) applied adversarial autoencoders with residual dense networks to estimate non-Gaussian hydraulic conductivities, highlighting the versatility of deep learning approaches in addressing complex subsurface characterization problems [21]. These successes in inland groundwater systems suggest that similar approaches could be adapted for more complex coastal aquifer environments, where density-driven flow and salinity gradients complicate model requirements.

In this study, we aim to develop a deep learning-based surrogate model to replace the forward simulation model (SEAWAT) within the iterative inversion process, thereby making the characterization of coastal aquifers more computationally efficient. We propose a novel integration of a CNN-based surrogate model with the ILUES to accomplish the joint inversion of contamination source parameters and hydraulic conductivity fields. The ILUES, a highly effective data assimilation technique, enables efficient parameter updating, especially for handling sparse and noisy observational data, which are typical challenges in aquifer characterization [11,24,25,26]. Utilizing a deep convolutional encoder–decoder architecture, our surrogate model closely approximates the input–output relationships of the SEAWAT model, delivering rapid predictions with acceptable accuracy. This integration not only significantly speeds up the inverse modeling process but also broadens the application scope of data assimilation techniques in real-time scenarios.

Our innovative approach holds particular significance for coastal aquifer systems that exhibit heterogeneous hydraulic conductivity, dynamic contamination plumes, and saltwater intrusion effects. Conventional methods often encounter difficulties in dealing with the high-dimensional and nonlinear nature of these challenges, frequently resulting in heavy computational burdens that hinder timely decision-making. The incorporation of CNN-based surrogate modeling within the ILUES framework offers a highly computationally efficient alternative, effectively addressing these issues and providing an innovative avenue for the efficient management of coastal groundwater resources.

2. Methodology

2.1. CNN Surrogate Model and ILUES-Based Data Assimilation Framework

In this proposed methodology, a deep convolutional neural network (CNN) surrogate model is integrated with the ILUES to jointly identify the contamination source parameters and the hydraulic conductivity field within a coastal aquifer system. The CNN is trained to act as a surrogate for the SEAWAT model, which is conventionally used to simulate the coupled groundwater flow and solute transport. The surrogate model is designed to alleviate the computational burden associated with the repeated simulations required for data assimilation.

The ILUES is an improved algorithm based on the Ensemble Smoother (ES), employing a Monte Carlo approximation of the Kalman filter [11]. This framework is structured to iteratively update both the model state and parameter ensembles using observational data. In this methodology, the traditional forward model (e.g., SEAWAT) is replaced by a CNN surrogate model, thereby facilitating a more efficient inversion process. The methodology encompasses the following key steps:

1. Initialization: Generate an ensemble of initial approximations for the hydraulic conductivity field and contamination source parameters, considering the inherent uncertainty in the aquifer properties and contamination scenario. Each ensemble member corresponds to a potential state of the system.

2. Prediction Step: At each iteration ($\mathrm{t}$), use the CNN surrogate model to forecast the state variables (such as hydraulic head and contaminant concentrations) for each ensemble member. Our surrogate model possesses an auto-regressive structure, utilizing the state at time ($\mathrm{t}$) as an input to predict the state at time ($\mathrm{t}+1$), thereby capturing the temporal dynamics.

3. Update Step: Incorporate observations at time ($\mathrm{t}$) into the model to update the ensemble. Compute the Kalman gain to adjust the ensemble predictions based on the differences between the predicted and observed states. Update the ensemble members using the calculated gain factor.

4. Constraint Handling: Enforce constraints during each update step to prevent parameters from deviating to unrealistic values. If any ensemble realization violates physical constraints (e.g., negative hydraulic conductivity), project it back to the boundary of the feasible space.

5. Iteration and Convergence: Iteratively execute the prediction and update steps until the convergence criteria are met, typically characterized by the stabilization of parameter estimates, a reduction in the ensemble spread, or reaching a predefined number of iterations.

2.2. AR-CNN-Based Surrogate Model

We adopted a deep convolutional encoder–decoder network proposed by Mo et al. (2019) [27], specifically AR-CNN, to develop a surrogate model for the SEAWAT numerical simulation. This network has demonstrated excellent reconstruction performance in high-dimensional pollution source inversion problems. A brief description is as follows.

The architecture of this network relies on two key components.

• Dense Block (depicted in Figure 1): Designed to extract essential features of the input data, specifically the information of the aquifer and contamination source. In a Dense Block containing $\mathrm{L}$ layers, each layer receives the concatenated outputs of all preceding layers as its input. This unique structure effectively reduces information loss and substantially enhances the training efficiency of the neural network. Each layer in the Dense Block performs three sequential operations: Batch Normalization (BN), ReLU activation, and convolution.

Figure 1

Schematic diagram of Dense Block structure (L = 3).

[Figure omitted. See PDF]

• Transition Layer (illustrated in Figure 2): Serves to smoothly connect two Dense Blocks and plays different roles in the encoding and decoding processes. In the encoding phase, the Transition Layer acts as an encoding layer, halving the number of feature maps and reducing the spatial dimensions through a combination of convolution and pooling operations. Conversely, in the decoding phase, it functions as a decoding layer, doubling the number of feature maps and increasing the spatial dimensions back to the original size using transposed convolution operations. This dual functionality is vital for ensuring the accurate reconstruction of output fields.

Figure 2

Schematic diagram of encoding and decoding layers.

[Figure omitted. See PDF]

The complete architecture of the AR-CNN-based surrogate model is illustrated in Figure 3. The inputs to the model consist of multipoint contamination source information, hydraulic conductivity fields, and the previous timestep’s pollution plume distribution. The model’s outputs include the current timestep pollutant concentration field, salinity, and hydraulic head.

To effectively capture the temporal dynamics of the coastal groundwater, an auto-regressive strategy was employed. This strategy enables the output at the current timestep to depend not only on the current inputs but also on the output from the previous timestep. Specifically, the input $\left(K,S_j,y_{j-1}\right)$ consists of the hydraulic conductivity field K; the source strength $S_j$ released at the time j; and the hydraulic head, contaminant concentration, and salinity at the time j − 1, denoted as $y_{j-1}=(C_S,C_C,{H)}_{j-1}$; the output $y_j$ is the hydraulic head, contaminant concentration, and salinity at the time j. Mathematically, this can be represented as:

(1)$y_j=\mathrm{f}\left(K,S,y_{j-1}\right)$

This approach allows the model to better handle the time-dependent nature of the groundwater processes.

During the training process, we employed the L1 loss function, given by:

(2)$ℒ\left(\widehat{y},y\right)={\sum }_{i=1}^n\left|\widehat{y_i}-y_i\right|$

The L1 loss function was chosen due to its robustness to outliers, which is particularly important given the inherent noise and variability in the groundwater observational data. The Adam optimizer was used to update the network weights.

By designing and training the AR-CNN-based surrogate model, we aimed to develop an accurate and computationally efficient approximation of the SEAWAT model. This surrogate model was integrated with the ILUES-based data assimilation framework to facilitate the inversion of coastal aquifer parameters in a computationally feasible manner.

2.3. ILUES Method

The data assimilation framework based on the ILUES was adopted as the inversion method for jointly identifying the contamination source parameters and the hydraulic conductivity field in coastal aquifers. A brief description of the ILUES method is provided herein. The relationship between the uncertain model parameters ($\mathrm{m}$) and the measurements ($\mathrm{d}$) is expressed as:

(3)$d=\mathrm{f}\left(m\right)+\epsilon$

where $\mathrm{f}\left(m\right)$ is the forward model (SEAWAT model, or AR-CNN-based surrogate model), and $\epsilon$ is the measurement error vector.

During each iteration of the ILUES process, the forward model plays a pivotal role. It takes in the ensemble members that represent the current state of the hydraulic conductivity and contamination source information, and predicts the outputs, which encompass the hydraulic head, contaminant concentration, and salinity. Subsequently, the ensemble is updated based on the discrepancies between the predicted outputs and the observed data. The update formula for the local set $\mathrm{J}\left(\mathrm{m}\right)$ consisting of ${\mathrm{N}}_{\mathrm{e}}$ samples of the uncertain groundwater model parameters ($\mathrm{m}$) is given by:

(4)$\mathrm{J}\left(m\right)={\mathrm{J}}_1\left(m\right)/{\mathrm{J}}_1^{\max }+{\mathrm{J}}_2\left(m\right)/{\mathrm{J}}_2^{\max }$

Here, ${\mathrm{J}}_1\left(m\right)={\left[\mathrm{f}\left(m\right)-d\right]}^{\mathrm{T}}{\mathrm{C}}_{\mathrm{D}}^{-1}\left[\mathrm{f}\left(m\right)-d\right]$ measures the distance between the model predictions $\mathrm{f}\left(m\right)$ and the observations $d$, and ${\mathrm{J}}_2\left(m\right)={\left[m-m_f\right]}^{\mathrm{T}}{\mathrm{C}}_{\mathrm{MM}}^{-1}\left[m-m_f\right]$ represents the distance between the model parameter $\mathrm{m}$ and the corresponding samples ${\mathrm{m}}_{\mathrm{f}}$. ${\mathrm{C}}_{\mathrm{D}}$ is the covariance matrix of observation errors, and ${\mathrm{C}}_{\mathrm{MM}}$ is the autocovariance matrix of model parameters. ${\mathrm{J}}_1^{\max }$ and ${\mathrm{J}}_2^{\max }$ are the maximum values of ${\mathrm{J}}_1\left(m\right)$ and ${\mathrm{J}}_2\left(m\right)$, respectively.

After calculating $\mathrm{J}\left(m\right)$ for each sample, the ILUES selects the $\mathrm{N}$ smallest $\mathrm{J}\left(m\right)$ values, and the corresponding sample parameter ${\mathrm{m}}_{\mathrm{f}}$ is defined as the local sample $m_{j,f}$ of the parameter. This parameter is updated using the following formula:

(5)$m_{j,i}^a=m_{j,i}^f+{\mathrm{C}}_{\mathrm{MD}}^{\mathrm{lf}}{\left({\mathrm{C}}_{\mathrm{DD}}^{\mathrm{lf}}+{\mathrm{C}}_{\mathrm{D}}\right)}^{-1}\left[d_j-\mathrm{f}\left(m_{j,i}^f\right)\right]$

where ${\mathrm{C}}_{\mathrm{MD}}^{\mathrm{lf}}$ is the cross-covariance matrix of $M_{j,f}^l=\left[m_{j,1}^f,m_{j,2}^f,\cdots ,m_{j,N}^f\right]$ and $D_{j,f}^l=\left[\mathrm{f}\left(m_{j,1}^f\right),\mathrm{f}\left(m_{j,2}^f\right),\cdots ,\mathrm{f}\left(m_{j,N}^f\right)\right],{\mathrm{C}}_{\mathrm{DD}}^{\mathrm{lf}}$ is the autocovariance matrix of $D_{j,f}^l$, and $d_j=d+\epsilon$ is the observed true value $\mathrm{d}$ after introducing errors $\epsilon$.

By using the surrogate model instead of the traditional SEAWAT forward model, significant computational advantages are achieved when solving inverse problems using the ILUES method. It allows for a larger ensemble size, thereby providing a more comprehensive representation of the possible states of the system. This is especially important for coastal aquifers, where the complex dynamics influenced by factors like tidal fluctuations necessitate a more detailed exploration of the parameter space. Furthermore, the faster prediction speed of the surrogate model allows for more iterations within a reasonable computational time, resulting in more accurate and stable results.

2.4. SEAWAT-ILUES and AR-CNN-ILUES Inversion Frameworks

We implemented two inversion frameworks—SEAWAT-ILUES and AR-CNN-ILUES—to address a coastal aquifer scenario characterized by multipoint contamination sources, geological heterogeneity, and observational uncertainties.

2.4.1. SEAWAT-ILUES Framework

The SEAWAT-ILUES framework employs the traditional SEAWAT model for simulating coupled groundwater flow and solute transport. SEAWAT, although robust and highly accurate, is computationally intensive, particularly for inversion tasks where multiple iterations are required for parameter updates. Within this framework, the ILUES assimilates observational data, such as hydraulic head, contaminant concentration, and salinity. It initiates by representing the hydraulic conductivity field through Karhunen–Loève (KL) expansion, followed by random sampling of the source prior distribution. The SEAWAT model forecasts the outputs, and the ILUES algorithm iteratively updates the ensemble based on observed values. This framework is effective when high accuracy is crucial but incurs significant computational costs due to the repeated forward model runs.

2.4.2. AR-CNN-ILUES Framework

The AR-CNN-ILUES framework tackles the computational limitations of SEAWAT-ILUES by substituting the SEAWAT forward model with a Convolutional Neural Network (AR-CNN) surrogate. Trained on synthetic data from SEAWAT, AR-CNN approximates the simulation output at lower computational cost, enabling faster ensemble updates for inversion. In this framework, the hydraulic conductivity field (represented by KL expansion coefficients) and contamination source parameters are fed into the AR-CNN surrogate model. The ILUES algorithm assimilates observational data to update hydraulic conductivity and contamination source parameters.

In both the SEAWAT-ILUES and AR-CNN-ILUES frameworks, the inversion of the hydraulic conductivity field is achieved through the KL expansion coefficients, which serve as a dimensionality-reduced representation of the field through Karhunen–Loève (KL) expansion. The KL expansion is used to parameterize the log-transformed hydraulic conductivity field $\mathrm{Y}\left(x\right)$, expressed as:

(6)$Y\left(x\right)\approx \overline{Y}\left(x\right)+{\sum }_{i=1}^{N_{KL}}\sqrt{{\mathsf{\tau }}_i}{\mathsf{\xi }}_i{s}_i\left(x\right)$

where $\overline{Y}\left(x\right)$ represents the mean component of the log-conductivity field, ${\tau }_i$ and $s_i\left(x\right)$ are the eigenvalues and eigenfunctions of the correlation function $C\left(x, x^{\prime }\right)$, and ${\xi }_i$ are independent standard Gaussian random variables. The number of terms $N_{KL}$ retained in the expansion determines the level of approximation. This method allows for an efficient reduction in dimensionality while preserving the dominant variability in the hydraulic conductivity field [28].

Both frameworks employ the ILUES algorithm for parameter iteration but differ in computational approaches. SEAWAT-ILUES provides detailed physical process simulations, whereas AR-CNN-ILUES uses machine learning for efficient approximation. Figure 4 illustrates the structural differences between the two frameworks. The AR-CNN-ILUES approach bypasses repeated SEAWAT runs by relying on the surrogate model for rapid predictions, thereby reducing computational complexity and enabling frequent updates—vital for inversion in dynamic coastal aquifers.

3. Case Study

3.1. Problem Description

To demonstrate the applicability of the proposed inversion framework introduced in Section 2, we utilized a two-dimensional model of groundwater flow through an isotropic, heterogeneous confined aquifer, as illustrated in Figure 5a. The groundwater flow is density-driven and transient, with a conservative pollutant whose migration is affected by the complex velocity field induced by density variations.

3.1.1. Scale Model Description

It is important to note that this study uses a scaled-down representation of a real coastal aquifer. Similar scaling approaches have been adopted in previous studies [7,29]. In this case, the following key simplifying assumptions were made:

• Unsaturated Flow: The model does not consider unsaturated flow. As discussed by Xin et al. [30], the impact of unsaturated flow on groundwater flow and solute transport is relatively small, except near the sediment–water interface.

• Non-Periodic Terms in Tidal and Wave Effects: Non-periodic terms in tidal and wave fluctuations (e.g., wind-driven waves and sea level variations) are neglected.

• Conservative Contaminants: The model assumes that the contaminants are conservative, ignoring adsorption and degradation processes.

• Lateral Flow Sources: Lateral flow sources, such as pumping wells or injection wells, are not considered. Therefore, a two-dimensional model is sufficient to describe the flow field and transport pathways in the coastal aquifer.

These assumptions allow the model to focus on the primary dynamics of seawater intrusion and contaminant transport, while maintaining computational efficiency and simplicity.

3.1.2. Model Domain and Boundary Conditions

The study area measures $80 \mathrm{cm}\times 20 \mathrm{cm}$ and was discretized using a $1 \mathrm{cm}\times 1 \mathrm{cm}$ finite-difference grid, resulting in $1600$ grid cells. The aquifer thickness was $20 \mathrm{cm}$. Note that this model represents a scaled-down version of a real coastal aquifer. Figure 5b shows the layout of the monitoring wells (W1–W10) and observation points.

At the left boundary of the aquifer (A–D), a constant inland hydraulic head boundary was applied, with the boundary head set at $22 \mathrm{cm}$. The top (D–C) and bottom (A–B) boundaries were no-flow boundaries. The seaward boundary (B–C) was subject to tidal influences, modeled as a time-varying sinusoidal tidal fluctuation:

(7)$H_{Tide}=H_{MSL}+A\sin \left({\displaystyle \frac{2\mathsf{\pi }}{T}}t\right)$

where $H_{Tide}\left[\mathrm{cm}\right]$ is the time-varying seawater level, $H_{MSL}\left[\mathrm{cm}\right]$ is the mean sea level (MSL), $A\left[\mathrm{cm}\right]$ is the tidal amplitude (set to $0.5 \mathrm{cm}$), $T\left[\min \right]$ is the tidal period (set to $20 \min$), and $t\left[\min \right]$ is the simulation time.

3.1.3. Hydraulic Conductivity Field

The hydraulic conductivity field was assumed to follow a log-normal distribution with mean ($\overline{\ln K}=4.0$) and variance (${\mathsf{\sigma }}_{\mathrm{lnK}}^2=1.0$). The correlation lengths in the x and z directions were 30 cm and 15 cm, respectively. The variogram type was exponential. The reference hydraulic conductivity field is shown in Figure 6.

In practical engineering applications, fully capturing the geological characteristics of a site is often challenging. Key parameters such as the mean, variance, and correlation lengths of the hydraulic conductivity field are typically only coarsely estimated based on limited observation well data. This introduces significant uncertainties that can affect the accuracy of source identification in coastal groundwater pollution studies. In this study, the initial survey provided approximate estimates: a hydraulic conductivity mean of 3.90 with a 2.5% error, a variance of 1.22 with a 22% error, and correlation lengths of ${\mathsf{\lambda }}_{\mathrm{x}}=25 \mathrm{cm}$ and ${\mathsf{\lambda }}_{\mathrm{z}}=12.5 \mathrm{cm}$, each with a 16.7% error. These uncertainties necessitate careful consideration of their impact on the characterization of coastal aquifer contamination and hydraulic conductivity field inversion. Table 1 lists the other main parameters of the coastal groundwater numerical model.

3.1.4. Pollution Sources and Simulation Parameters

The total simulation duration was 100 min, divided into 10 stress periods. In practical groundwater pollution scenarios, while the suspected locations of pollution sources can be identified through preliminary field investigations, the temporal release history of contaminants from these sources often remains unknown. In this study, three pollution sources were considered: $S_1$ and $S_2$ were active sources continuously releasing pollutants, while $S_3$ was inactive and did not release any contaminants. The sources were located at $S_1\left(12,19\right)$, $S_2\left(27,18\right)$, and $S_3\left(17,8\right)$. To characterize the release intensity of each source across the 10 stress periods, a total of 30 source parameters were introduced, denoted as $S_s=\left\{S_{si}^j\right\}$, where $\left(i=1, 2, 3\right)$ represents the source index, and $\left(j=1, 2, \dots , 10\right)$ denotes the stress period. The prior distribution of these parameters was assumed to be uniform, with the prior range and reference values detailed in Table 2.

In this study, 100 Karhunen–Loève (KL) terms were selected to characterize the heterogeneous coastal aquifer, capturing approximately 95.08% of the variance in the hydraulic conductivity field. Consequently, a total of 130 parameters were estimated for groundwater pollution source identification, comprising the release intensity of three pollution sources across 10 stress periods and the 100 KL terms.

To achieve accurate identification of both the pollution sources and the hydraulic conductivity field, observational data were collected, including salinity, pollutant concentration, and hydraulic head measurements at 70 observation points over 10 time intervals. This resulted in a dataset of 2100 observations: 700 for seawater intrusion salinity, 700 for pollutant concentration, and 700 for transient hydraulic head. The observational data were subject to a noise level of ${\mathrm{L}}_{\mathrm{noi}}=5\%$.

3.2. Assessment Criteria

• (1). Surrogate Model Evaluation Metrics

To evaluate the predictive accuracy of the deep convolutional network surrogate model, two evaluation metrics were adopted: the coefficient of determination ($R^2$) and the root mean square error (RMSE). A higher $R^2$ value indicates that the surrogate model’s predictions closely match the observed values. Meanwhile, a lower RMSE signifies smaller deviations between the predicted and observed values, reflecting the higher predictive accuracy of the surrogate model.

The coefficient of determination is calculated as:

(8)$R^2=1-{\displaystyle \frac{{\sum }_{i=1}^{N_{e_{test}}}{\left(y_i-\widehat{y_i}\right)}^2}{{\sum }_{i=1}^{N_{e_{test}}}{\left(y_i-\overline{y}\right)}^2}}$

where $\overline{y}$ is the mean of all outputs in the test dataset, calculated as:

(9)$\overline{y}={\displaystyle \frac{1}{N_{e_{test}}}}{\sum }_{i=1}^{N_{e_{test}}}{y}_i$

The root mean square error is calculated as:

(10)$RMSE=\sqrt{{\displaystyle \frac{{\sum }_{i=1}^{N_{e_{test}}}{\left(y_i-\widehat{y_i}\right)}^2}{N_{e_{test}}}}}$

In these equations, $N_{e_{test}}$ is the number of test samples, y represents the outputs from the original model (SEAWAT code), and $\widehat{y}$ denotes the corresponding outputs from the surrogate model.

• (2). Parameter Inversion Evaluation Metrics

To quantitatively assess the performance of the ILUES inversion framework in groundwater pollution source identification, two evaluation metrics were employed: Mean Absolute Error (MAE) and Average Ensemble Spread (AES). MAE quantifies the deviation between the estimated values of unknown parameters and their true values. A lower MAE indicates higher accuracy in parameter inversion, thus reflecting the effectiveness of the inversion method. The MAE is calculated as follows:

(11)$MAE={\displaystyle \frac{1}{N_x}}{\sum }_{i=1}^{N_x}\left|{\left(V_{est}\right)}_i-{\left(V_{true}\right)}_i\right|$

AES, as an additional evaluation metric, measures the dispersion of the parameter ensemble relative to their estimated means. Specifically, a smaller AES value indicates that the samples in the ensemble are closer to the mean, exhibiting lower dispersion. This metric reflects the degree of uncertainty in the inversion results and is crucial for evaluating their stability and reliability. Effectively reducing the AES value can further optimize the inversion algorithm. The AES is calculated as:

(12)$AES=\sqrt{{\displaystyle \frac{1}{N_x}} [{{\displaystyle \sum }}_{i=1}^{N_x}\left\{{{\displaystyle \sum }}_{j=1}^{N_e} [{(V_{est})}_{ij}-{\left(V_{est}{)}_i\right]}^2\right\}}]$

In these equations, ${(V_{true})}_i$ and ${(V_{est})}_i$ represent the true value and the estimated mean of the i-th parameter, respectively; ${(V_{est})}_{ij}$ is the estimated value in the j-th ensemble sample for the i-th parameter; $N_x$ is the total number of parameters to be estimated; and $N_e$ is the number of ensemble samples.

3.3. Network Structure Design

Table 3 presents the architecture of the AR-CNN utilized in this study. The network begins with an initial convolutional layer, followed by three Dense Blocks with internal layer counts L set to [5, 10, 5], respectively. Between each pair of Dense Blocks, encoding and decoding layers are employed to facilitate smooth transitions. The final layer is a decoding layer that outputs images with dimensions $N_f\times H_f\times W_f$ = 3$\times 20\times 80$.

To ensure that the predicted concentration and hydraulic head values were non-negative, the softplus activation function was applied in the last layer of the network. During training, the L1 loss function was chosen as the objective function due to its robustness to outliers. The network was trained using the Adam optimizer with a learning rate of 0.005, and the number of training epochs was set to 200 to ensure convergence.

The computational experiments were performed on a workstation equipped with an Intel Core i7-7700K CPU (4.20 GHz, 4 cores) (Intel, Santa Clara, CA, USA), 16 GB DDR4 RAM, an NVIDIA GeForce GTX 1050 Ti GPU (NVIDIA, Santa Clara, CA, USA), and 500 GB SSD storage.

4. Results and Discussion

4.1. Performance of the Surrogate Model

In this study, an AR-CNN was utilized to construct the surrogate model, employing a training sample size of 2000 and a test sample set consisting of 500 randomly generated samples. The network training time was approximately 4.84 h. (A discussion on the impact of training sample size on surrogate model performance is detailed in Section 4.2.)

4.1.1. Training Sample Generation

The 2000 training samples and 500 test samples used in this study were generated through synthetic simulations using the SEAWAT model, which simulates coupled groundwater flow and solute transport in coastal aquifers. The generation of these samples was based on the groundwater model described in Section 3.1, where the aquifer’s boundary conditions, hydraulic conductivity field, and pollution source characteristics were defined.

Specifically, the training and test samples were generated by varying the key parameters of the aquifer system, including the hydraulic conductivity field (represented by Karhunen–Loève expansion coefficients) and the pollution source strengths across different stress periods. Each sample corresponds to a unique combination of these parameters, resulting in a diverse set of scenarios that reflect the inherent uncertainty and heterogeneity of coastal aquifers. The outputs of each simulation include the spatial distributions of salinity, contaminant concentration, and hydraulic head at multiple timesteps, which were then used as the target outputs for training and evaluating the AR-CNN surrogate model.

Notably, the test samples were generated independently from the training samples, using a separate set of parameter combinations that were not included in the training dataset. This ensured that the test samples were completely unseen by the model during training, providing an unbiased evaluation of the surrogate model’s generalization capability.

4.1.2. Surrogate Model Evaluation

The surrogate model was used to predict the seawater intrusion salinity distribution, contaminant plume distribution, and hydraulic head field at different times. Comparisons were made between the AR-CNN-predicted fields ($\widehat{y}$) and the true fields (y) obtained from the test samples generated by the SEAWAT model. The results are presented in Figures S1–S3 in the Supplementary Materials. Specifically:

• Figure S1 shows the predicted, true, and absolute error fields for seawater intrusion salinity at different times.

• Figure S2 illustrates the predicted, true, and absolute error fields for the contaminant plume distribution.

• Figure S3 depicts the predicted, true, and absolute error fields for the hydraulic head field.

From the results shown in Figures S1–S3, it can be observed that the surrogate model developed using AR-CNN was able to reasonably replicate the true results of the SEAWAT model. There was good agreement in the salinity, contaminant concentration, and hydraulic head fields between the predictions and the true values. However, notable errors were observed in the predicted results near the pollution sources and the saline wedge. This discrepancy is possible due to significant gradients in concentration and salinity in these areas. Such steep gradients present challenges for the surrogate model, making it difficult to accurately capture localized dynamic changes.

To further evaluate the accuracy of the surrogate model quantitatively, the root mean square error (RMSE) and coefficient of determination ($R^2$) were calculated for both the training and testing samples. The results showed that the RMSE during training (${\mathrm{RMSE}}_{train}$) was 0.0437, while the RMSE during testing (${\mathrm{RMSE}}_{test}$) was 0.0442. The value for training ($R_{train}^2$) was 0.9994, and for testing ($R_{test}^2$) it was 0.9993. These results indicate that the surrogate model has excellent predictive accuracy and a stable generalization capability. The minimal difference between the training and testing metrics suggests that the model is not overfitting and can effectively predict unseen data.

4.2. Influence of Training Sample Size

To explore the impact of the training sample size ${\mathrm{N}}_{\mathrm{e}}$ on the predictive accuracy of the surrogate model, we considered four different sample sizes: 500, 1000, 2000, and 3000 training samples. The training dataset included seawater intrusion salinity distribution, contaminant plume distribution, and hydraulic head field distribution at 10 specific times $\left[55, 60, 65, 70, 75, 80, 85, 90, 95, 100\right]$ [min].

As shown in Figure 7, the predictive error of the surrogate model began to decrease significantly after approximately 80 training epochs, with further improvements observed as the number of epochs increased. When the number of training samples was increased from 500 to 1000, the RMSE of the model decreased significantly from 0.6693 to 0.0567, representing a reduction of 91.53%. This indicates that increasing the number of training samples can effectively improve the model’s predictive performance. When the training sample size was further increased to 2000, the RMSE dropped from 0.0567 to 0.0437 (a reduction of approximately 22.93%), indicating a noticeable yet smaller improvement in accuracy. However, when increasing the training sample size from 2000 to 3000, the decrease in RMSE tended to plateau, reducing only from 0.0437 to 0.0396 (a reduction of approximately 9.38%).

Additionally, $R^2$ values were calculated for both the training and testing samples. With the training sample size of 500, the $R^2$ values for the training and testing samples were 0.8630 and 0.8644, respectively. These relatively low $R^2$ values indicate that the surrogate model’s accuracy and generalization capability were not yet at an ideal level. As the training sample size increased to 2000, the $R^2$ values for the training and testing samples increased significantly to 0.9994 and 0.9993, respectively, demonstrating that the surrogate model could accurately capture the complex relationships between input and output data. As the number of training samples continued to increase from 2000 to 3000, the improvement in $R^2$ value became minimal, with the $R^2$ for the training set increasing to 0.9995 (an improvement of about 0.01%) and the $R^2$ for the test set to 0.9997 (an improvement of about 0.04%). This small increase indicates that after reaching a certain training sample size, the performance improvement from further increasing the sample size becomes limited, while the computational cost increases significantly (from 4.84 h of training time for 2000 samples to 7.88 h for 3000 samples).

Therefore, considering the trade-off between surrogate model accuracy and computational cost, the AR-CNN-based surrogate model with 2000 training samples was deemed the optimal choice. This configuration ensures high accuracy and good generalization capability while maintaining training costs and time.

4.3. Parameter Inversion of Coastal Groundwater Model Based on AR-CNN-ILUES Framework

Using the AR-CNN surrogate model with 2000 training samples, we applied the AR-CNN-ILUES framework to identify the parameters of the coastal groundwater model, successfully estimating the release history of the three pollution sources ${\mathrm{S}}_i^j$.

Figure 8 presents selected results for ${\mathrm{S}}_i^j$ where $i=1, 2, 3$ and $j=1, 3, 5, 7, 9$. By examining the prior and posterior distributions of the pollution source intensities, it can be observed that the AR-CNN-ILUES framework accurately identified the pollution source strengths.

Figure 9 illustrates the characterization results of the hydraulic conductivity field. In Figure 9, the hydraulic conductivity field is represented by the parameter lnK, which denotes the natural logarithm of the hydraulic conductivity K. The use of lnK is motivated by the fact that hydraulic conductivity in natural aquifers typically follows a log-normal distribution, as shown in Figure 6. This log-normal distribution reflects the inherent heterogeneity of aquifer materials, where hydraulic conductivity values can span several orders of magnitude.

Figure 9b–d illustrate three randomly selected posterior estimate fields of the hydraulic conductivity obtained through the AR-CNN-ILUES framework. Each figure represents a distinct realization of the hydraulic conductivity field, reflecting the variability and uncertainty in the inversion results. The key observations from these figures are as follows: (1) The three posterior estimate fields show a close resemblance to the reference field (Figure 9a) in many regions, particularly in areas with dense observational data. (2) Despite the overall similarity, there are noticeable variations between the three posterior estimate fields. These variations reflect the inherent uncertainty in the inversion process, particularly in regions with sparse observational data or complex geological structures. For example, the left edge of the aquifer, where monitoring wells are sparse, shows more variability compared to the right side, where data are more abundant.

From the variance field in Figure 9f, it can be seen that after seven iterations, the uncertainty in most areas decreased to a low level, although some uncertainty remained along the left edge where observational data are sparse. However, comparing the posterior mean field in Figure 9e to the reference field in Figure 9a, it is evident that the identification accuracy of the highly permeable region on the left side of the aquifer (within the dashed box) was relatively low. The main reasons for this are as follows:

(1) The active pollution sources in this case were located near the upper portion of the aquifer. Consequently, the pollutant plume distribution was primarily concentrated in the upper and middle parts of the aquifer, making it difficult for monitoring wells in the dashed box area to capture pollutant concentration information, thus affecting the inversion results.

(2) Compared to the densely distributed monitoring wells near the saline wedge on the right side, the monitoring wells within the dashed box area were sparse. This lack of monitoring wells might have contributed to the poor inversion performance in that region. Increasing the number of monitoring wells may improve the inversion accuracy in this area.

(3) While the surrogate model achieved high accuracy in single runs, the ILUES algorithm requires tens of thousands of model calls, which could lead to the accumulation and amplification of errors over iterations.

4.4. Comparison with Traditional SEAWAT-ILUES Framework

The primary objective of this section is to quantify the impact of using a surrogate model (AR-CNN-ILUES) in place of the SEAWAT model on parameter identification accuracy during coastal groundwater parameter inversion. For this purpose, we compared the inversion results based on the original SEAWAT-ILUES framework and the AR-CNN-ILUES framework.

Figure 10 and Table 4 show the estimated values of pollution source intensities $S_i^j$(i = 1, 2, 3; j = 1, 2, 3, 10) obtained from both the SEAWAT-ILUES and AR-CNN-ILUES frameworks. Although the parameter identification results from the original SEAWAT model were closer to the true release intensities of the pollution sources, the data assimilation results using the AR-CNN surrogate model were also effective in estimating the unknown parameters, with significant advantages in computational efficiency.

Table 5 presents the accuracy of the parameter inversion results for both the original SEAWAT model (denoted as “O”) and the surrogate model (denoted as “S”). Notably, under the SEAWAT-ILUES framework, the uncertainty of parameters is significantly reduced after iterative updates, leading to a more precise estimation of the release intensities of multipoint pollution sources. This is specifically reflected in a lower Mean Absolute Error (${\mathrm{MAE}}_s$ = 0.4804) and a reduced Average Ensemble Spread (${\mathrm{AES}}_s$ = 0.4753). In contrast, the inversion results obtained using the AR-CNN-ILUES framework demonstrate higher ensemble spread, with an ${\mathrm{MAE}}_s$ of 1.1395 and an ${\mathrm{AES}}_s$ of 2.3852. This discrepancy is likely due to the limitations of surrogate models when handling highly complex dynamic systems, particularly in scenarios requiring the accurate capture of temporal dynamics and intricate geological variability.

Figure 11 compares the reference hydraulic conductivity field with the inversion results under both frameworks. The inversion results based on the original SEAWAT model accurately captured the general contours of high- and low-permeability areas, as evidenced by the lower MAE (${\mathrm{MAE}}_k$ = 0.6059) and AES (${\mathrm{AES}}_k$ = 0.3072). However, outside the saline wedge region where monitoring points are densely distributed, particularly in the left-side boxed area, the inversion accuracy decreases significantly due to the lack of pollutant plume observation data and sparse monitoring wells. This limitation is more pronounced in the results from the surrogate model framework, where ${\mathrm{MAE}}_k$ increased to 0.7610 and ${\mathrm{AES}}_k$ slightly increased to 0.3222. The surrogate-based ILUES method struggles in this region because of the limited observational data and steep concentration gradients, which reduce the ability of the surrogate model to constrain the inversion effectively.

In summary, the surrogate model exhibits certain limitations when dealing with highly complex dynamic systems, particularly in regions with strong nonlinearity in geological structures and dynamic changes. Nevertheless, the inversion errors under the AR-CNN-ILUES framework remain within an acceptable range, maintaining a satisfactory level of accuracy, and ensuring the reliability of the estimated values.

4.5. Computational Efficiency

While the AR-CNN-ILUES framework introduces some errors compared to the SEAWAT-ILUES framework, it significantly improves computational efficiency. To quantify the efficiency of parameter inversion in both frameworks, this section compares their computational performance.

When using the SEAWAT-ILUES algorithm for coastal groundwater parameter identification, the main computational cost comes from the forward simulations of the groundwater numerical model. This cost is especially high because the ILUES data assimilation process requires multiple iterations with the original SEAWAT model. The total time required for the SEAWAT-ILUES, denoted as $T_{ILUES}$, is formulated as follows:

(13)$T_{ILUES}=\left(1+N_i\right)\times N_e\times t_c$

Here, $N_i$ is the number of ILUES algorithm iterations, $N_e$ is the number of ensemble samples, and $t_c$ is the runtime for a single SEAWAT model execution.

Specifically, in our pollutant source identification problem involving three pollution sources and geological survey uncertainty, we set $N_e=2000$ and $N_i=7$. The runtime for a single model execution was $t_c\approx 15.843 \mathrm{s}$. Therefore, the total inversion time using the SEAWAT-ILUES was approximately $T_{ILUES}\approx 2.535\times {10}^5$ s (approximately 70.4 h).

In contrast, the AR-CNN-ILUES framework requires time for generating training samples, training the surrogate model, and performing inversion calculations after training. The total time required, $T_{CNN}$, can be expressed as:

(14)$T_{CNN}=N_e\times t_c+t_s+t_v$

where $t_s$ is the time to train the surrogate model, and $t_v$ is the time for ILUES inversion calculations after model training.

For our case, $N_e\times t_c=8.8 \mathrm{h}$, $t_s\approx 6.2 \mathrm{h}$, and $t_v=1.22 \mathrm{h}$. Thus, the total inversion time using the AR-CNN-ILUES framework, ${\mathrm{T}}_{\mathrm{CNN}}$, is approximately $5.84\times {10}^4$ s (16.2 h). Compared to the SEAWAT-ILUES framework, this represents a computational efficiency improvement of approximately 77%. Notably, after training the surrogate model, the efficiency improvement for each iteration reaches nearly 99%.

For more complex estuarine groundwater pollution identification problems—such as those involving slopes, impermeable walls, or chemical reactions—the runtime for a single SEAWAT model execution could extend to several hours. Additionally, achieving more reliable inversion results might require increasing both the ensemble sample count $N_e$ and the number of iterations $N_{iter}$ of the ILUES algorithm. In such cases, using a surrogate model for high-dimensional inversion problems would result in even more significant time savings.

4.6. Why Use AR-CNN?

In this study, we selected AR-CNN as the surrogate model for the joint inversion of hydraulic conductivity fields and contamination source parameters. The choice of AR-CNN was motivated by its unique advantages in handling the complex dynamics and spatial heterogeneity of coastal aquifer systems. Below, we provide a detailed explanation of why AR-CNN is particularly well suited for this task:

• (1). Spatial Feature Extraction

The convolutional architecture of AR-CNN is highly effective at capturing spatial heterogeneity in hydraulic conductivity fields. Coastal aquifers are characterized by complex geological structures and significant spatial variability, which are critical to accurately model. The convolutional layers in AR-CNN excel at extracting spatial features from input data, enabling the model to effectively represent the high-dimensional and nonlinear relationships inherent in aquifer systems.

• (2). Temporal Dynamics

The auto-regressive strategy in AR-CNN allows it to handle temporal dynamics, making it suitable for modeling time-dependent processes such as contaminant transport and saltwater intrusion. By incorporating the output from the previous timestep as part of the input for the current step, AR-CNN can capture the temporal evolution of groundwater flow and solute transport. This capability is essential for accurately simulating the transient behavior of coastal aquifers.

• (3). Computational Efficiency

Compared to other deep learning methods like Generative Adversarial Networks (GANs) or Recurrent Neural Networks (RNNs), AR-CNN is more computationally efficient and stable. GANs, while powerful for data generation tasks, often require extensive training data and computational resources, and their training process can be unstable. RNNs, particularly Long Short-Term Memory (LSTM) networks, are effective for sequential data but struggle with spatial feature extraction and are computationally intensive for high-dimensional problems. In contrast, AR-CNN strikes a balance between accuracy and efficiency, making it ideal for iterative inversion processes like the Iterative Local Updating Ensemble Smoother (ILUES).

• (4). Suitability for Joint Inversion

The combination of spatial feature extraction, temporal dynamics handling, and computational efficiency makes AR-CNN particularly well suited for the joint inversion of hydraulic conductivity fields and contamination source parameters. These capabilities enable AR-CNN to accurately approximate the complex input–output relationships of coastal aquifer systems while significantly reducing the computational burden compared to traditional numerical models like SEAWAT.

5. Conclusions

In this study, we developed a deep learning-based approach to address the challenges of characterizing coastal aquifer contamination and performing hydraulic conductivity field inversion. By integrating an AR-CNN-based surrogate model with the ILUES framework, we aimed to enhance computational efficiency without significantly compromising accuracy. The key findings of this study are summarized as follows:

1. The AR-CNN surrogate model demonstrated excellent predictive accuracy, with the training and testing $R^2$ values exceeding 0.999 and the RMSE values below 0.05. The model effectively replicated the behavior of the SEAWAT model in simulating seawater intrusion salinity distribution, contaminant plume migration, and hydraulic head fields. However, minor discrepancies were observed near pollution sources and saline wedges, likely due to steep concentration gradients that are challenging for the surrogate model to capture fully.

2. The AR-CNN-ILUES framework successfully estimated the release history of multipoint pollution sources and characterized the heterogeneous hydraulic conductivity field of the coastal aquifer. The inversion results showed that the framework could accurately identify pollution source strengths and capture the general contours of high- and low-permeability areas. However, limitations were noted in regions with sparse monitoring points and complex geological structures, where inversion accuracy decreased.

3. The AR-CNN-ILUES framework significantly improved computational efficiency, reducing the total inversion time from approximately 70.4 h to 16.2 h, an improvement of about 77%. The efficiency gains were even more substantial in the iterative inversion process after surrogate model training, with improvements reaching nearly 99%.

The results of this study highlight the potential of integrating deep learning techniques with data assimilation frameworks to overcome computational bottlenecks in high-dimensional, nonlinear inverse problems associated with coastal aquifers. The AR-CNN-ILUES framework demonstrates a promising avenue for efficient and accurate parameter inversion, which is crucial for predicting contaminant transport and assessing seawater intrusion. However, several limitations should be noted:

1. Sparse Monitoring Points: The accuracy of the inversion results is influenced by the distribution of monitoring points. In regions with sparse monitoring data, the inversion performance is less reliable. Future studies could explore strategies such as optimizing sensor placement or incorporating remote sensing data to address this issue.

2. Steep Gradients and Complex Structures: The surrogate model faces challenges in accurately capturing steep concentration gradients and complex geological structures. Alternative approaches, such as incorporating domain-specific knowledge or refining the network architecture, could improve model performance in these areas.

3. Synthetic Data Limitations: The current framework relies on synthetic data, which may not fully capture the complexities of real-world scenarios. Future work will focus on integrating field measurements to validate the framework in practical settings.

Author Contributions

Conceptualization, X.C. and S.J.; Formal analysis, X.C. and Y.C.; Methodology, S.J.; Project administration, X.C. and S.J.; Resources, X.C. and Y.C.; Software, C.W.; Supervision, S.J.; Validation, C.T. and D.L.; Writing—original draft, X.C. and C.W.; Writing—review and editing, S.J. All authors have read and agreed to the published version of the manuscript.

Data Availability Statement

The code supporting this study is currently undergoing processing for archiving and will be made publicly available upon completion.

Conflicts of Interest

Chao Wu was employed by the company Yangtze Ecology and Environment Co., Ltd. The remaining authors declare that this research was conducted in the absence of any commercial or financial relationships that could be construed as potential conflicts of interest.

Footnotes

Disclaimer/Publisher’s Note: The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.

Figures and Tables

Figure 3. Schematic diagram of AR-CNN-based surrogate model structure.

Figure 4. Comparison of SEAWAT-ILUES and AR-CNN-ILUES inversion frameworks.

Figure 5. (a) shows a cross-section of the coastal confined aquifer, including the model domain, pollutant source locations, and boundary conditions, while (b) shows the model grid, observation wells (W1–W10), and observation points.

Figure 6. The reference hydraulic conductivity field in the coastal aquifer cross-section.

Figure 7. The trend in RMSE values for different training sample sizes as the number of training epochs increases.

Figure 8. Estimated values of contaminant source strength [Forumla omitted. See PDF.] based on the AR-CNN-ILUES framework.

Figure 9. Characterization results of the hydraulic conductivity field based on the AR-CNN-ILUES framework: (a) reference hydraulic conductivity field; (b–d) three posterior estimate fields; (e) posterior mean field of hydraulic conductivity; (f) variance field of posterior distribution.

Figure 10. Inversion results of pollution source intensities using SEAWAT-ILUES and AR-CNN-ILUES frameworks.

Figure 11. Comparison of the reference hydraulic conductivity field and the inversion results: (a) reference hydraulic conductivity field; (b) inversion results under the SEAWAT-ILUES framework; (c) inversion results under the AR-CNN-ILUES framework.

Table 1

Main parameters of the coastal groundwater numerical model.

Parameter	Symbol	Unit	Value
Hydraulic conductivity	K	cm/min	KL expansion
Molecular diffusion	D	cm2/s	0.066
Specific storage	Ss	1/cm	1.0 × 10−5
Longitudinal dispersion	${\alpha }_L$	cm	0.50
Transverse dispersion	${\alpha }_T$	cm	0.05
Density of the freshwater	${\rho }_f$	g/cm3	1.000
Density of the seawater	$\rho$	g/cm3	1.025
Freshwater salinity	$C_f$	g/L	0.0
Seawater salinity	$C_{sw}$	g/L	35.0
Aquifer porosity	$\theta$	-	0.30

Table 2

The reference values and prior range of pollution source strength.

Stress Period	Unit	Prior Range	Reference Value
			S1	S2	S3
SP1	g/L	u [0,50]	20.00	22.50	0.00
SP2	g/L	u [0,50]	25.00	25.00	0.00
SP3	g/L	u [0,50]	28.00	27.60	0.00
SP4	g/L	u [0,50]	22.00	34.20	0.00
SP5	g/L	u [0,50]	30.00	22.00	0.00
SP6	g/L	u [0,50]	32.00	28.50	0.00
SP7	g/L	u [0,50]	24.00	30.80	0.00
SP8	g/L	u [0,50]	26.00	32.40	0.00
SP9	g/L	u [0,50]	34.00	33.60	0.00
SP10	g/L	u [0,50]	27.00	35.00	0.00

Table 3

AR-CNN structure.

Structure	Convolution Kernel (k,s,p)	${\mathit{N}}_{\mathit{f}}\times {\mathit{H}}_{\mathit{f}}\times {\mathit{W}}_{\mathit{f}}$
Input Layer	-	3$\times 20\times 80$
Initial Convolutional Layer	(7,2,3)	48$\times 10\times 40$
Dense Block 1 (5 layers)	(3,1,1)	248$\times 10\times 40$
Encoding Layer 1 (Convolution)	(1,1,0), (3,2,1)	124$\times 5\times 20$
Dense Block 2 (10 layers)	(3,1,1)	524$\times 5\times 20$
Decoding Layer 1 (Deconvolution)	(1,1,0), (3,2,1)	262×10$\times 40$
Dense Block 3 (5 layers)	(3,1,1)	462$\times$10$\times 40$
Output Layer (Decoding Layer 2, Deconvolution)	(1,1,0), (6,2,1)	3$\times$20$\times 80$

Note: $N_f$ represents the number of feature maps, and $H_f\times W_f$ denotes the size of the feature maps. In the convolution kernel column, k indicates the kernel size, s is the stride, and p is the amount of zero padding applied.

Table 4

Comparison of pollution source intensity estimation results between SEAWAT-ILUES (O) and AR-CNN-ILUES (S) frameworks.

Parameter	Reference	Mean (O)	Mean (S)	Std (O)	Std (S)	95% CI (O)	95% CI (S)
$S_1^1$	20.00	20.14	20.37	1.81	2.04	(16.59, 23.69)	(16.37, 24.37)
$S_1^2$	25.00	24.73	24.54	1.83	2.12	(21.14, 28.32)	(20.38, 28.70)
$S_1^3$	28.00	28.69	26.61	1.85	2.64	(25.06, 32.32)	(21.44, 31.78)
$S_1^{10}$	27.00	26.78	25.53	1.94	2.56	(22.98, 30.58)	(20.51, 30.55)
$S_2^1$	22.50	22.73	23.32	2.03	1.64	(18.75, 26.71)	(20.11, 26.53)
$S_2^2$	25.00	24.97	24.63	1.86	2.66	(21.32, 28.62)	(19.42, 29.84)
$S_2^3$	27.60	25.87	25.48	2.31	2.78	(21.34, 30.40)	(20.03, 30.93)
$S_2^{10}$	35.00	34.57	35.62	2.34	2.98	(29.98, 39.16)	(29.78, 41.46)
$S_3^1$	0.00	0.23	2.41	1.94	1.78	(−4.04, 4.50)	(−1.39, 6.21)
$S_3^2$	0.00	0.36	1.78	1.78	1.78	(−5.01, 5.73)	(−1.77, 5.27)
$S_3^3$	0.00	0.27	2.61	2.27	2.27	(−5.08, 5.62)	(−1.84, 7.06)
$S_3^{10}$	0.00	0.22	1.59	1.86	1.86	(−5.15, 5.59)	(−2.06, 5.24)

Note: Reference: Reference values for pollution source intensities. Mean (O): Mean values from the original SEAWAT-ILUES framework. Mean (S): Mean values from the AR-CNN-ILUES framework. Std (O): Standard deviation from the original SEAWAT-ILUES framework. Std (S): Standard deviation from the AR-CNN-ILUES framework. 95% CI (O): 95% confidence interval from the original SEAWAT-ILUES framework. 95% CI (S): 95% confidence interval from the AR-CNN-ILUES framework.

Table 5

Parameter inversion accuracy of the original model (O) and surrogate model (S).

	$\mathbf{M}\mathbf{A}{\mathbf{E}}_{\mathit{k}}$	$\mathbf{M}\mathbf{A}{\mathbf{E}}_{\mathit{s}}$	$\mathbf{A}\mathbf{E}{\mathbf{S}}_{\mathit{k}}$	$\mathbf{A}\mathbf{E}{\mathbf{S}}_{\mathit{s}}$
Original Model (O)	0.6059	0.4804	0.3072	0.4753
Surrogate Model (S)	0.7610	1.1395	0.3222	2.3852

Supplementary Materials

The following supporting information can be downloaded at: https://www.mdpi.com/article/10.3390/w17020255/s1, Figure S1: the predicted, true, and absolute error fields for seawater intrusion salinity at different times; Figure S2: the predicted, true, and absolute error fields for the contaminant plume distribution; Figure S3: the predicted, true, and absolute error fields for the hydraulic head field.

References

1. Karthikeyan, S.; Kulandaisamy, P.; Senapathi, V.; Thangaraj, K.; Arumugam, M.; Sekar, S.; Ranganathan, P.C. Issues of Coastal Groundwater Contamination. Groundwater Contamination in Coastal Aquifers: Assessment and Management; Elsevier: Amsterdam, The Netherlands, 2022; pp. 9-18.

2. Stein, S.; Shalev, E.; Sivan, O.; Yechieli, Y. Challenges and Approaches for Management of Seawater Intrusion in Coastal Aquifers. Hydrogeol. J.; 2023; 31, pp. 19-22. [DOI: https://dx.doi.org/10.1007/s10040-022-02575-5]

3. Re, V.; Sacchi, E.; Mas-Pla, J.; Menció, A.; El Amrani, N. Identifying the Effects of Human Pressure on Groundwater Quality to Support Water Management Strategies in Coastal Regions: A Multi-Tracer and Statistical Approach (Bou-Areg Region, Morocco). Sci. Total Environ.; 2014; 500–501, pp. 211-223. [DOI: https://dx.doi.org/10.1016/j.scitotenv.2014.08.115]

4. Basack, S.; Loganathan, M.K.; Goswami, G.; Khabbaz, H. Saltwater Intrusion into Coastal Aquifers and Associated Risk Management: Critical Review and Research Directives. J. Coast. Res.; 2022; 38, pp. 654-672. [DOI: https://dx.doi.org/10.2112/JCOASTRES-D-21-00116.1]

5. Hussain, M.S.; Abd-Elhamid, H.F.; Javadi, A.A.; Sherif, M.M. Management of Seawater Intrusion in Coastal Aquifers: A Review. Water; 2019; 11, 2467. [DOI: https://dx.doi.org/10.3390/w11122467]

6. Dey, S.; Prakash, O. Coupled Sharp-Interface and Density-Dependent Model for Simultaneous Optimization of Production Well Locations and Pumping in Coastal Aquifer. Water Resour. Manag.; 2022; 36, pp. 2327-2341. [DOI: https://dx.doi.org/10.1007/s11269-022-03145-w]

7. Wu, C.; Jiang, S.; Xia, X.; Sun, Y.; Li, Z.; Ju, M.; Li, S.; Liu, S. Estimation of Pollution Sources and Hydraulic Conductivity Field in a Coastal Aquifer under Tidal Effects. Mar. Georesources Geotechnol.; 2024; 42, pp. 1415-1427. [DOI: https://dx.doi.org/10.1080/1064119X.2023.2280636]

8. Dodangeh, A.; Rajabi, M.M.; Carrera, J.; Fahs, M. Joint Identification of Contaminant Source Characteristics and Hydraulic Conductivity in a Tide-Influenced Coastal Aquifer. J. Contam. Hydrol.; 2022; 247, 103980. [DOI: https://dx.doi.org/10.1016/j.jconhyd.2022.103980] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/35245819]

9. Singh, A.; Minsker, B.S.; Valocchi, A.J. An Interactive Multi-Objective Optimization Framework for Groundwater Inverse Modeling. Adv. Water Resour.; 2008; 31, pp. 1269-1283. [DOI: https://dx.doi.org/10.1016/j.advwatres.2008.05.005]

10. Markovich, K.H.; White, J.T.; Knowling, M.J. Sequential and Batch Data Assimilation Approaches to Cope with Groundwater Model Error: An Empirical Evaluation. Environ. Model. Softw.; 2022; 156, 105498. [DOI: https://dx.doi.org/10.1016/j.envsoft.2022.105498]

11. Zhang, J.; Lin, G.; Li, W.; Wu, L.; Zeng, L. An Iterative Local Updating Ensemble Smoother for Estimation and Uncertainty Assessment of Hydrologic Model Parameters With Multimodal Distributions. Water Resour. Res.; 2018; 54, pp. 1716-1733. [DOI: https://dx.doi.org/10.1002/2017WR020906]

12. Steklova, K.; Haber, E. Joint Hydrogeophysical Inversion: State Estimation for Seawater Intrusion Models in 3D. Comput. Geosci.; 2017; 21, pp. 75-94. [DOI: https://dx.doi.org/10.1007/s10596-016-9595-y]

13. Rao, S.; Sreenivasulu, V.; Bhallamudi, S.M.; Thandaveswara, B.; Sudheer, K. Planning Groundwater Development in Coastal Aquifers/Planification Du Développement de La Ressource En Eau Souterraine Des Aquifères Côtiers. Hydrol. Sci. J.; 2004; 49, pp. 155-170. [DOI: https://dx.doi.org/10.1623/hysj.49.1.155.53999]

14. Miao, T.; Guo, J. Application of Artificial Intelligence Deep Learning in Numerical Simulation of Seawater Intrusion. Environ. Sci. Pollut. Res.; 2021; 28, pp. 54096-54104. [DOI: https://dx.doi.org/10.1007/s11356-021-13680-5] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/34046828]

15. Dai, Z.; Zhan, C.; Soltanian, M.R.; Ritzi, R.W.; Zhang, X. Identifying Spatial Correlation Structure of Multimodal Permeability in Hierarchical Media with Markov Chain Approach. J. Hydrol.; 2019; 568, pp. 703-715. [DOI: https://dx.doi.org/10.1016/j.jhydrol.2018.11.032]

16. Zhan, C.; Dai, Z.; Samper, J.; Yin, S.; Ershadnia, R.; Zhang, X.; Wang, Y.; Yang, Z.; Luan, X.; Soltanian, M.R. An Integrated Inversion Framework for Heterogeneous Aquifer Structure Identification with Single-Sample Generative Adversarial Network. J. Hydrol.; 2022; 610, 127844. [DOI: https://dx.doi.org/10.1016/j.jhydrol.2022.127844]

17. Tartakovsky, A.M.; Marrero, C.O.; Perdikaris, P.; Tartakovsky, G.D.; Barajas-Solano, D. Physics-Informed Deep Neural Networks for Learning Parameters and Constitutive Relationships in Subsurface Flow Problems. Water Resour. Res.; 2020; 56, e2019WR026731. [DOI: https://dx.doi.org/10.1029/2019WR026731]

18. Zheng, N.; Jiang, S.; Xia, X.; Kong, W.; Li, Z.; Gu, S.; Wu, Z. Efficient Estimation of Groundwater Contaminant Source and Hydraulic Conductivity by an ILUES Framework Combining GAN and CNN. J. Hydrol.; 2023; 621, 129677. [DOI: https://dx.doi.org/10.1016/j.jhydrol.2023.129677]

19. Zhou, Z.; Zabaras, N.; Tartakovsky, D.M. Deep Learning for Simultaneous Inference of Hydraulic and Transport Properties. Water Resour. Res.; 2022; 58, e2021WR031438. [DOI: https://dx.doi.org/10.1029/2021WR031438]

20. Jardani, A.; Vu, T.; Fischer, P. Use of Convolutional Neural Networks with Encoder-Decoder Structure for Predicting the Inverse Operator in Hydraulic Tomography. J. Hydrol.; 2022; 604, 127233. [DOI: https://dx.doi.org/10.1016/j.jhydrol.2021.127233]

21. Mo, S.; Zabaras, N.; Shi, X.; Wu, J. Integration of Adversarial Autoencoders with Residual Dense Convolutional Networks for Estimation of Non-Gaussian Hydraulic Conductivities. Water Resour. Res.; 2020; 56, e2019WR026082. [DOI: https://dx.doi.org/10.1029/2019WR026082]

22. He, T.; Wang, N.; Zhang, D. Theory-guided full convolutional neural network: An efficient surrogate model for inverse problems in subsurface contaminant transport. Adv. Water Resour.; 2021; 157, 104051. [DOI: https://dx.doi.org/10.1016/j.advwatres.2021.104051]

23. Mo, S.; Zabaras, N.; Shi, X.; Wu, J. Deep autoregressive neural networks for high-dimensional inverse problems in groundwater contaminant source identification. Water Resour. Res.; 2020; 55, pp. 3856-3881. [DOI: https://dx.doi.org/10.1029/2018WR024638]

24. Zhang, R.; Zhou, N.; Xia, X.; Zhao, G.; Jiang, S. Joint Estimation of Hydraulic and Biochemical Parameters for Reactive Transport Modelling with a Modified ILUES Algorithm. Water; 2020; 12, 2161. [DOI: https://dx.doi.org/10.3390/w12082161]

25. Cao, C.; Zhang, J.; Gan, W.; Nan, T.; Lu, C. A Deep Learning-Based Data Assimilation Approach to Characterizing Coastal Aquifers Amid Non-Linearity and Non-Gaussianity Challenges. Water Resour. Res.; 2024; 60, e2023WR036899. [DOI: https://dx.doi.org/10.1029/2023WR036899]

26. Nan, T.; Zhang, J.; Xie, Y.; Cao, C.; Wu, J.; Lu, C. Effective Characterization of Fractured Media With PEDL: A Deep Learning-Based Data Assimilation Approach. Water Resour. Res.; 2024; 60, e2023WR036673. [DOI: https://dx.doi.org/10.1029/2023WR036673]

27. Mo, S.; Zhu, Y.; Zabaras, N.; Shi, X.; Wu, J. Deep Convolutional Encoder-Decoder Networks for Uncertainty Quantification of Dynamic Multiphase Flow in Heterogeneous Media. Water Resour. Res.; 2019; 55, pp. 703-728. [DOI: https://dx.doi.org/10.1029/2018WR023528]

28. Zhang, D.; Lu, Z. An efficient, high-order perturbation approach for flow in random porous media via Karhunen-Loeve and poly-nomial expansions. J. Comput. Phys.; 2004; 194, pp. 773-794. [DOI: https://dx.doi.org/10.1016/j.jcp.2003.09.015]

29. Panjehfouladgaran, A.; Rajabi, M. Contaminant source characterization in a coastal aquifer influenced by tidal forces and density-driven flow. J. Hydrol.; 2022; 610, 127807. [DOI: https://dx.doi.org/10.1016/j.jhydrol.2022.127807]

30. Xin, P.; Robinson, C.; Li, L.; Barry, D.; Bakhtyar, R. Effects of wave forcing on a subterranean estuary. Water Resour. Res.; 2010; 46, e2010WR009632. [DOI: https://dx.doi.org/10.1029/2010WR009632]