Physics-Informed DeepONet Coupled with FEM for Convective Transport in Porous Media with Sharp Gaussian Sources

TL;DR

This paper introduces a hybrid FEM and physics-informed DeepONet framework for modeling convective transport in porous media with sharp Gaussian sources, achieving accurate and fast predictions.

Contribution

It combines FEM with DeepONet to efficiently model fluid transport, incorporating adaptive sampling for steep gradients, and demonstrates significant speedups over traditional methods.

Findings

Good agreement with reference solutions

Orders of magnitude speedup

Effective handling of sharp source gradients

Abstract

We present a hybrid framework that couples finite element methods (FEM) with physics-informed DeepONet to model fluid transport in porous media from sharp, localized Gaussian sources. The governing system consists of a steady-state Darcy flow equation and a time-dependent convection-diffusion equation. Our approach solves the Darcy system using FEM and transfers the resulting velocity field to a physics-informed DeepONet, which learns the mapping from source functions to solute concentration profiles. This modular strategy preserves FEM-level accuracy in the flow field while enabling fast inference for transport dynamics. To handle steep gradients induced by sharp sources, we introduce an adaptive sampling strategy for trunk collocation points. Numerical experiments demonstrate that our method is in good agreement with the reference solutions while offering orders of magnitude speedups over traditional solvers, making it suitable for practical applications in relevant scenarios. Implementation of our proposed method is available at https://github.com/erkara/fem-pi-deeponet.

Tables

Table 1. Table 1 : Physical parameters used in the Darcy and convection–diffusion PDE system.

Parameter	Value	Unit	Explanation
$\Omega$	$[0,10]\times [0,10]$	mm	Spatial domain
K	${10}^{-9}$	mm2	Intrinsic permeability
$\mu$	$9\times {10}^{-4}$	Pa $\cdot$ s	Viscosity of water
$\alpha$	${10}^{-2}$	-	Sink term
D	$4\times {10}^{-6}$	mm2/s	Diffusion coefficient
${\beta }_1$	$\frac{5}{60}$	mm2/s	Fluid addition rate
${\beta }_2$	$\frac{5}{240}$	mmol/s	Concentration addition rate

Table 2. Table 2 : Sampling configurations used in the PI-DeepONet training pipeline.

Parameter	Value	Explanation
$N$	2000	Number of Gaussian input functions sampled
$m$	$30\times 30$	Number of branch sensor points (grid)
$P_{\mathrm{bcs}}$	100	Number of boundary condition (BC) collocation points
$P_{\mathrm{ics}}$	5	Number of initial condition (IC) collocation points
$n_r$	30	Radial points around the source ($r\le 3\sigma$)
$n_{\theta }$	30	Angular points around the source
$n_{\mathrm{rand}}$	300	Number of random collocation points
$N_{test}$	30	Number of input functions used for model testing
$k$	$80\times 80$	dimensions of the testing grid

Table 3. Table 3 : Relative L 2 errors for single-source test cases. Here, ℰ full \mathcal{E}_{\text{full}} denotes the average error across all time steps, and ℰ T \mathcal{E}_{T} denotes the error at final time t = T t=T .

Final Time	${ℰ}_{\text{full}}$	${ℰ}_T$
$T=50$	9.99%	1.97%
$T=250$	10.06%	2.72%
$T=500$	10.08%	2.86%

Table 4. Table 4 : Relative L 2 errors for the case where we randomly sample one or two sources with fixed width. Here, ℰ full \mathcal{E}_{\text{full}} denotes the average error across all time steps, and ℰ T \mathcal{E}_{T} denotes the error at final time t = T t=T .

Final Time	${ℰ}_{\text{full}}$	${ℰ}_T$
$T=50$	9.73%	1.90%
$T=250$	8.86%	1.93%
$T=500$	10.06%	3.19%

Table 5. Table 5 : Relative L 2 errors for the case where we have triple source fixed width. Here, ℰ full \mathcal{E}_{\text{full}} denotes the average error across all time steps, and ℰ T \mathcal{E}_{T} denotes the error at final time t = T t=T .

Final Time	${ℰ}_{\text{full}}$	${ℰ}_T$
$T=50$	12.74%	8.92%
$T=250$	12.48%	9.44%
$T=500$	14.22%	10.52%

Table 6. Table 6 : Relative L 2 errors for the case where we have triple source with fixed centers and varying width. Here, ℰ full \mathcal{E}_{\text{full}} denotes the average error across all time steps, and ℰ T \mathcal{E}_{T} denotes the error at final time t = T t=T .

Final Time	${ℰ}_{\text{full}}$	${ℰ}_T$
$T=50$	9.81%	2.19%
$T=250$	9.91%	2.81%
$T=500$	10.05%	3.12%

Table 7. Table 7 : Relative L 2 errors for two difference experiments with permeability K = 10 − 8 K=10^{-8} . ℰ full \mathcal{E}_{\text{full}} denotes the average error across all time steps, and ℰ T \mathcal{E}_{T} denotes the error at final time t = 500 t=500 .

Experiment (T=500)	${𝓔}_{\text{full}}$	${𝓔}_{𝑻}$
Single Source with Varying Center and Width	10.84	5.70
Triple Source with Fixed Centers and Varying Width	10.56	4.04

Table 8. Table 8 : Ablation study evaluating the impact of residual point sampling strategy on model performance. For each experiment, the top row shows results using hybrid sampling (structured polar grid around the source plus random points), while the bottom row uses only randomly sampled residual points. The total number of residual points is kept fixed across both configurations.

Experiment	$n_r$	$n_{\theta }$	$n_{\mathrm{rand}}$	${ℰ}_{\text{full}}$	${ℰ}_T$
Single Source with Varying Center and Width	30	30	300	10.70%	4.90%
	0	0	1200	14.76%	14.10%
Triple Source with Fixed Width	30	30	300	12.59%	8.92%
	0	0	3000	31.11%	29.10%
Triple Source with Fixed Centers and Varying Width	30	30	300	10.33%	3.73%
	0	0	3000	36.53%	38.32%

Table 9. Table 9 : Average execution times (in seconds) for FEM and PI-DeepONet over 30 test cases per experiment.

Experiment (T=500)	FEM (s)	Training Time (min)	Model (s)	Speed-up
Single Source with Varying Center and Width	5.265	84.53	0.038	138.399
Random Sampling of One or Two Sources with Fixed Width	6.564	83.72	0.039	168.713
Triple Source with Fixed Width	9.081	102.13	0.031	289.739
Triple Source with Fixed Centers and Varying Width	11.809	102.04	0.036	328.153

Equations

$$\mathbf{v}$$

$$\nabla \cdot \mathbf{v}$$

$$\frac{\partial c}{\partial t}$$

$$p$$

$$\mathbf{n}\cdot \nabla c$$

$$f(\mathbf{x})=\frac{1}{\mathcal{N}}\exp \left(-\frac{\Vert \mathbf{x}-{\mathbf{x}}_0{\Vert }^2}{2{\sigma }^2}\right),\mathcal{N}={\int }_{\Omega }\exp \left(-\frac{\Vert \mathbf{x}-{\mathbf{x}}_0{\Vert }^2}{2{\sigma }^2}\right)d\mathbf{x}.$$

$$\mathbf{V}$$

$$Q$$

$$(\mathbf{v},\mathbf{w})+\left(p,\nabla \cdot \mathbf{w}\right)$$

$$(\nabla \cdot \mathbf{v},q)+(\alpha p,q)$$

$$\left(\frac{c^{n+1}-c^n}{\Delta t},\varphi \right)-(D\nabla c^{n+1},\nabla \varphi )+(\mathbf{v}{c}^{n+1},\nabla \varphi )+\tau (\mathbf{v}\cdot \nabla \varphi ,R^{n+1}(c))=({\beta }_{2}f,\varphi ),$$

$$R^{n+1}(c):=\frac{c^{n+1}-c^n}{\Delta t}+\nabla \cdot (\mathbf{v}{c}^{n+1})-\nabla \cdot (D\nabla c^{n+1})-{\beta }_{2}f,$$

$$r_i=\max \left(k_i{\sigma }_i,0.1\cdot \min (L_x,L_y)\right),{\delta }_i=\frac{r_i}{2},$$

$$k_i=k_0+k_1\cdot \frac{{\sigma }_{\text{ref}}-{\sigma }_i}{{\sigma }_{\text{ref}}},$$

$${\psi }_i(x,y):=\frac{1}{2}\left(1+\tanh \left(\frac{\Vert (x,y)-{\mathbf{x}}_i\Vert -r_i}{{\delta }_i}\right)\right),$$

$$h(x,y)=h_{\min }+(h_{\max }-h_{\min })\cdot \prod _{i=1}^N{\psi }_i(x,y),$$

$$f\mapsto c,$$

$$\mathcal{G}:\mathcal{A}\to \mathcal{B},c(\mathbf{x},t)=\mathcal{G}(f)(\mathbf{x},t),(\mathbf{x},t)\in \Omega \times [0,T],$$

$${\mathcal{G}}_{\theta }(f)(\mathbf{x},t)=⟨B_{{\theta }_b}(f(z_1),\dots ,f(z_m)),T_{{\theta }_t}(\mathbf{x},t)⟩,$$

$${\mathcal{Y}}_r={\mathcal{Y}}_{\text{source}}\cup {\mathcal{Y}}_{\text{rand}}\cup {\mathcal{Y}}_{\text{bcs}}\cup {\mathcal{Y}}_{\text{ics}},$$

$$\mathcal{Y}_{\text{source}}=\left\{(\mathbf{x}_{i},t_{i})\,\big{|}\,\mathbf{x}_{i}=\mathbf{x}_{0}+r_{j}(\cos\theta_{\ell},\sin\theta_{\ell}),\;r_{j}\in[0,3\sigma],\;\theta_{\ell}\in[0,2\pi),\;t_{i}\sim\mathcal{U}[0,T]\right\},$$

$${\mathcal{Y}}_{\mathrm{rand}}=\{({\mathbf{x}}_i,t_i)\mid {\mathbf{x}}_i\sim \mathcal{U}(\Omega ),t_i\sim \mathcal{U}[0,T]\},|{\mathcal{Y}}_{\mathrm{rand}}|=n_{\mathrm{rand}}.$$

$${\mathcal{Y}}_{\mathrm{bcs}}=\{({\mathbf{x}}_j,t_j)\mid {\mathbf{x}}_j\sim \mathcal{U}(\partial \Omega ),t_j\sim \mathcal{U}[0,T]\},|{\mathcal{Y}}_{\mathrm{bcs}}|=P_{\mathrm{bcs}},$$

$${\mathcal{Y}}_{\mathrm{ics}}=\{({\mathbf{x}}_j,0)\mid {\mathbf{x}}_j\sim \mathcal{U}(\Omega )\},|{\mathcal{Y}}_{\mathrm{ics}}|=P_{\mathrm{ics}}.$$

$${\mathbf{b}}^{(n)}=\left(f^{(n)}(z_1),\dots ,f^{(n)}(z_{m^2})\right)\in {\mathbb{R}}^{m^2}.$$

$${\mathcal{D}}_{\text{train}}={\left\{\left({\mathbf{b}}^{(n)},{\mathcal{Y}}_r^{(n)}\right)\right\}}_{n=1}^N.$$

$${\mathcal{L}}_{\text{total}}(\theta )$$

$${\mathcal{L}}_{\text{res}}(\theta )=\frac{1}{|{\mathcal{Y}}_{\text{res}}|}\sum _{(\mathbf{x},t)\in {\mathcal{Y}}_{\text{res}}}{\mathcal{R}}_{\theta }(\mathbf{x},t{)}^2,$$

$${\mathcal{R}}_{\theta }(\mathbf{x},t)=\frac{\partial c_{\theta }}{\partial t}(\mathbf{x},t)-\nabla \cdot (D\nabla c_{\theta }(\mathbf{x},t))+\nabla \cdot (v(\mathbf{x})c_{\theta }(\mathbf{x},t))-{\beta }_{2}f(\mathbf{x})$$

$${\mathcal{L}}_{\text{bcs}}(\theta )=\frac{1}{|{\mathcal{Y}}_{\text{bcs}}|}\sum _{(\mathbf{x},t)\in {\mathcal{Y}}_{\text{bcs}}}{\left(\mathbf{n}\cdot \nabla c_{\theta }(\mathbf{x},t)\right)}^2,$$

$${\mathcal{L}}_{\text{ics}}(\theta )=\frac{1}{|{\mathcal{Y}}_{\text{ics}}|}\sum _{\mathbf{x}\in {\mathcal{Y}}_{\text{ics}}}{c}_{\theta }(\mathbf{x},0{)}^2.$$

$${\mathcal{E}}_{\text{full}}$$

$${\mathcal{E}}_T$$

Full text

1 Introduction

Fluid flow and solute transport in porous media play a critical role in many scientific and engineering applications, particularly in areas such as subsurface energy storage, carbon sequestration, geothermal systems, and groundwater contamination [1, 2, 3, 4, 5, 6]. These processes are often modeled using coupled systems of partial differential equations (PDEs) that describe fluid pressure, velocity, solute concentration and other variables of interest over time and space. In practical scenarios, these governing equations depend on many parameters, including source characteristics, medium properties, and boundary conditions, which can vary significantly across applications [7, 8].

Classical numerical methods such as FEM has been the gold standard in the qualitative and numerical study of the underlying governing PDEs [8, 9, 10]. By construction, these PDEs often naturally involve several parameters such as problem geometry, initial and boundary conditions as well as physical parameters representing permeability, porosity, diffusivity and so on. One of the fundamental limitation of conventional methods is that they usually solve one instance of the underlying PDE, thus limiting their practical use in applications requiring fast inference since every parameter configuration requires to run the solver again.

In recent years, physics-informed machine learning has emerged as a promising avenue to overcome these limitations by constructing fast surrogate models for PDEs. Physics-Informed Neural Networks (PINNs) were introduced as deep learning frameworks that embed physical laws (e.g. underlying PDEs) into the training objective [11]. By penalizing the PDE residual in the loss function, PINNs can learn the solution of a given differential equation without requiring extensive labeled data, yielding a surrogate that obeys the underlying physics. This approach has been successfully applied to forward and inverse problems for fluid flow and other nonlinear PDE systems [12, 13]. However, a standard PINN is typically trained for one specific problem setup at a time (with fixed domain and parameters), which can limit its reusability across varying inputs.

To address this limitation, neural operator learning has emerged as a powerful alternative, enabling the learning of solution mappings across varying domain geometries, input parameters, and boundary conditions. Unlike PINNs, which approximate a single solution, neural operators learn mappings between function spaces, effectively capturing the solution operator of an entire parametric family of PDEs [14, 15, 16]. Although they require upfront training, neural operators allow fast inference across diverse scenarios without repeatedly solving the PDE from scratch, with ideally minimal loss of accuracy.

Among these approaches, one of the most well-known and successful one is DeepOnet framework Lu et al [17]. DeepONet uses a dual-network architecture (with a “branch” net encoding input functions, and a “trunk” net encoding the output coordinates) to approximate nonlinear operators, supported by a universal approximation theorem for operators [18]. This enables learning an entire family of PDE solutions: for instance, mapping different source terms or boundary conditions to the corresponding solution fields [19, 15]. In fact, these surrogates have been reported to produce solutions up to 1000x faster than standard FEM-based PDE solvers while maintaining good accuracy. Literature is rich with several extensions of the original DeepOnet, such Bayesian DeepOnet [20], multi-fidelity DeepOnet [21], multi input-output DeepOnet [19].

In its original form, DeepONet is a supervised machine learning model and therefore requires input–output pairs. This implies that one must obtain such solutions from the underlying PDEs. Even for moderately sized problems, this involves collecting a large amount of high-fidelity data from numerical solvers, which can be very expensive. To address these challenges, Wang et al. incorporated physical constraints into the training process, giving rise to physics-informed DeepONets (PI-DeepONets), which use PDE residuals instead of supervised solution labels [22]. The fundamental idea relies on the fact that DeepONet outputs are differentiable with respect to input parameters (i.e., trunk and branch inputs), allowing for PINN-style physics enforcement in DeepONet via automatic differentiation [23]. The implication is that the network can learn operator mappings without requiring paired input–output datasets. The same team also introduced a modified DeepONet architecture designed to mitigate vanishing gradient issues and improve convergence in physics-informed operator learning via NTK theory [24]. In this architecture, authors also augmented the conventional branch and trunk networks of original DeepOnet with encoder layers that embed the input function and evaluation coordinates, respectively, into higher-dimensional latent spaces. These embeddings are then injected into each hidden layer via pointwise operations, enabling hierarchical interaction between the networks throughout the forward pass. This architectural modification facilitates more stable signal propagation and has been empirically shown to outperform conventional DeepONets, particularly in the physics-informed setting without paired training data. This work forms the foundation of our work in this paper that will be detailed in the subsequent sections.

In the study of porous media flow, several data-driven studies have been proposed in recent years. A recent review by Yang et al. provides a broad overview of such methods [1]. In the context of physics-guided machine learning, there is a relatively rich literature applying PINNs to porous media problems, motivated by the need for fast surrogates in real-world applications.

Tartakovsky et al. developed a PINN framework for inferring spatially varying parameters and unknown constitutive relations in subsurface flow, demonstrating accurate recovery from sparse data under Darcy and Richards models [25]. Jagtap and Karniadakis introduced XPINNs, a space–time domain decomposition approach where separate PINNs are deployed in irregular subdomains to improve parallelism and flexibility for complex PDEs [26]. Yan et al. proposed a gradient-based PINN (GDNN) that regularizes training through physically motivated differential operators; applied to CO2 injection scenarios, the model captures nonlinear multiphase flow and extends naturally to convolutional architectures [27]. Zhang et al. designed a physics-informed deep convolutional network (PIDCNN) for transient Darcy flow, combining finite volume loss enforcement with a well model to capture steep pressure gradients [28]. Cai et al. developed a PINN framework for Biot’s poroelasticity equations using a fixed-stress splitting strategy to decouple pressure and displacement, with Darcy flow as the governing mechanism for fluid motion [29]. For additional PINN applications in porous media, see also [30, 31, 32].

Compared to PINNs, operator learning studies in porous media are more limited. Huang et al. introduced the Porous-DeepONet framework, which augments DeepONet with convolutional networks to encode image-based porous structures [33]. They also proposed Porous-PI-DeepONet for physics-informed training and Porous-DeepM&Mnet for multi-physics coupling, demonstrating strong generalization across various reactive transport problems. Yan et al. employed Fourier Neural Operators to model multiphase flow under heterogeneous permeability, showing accurate predictions and speedups in CO2 sequestration scenarios [34]. Kumar et al. introduced Multi-task DeepONet (MT-DeepONet), which solves diverse PDEs across varying source terms and geometries using a modified branch network and binary masking [35]. Mandl et al. proposed Physics-Informed Separable DeepONet (PI-Sep-DeepONet), which factorizes the trunk network for efficient operator learning in high-dimensional problems, and applied it to Biot’s consolidation system, a Darcy-type flow setting without transport [36].

In this paper, we study the problem of fluid infusion into a porous medium from a sharp, localized Gaussian source. The governing model consists of three coupled PDEs: a steady-state Darcy system that determines the pressure and velocity fields, and a time-dependent convection–diffusion equation that models the transport of the infused fluid. Our objective is to learn the mapping from the source function to the resulting concentration profile. While operator learning has been applied to porous media problems in various forms, our work introduces several key innovations.

First, we propose a modular solution strategy in which the Darcy flow problem is solved using a finite element method (FEM), and only the transport component is learned using a physics-informed DeepONet. Unlike most existing approaches that treat the entire PDE system as a single learning task, we couple high-fidelity numerical solvers and operator learning in a sequential and physically consistent way. This allows us to preserve FEM-level accuracy in the velocity field while benefiting from the fast inference and generalization properties of DeepONet in the transport stage. Notably, despite the growing interest in operator learning, the use of PI-DeepONet, particularly in porous media applications, remains extremely rare in the current literature.

Second, to handle the sharp spatial localization of the Gaussian sources, we introduce a simple yet effective adaptive sampling strategy for trunk collocation points. This improves the resolution of steep gradients in the solution, which would otherwise challenge standard PINN or DeepONet training pipelines. To the best of our knowledge, this work is one of the few that applies physics-informed DeepONet to porous media transport problems and the first to combine it in a tightly coupled manner with FEM solvers.

2 Problem Description

2.1 Convective Fluid Transport in Porous Media

In this section, we describe the coupled systems of partial differential equations(PDE) to model the delivery and transport of an externally infused fluid carrying a dissolved solute into a porous medium. The system consists of two components: a steady-state Darcy flow and mass balance model that determines the pressure and velocity fields in the porous medium, and a time-dependent convection–diffusion equation that governs the evolution of the solute concentration.

[TABLE]

where $\mathbf{v}$ is the Darcy velocity, $p$ is the fluid pressure, and $c$ is the solute concentration. Here, $\mathbf{K}$ is the symmetric, positive-definite permeability tensor of the medium, $\mu$ is the dynamic viscosity, $D$ is the diffusion tensor while $f(x)$ is the spatially distributed source intensity, scaled by $\beta_{1}$ and $\beta_{2}$ for the fluid and solute, respectively.

We impose the following boundary conditions on the domain boundary $\partial\Omega$:

[TABLE]

where $\mathbf{n}$ denotes the outward unit normal vector.

The first equation in the system (1a) is Darcy’s law that models fluid motion at low Reynolds number, where viscous forces dominate and inertial effects can be neglected. Darcy’s law expresses that the fluid flux is proportional to the pressure gradient, modulated by the local permeability and viscosity. Second equation (1b) expresses the local mass conservation which accounts for two effects: fluid loss through pressure-dependent sinks, such as leakage or matrix absorption ($-\alpha p$), and external fluid infusion through the source term $\beta_{1}f(x)$. While the sink term is general to porous media flow, we note that a similar formulation appears in biology, where Starling’s law describes fluid exchange across semipermeable membranes as a balance of hydraulic and osmotic forces. Together, these equations form the system, which provides the stationary flow field and pressure distribution induced by the external infusion and internal leakage mechanisms.

The third equation (1c), which is of the standard convection-diffusion form, governs the transport of the dissolved solute. The first term on its right-hand side models diffusive or dispersive spreading of the solute, while the second term accounts for its advective transport by the fluid flow $\mathbf{v}$. The source term $\beta_{2}f(x)$ models direct injection of solute into the domain, coinciding spatially with the fluid infusion.

The boundary conditions are chosen to reflect physically meaningful constraints. The zero-pressure condition, $p=0$, assumes that the pressure has equilibrated at the domain boundary, meaning that the infused fluid has relaxed back to its environmental baseline pressure at the outer edge of the modeled region. This models an open system where pressure disturbances dissipate away from the injection site, and no artificial buildup of pressure occurs at the boundary. The no-flux condition on the velocity, $\mathbf{v}\cdot\mathbf{n}=0$, reflects the presence of impermeable boundaries that do not permit net fluid flow across the boundary, ensuring that the fluid exchange is entirely governed by the internal sources and sinks. For the solute, the homogeneous Neumann condition, $\mathbf{n}\cdot\nabla c=0$, enforces zero diffusive flux across the boundary, preventing solute from escaping the domain or being artificially introduced through the boundary.

In this work, we are primarily interested in fluid infusion driven by a localized source. To model this, we represent the source term $f(\mathbf{x})$ as a normalized Gaussian profile centered at a prescribed location $\mathbf{x_{0}}$

[TABLE]

where $\sigma$ characterizes the spatial width of the source and determines how sharply localized the injection is around $\mathbf{x_{0}}$. This formulation is relevant in many practical applications where infusion occurs through a focused delivery mechanism, such as a point source in contaminant transport modeling [37, 38, 39], hydraulic fracturing [40, 41] and targeted drug delivery in biomedical applications [42, 43, 44]. The Gaussian profile captures the physical intuition that the intensity of injection is highest at the source center and decays smoothly away from it in space.

2.2 Finite Element Formulation of the Darcy-Transport System

In this section, we will describe the numerical solution of the coupled PDE system. Here, we adopt a sequential strategy that takes advantage of the decoupled structure of the problem. Specifically, we first solve the steady-state Darcy flow system to compute the pressure and velocity fields, and then use the computed velocity as input to the time-dependent convection–diffusion equation governing solute transport. The finite element method (FEM) implementation is developed using the Firedrake framework [45]. Since the source term $f(\mathbf{x})$ can become highly localized for small values of $\sigma$, we employ an adaptive mesh refinement strategy to ensure sufficient spatial resolution near the source region.

We begin by formulating the steady-state Darcy system using a mixed finite element approach for the governing equations (1a)-(1b) system. Let us introduce two distinct trial and test spaces for $\mathbf{v}$ and $p$. Let

[TABLE]

where $H(\text{div};\Omega)$ denotes the space of square-integrable vector fields with square-integrable divergence. The homogeneous Neumann condition on the velocity field is imposed essentially via the definition of $\mathbf{V}$, while the Dirichlet condition on pressure is imposed naturally in the variational formulation.

The weak form then reads; find $(\mathbf{v},p)\in\mathbf{V}\times Q$ such that:

[TABLE]

To discretize the system, we select conforming finite-dimensional subspaces $\mathbf{V}_{h}\subset\mathbf{V}$ and $Q_{h}\subset Q$ such that the pair satisfies the Ladyzhenskaya–Babuška–Brezzi (LBB) condition. In our implementation, we use second-order Brezzi–Douglas–Marini (BDM2) elements for velocity, which provide $H(\text{div})$-conforming approximations with continuous normal components across element interfaces, and discontinuous piecewise linear (DG1) elements for pressure. This BDM2–DG1 pairing is a classical choice for mixed Darcy problems, known to ensure stability and optimal convergence [46, 47].

To generate reference solutions to test the proposed model, we discretize the transport equation defined in Equation-(1c) using a stabilized Galerkin finite element method [48, 49]. The Darcy velocity field $\mathbf{v}(x)$ is computed from the preceding system and treated as fixed input. We use the same adaptive mesh as that of the Darcy solver. Details of this mesh generation will be explained shortly.

Regarding the transport equation (1c), let $\{\phi\}\subset V$ denote the set of admissible test functions for $V:=H^{1}(\Omega)$. We consider a uniform temporal discretization with time step $\Delta t$, and approximate the solution $c(t,x)$ by $c^{n}(x)\approx c(t_{n},x)$ at discrete tims $\{t_{n}\}$ using the backward Euler method. We seek $c^{n+1}\in V$ such that for all $\phi\in V$,

[TABLE]

where the streamline residual is defined as

[TABLE]

Stabilization parameter is computed elementwise as $\tau=\frac{h}{2\|\mathbf{v}\|},$ where $h$ denotes the local cell diameter and $\|\mathbf{v}\|$ the magnitude of the advective velocity. The spatial discretization uses continuous piecewise linear (P1) finite elements over the shared adaptive mesh.

Lastly we will describe the adaptive mesh generation procedure. To accurately resolve the steep gradients induced by the localized source term $f(x)$, we generate an adaptively refined mesh that concentrates resolution near the centers of the Gaussian profiles. This is particularly important when the source width $\sigma$ is small, leading to highly localized source regions. The refinement procedure is implemented using Gmsh and designed to scale with the source configuration parameters [50].

Let $\{\mathbf{x}_{i}\}_{i=1}^{N}$ denote the centers of the Gaussian source components, with corresponding spatial widths $\{\sigma_{i}\}_{i=1}^{N}$. For each Gaussian, we define a target refinement radius $r_{i}$ and a transition width $\delta_{i}$ as:

[TABLE]

where $k_{i}$ is a scaling factor that increases as $\sigma_{i}$ decreases, thus enforcing finer resolution around sharper sources. Specifically, $k_{i}$ is defined in the form of:

[TABLE]

with $\sigma_{\text{ref}}$ a nominal reference scale and $k_{0},k_{1}$ tunable constants. We set $\sigma_{\text{ref}}=1e-3,k_{0}=k_{1}=10$.

To smoothly transition between fine and coarse mesh regions, we define a background mesh size field $h(x,y)$ over the domain using a product of smooth decay functions centered at each source. For each source, we define:

[TABLE]

and assemble the full mesh size field as:

[TABLE]

where $h_{\min}$ is a user-defined minimum mesh size and $h_{\max}$ is the maximum allowable element size. Our aim with this construction is to have $h(x,y)\approx h_{\min}$ near each Gaussian center and $h(x,y)\approx h_{\max}$ far from all sources while having smooth variation in between, avoiding abrupt transitions. In our implementation, we set $h_{\min}=\min_{i}\frac{\sigma_{i}}{6}$ and $h_{max}=0.04\max(L_{x},L_{y})$.

Figure 1 illustrates the complete FEM pipeline, from adaptive mesh generation to the solution of the coupled Darcy–Transport system (1a), (1b), (1c), for representative single- and triple-source cases.

3 Methodologies

3.1 Proposed Method

In this section, we outline our FEM–DeepONet coupling strategy and explain the mathematical and computational motivations behind its design. In principle, one could formulate the learning task as the operator mapping $f\mapsto(\mathbf{v},p,c)$, seeking to predict the full velocity, pressure, and concentration fields from the source term. However, this approach presents both mathematical and computational challenges. First, the velocity and pressure fields are governed by the steady-state Darcy system, which is elliptic and typically well-conditioned. These equations can be solved efficiently and with high accuracy using a suitable numerical method(e.g. finite element method), as they converge rapidly and incur relatively low computational cost compared to the time-dependent transport problem. Consequently, attempting to learn the mapping $f\mapsto(\mathbf{v},p)$ offers limited computational benefit.

More critically, the system exhibits a natural decoupling: the velocity field $\mathbf{v}$ enters the transport equation as an advective term, meaning that even small perturbations in the predicted $\mathbf{v}$ or $p$ fields can propagate nonlinearly into the concentration $c$, significantly amplifying errors. Learning the full coupled map $f\mapsto(\mathbf{v},p,c)$ therefore becomes a considerably harder problem, both in terms of expressivity and stability. To address these issues, we propose a hybrid framework: we first solve the steady-state Darcy system using FEM to obtain accurate estimates of $\mathbf{v}$ and $p$, and then use this computed velocity field to drive the time-dependent convection–diffusion equation. We train DeepONet to learn the reduced mapping $f\mapsto c$, where the input incorporates both the source term $f$ and the advective field $\mathbf{v}$. This approach is particularly effective because the time-dependent transport equation constitutes the main computational bottleneck, while the steady Darcy system is comparatively cheap to solve numerically.

To our knowledge, such a hybrid coupling of FEM solvers with operator-learning models like DeepONet has received little attention in the current literature. Beyond the clear computational advantages, this coupling offers several important benefits. First, it provides a modular framework, where the Darcy system can be solved with a high-fidelity, physics-based solver and seamlessly combined with a data-driven model for the more demanding transport component. This modularity allows the method to generalize more easily across heterogeneous media or varying flow conditions since these will be handled by FEM solver. Second, by explicitly solving the steady-state flow problem, the approach retains physical interpretability and numerical error control, avoiding the black-box nature of end-to-end learning and ensuring that discretization errors in the flow field can be systematically monitored. Third, the learning task becomes significantly more data-efficient, as the network no longer needs to implicitly learn the physics of fluid flow but can focus its capacity on approximating the transport dynamics. Finally, framework offers transferability across applications: by adapting only the flow solver component, the overall system can be applied to a wide range of scientific and engineering problems.

The second key aspect of our study concerns the choice of spatial sampling strategy used during training, particularly on the trunk network side of the DeepONet architecture. In much of the recent operator learning literature, the training inputs are drawn from smooth function spaces, often using Gaussian random field (GRF) priors, and the solution is learned over a set of spatial points that are typically randomly sampled from a uniform distribution. This approach has been successful in many applications because the input functions are smooth and the solution varies relatively smoothly across the domain, allowing standard sampling from a uniform distribution to capture the main features of the solution.

However, in our setting, the source term $f(\mathbf{x})$ is modeled as a sharply localized Gaussian, representing a focused infusion into a porous medium. This leads to solution fields, particularly the concentration $c(\mathbf{x},t)$, that exhibit steep gradients near the source, with fine-scale features that are poorly captured by uniform trunk sampling. As we will demonstrate in later sections, the classical random uniform sampling approach fails to resolve these sharp features. To address this challenge, we propose a simple yet effective adaptive trunk sampling strategy, which places a denser set of evaluation points near the regions of high variation induced by the Gaussian source. This ensures that the DeepONet accurately captures the local solution behavior where it is most sensitive, while avoiding unnecessary oversampling in smooth regions. In addition, we emphasize that unlike many prior studies, which rely on discrete input samples interpolated onto a fixed branch sensor grid, our Gaussian source functions are known analytically and can be evaluated directly on any desired grid. This eliminates the need for a separate mapping step on the branch side and allows us to focus our adaptive sampling efforts entirely on the trunk network.

3.2 FEM-PI-DeepOnet Coupling with Adaptive Sampling

We now describe the complete workflow that integrates finite element computations with physics-informed DeepONet to construct the operator mapping

[TABLE]

where $f$ is the source function and $c$ is the resulting solute concentration field. As outline above, The steady-state Darcy system is solved using FEM to obtain the velocity field $v$, while the time-dependent convection–diffusion equation is learned by a PI-DeepONet model.

Let $\mathcal{A}$ and $\mathcal{B}$ be Banach spaces of input and output functions, respectively, and let

[TABLE]

denote the operator mapping of interest. DeepONet approximates $\mathcal{G}$ using a neural network of the form:

[TABLE]

where:

• •

$\{z_{i}\}_{i=1}^{m}\subset\Omega$ are fixed sensor points for evaluating the input function $f$,

• •

$B_{\theta_{b}}:\mathbb{R}^{m}\to\mathbb{R}^{q}$ is the branch network,

• •

$T_{\theta_{t}}:\Omega\times[0,T]\to\mathbb{R}^{q}$ is the trunk network,

• •

$\langle\cdot,\cdot\rangle$ denotes the Euclidean inner product in $\mathbb{R}^{q}$,

• •

$\theta=(\theta_{b},\theta_{t})$ denotes all trainable parameters.

In its original supervised form, DeepONet is trained using labeled pairs $\{(f^{(n)},c^{(n)}(x))\}_{n=1}^{N}$. However, in the physics-informed setting, we exploit the differentiability of DeepONet outputs with respect to trunk inputs to impose PDE constraints directly, allowing unsupervised training without solution labels. This yields a Physics-Informed DeepONet (PI-DeepONet) that minimizes a composite residual loss over a collocation set in space–time [22]. We will describe this next.

To evaluate the PDE residuals during training, we construct a specialized set of trunk collocation points in space–time. The total collocation set is denoted by

[TABLE]

where each subset serves a specific role.

Source-Centered Residual Points: To resolve the steep gradients induced near the source, we place residual trunk points more densely in its vicinity using a polar sampling scheme:

[TABLE]

This structured sampling yields $|\mathcal{Y}_{\mathrm{source}}|=n_{r}\cdot n_{\theta}$ where $n_{r}$ and $n_{\theta}$ denote the number of radial and angular subdivisions

Random Residual Points: To ensure broader spatial-temporal coverage, we include a uniformly random set over the domain:

[TABLE]

Boundary and Initial Condition Points: In addition to interior residual points, we generate collocation points to enforce Neumann boundary and initial conditions:

[TABLE]

[TABLE]

Figure-2 illustrates these three subsets for a representative Gaussian sampled.

The training process proceeds as follows. For each training instance:

A Gaussian source function $f$ is sampled by drawing the center $\mathbf{x_{0}}$ and width $\sigma$. 2. 2.

The steady-state Darcy system is solved using FEM on an adaptive mesh refined around $x_{0}$, yielding the velocity field $v(\mathbf{x})$. 3. 3.

This velocity field and its derivatives are interpolated onto the union $\mathcal{Y}_{\text{source}}\cup\mathcal{Y}_{\text{rand}}$, and used in evaluating the PDE residuals. 4. 4.

The source function is evaluated at a fixed set of branch sensor points $\{z_{i}\}_{i=1}^{m^{2}}$ to construct the input vector:

[TABLE] 5. 5.

Each point $(\mathbf{x}_{j},t_{j})\in\mathcal{Y}_{r}^{(n)}$ is paired with $\mathbf{b}^{(n)}$ to evaluate the network output $c_{\theta}(\mathbf{x}_{j},t_{j})$ and its derivatives.

Here, we used $c_{\theta}(\mathbf{x},t):=\mathcal{G}_{\theta}(f)(\mathbf{x},t)$ for simplicity. With this set-up, each training instance is thus a tuple $\left(\mathbf{b}^{(n)},\mathcal{Y}_{r}^{(n)}\right)$, and the full training dataset can we written as:

[TABLE]

PI-DeepONet is constructed in a similar manner to PINNs. In this sense, a composite loss function can be defined as follow

[TABLE]

with fixed weights $\lambda_{\text{res}},\lambda_{\text{bcs}},\lambda_{\text{ics}}$ chosen empirically.

The first term corresponds to the residual of the governing PDE:

[TABLE]

where

[TABLE]

is the residual associated with the convection–diffusion equation (1c).

The second term penalizes violations of the homogeneous Neumann boundary condition:

[TABLE]

where $\mathbf{n}$ denotes the outward unit normal at the boundary point $\mathbf{x}$, computed analytically from the rectangular domain geometry (e.g., $\mathbf{n}=(1,0)$ on the right edge, $\mathbf{n}=(0,-1)$ on the bottom, etc.).

The third and final component enforces the initial condition:

[TABLE]

We should note that, in [24], authors incorporated a weighting strategy inspired by Neural Tangent Kernel (NTK) theory to balance the magnitude of gradients across residual, boundary, and initial loss terms. This adaptive scheme is shown to accelerate convergence and improve predictive accuracy, particularly in fully self-supervised settings in expense of longer training times. In our case, however, we impose Neumann boundary conditions on the concentration field, and we observed that the directional gradients required by NTK-based weighting become ill-defined or unstable near the domain edges. As a result, we do not adopt the NTK reweighting strategy. Instead, we use a fixed weighting scheme for the composite loss function as defined in Equation (4), with $\lambda_{\mathrm{res}}=10.0$, $\lambda_{\mathrm{bcs}}=10^{-3}$, and $\lambda_{\mathrm{ics}}=1.0$.

Fig-3 illustrates end-to-end training pipeline for the proposed method.

Testing Procedure: To evaluate the model’s predictive accuracy, we compare the PI-DeepONet outputs against finite element (FEM) solutions computed on an adaptive mesh. Both predictions and reference solutions are mapped onto a uniform $k\times k$ Cartesian grid to ensure consistency. For each test source, the FEM pipeline generates a high-fidelity concentration field by solving the Darcy and transport equations sequentially, followed by interpolation onto the evaluation grid. In parallel, the PI-DeepONet model uses the same source and FEM-computed velocity field to generate its prediction directly on the same grid. We illustrate the full testing pipeline in Figure-4.

3.3 Model Architecture and Baseline Parameters

In this work, we adopt the modified DeepONet architecture introduced in [24], which we discussed in the introduction. Unless otherwise stated, both the branch and trunk networks consist of 4 hidden layers with 128 neurons per layer. All weights are initialized using the Glorot normal scheme, and hyperbolic tangent (tanh) activations are used. Training is conducted using mini-batch stochastic gradient descent via the Adam optimizer, with exponential learning rate decay. We set the initial learning rate to $10^{-3}$, decay steps to 5000, and decay rate to 0.95. All models are trained for $3\times 10^{5}$ iterations with a batch size of 200. While preliminary tests showed that fewer epochs (e.g., 300k–600k) were sufficient for simpler cases, we adopted a uniform upper-bound setting to ensure convergence across all configurations and to simplify comparison.

Table 1 summarizes the physical parameters used in the governing PDE system, while Table 2 lists the key modeling and architectural settings used in the PI-DeepONet experiments.

We should note that to keep the training pipeline lightweight and modular, we use set of fixed sampling-related hyperparameters across all experiments. The values provided in Table-2 were found to provide a good trade-off between accuracy and computational cost in both single and multi-source settings, which will be detailed in the subsequent section.

Error Metrics: All evaluation metrics are computed on a uniform $80\times 80$ Cartesian spatial grid over the domain $\Omega=[0,10]\times[0,10]$, and at $N_{t}=10$ time steps in the interval $(0,T]$. The test set consists of $N_{\text{test}}=30$ independently sampled source functions.

Let $c(\mathbf{x},t)$ denote the ground truth concentration field obtained from the FEM solver, and $c_{\theta}(\mathbf{x},t)$ denote the prediction from the trained PI-DeepONet model with parameters $\theta$, as defined earlier. We define two relative error metrics computed over all test instances of interest:

[TABLE]

where $\|\cdot\|_{2}$ denotes the standard Euclidean norm over the relevant grid. The metric $\mathcal{E}_{\text{full}}$ captures the model’s overall fidelity in approximating the full spatiotemporal evolution of the solution. Since our primary interest lies in accurately capturing the final concentration profile at $t=T$, we also report $\mathcal{E}_{T}$, which quantifies the model’s accuracy at the final time step. A detailed discussion of this design choice appears in Section-5.

4 Results

In this section, we consider several configurations for evaluating the operator learning capabilities of the proposed PI-DeepONet pipeline. Our primary focus is the mapping from a single Gaussian source to the concentration field at the final time, with the center and width of the Gaussian are randomly sampled. Building on this foundation, we extend the setup to more complex source terms defined as mixtures of multiple Gaussians. These multi-source functions are expressed as normalized sums of individual Gaussians, each with its own center and width:

[TABLE]

where $\mathbf{x}_{i}\sim\mathcal{U}(\Omega_{\text{src}})$ and $\sigma_{i}\sim\mathcal{U}(\sigma_{\min},\sigma_{\max})$ denote the center and spatial width of the $i$-th Gaussian, drawn independently from uniform distributions. Here, $\Omega_{\text{src}}\subset\Omega$ specifies a source-specific region within the domain, and $\mathcal{N}$ is a normalization constant. However, jointly varying both the center and the width across multiple sources significantly increases the complexity of the learning task. Specifically, the number of possible source configurations grows exponentially with the number of Gaussians, and the resulting input functions become increasingly difficult to represent due to overlapping supports, varying spatial scales, and nontrivial interaction patterns.

To explore what can be achieved under the same sampling and training constraints used in the single-source case, we focus on three practical subcases where a constraint is introduced into the source configuration. In the first, each training instance contains either one or two Gaussians, randomly selected, allowing the model to experience a range of source complexities. In the second, we consider the case of exactly three fixed-width Gaussians per sample. In the third subcase, we fix the positions of the three Gaussians and instead allow their widths to vary. From a practical standpoint, fixing the Gaussian width corresponds to using a consistent injection mechanism or fluid dispersal device with known spatial characteristics while fixing the centers reflects repeatable positioning of hardware or predefined injection zones.

4.1 Single Source with Varying Center and Width

We begin by evaluating the model’s performance in the baseline setting, where the source term consists of a single Gaussian profile with randomized spatial center and width. The source function is given by:

[TABLE]

where the center $\mathbf{x}_{0}\sim\mathcal{U}(\Omega_{\text{src}})$ is drawn uniformly from a square subregion $\Omega_{\text{src}}=[3,7]\times[3,7]\subset\Omega$, and the spatial width $\sigma\sim\mathcal{U}(0.25,0.60)$. This means its effective support (within $3\sigma$) covers on average only about 3.7% of the total domain area.

PI-DeepONet is able to recover the final-time concentration field with high accuracy. As shown in Table-3, the mean relative error across all time steps remains around $10\%$, while the final-time error drops is around 4% for moderate and long time horizons. Figure-5 illustrates representative comparisons between the predicted and FEM-simulated concentration fields for several test cases, showing close agreement in both shape and peak location.

4.2 Random Sampling of One or Two Sources with Fixed Width

We next consider a controlled extension of the baseline setting by allowing each input function to consist of either one or two Gaussian sources, randomly selected for each training and test instance. In this setup, the spatial width of each Gaussian is held fixed at a representative value $\sigma=0.45$, while the centers are drawn independently from the same region as before: $\mathbf{x}_{i}\sim\mathcal{U}(\Omega_{\text{src}})$, with $\Omega_{\text{src}}=[3,7]\times[3,7]$. The source function is defined as:

[TABLE]

where $N_{\text{source}}\in\{1,2\}$ is sampled uniformly, and $\mathcal{N}$ denotes a normalization factor ensuring consistent amplitude across inputs.

As summarized in Table-4, the model performs remarkably well in this setting, achieving final-time prediction accuracy comparable to the single-source case. The PI-DeepONet learns to generalize across a diverse set of overlapping and spatially separated source profiles without requiring a larger model or increased collocation budget. This formulation is also preferable to training a separate model with exactly two fixed sources, as it introduces less structural rigidity and better reflects real-world variability. In practical terms, it corresponds to scenarios where the number of injection sites may vary between deployments or trials, e.g., due to experimental error, adaptive control, or operational decisions

Figure-6 provides representative examples across final times $T=50,250,500$ showing model predictions for both one-source and two-source cases sampled at random. For example, in the second row (panel b2), the two Gaussians are located very close to one another and exhibit substantial overlap. Despite this configuration, the model still accurately reconstructs the final concentration field, capturing both the peak structure and spatial extent of the true solution.

4.3 Triple Source with Fixed Width

In this configuration, we consider a more structured but challenging multi-source scenario in which each input sample contains exactly three Gaussian sources with fixed spatial width. As before, the centers are drawn independently from the source subdomain $\Omega_{\text{src}}=[3,7]\times[3,7]$, and the spatial width is set to $\sigma=0.45$ for all components. The input function is defined as:

[TABLE]

where $\mathcal{N}$ is a normalization constant. This case probes the model’s ability to learn complex source–solution mappings under compounded spatial interactions. Compared to the mixed one/two-source case, the triple-source setting results in denser and often overlapping source profiles, with higher peak magnitudes and more intricate downstream transport behavior.

To accommodate the increased complexity of the input function class in this triple-source setting, we made modest adjustments to the model architecture and training schedule. Specifically, we use a deeper network with 4 hidden layers and 256 neurons per layer, and set the learning rate decay schedule with a step size of 10,000 and decay rate of 0.97. The model is trained for $7\times 10^{5}$ epochs using the same number of collocation and residual points as in the single-source case.

The results are summarized in Table-5. While the relative errors are slightly higher than those reported in earlier cases, as expected due to the added spatial complexity, the model still achieves strong predictive performance across all final times. Notably, the error at the final time step remains under 9% for all three test horizons.

4.4 Triple Source with Fixed Centers and Varying Width

As a complementary experiment to the previous section, we now consider the inverse configuration: the centers of the three Gaussian sources are held fixed, while their widths are independently varied. Each sample contains three Gaussians centered at fixed locations and the spatial width of each Gaussian is drawn from the uniform distribution $\sigma_{i}\sim\mathcal{U}(0.25,0.60)$. The resulting source term is given by:

[TABLE]

where $\mathcal{N}$ is a normalization constant, and the centers $\{\mathbf{x}_{i}\}$ are fixed to coordinates $(3,3)$, $(4,6)$, and $(7,4)$ across all training and test instances. The network architecture and training schedule are kept identical to those used in Section 3.2.2.

As summarized in Table-6, the model exhibits lower error compared to the fully variable triple-source case. Across all three final times, the relative error at the final time step drops to the 2–4% range, suggesting that the task is more learnable when the spatial layout of the sources remains fixed.

Figure-8 presents representative model outputs for final times $T=50$, $250$, and $500$. The concentration fields exhibit close agreement with the FEM ground truth, including accurate peak heights and spatial spread. In the last row, some minor visual artifacts are present in the PI-DeepONet prediction, likely due to subtle convergence irregularities, but the overall shape and magnitude remain well-preserved.

4.5 Performance in a Higher-Péclet Regime with Increased Permeability

Throughout our previous experiments, the intrinsic permeability was fixed at $K=10^{-9}$, corresponding to a maximum Péclet number of approximately $0.48$, placing the system in a diffusion-dominated regime. To assess the performance of the PI-DeepONet framework under more convection-dominated conditions, we increased the permeability by one order of magnitude to $K=10^{-8}$, which raises the maximum Péclet number to about $4.8$.

We repeated the experiments described in Sections 4.1 and 4.4 under this new setting. The results are summarized in Table-7.

In both experiments, the model maintains accuracy comparable to that observed in the low-convective regime. Figure-9 shows representative test cases for each setting. Due to the increased convective transport, peak concentration values are reduced: from approximately $12\text{mmol}$ to $7\text{mmol}$ in the single-source case, and from about $3.1$ to $2.5$ in the triple-source case, as the injected solute is advected away from the source region more rapidly.

4.6 Effect of Residual Point Sampling Strategy

To better understand the role of our spatial residual point sampling strategy, we conducted an ablation study comparing the standard hybrid sampling scheme to a purely random baseline. In our standard setup, residual collocation points are composed of a structured polar grid, using $n_{r}=30$ and $n_{\theta}=30$, centered at the Gaussian source, combined with $n_{\mathrm{rand}}=300$ additional points randomly distributed over the domain. This totals to 1200 and 3000 residual points per instance for single source and triple source cases respectively.

To isolate the effect of the structured sampling, we reran each experiment with the same total number of residual points but removed the polar grid entirely, setting $n_{r}=n_{\theta}=0$ and allocating all points to random sampling. This comparison was repeated for three representative experiments in Sec-4.1,4.3 and 4.4. All other model and training parameters were kept constant, and results were reported for the final-time case $T=500$.

The results are summarized in Table-8. Across all cases, replacing the dense grid with random points leads to significantly worse performance. In particular, the final-time error $\mathcal{E}_{T}$ increases by a significant margin across all experiments. This confirms that dense sampling around the source locations is critical for guiding the model to resolve steep gradients and local transport effects.

4.7 Runtime Comparison: FEM vs. PI-DeepONet

To evaluate the computational benefits of the learned PI-DeepONet model, we compare its inference time against the traditional FEM solver across four representative test cases at final time $T=500$. For a fair comparison, all file I/O operations on the FEM side were disabled, and timings were averaged over 30 test samples per experiment. As shown in Table 9, the PI-DeepONet model consistently achieves two to three orders of magnitude speed-up over the FEM baseline. The speed-up becomes especially pronounced in the more complex triple-source cases, reaching nearly $\times 300$ in multi-source setting.

5 Discussion and Future Directions

Our results demonstrate that the proposed FEM–PI-DeepONet pipeline can accurately learn the map $f\mapsto c$ even in the presence of sharp, spatially localized Gaussian sources and a coupled Darcy–Transport system. Unlike standard DeepONet demonstrations, which often focus on isolated PDEs with smooth input functions drawn from Gaussian random fields or low-frequency bases, our setup involves a fully physics-informed training objective, where the target equation is a time-dependent convection–diffusion system governed by a precomputed Darcy velocity field. While the Gaussian source functions used here are smooth, they are highly localized and sharply peaked, leading to steep gradients and requiring fine spatial resolution. Despite these challenges, we show that accurate learning is possible with as few as 2000 training functions and moderate computational resources, a consumer laptop with an Intel i7-11800H CPU and NVIDIA RTX 3080 GPU was used for all experiments.

The resulting model achieves up to $\times 300$ speed-up compared to full FEM solvers while maintaining strong agreement in final-time solute concentration profiles. This acceleration is especially valuable for settings that require repeated queries, such as uncertainty quantification, inverse design, or real-time control, where traditional FEM approaches become prohibitively expensive. A key advantage of our approach is its modularity: the FEM solver is used solely to compute the velocity field, while the operator network independently learns the mapping from the source function to concentration. The meshing, Darcy solver, and PI-DeepONet architecture are fully decoupled, allowing one to easily swap out components. For example, one can retrain the operator network using a different residual grid or with updated flow inputs, without modifying the FEM solver infrastructure. Since the velocity field is precomputed and mapped to the collocation domain before training, the FEM and PI-DeepONet components remain modular and independently configurable.

Our model performs well within the convective regime considered in this work; however, pushing to more strongly advection-dominated flows becomes increasingly difficult due to the combined effect of high Péclet numbers and sharply localized Gaussian sources. Moreover, while the Péclet number can be increased further by raising $K$, the combination of higher advection and sharply localized Gaussians leads to steep gradients that are difficult even for FEM to resolve, making such scenarios especially demanding for PI-DeepONet. Another important consideration is that, although trunk network collocation points are adaptively sampled, the branch network sensor grid is fixed by design. Recall that, in our experiments, the Gaussian width is drawn from $\sigma\sim U(0.25,0.60)$, which corresponds to an average coverage of only about $3.7\%$ of the domain area. Using a fixed $30\times 30$ branch grid is adequate in this range, but pushing $\sigma$ to smaller values would require a finer grid to capture the source accurately. This, in turn, would substantially increase the dimensionality of the branch network input layer, and hence the total number of trainable parameters, without necessarily adding proportional expressive power. This sensitivity to sensor resolution is a general limitation of DeepONet-based models and several recent works attempt to address this limitation [51, 52, 53].

From a physical modeling standpoint, an important limitation arises when the source term in the transport equation becomes extremely weak. Recall that the key parameter governing this term is the source coefficient $\beta_{2}$, which represents the rate of solute mass injection per unit volume. Physically, $\beta_{2}$ is the product of the volumetric infusion rate, which we already defined as $\beta_{1}$, and the solute concentration of the infused fluid $c_{\inf}$, so that $\beta_{2}=\beta_{1}c_{\inf}$. When $c_{\inf}$ is small, as is typically the case in realistic applications, the effective source term $\beta_{2}f$ becomes negligible. As a result, the solute concentration field is primarily governed by advection, and the mapping $f\mapsto c$ becomes nearly unidentifiable. To retain a learnable signal in this study, we used $c_{\inf}=1/4$, but this is significantly higher than what is physically realistic. In future work, we plan to address this limitation by learning the full mapping $(\mathbf{v},f)\mapsto c$, which more accurately reflects the underlying physics and may enable learning even when the source signal is weak. This extension introduces additional challenges related to velocity field representation and input dimensionality, but it is a necessary step toward modeling more realistic systems.

While the model performs very well at final times, its accuracy is more limited during earlier stages of the simulation. For example, in the case with final time $T=500$, the final-time relative $L^{2}$ error is approximately 5%, while the average error across all time steps is around 10%. This higher error is mainly due to the model not converging well at early times, which is expected for two reasons. First, the framework lacks an explicit notion of temporal causality which train over successive time steps to enforce sequential consistency, with significantly higher computational cost [54, 55]. Second, the model is subject to spectral bias: the localized Gaussian sources typically span only a small portion of the domain, which means that early-time signals are weak and spatially confined, making them harder to resolve. Despite these challenges, the model successfully captures the full concentration profile and achieves high accuracy at the final time, which is the primary target of interest in this study.

The training process also reveals that progressive learning, from single to multi-source configurations, improves stability and generalization. Although not explored in detail in this paper, initial experiments indicate that curriculum-based training strategies may be highly beneficial for operator learning in multiscale PDE systems. This is consistent with broader trends in curriculum learning literature, though applying such ideas to DeepONet remains challenging due to catastrophic forgetting and a lack of robust scheduling strategies [56, 57, 58]. .

In our pipeline, model training is performed in JAX [59], while data loading and batching are handled using PyTorch utilities [60] to enable efficient preprocessing and integration with standard ML workflows. All FEM-related computations are carried out using Firedrake [45]. We should note that a substantial software engineering effort is required to integrate these tools into a cohesive and performant workflow. This combination is powerful but technically demanding, often prone to subtle bugs and compatibility issues. Our implementation abstracts these complexities away and will be released in a Dockerized environment, allowing users to reproduce our experiments with no installation burden. We believe this kind of reproducibility and accessibility is essential for pushing the field forward. This effort aligns with a growing body of work aiming to bridge the gap between scientific computing and machine learning through differentiable programming and unified toolchains [61, 62, 63]. Projects like JAX-Fluids [61], JAX-MD [62], and JAX-FEM [63] are some examples of this trend trying to incorporate differentiable simulations with modern ML libraries. Our contribution follows this trajectory, demonstrating that tightly coupled PDE solvers and neural networks can be utilized in a unified, reproducible manner in ML-driven scientific scientific computing.

Acknowledgments

This project was completed with support from the U.S. Department of Energy, Advanced Scientific Computing Research program, under the ”Uncertainty Quantification for Multifidelity Operator Learning (MOLUcQ)” project (Project No. 81739). The research was also funded by the Simons Foundation through a grant made to Spelman College. Pacific Northwest National Laboratory (PNNL) is a multi-program national laboratory operated for the U.S. Department of Energy (DOE) by Battelle Memorial Institute under Contract No. DE-AC05-76RL01830.