GeoFUSE: A High-Efficiency Surrogate Model for Seawater Intrusion Prediction and Uncertainty Reduction

The Beaver Creek model from Yabusaki et al. (2020) was developed using PFLOTRAN Hammond et al. (2014), a computational tool designed to simulate subsurface flow and reactive transport processes in variably saturated porous media. Given these capabilities, PFLOTRAN is particularly useful for studying seawater intrusion as it effectively captures hydrological and biogeochemical dynamics in complex environments such as aquifers and floodplains. By employing advanced numerical techniques to solve flow and transport equations, PFLOTRAN enables researchers to simulate water movement, salinity intrusion, and contaminant transport scenarios.

Yabusaki et al. (2020) implemented a single-phase solute transport model in variably saturated porous media to simulate the Beaver Creek system. The 2D flow and salinity transport simulations were conducted with the PFLOTRAN simulator Hammond et al. (2014). Darcy’s and Richard’s equations were applied to describe 2D flow and transport process, making them well-suited for capturing the dynamics of groundwater flow and salinity intrusion. The governing mass conservation equation is written as

where $\phi$ is the porosity [-], $S$ is saturation [-], $\eta$ is the water molar density [kmol/m³], $\mathbf{q}$ represents the water flux [m/s], and $Q_{w}$ represents the source/sink term [kmol/m³/s]. The Darcy flux $\mathbf{q}$ is defined by

where $\mathbf{k}$ denotes the intrinsic permeability tensor [m²], $k_{r}(S)$ denotes the relative permeability [-] (the relationship between $k_{r}$ and $S$ is derived on laboratory experiments), $\mu$ is the viscosity [Pa$\cdot$s], $P$ is the pressure [Pa], $\rho$ is the water density [kg/m³], $g$ is the gravitational field [m/s²], and $z$ is the vertical position [m]. The van Genuchten function is used for the water retention and relative permeability curve, with specific parameters such as a residual saturation of 0.15, alpha of $2\times 10^{-4}$ Pa, and m of 0.2908.

For the seawater intrusion model, the reactive transport equation for salt mass considering advection and diffusion is written as

where $c$ is the contaminant concentration [kg/m³], $\mathbf{D}$ denotes the diffusion coefficient [m²/s], and $Q_{r}$ is the source/sink term of salinity [kg/m³/s].

The system is a 2D cross-section from the upland to the stream with the same dimensions and elevation as the cross-section described in Yabusaki et al. (2020). The flow model has dimensions of 84 m $\times$ 4 m and is defined on a grid with 28 $\times$ 40 cells. As shown in Fig. 1, the flow system is driven by Dirichlet pressure boundary conditions on the upland and stream sides. In contrast, a seepage boundary condition is applied to the ground surface and a no-flow boundary condition at the bottom. The source of salinity is the tidal stream, with the boundary conditions set based on one-year observations from hillslope wells and stream data Yabusaki et al. (2020). For more details on the model, please refer to Yabusaki et al. (2020).

Compared to the homogeneous case in Yabusaki et al. (2020), we simulate a more complicated heterogeneous field in this work. Figure 2 present three random realizations of log-permeability. These Gaussian realizations are generated using sequential Gaussian simulation with SGeMS Remy et al. (2009). The correlation lengths $l_{x}$ and $l_{z}$ are 5 and 5 grids, respectively. The mean and standard deviation of log-permeability (md) are 4.5 and 1, while the mean and standard deviation of porosity are 0.5 and 0.05. The same Gaussian field is used for both permeability and porosity, which serves as the input for the surrogate model. The other parameters are set based on Yabusaki et al. (2020). The model in Fig. 2(a) is used as the “true” model for data assimilation in Section 3.

The model is initially simulated for four years using annual observations to spin up the system to a dynamic steady state of groundwater levels, saturation, and total mass. After the spin-up, the simulation represents the initial freshwater conditions before seawater intrusion, which is then modeled over 7260 days ($\sim$20 years). Monitoring wells, as shown in Fig. 1, are placed in the same locations as in Yabusaki et al. (2020) for model calibration and validation.

To efficiently address the uncertainty quantification and decision-making challenges in subsurface flow systems, it is necessary to discretize and solve equations 1 to 3 thousands to millions of times, which can be computationally prohibitive. To overcome this challenge, we develop a deep-learning-based surrogate model, U-FNO, to approximate the flow and transport processes for various permeability and porosity fields under consistent boundary conditions. The flow simulation process is expressed as

where $\mathbf{m}\in\mathbb{R}^{n_{c}}$ denotes the geological model parameters (the permeability values in each cell, with $n_{c}=n_{x}\times n_{z}$ as the total number of cells, and $n_{x}$ and $n_{z}$ denoting the number of cells in the 2D x-z domain), $f$ is the flow and transport simulation process, and $\mathbf{d}$ represents the dynamic states of pressure $\mathbf{P}\in\mathbb{R}^{n_{c}\times n_{t}}$ and salinity $\mathbf{c}\in\mathbb{R}^{n_{c}\times n_{t}}$ in each cell over $n_{t}$ time steps.

The forward process can be approximated by a surrogate model as

where $\hat{\mathbf{d}}\in\mathbb{R}^{2n_{c}\times n_{t}}$ represents the approximate flow responses, $\hat{f}$ denotes the surrogate model, and $\boldsymbol{\theta}$ represents the parameters of the neural networks.

In this work, we apply the U-FNO framework proposed by Wen et al. (2022). The U-FNO architecture is modified based on the FNO Li et al. (2020) with modifications that append a U-Net to the Fourier layer to improve prediction accuracy. U-FNO is designed to conduct an infinite-dimensional-space input output mapping. The architecture of U-FNO used in this work is presented in Fig. 3. The geomodel input $\mathbf{m}$ is discretized on the 2D domain $\mathbf{x}$ and is repeated $n_{t}$ times as the function input $m(\mathbf{x},t)$ for the networks. The input $a$ is then lifted to a higher-dimensional output $v_{0}(\mathbf{x},t)$ through fully connected neural networks $P(\cdot)$. Multiple iterative Fourier layers and U-Fourier layers are then applied to capture the spatio-temporal features of flow states. Another fully connected layer $Q(\cdot)$ is applied in the end to map the U-Fourier layer output to surrogate model output $d(\mathbf{x},t)$.

In each Fourier layer $i$, $i=1,\ldots,L$, the Fast Fourier Transform (FFT) is applied to the previous layer output $v_{i-1}(\mathbf{x},t)$ to generate Fourier domain $\mathcal{F}(v_{i-1})$ and capture the latent features. A linear transformation $\mathcal{R}(\cdot)$ is then applied, followed by the inverse Fourier transform $\mathcal{F}^{-1}(\cdot)$ to map the Fourier domain back to the spatio-temporal domain. Additionally, a linear transform $W(\cdot)$ is applied directly to the function space $v_{i-1}(\mathbf{x},t)$ and added to the output. The overall transformation process is written as

where $\sigma$ denotes the activation function, $L$ is the number of Fourier layer used in U-FNO. K U-Fourier layers are applied after $L$ Fourier layers. U-Net $U(\cdot)$ is attached in each Fourier layer to improve the surrogate model performance. The output $v_{i}(\mathbf{x},t)$ of U-Fourier layer is expressed as

for $i=L+1,\cdots,L+K$. The U-FNO architecture is written as

Please refer to Wen et al. (2022) for more details on U-FNO.

The training objective is to minimize the difference between the flow simulation output with the surrogate predictions from U-FNO. The minimization process of U-FNO parameters $\boldsymbol{\theta}$ is described as

where $n_{smp}$ is the number of training samples, $\lambda$ is the extra weight applied to the data in the active cells $\mathbf{d}^{active}$. These active cells are determined based on the digital elevation model for the 2D vertical cross-section. The adaptive moment estimation (ADAM) in Kingma & Ba (2014) is applied to conduct the minimization process.

To improve flow prediction accuracy, we apply data assimilation to calibrate the geological model, employing the GeoFUSE framework (Geological Forecasting using Surrogates and ESMDA). GeoFUSE integrates PCA for dimensionality reduction, the U-FNO surrogate model for high-fidelity flow prediction, and the ESMDA for iterative parameter calibration. This combination enables efficient model calibration, enhancing prediction accuracy while reducing uncertainty. A Bayesian framework is employed within the data assimilation process to generate posterior probability distributions.

In practical applications, monitoring data is collected from wells or geophysical observations. For synthetic cases like ours, the observed data is generated by adding random noise to a synthetic “true” model (Fig. 2(a)) as follows

where $\boldsymbol{\epsilon}$ is the measurement error. When no structural bias is considered, the measurement error $\boldsymbol{\epsilon}$ is assumed to be sampled from Gaussian distribution $N(\mathbf{0},C_{D})$, where $C_{D}$ denotes the error covariance matrix. The posterior probability density function of model parameter $\mathbf{m}$ conditioned to the observation $\mathbf{d}_{\text{obs}}$ using Bayesian framework is given by

where $p(\mathbf{m})$ is the prior distribution of model parameters, $p(\mathbf{d}_{\text{obs}}|\mathbf{m})$ is the likelihood function, and $p(\mathbf{d}_{\text{obs}})$ is a normalization factor that can be ignored.

Within the GeoFUSE framework as shown in Fig. 4, ESMDA Emerick & Reynolds (2013) is employed to sample posterior models, as illustrated in Eq. 11. Parameterization is generally needed to reduce the dimension of high-dimensional geological models and preserve the geological features. In this work, we consider a 2D Gaussian geomodel, and PCA is applied. The 2D geomodel can be generated from $\mathbf{m}\approx\mathbf{m}_{\text{PCA}}(\boldsymbol{\xi})=\Phi\boldsymbol{\xi}+\bar{\mathbf{m}}$, where $\boldsymbol{\xi}\in\mathbb{R}^{n_{l}}$ denotes the latent variables, $\Phi$ denotes the basis matrix generated from singular value decomposition, and $\bar{\mathbf{m}}$ is the average value of model parameters. Using U-FNO as a surrogate model, the flow responses can be generated from $\hat{f}(\mathbf{m}_{\text{PCA}}(\boldsymbol{\xi}))$.

In the GeoFUSE data assimilation framework, the latent variables $\boldsymbol{\xi}$ are updated using

for $k=1,\ldots,N_{a}$ and $i=1,\ldots,N_{r}$, where $N_{a}$ is the number of data assimilation steps, and $N_{r}$ is the number of ensemble realizations. In this work, we use $N_{a}=4$, and the inflation coefficients $\alpha_{k}$ are set to 4. And $C_{\xi\mathbf{d}}$ is the cross-covariance of latent variables $\boldsymbol{\xi}$ and data variables $\mathbf{d}$, and $C_{\mathbf{d}}$ is the data covariance matrix. These covariance matrices are updated at each assimilation step, corresponding to the updated ensemble. Random noise $\mathbf{e}$, sampled from a Gaussian distribution $N(\mathbf{0},C_{D})$, is added in each iteration for mismatch calculation. After all assimilation steps are completed, the posterior models $\mathbf{m}_{\text{post}}$ are generated using PCA on the posterior latent variables $\boldsymbol{\xi}_{\text{post}}$, and the posterior predictions $\mathbf{d}_{\text{post}}$ are generated using U-FNO on the posterior models.

Figure 4 presents the GeoFUSE framework for data assimilation in the seawater intrusion problem. First, latent variables $\boldsymbol{\xi}$ are sampled from a Gaussian distribution $N(\mathbf{0},I)$ to generate a prior ensemble with $N_{r}$ samples. PCA is then applied to the latent variables $\boldsymbol{\xi}$ to generate prior geomodels $\mathbf{m}$. The surrogate model is then applied to generate flow responses $\mathbf{d}$. ESMDA is performed, conditioned on the observations $\mathbf{d}_{\text{obs}}$, to calibrate the latent variables $\boldsymbol{\xi}$. In each iteration, PCA and U-FNO are applied to generate pressure and salinity data that match the observations. After $N_{a}$ iterations, we generate posterior models $\mathbf{m}_{\text{post}}$ and posterior predictions $\mathbf{d}_{\text{post}}$.

We simulate 1600 prior realizations to collect pressure and salinity data at 243 time steps, with measurements taken every 30 days. Of these, 1500 realizations are used to train the U-FNO-based surrogate model, while the remaining 100 are reserved for evaluating its predictive performance. The U-FNO architecture, comprising over 126 million parameters, is employed to train separate models for pressure and salinity outputs (Table 1).

The training process, conducted on a single Nvidia Tesla A100 GPU, takes 6.5 hours over 200 epochs for each network. In comparison, a full PFLOTRAN forward simulation required approximately 1.5 hours with one CPU. However, the U-FNO surrogate model predictions are computed in just 0.015 seconds, resulting in a computational speedup of around 360,000 times.

To assess the predictive accuracy of the U-FNO surrogate model, we compare its results with the PFLOTRAN simulations. Figures 5 and 6 demonstrate the close agreement between the predicted and simulated salinity distributions over time and across multiple realizations in Fig. 2. The top row presents PFLOTRAN results, the middle row shows U-FNO predictions, and the bottom row highlights the absolute errors between the two. The comparisons are made for three key time points: 1050 days, representing the early-stage (approximately 3 years), 3570 days, marking the mid-stage (approximately 10 years), and 7170 days, indicating the late-stage (approximately 20 years), clearly illustrating the salinity dynamics over the 20-year period. As shown in Fig. 5, the U-FNO model effectively captures both the temporal and spatial dynamics of salinity. The absolute error for this realization falls in the middle range of the 100 test samples, demonstrating the model’s robustness in approximating complex flow and transport processes.

Figure 6 presents the variability in salinity plume size and shape across three different permeability realizations. This variability highlights the challenges of modeling heterogeneous systems, yet U-FNO successfully replicates these dynamics. The surrogate model consistently delivers high accuracy in its predictions, with the errors between U-FNO and PFLOTRAN being notably small, as shown in the bottom rows showing errors.

In addition to salinity distribution, Fig.  7 compares salinity accumulation between the U-FNO surrogate predictions and PFLOTRAN simulations for 100 test realizations. The results show the 10^th, 50^th, and 90^th percentiles (P₁₀, P₅₀, and P₉₀) over a period of 7,260 days. The black curves represent the PFLOTRAN results, while the blue dashed curves represent the U-FNO predictions. Over the 20-year period, the consistent increase in salinity accumulation and the close agreement between the surrogate model and full-scale simulations demonstrate the accuracy of U-FNO in approximating these long-term cumulative effects.

The development of the U-FNO model represents a significant advancement in computational efficiency and accuracy for subsurface flow simulations. It achieves a remarkable speedup of approximately 360,000 times over traditional simulation models like PFLOTRAN, offering substantial benefits for large-scale simulations where computational resources and time are limited Wen et al. (2022). This efficiency is especially critical for scenarios requiring extensive Monte Carlo sampling for uncertainty quantification. As noted by Zhong et al. (2019), their cDC-GAN model similarly reduced simulation time, enabling numerous simulations for uncertainty analysis. However, U-FNO faces challenges in accurately predicting saltwater front movements, where its accuracy diminishes. This issue is not unique to U-FNO, as convolutional neural networks (CNN)-based methods Mo et al. (2019); Tang et al. (2020) have exhibited similar discrepancies when approximating multiphase flow behaviors in heterogeneous systems. Despite these challenges, U-FNO’s overall predictive accuracy remains impressive, making it a transformative tool for large-scale, data-intensive studies.

The GeoFUSE framework substantially accelerates the subsurface flow modeling and data assimilation process while significantly reducing the uncertainty in geological models and long-term predictions, particularly for seawater intrusion and groundwater management in coastal areas. This approach aligns with previous studies that emphasize the importance of addressing geological uncertainty to enhance predictive capabilities Sonnenborg et al. (2015); He et al. (2015). Within the GeoFUSE framework, PCA is used to reduce the dimension of geological models while preserving essential features, enabling more efficient uncertainty quantification for large-scale Gaussian models. Additionally, the ability of ESMDA to assimilate sparse and noisy data makes the framework both simple and robust for subsurface flow predictions. The results, which show improved salinity distribution predictions over a 20-year simulation period, further highlight the framework’s effectiveness in model calibration and uncertainty reduction. Overall, GeoFUSE provides a transformative approach to subsurface flow simulations and uncertainty quantification in complex coastal environments by accelerating flow and transport simulations, efficiently calibrating geological models, and reducing the uncertainty of future predictions.

GeoFUSE offers a powerful tool for decision-makers involved in floodplain management, tidal restoration, and salinity control by providing rapid and accurate predictions and model calibrations. This capability enhances water resource management in the face of rising sea levels and climate change. By offering predictive modeling for salinity intrusion in vulnerable coastal regions, GeoFUSE plays a critical role in broader coastal management efforts. Its ability to deliver timely, accurate insights makes it essential for adaptive management, especially as seawater intrusion intensifies due to climate change. This allows decision-makers to respond swiftly to evolving conditions in coastal aquifers.

The integration of monitoring data into predictive models, as demonstrated in this study, significantly reduces model uncertainty and improves the management of coastal aquifers. GeoFUSE’s use of a surrogate model reduces computational costs, enabling rapid predictions for uncertainty quantification and decision-making. By assimilating monitoring data, the framework improves salinity forecasts by reducing uncertainty in the geological models. The model’s performance, validated against PFLOTRAN simulation results and prior statistics, confirms its predictive accuracy and robustness in model calibration.

GeoFUSE has demonstrated substantial potential in subsurface flow modeling and uncertainty quantification, but there is room for further improvement and scalability. One area of focus is refining the U-FNO model, which has effectively captured complex multiphase flow dynamics Wen et al. (2022). However, further adjustments could enhance its ability to account for fine-scale geological variations, which are crucial in heterogeneous environments. Additionally, incorporating more complex factors, such as boundary conditions influenced by climate change, could increase the model’s robustness. Climate-induced variability significantly impacts hydrological cycles and groundwater dynamics He et al. (2015), and addressing these complexities would enable GeoFUSE to more accurately simulate the intricate interactions between geological formations and hydrological processes across diverse environments.

Future research could also focus on extending GeoFUSE to three-dimensional (3D) models, which are essential for accurately simulating subsurface flow in coastal regions increasingly affected by seawater intrusion. Recent studies highlight the importance of 3D geological modeling in understanding fluid dynamics and permeability in geothermal systems Siler et al. (2019); Torresan et al. (2020). A 3D framework would better integrate geological, geophysical, and geochemical data, providing a more comprehensive approach to evaluating both anthropogenic and natural influences on subsurface flow Strati et al. (2017); Dwivedi et al. (2018). By improving scalability and addressing these key areas, GeoFUSE could evolve in a versatile tool for groundwater management in diverse environments, particularly in coastal regions facing challenges such as groundwater depletion and climate change Messier et al. (2015); Ebong et al. (2020).