{'fraction_non_alphanumeric': 0.041062986121702, 'fraction_numerical': 0.013458898886685984, 'mean_word_length': 4.66012084592145, 'pattern_counts': {'":': 0, '<': 1, '<?xml version=': 0, '>': 1, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 3, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Physics-Informed Neural Networks for Multiphysics Data Assimilation with Application to Subsurface Transport

Data assimilation for parameter and state estimation in subsurface transport problems remains a significant challenge due to the sparsity of measurements, the heterogeneity of porous media, and the high computational cost of forward numerical models. We present a physics-informed deep neural networks (DNNs) machine learning method for estimating space-dependent hydraulic conductivity, hydraulic head, and concentration fields from sparse measurements. In this approach, we employ individual DNNs to approximate the unknown parameters (e.g., hydraulic conductivity) and states (e.g., hydraulic head and concentration) of a physical system, and jointly train these DNNs by minimizing the loss function that consists of the governing equations residuals in addition to the error with respect to measurement data. We apply this approach to assimilate conductivity, hydraulic head, and concentration measurements for joint inversion of the conductivity, hydraulic head, and concentration fields in a steady-state advection--dispersion problem. We study the accuracy of the physics-informed DNN approach with respect to data size, number of variables (conductivity and head versus conductivity, head, and concentration), DNNs size, and DNN initialization during training. We demonstrate that the physics-informed DNNs are significantly more accurate than standard data-driven DNNs when the training set consists of sparse data. We also show that the accuracy of parameter estimation increases as additional variables are inverted jointly.


Introduction
Modeling of transport in heterogeneous porous media is a part of many environmental and engineering applications, including hydrocarbon recovery [1], hydraulic fracking [2], exploitation of geothermal energy [3], geologic disposal of radioactive waste [4], and groundwater contamination assessment [5]. Numerical models of transport in porous media require deterministic or statistical knowledge of the subsurface properties (e.g., hydraulic conductivity) and initial and boundary conditions [6,7]. However, because of heterogeneity and data sparsity, the parameters of natural systems are often unknown. Despite significant research in inverse methods [8,9,5,10], parameter estimation at a resolution required for accurate modeling of transport processes remains a challenge.
Parameter estimation is complicated by the fact that the parameters of interest (e.g., hydraulic conductivity) are difficult to measure directly. Most inverse (parameter estimation) methods use indirect measurements (in addition to direct measurements) to estimate parameters. Data assimilation (or model-data integration) has been well recognized as an effective technique to reduce predictive uncertainties and improve model accuracy. Data assimilation is a process where model parameters and system states are updated using measurements and governing equations [5,11,12]. Data assimilation has been used in many fields, including atmospheric and oceanic sciences [13,14], hydrology [15,9], subsurface transport [5,10,16], and uncertainty quantification [17,18]. What makes subsurface flow and transport applications challenging is that the inverse problems for the subsurface flow and transport equations are highly nonlinear and the states and parameters are non-Gaussian [16]. The non-linearity presents a challenge for both the direct inverse methods and Bayesian parameter estimation methods.
Recent advances in machine learning (ML) methods and automatic differentiation (including software infrastructures such as TensorFlow [19]) has made them potentially powerful tools for parameter estimation and data assimilation.
For example, Schmidt and Lipson [20] applied symbolic regression to learn conservation laws, and Brunton et al. [21] used sparse regression to discover equations of nonlinear dynamics directly from data. Physics-informed deep neural networks (PINNs) were used to learn solutions and parameters in partial and ordinary differential equations [22,23,24,22]. Recently, PINNs were extended for inverse problems associated with partial differential equations with random coefficients (e.g., to estimate the space-dependent hydraulic conductivity using sparse measurements of conductivity and hydraulic head) [25]. Deep neural networks (DNNs) were also combined with the finite element method (FEM) for estimating parameters in partial differential equations (PDEs) given that the PDE states are fully known [26].
In this study, we extend the PINN-based parameter estimation method of [25] to assimilate measurements of hydraulic conductivity, hydraulic head, and solute concentration. In the PINN method, we use the Darcy and advection-dispersion equations jointly with data to train DNNs that represent space-dependent conductivity, head, and concentration fields. During training of the DNNs, the governing equations and the associated boundary conditions are enforced at the collocation points over the domain. We demonstrate that for sparse data, PINNs significantly improve the accuracy of parameter and state estimation as compared to standard DNNs trained with data only. This paper is organized as follows. In Section 2, we briefly review the feedforward fully-connected deep networks and the PINN approach. The extension of PINN to the assimilation of multiphysics data is presented in Section 3.
Performance of PINNs, including the dependence of the estimation errors on the number of measurements, is given in Section 4. The effect of the neural network size and the conductivity field complexity on the accuracy of the parameter estimation is discussed in Section 5. Conclusions are given in Section 6.

Physics-informed neural networks
In this section, we first give a brief overview of feed-forward fully-connected DNNs as a function approximation and regression tool [27]. Next, we describe the PINN framework that enforces physics in the regression and enables assimilation of different types of data.

Hidden layers
Output layer ( ) We employ DNNs to approximate and regress a function u(x) using sparse measurements of u. We use a fully-connected feed-forward network architecture known as multilayer perceptrons, where the basic computing units (neurons) are stacked in layers, as shown in Figure 1. The first layer is called the input layer, and the last layer is the output layer, while all the intermediate layers are known as hidden layers. The DNN approximation of u(x) is given as: where(·) denotes a DNN approximation, and ...
y n l +1 (y n l ) = σ(W n l y n l + b n l ). ( Here, n l denotes the number of hidden layers, σ is the activation function, x ∈ R d denotes the input (d is the number of spatial dimensions), y n l +1 is the output vector, and θ is the vector of the weights and biases parameters: In the "data-driven" approach, θ is estimated by minimizing the loss function where T u = {x 1 , x 2 , ..., x |Tu| } ⊂ Ω denotes a set of measurement locations, Ω is the domain of function u, and u * (x), x ∈ T u are the measured values of u at these locations.
Some of the commonly used activation functions include logistic sigmoid, hyperbolic tangent, ReLu, and leaky ReLu. As the objective of this study is to approximate differentiable functions (space-dependent parameters and state variables in PDEs), we adopt the hyperbolic tangent activation function σ(x) = tanh(x), which is infinitely differentiable.

Enforcing physics constraints
Given a sufficiently large number of hidden layers, DNNs have excellent representative properties but require a lot of data to train them. This creates a challenge in applying DNNs to subsurface problems where measurements are usually sparse. For the purpose of this work, we define sparse measurements as measurements that do not sufficiently cover the computational domain to accurately determine θ from the unconstrained minimization problem (4). In [25], we demonstrated that the Darcy law can be used as a constraint for solving the minimization problem (4) and learning hydraulic conductivity from sparse measurements of conductivity and hydraulic head using the PINN approach [22].
In the rest of this section, we extend the PINN parameter estimation method of [25] to a data assimilation problem where different types of measurements (e.g., conductivity, hydraulic head, and concentration measurements) are used to estimate parameters (i.e., conductivity) and states (i.e., head and concentration fields). Consider a system of PDEs forming the boundary value problem defined on the domain Ω ⊂ R d with the boundary ∂Ω: where u is the (unknown) solution vector, p is the (unknown) system parameter vector, L denotes the known (nonlinear) differential operator, and the operator B expresses arbitrary boundary conditions associated with the problem. The boundary conditions can be of the Dirichlet and Neumann types applied on ∂ D Ω and ∂ N Ω, respectively, such that ∂ D Ω ∪ ∂ N Ω = ∂Ω and ∂ D Ω ∩ ∂ N Ω = ∅.
Following [25], the DNNs are used to approximate unknown states and parameters, u(x) ≈û(x; θ) and p(x) ≈p(x; γ), x ∈ Ω. To determine these parameters, we minimize the loss function J(θ, γ) with physics-informed penalty terms: where Here, J d (θ, γ) is the loss due to a mismatch with the data (i.e., the measurements of u and p): J f (θ, γ) is the loss due to mismatch with the governing PDEs L(u(x); p(x)) = 0: and J b (θ, γ) is the loss due to mismatch with the boundary conditions B(u(x); p(x)) = 0: In (7), ω f , and ω b are weights that determine how strongly mismatch with the governing PDEs and boundary conditions are penalized relative to data mismatch. In this work, we assume that the measurements and physics model are exact and set ω f = ω b = 1. The sets T u = {x 1 , x 2 , ..., x |Tu| } ⊂ Ω and ., x |Tp| } ⊂ Ω denote the measurement locations of u and p, respectively, and u * (x), x ∈ T u and p * (x), x ∈ T p are the measured values of u and p at these locations. The sets and J b (θ, γ) are, respectively, minimized. The penalty terms J f (θ, γ) and J b (θ, γ) force the DNN approximations of u and p to satisfy the governing equation (5) at the residual points. Note that while it is preferable to enforce physics over the whole domain, the computational cost of estimating and minimizing the loss function (7) increases with the number of residual points. In this work, we demonstrate convergence of the solution of (6) with an increasing number of residual points, meaning that the DNNsû(x; θ) andp(x; γ) can be accurately trained using a finite number of residual points. Similar convergence results for solving partial differential equations with the PINN method were also observed in [22,28,25,23,29].
The loss J f (θ, γ) is evaluated by computing spatial derivatives ofû(x; θ) and p(x; γ) using automatic differentiation (AD) [30]. AD is also used to evaluate the normal derivative n · ∇ in the Neumann boundary condition in the loss J b (θ, γ) (see details in Section 3). AD is implemented in most ML libraries, including TensorFlow and Pytorch [31], where it is mainly used to compute derivatives with respect to the DNN weights (i.e., θ and γ). In the PINN method, AD allows the implementation of any PDE and boundary condition constraints without numerically discretizing and solving the PDEs.
Another benefit of enforcing PDE constraints via the penalty term J f (θ, γ) is that we can use the corresponding weight ω f to account for the "fidelity" of the PDE model. For example, we can assign a smaller weight to a lowfidelity PDE model. In general, the number of unknown parameters in θ and γ is much larger than the number of measurements, and training of the DNNs requires regularization. One can consider the losses J b (θ, γ) and J f (θ, γ) in the minimization problem (6) as a physics-informed regularization terms [25,32]. at the x 1 = 0 boundary. This data includes the measurements of conductivity , respectively, where N K , N h , and N C are the number of measurements of each variable. Our objective is to learn the conductivity, head, and concentration fields based on these measurements.
We further assume that the concentration, hydraulic head, and conductivity data can be accurately modeled by the steady-state Darcy flow: and advection-dispersion equation: where φ is the effective porosity of the medium, v is the average pore velocity, and D is the dispersion coefficient: Here, I is the identity tensor, D w is the diffusion coefficient, τ is the tortuosity of the medium, and α is the dispersivity tensor with the diagonal components α L and α T . The conductivity K(x) is assumed to be unknown except at the In the following simulations, we set the parameters as: L 1 = 1 m, L 2 = 0.5 m, We start by defining the DNN representations of K(x), h(x), and C(x) as: where θ K , θ h , and θ C are the vectors of parameters associated with each neural network. Next, the residuals of equations (11) and (12) are expressed in terms of θ K , θ h , and θ C as: To enforce the Neumann boundary conditions for equations (11) and (12), we define DNNs approximating fluxes at the boundaries: and The loss function is then defined as: where the loss due to mismatch with data is and the losses due to PDE constraints and boundary conditions are: In Eq. (20), PDEs (11) and (12) are enforced at the residual points given by the In this work, we compare three approaches for trainingK: jointly training the DNNsK,ĥ, andĈ by minimizing the loss function (18) that we refer to as "multiphysics-informed neural networks" (MPINN); jointly trainingK,ĥ by only enforcing the Darcy equation and boundary conditions (11) that we refer to as "PINN-Darcy;" and separately trainingK,ĥ, andĈ using only data that we refer to as "data-driven DNN." The performance of these three approaches is investigated and compared in sections 4 and 5.
Considering that the loss function is highly nonlinear and non-convex with respect to the network parameters θ K , θ h , and θ C , we use the gradient descent minimization algorithms, including the Adam [33], and L-BFGS-B [34] algorithms.
As suggested in [22,25,35,36], the L-BFGS-B, a quasi-Newtown method, shows superior performance with a better rate of convergence, lower gradient vanishing, and a lower computational cost for problems with a relatively small amount of training data. Therefore, in the following numerical examples, we primarily adopt the L-BFGS-B method with the default settings from Scipy [37] and the Xavier initialization scheme [38].
Adding physics constraints complicates the loss function landscape and makes the minimization process more challenging. We demonstrate that the parameter and state estimation could be improved by first pre-trainingK andĥ with the PINN-Darcy approach and then retrainingK andĥ withĈ using the MPINN approach. We denote this approach as the "sequential training." The approach where all three DNNs are initialized and trained jointly we will refer to as the "simultaneous training."

Numerical example 1: periodic conductivity field
In this study, we test the proposed physics-informed neural network methods on synthetic data obtained by sampling numerical solutions of equations (11)- (12).
These equations are solved using the finite-volume Subsurface Transport Over Multiple Phase (STOMP) code [7] with 256 × 128 cells. We obtain numerical solutions for several conductivity fields and refer to these conductivity fields and the corresponding numerical solutions as reference fields and use them to test the accuracy of the proposed methods.
We first consider the reference conductivity field K(x) = 0.5 sin(4πx 1 ) sin(4πx 2 )+ 1, shown in Figure 3  We quantify the error between the estimated and reference fields in terms of the relative L 2 errors, defined as: where K(x), h(x), and C(x) are the reference fields. We first test the accuracy of the data-driven DNNs for learning K, h, and C fields from their corresponding measurements. The DNN parameters θ K , θ h , and

MPINN
In this section, we jointly trainK(x; θ K ),ĥ(x; θ h ), andĈ(x; θ C ) for estimating K, h, and C using the MPINN method with the proposed sequential training scheme. Figure 6 shows the mean error in the MPINN-estimated fields  C(x; θ C ). The shallow DNNs perform the worst because of their limited representation ability. The medium-size DNNs perform the best with a smaller number of K measurements, and the deep DNNs have the best performance with a larger number of measurements. This is because deeper DNNs are more representative but need more data for training. We note that the PDE constraints allow deeper DNNs to be trained with sparse data; see also our results in Section 4.2.1.
It is also worth noting that adding physics constraints makes it more difficult and expensive to minimize the loss function. The results in Figure 6 are obtained using sequential training of DNNs in MPINN. The sequential training is designed to simplify minimization of the loss function with multiple PDE constraints.
However, the sequential training also has the potential to introduce errors due to "decoupling" of the coupled physical processes. In Figure 7, we show the mean relative error in the estimated K, h, and C when the three DNNs are trained    chosen DNN architecture is given in Table 1. The conductivity fields in Figure 8 are generated on the domain Ω = [0, 1] × [0, 0.5] on a 256 × 128 grid with 32,768 grid points. Here, we use values of K at 20,000 grid points to trainK(x, θ K ) (without any physics constraints) and use K values at all 32,768 grid points to evaluate theK(x, θ K ) accuracy.
For this large number of measurements, we find that the L-BFGS-B algorithm is not efficient for minimizing the loss function, especially for the field with λ = 0.2. Therefore, here we adopt the Adam stochastic gradient descent (SGD) method [33] with the experimentally determined initial learning rate of 0.0002 and the batch size of 1000. We find that 4 × 10 5 iterations are needed to train K(x, θ K ) for λ = 0.2, 3 × 10 5 for λ = 0.5, and 2 × 10 5 for λ = 1.0 to achieve a sufficiently low training error. Our results show that less iterations are needed to train DNNs for smoother conductivity fields with larger correlation lengths. size increases because of the increasing DNN representation ability. We can see that a smaller DNN is sufficient to represent a smoother field with larger correlation length. We also see that for fields with the correlation length 0.5 and 1, the approximation error increases due to overfitting once the DNN size exceeds the "optimal size." For example, for the K field with λ = 0.5 (see Figure 9 (b)), the smallest relative error of 0.52% is reached at m h = 60. DNNs with m h < 60 are not representative enough, and DNNs with m h > 60 cause overfitting. Therefore, we postulate that the DNN with three hidden layers and m h = 60 is the optimal-size DNN for this K field. For the K fields with λ = 0.2 and λ = 1.0, the optimal DNN size is reached at m h = 90 and m h = 40, respectively. For the K field with λ = 1.0, the minimum of the mean error function (≈ 0.3%) is very shallow, as shown in Figure 9 (c). Hence, we select m h = 40 as the optimal DNN width because it results in the smallest error standard deviation.  Figure 10 shows an approximately power-law dependence of the optimal DNN size (in terms of the total number of DNN parameters) as a function of the correlation length of the modeled field. It is important to note that in addition to the correlation length of the modeled field, the optimal DNN size depends on many other factors, including the type of activation function and the number of hidden layers. In this study, we fix the number of hidden layers and the activation function. Therefore, the results in Figure 10 might not apply to other DNN architectures. forK, as this network size produces a reasonable (but not optimal) fit for these three correlation lengths. To optimize the loss function, we adopt the Adam algorithm with the learning rate of 0.0002, then switch to the L-BFGS method until the training loss reaches a low threshold of 0.0005.  For all three correlation lengths, the approximation error decreases with increasing N K . We also see that adding physics constraints improves the accuracy of the DNN approximation of the K field. The biggest reduction in the  Figure 11 (a)), the relative error of parameter estimation drops from more than 180% when using the data-driven DNN method to 66.8% in the PINN-Darcy estimate to 57.7% in the MPINN estimate.

Data assimilation
MPINN also provides a significantly improved hydraulic head and concentration estimations as compared to the data-driven DNN method, as evident in the comparison of the mean errors h and C in Table 2. We note that the number of measurements N h and N C is fixed, and the data-driven DNN hydraulic head and concentration estimates do not depend on N K . On the other hand, the PINN-Darcy and MPINN estimations of h and C improve with increasing N K .
This demonstrates the capability of the physics-informed DNNs to learn from indirect measurements.
The improvements due to using physics constraints and indirect measurements is particularly pronounced for estimating (highly nonlinear) C(x) can be seen in the comparison of the data-driven and MPINN C estimation errors. For example, for λ = 1 and N C = 100, C decreases from 16.72% in the data-driven DNN to 1.35%.
As we previously observed in Section 4.   Figure 12 with the reference K field in Figure 8 (b) shows that PINN-Darcy significantly improves the data-driven DNN prediction.
The MPINN further improves the K estimation as indicated by the smaller K .
The data-driven DNN approximation near the upper left corner significantly differs from the ground truth K field due to the lack of measurements in this region. However, the approximation error around this area is greatly reduced in  (c) K = 6.62% Although a good approximation of h can be obtained with all three methods (see Figure 13), we still observe that PINN and MPINN improve the data-driven DNN h estimation as indicated by the smaller mean h errors. For the highly- nonlinear C field, the data-driven DNN estimate is significantly less accurate than MPINN, both in terms of the concentration distribution and the mean C error, as seen in Figure 14. Notably, MPINN is able to accurately describe the eye of the concentration plume with very few direct measurements near this region. Once again, this demonstrates that MPINN can use sparse direct and indirect measurements in combination with PDEs to capture local features that otherwise cannot be described with only direct (sparse) measurements.

The effect of DNN size on the estimation error
Finally, we investigate whether using the optimal sizeK, as determined in Section 5.1, would reduce the error in the estimated K, h, and C fields. As an example, we choose the case with λ = 0.2. According to Figure 9(a), the optimalK size to represent a K field with this correlation length is m h = 90. In Section 5.2, we usedK with m h = 60 to estimate K, h, and C. Here, we repeat this study with m h = 90 and, for comparison, with m h = 120. The resulting estimation errors for the threeK DNN sizes are given in Table 3. In this comparison study, we fix theĥ andĈ DNNs' size at m h = 60 and use N K = 80, N h = 40, and N C = 100 measurements and N h f = 1000 and N C f = 1000 residual points. We can see that the optimal-sizeK produces the smallest estimation errors not only for the K field but also for the h field in all the data-driven DNN, PINN-Darcy, and MPINN methods. For the C field, the smallest error is achieved with m h = 60 in theK DNN. This shows that a smaller estimation error in K and h does not always translate into a smaller error in C.

Conclusion
In this study, we presented a physics-informed neural network approach for data and model assimilation for parameter and state estimation in multiphysics problems with application to subsurface transport problems. In this approach, all unknown parameter fields and states are modeled with DNNs, which are jointly trained by minimizing the loss function containing the multiphysics data (e.g., conductivity, hydraulic head, and concentration measurements ) and the associated physical models constraints. As a result, the DNNs can be trained using indirect measurements, which is important when the data is sparse.
By representing the space-dependent parameter and state variables with DNNs, the physics-informed approach offers a flexible and unified framework to deal with sparse and multiphysics data. In this study, we have compared three approaches: 1) the pure data-driven DNN approach, which only uses data to train DNNs; 2) the physics-informed DNN approach, called "PINN- We investigated the effect of the neural network size on the accuracy of parameter and state estimation as a function of the correlation length of the modeled K field. We demonstrated that small and large networks might result in poor representability or overfitting of DNNs and identified an optimal DNN size that has a power-law dependence on the correlation length. The physics constraints and added measurements reduce dependence of the DNN prediction on the DNN size given that the DNN is large (representative) enough. However, for a small number of measurements, we demonstrated that an optimal-size DNN outperforms the larger and smaller DNNs.
In subsurface applications, data is usually sparse and is often indirect. Therefore, the physics-informed DNNs provide an attractive alternative to standard data-driven machine learning methods. Since the proposed method involves training DNNs by minimizing the loss function, the performance of training algorithms is crucial. In our study, we found that introducing nonlinear PDE constraints into the loss function increases the computational cost of training.
Application of the physics-informed DNNs to large-scale problems will require access to multi-GPU computers and scalable training algorithms. The selection of training algorithms and hyperparameters (width and depth of DNNs) should also be studied in more details.