\section*{Introduction}

The volumetric displacement of a porous medium caused by the changes in fluid pressure inside the pore spaces is essential for many applications, including groundwater flow, underground heat mining, fossil fuel production, earthquake mechanics, and biomedical engineering 
\cite{nick2013reactive,bisdom2016impact,lee2016pressure,juanes2016were,vinje2018fluid}. Such volumetric deformation may impact the hydraulic storability and permeability of porous material, which influences the fluid flow behavior. This multiphysics problem, i.e., the coupling between fluid flow and solid deformation, can be captured through the application of Biot's equations of poroelasticity \cite{biot1941general, biot1957elastic}. \rev{The Biot's equations can be solved by analytical solutions for simple cases \cite{terzaghi1951theoretical,wang2017theory}. In complex cases, the finite difference approximation \cite{masson2010finite, wenzlau2009finite}, finite volume discretization \cite{nordbotten2014cell, sokolova2019multiscale}, or more commonly finite element methods such as the mixed formulation or discontinuous/enriched Galerkin, \cite{Haga2012,murad2013new, wheeler2014coupling,bouklas2015nonlinear,bouklas2015effect,choo2018enriched, liu2018lowest,liu2009coupled}, can be used. These numerical methods, however, require significant computational resources, and the accuracy of the solution depends heavily on the quality of the generated mesh. Hence, these methods may not be suitable to handle an inverse problem \cite{hansen2010discrete,hesthaven2016certified} or a forward problem with complex geometries \cite{wang2017physics,raissi2018deep,raissi2019physics}.}







\rev{Recent proposals have speculated that neural-network-based approaches such as deep learning might be an appealing alternative in solving physical problems, which are governed by partial differential equations since they operate without a mesh and are scalable on multithreaded systems \cite{kutz2017deep,wang2018model,zhang2019application,raissi2018deep, raissi2019physics}.} Moreover, deep learning has been successfully applied to many applications \cite{krizhevsky2012imagenet, alipanahi2015predicting, shen2017deep} because of its capability to handle highly nonlinear problems \cite{lecun2015deep}. \rev{As this technique, in general, requires a significantly large data set to reach a reasonable accuracy \cite{ahmed2015improved}, its applicability to many scientific and industrial problems can be a challenge \cite{raissi2017inferring}.} Use of a priori knowledge, however, can reduce the amount of examples needed. \rev{For instance, the idea of encoding physical information into the architectures and loss functions of the deep neural networks has been successfully applied to computational fluid dynamics problems \cite{xiao2016quantifying, wang2017physics, raissi2019physics,kutz2017deep,wang2018model,zhang2019application}.} The published results illustrate that by incorporating the physical information \rev{in the form of regularization terms of the loss functions}, the neural networks \rev{ can be trained to} provide good accuracy with a reasonably sized data set. 



Since the coupled fluid and solid mechanics process is highly nonlinear and generally involves complex geometries \cite{ruiz2019modelling, choo2018enriched, Kadeethum2019}, it seems to fit well into the context of physics-informed neural networks (PINN). For this reason, we propose to apply PINN for solving the nonlinear diffusivity and Biot's equations, both concerning forward and inverse modeling. The rest of the paper is organized as follows. The governing equations of the coupled solid and fluid mechanics are presented in the methodology section. Subsequently, the PINN architecture and its loss functions are defined. We then present forward and inverse modeling results for both the nonlinear diffusivity equation as well as Biot's equations. \rev{We also study the impact of the stochastic variations of the training procedures on the accuracy of the predicted primary variables and estimated parameters.} Finally, we conclude the findings and describe the possibilities that can enhance the capability of this model, which should be addressed in future works.

\section*{Methodology}

\subsection*{Governing equations}

We are interested in solving the nonlinear Biot's equations on the closed domain \rev{($\Omega$)} which amounts to a time-dependent multiphysics problem coupling solid deformation with fluid flow. Let $\Omega \subset \mathbb{R}^d$ be the domain of interest in $d$-dimensional space where $d = 1,$ $2,$ or $3$ and bounded by boundary, $\partial \Omega$. $\partial \Omega$ can be decomposed into displacement and traction boundaries, $\partial \Omega_u$ and $\partial \Omega_{\sigma}$, respectively, for the solid deformation problem. For the fluid flow problem, $\partial \Omega$ is decomposed into pressure and flux boundaries, $\partial \Omega_p$ and $\partial \Omega_q$, respectively. In short, $\partial \Omega_u$ and $\partial \Omega_p$ represent the first-kind boundary condition or Dirichlet boundary condition ($\partial \Omega_D$). The $\partial \Omega_{\sigma}$ and $\partial \Omega_q$, on the other hand, represent the second-kind boundary condition or Neumann boundary condition ($\partial \Omega_N$). The time domain is denoted by $\mathbb{T} = \left(0,\tau\right]$ with $\tau>0$.  \par

As just stated, the coupling between the fluid flow and solid deformation can be captured through the application of Biot's equations of poroelasticity, which is composed of linear momentum and mass balance equations \cite{biot1941general}. The linear momentum balance equation can be written as follows:  \par

\begin{equation}
\nabla \cdot \bm{\sigma}(\bm{u},p) =\bm{f},
\end{equation}

\noindent
where $\bm{u}$ is displacement, $p$ is fluid pressure, $\bm{f}$ is body force. The bold-face letters or symbols denote tensors and vectors, and the normal letters or symbols denote scalar values. Here, $\bm{\sigma}$ is \rev{the} total stress, which is defined as:

\begin{equation}
\bm{\sigma}:=\bm{\sigma^{\prime}}(\bm{u}) - \alpha p \bm{I},
\end{equation}

\noindent
where $\bm{I}$ is the identity tensor and $\alpha$ is Biot's coefficient defined as \cite{Jaeger2010}:

\begin{equation} \label{eq:biot_coeff}
\alpha:=1-\frac{K}{K_{{s}}},
\end{equation}

\noindent
with the bulk modulus of a rock matrix $K$ and the solid grains modulus $K_s$. In addition, $\bm{\sigma^{\prime}}$ is an effective stress defined as: 

\begin{equation}
\bm{\sigma^{\prime}}(\bm{u}):=2 \mu_{l} \bm{\varepsilon}(\bm{u})+\lambda_{l}  \boldsymbol{u} \bm{I},
\end{equation}

\noindent
where $\lambda_{l}$ and $\mu_{l}$ are Lam\'e constants, $\bm{\varepsilon}(\bm{u})$ is strain assuming infinitesimal \rev{displacements defined as:}

\begin{equation}
\bm{\varepsilon}(\bm{u}) :=\frac{1}{2}\left(\nabla \boldsymbol{u}+\nabla^{T} \boldsymbol{u}\right).
\end{equation}

\noindent
We can write the linear momentum balance and its boundary conditions as:

\begin{equation}\label{eq:linear_balance}
\nabla \cdot \bm{\sigma^{\prime}}(\bm{u}) -\alpha \nabla \cdot p \bm{I} = \bm{f} \text { in } \Omega \times \mathbb{T},
\end{equation}

\begin{equation}\label{eq:linear_balance_D}
\boldsymbol{u}=\boldsymbol{u}_{D} \text { on } \partial \Omega_{u} \times \mathbb{T},
\end{equation}

\begin{equation}\label{eq:linear_balance_Omega_t}
\boldsymbol{\sigma} \cdot \bm{n}=\bm{\sigma_{D}} \text { on } \partial \Omega_{\sigma} \times \mathbb{T},
\end{equation}

\begin{equation}
\bm{u}=\bm{u}_{0} \text { in } \Omega \text { at } t = 0,
\end{equation}

\noindent


\noindent
where $\bm{u_D}$ and ${\bm{\sigma_D}}$ are prescribed displacement and traction at boundaries, respectively, $\bm{n}$ is a normal unit vector, and $t$ is time.

The mass balance equation is written as \cite{coussy2004poromechanics, Kadeethum2019}:

\begin{equation} \label{eq:mass_balance}
\left(\phi c_{f}+\frac{\alpha-\phi}{K_{s}}\right) \frac{\partial p}{\partial t}+ \alpha \frac{\partial \nabla \cdot \boldsymbol{u}}{\partial t}-\nabla \cdot \mathcal{N}[\bm{\kappa}](\nabla p-\rho \mathbf{g})=g \text { in } \Omega \times \mathbb{T},
\end{equation}

\begin{equation} \label{eq:mass_balance_D}
p=p_{D} \text { on } \partial \Omega_{p} \times \mathbb{T},
\end{equation}

\begin{equation} \label{eq:mass_balance_N}
-\mathcal{N}[\bm{\kappa}](\nabla p-\rho \mathbf{g}) \cdot \boldsymbol{n}=q_{D} \text { on } \partial \Omega_{q} \times \mathbb{T},
\end{equation}

\begin{equation} \label{eq:mass_balance_I}
p=p_{0} \text { in } \Omega \text { at } t = 0,
\end{equation}

\noindent
where $\rho$ is fluid density, $\phi$ is initial porosity and remains constant throughout the simulation (the volumetric deformation is represented by $\partial \nabla \cdot \boldsymbol{u} / \partial t$), $c_f$ is fluid compressibility, $\mathbf{g}$ is a gravitational vector, $g$ is sink/source, $p_D$ and $q_D$ are specified pressure and flux, respectively, $\mathcal{N}[\cdot]$ represents a nonlinear
operator, and $\bm{\kappa}$ is hydraulic conductivity defined as:

\begin{equation}  \label{eq:permeability_matrix}
\bm{\kappa}:=\left[ \begin{array}{lll}{\kappa^{xx}} & {\kappa^{xy}} & {\kappa^{xz}} \\ {\kappa^{yx}} & {\kappa^{yy}} & {\kappa^{yz}} \\ {\kappa^{zx}} & {\kappa^{zy}} & \kappa^{zz}\end{array}\right],
\end{equation}


\noindent
where the tensor components characterize the transformation of the components of the gradient of fluid pressure into the components of the velocity vector. The $\kappa^{xx}$, $\kappa^{yy}$, and $\kappa^{zz}$ represent the matrix permeability in x-, y-, and z-direction, respectively. \rev{In this study, all off-diagonal terms are zero because we assume that a porous media is isotropic and homogeneous \cite{durlofsky1991numerical,holden1992tensor}.} Note that the diffusivity equation is a specific case of the Biot's equations since Eqs~(\ref{eq:linear_balance}) and (\ref{eq:mass_balance}) decouple when $\alpha = 0$. The details of this equation are presented in \rev{the results and discussion section.} 

\subsection*{Physics-informed neural networks model}

\subsubsection*{Neural network architecture}
The neural network architecture used in this study is presented in Fig~\ref{fig:gerneral_nn} \cite{rumelhart1986learning,lecun2015deep,hinton2006reducing}. The number of input and output nodes in the neural networks are determined from the problem formulation; for example, if the problem is time-dependent poroelasticity (as discussed in the previous section) and bounded by $\Omega = \left[ 0,1\right]^1$, we have two input nodes ($x$ and $t$) and two output nodes ($u$ and $p$) where $x$ is coordinate in x-direction, $t$ is time, $u$ is the displacement in x-direction, and $p$ is fluid pressure. The number of hidden layers ($N_{hl}$) and the number of neurons ($N_n$) act as so-called  hyperparameters \cite{goodfellow2016deep}. \rev{Each neuron (e.g., $\mathrm{H_{1, 1}}$ $...$ $\mathrm{H_{1, N_n}}$) is connected to the nodes of the previous layer with adjustable weights and also has an adjustable bias. We denote the set of weights and biases as ($\mathrm{W}$) and ($\mathrm{b}$), respectively. These variables are learned during a training phase \cite{hinton2006reducing, goodfellow2016deep}. In this paper, we define all hidden layers to have the same number of neurons.}

\rev{Physics-informed neural networks encode the information given by the differential operators as specific regularizing terms of the loss functions used when training the networks (see the section below). Because the training examples that are used to evaluate these extra regularizing terms, in general, are different from those used to train the network shown in Fig~\ref{fig:gerneral_nn}, one conceptually introduces an additional neural network - denoted the physical informed neural network \cite{raissi2019physics}. This additional neural network is dependent on all the $\mathrm{W}$ and $\mathrm{b}$ of the first neural network, but it introduces some extra variables to be learned for the inverse modeling case. This is elaborated in more detail in the below section on training the PINN. Note that while the neural network shown in Fig~\ref{fig:gerneral_nn} can be seen as point-to-point learning, the physics-informed regularization indirectly contributes to the local interactions (the stencil). The neural networks are built on the Tensorflow platform \cite{tensorflow2015-whitepaper}.} The results produced using either the rectified linear unit (ReLU) or the hyperbolic tangent ($\tanh$) were comparable. Hence, we only present the results using the $\tanh$ activation function in this paper.

\begin{figure}[!ht]
   \centering
        \includegraphics[width=10.0cm, height=7.0cm]{pictures/general_nn.pdf}
   \caption{
{\bf General neural network architecture used in this study \cite{rumelhart1986learning,lecun2015deep,hinton2006reducing}.} The input layer contains up to $i$ input nodes, and the output layer is composed of $1,...,k$ output nodes. $N_{hl}$ refes to the number of hidden layers, and each hidden layer is composed of $N_n$ neurons. Each neuron (e.g., $\mathrm{H_{1, 1}}$ $...$ $\mathrm{H_{1, N_n}}$) is connected to the nodes of the previous layer with adjustable weights   and also has an adjustable bias.}
   \label{fig:gerneral_nn}
\end{figure}

\subsubsection*{Physics-informed function}

We encode the underlying physical information to the neural networks through the \rev{so-called} physics-informed function ($\Pi$), acting as additional regularizing terms in the loss function defined below.

For the linear momentum balance equation Eq~(\ref{eq:linear_balance}), we define $\Pi_{\bm{u}}$ as follows:

\begin{equation}\label{eq:linear_balance_pinn}
\Pi_{\bm{u}} :=\nabla \cdot \bm{\sigma^{\prime}}(\bm{u}) -\alpha \nabla \cdot p \bm{I} - \bm{f} \text { in } \Omega \times \mathbb{T},
\end{equation}

\noindent
and with reference to Eq~(\ref{eq:linear_balance_Omega_t}) for its $\partial \Omega_{\sigma}$:

\begin{equation}\label{eq:linear_balance_N_pinn}
\Pi_{\bm{u_{\sigma}}} := \boldsymbol{\sigma} \cdot \bm{n} - \bm{\sigma_{D}} \text { on } \partial \Omega_{\sigma} \times \mathbb{T},
\end{equation}

\noindent
and for the mass balance equation Eq~(\ref{eq:mass_balance}), we define the $\Pi_p$ as:

\begin{equation} \label{eq:mass_balance_pinn}
\Pi_p := \left(\phi c_{f}+\frac{\alpha-\phi}{K_{s}}\right) \frac{\partial p}{\partial t}+ \alpha \frac{\partial \nabla \cdot \boldsymbol{u}}{\partial t}-\nabla \cdot \mathcal{N}[\bm{\kappa}](\nabla p-\rho \mathbf{g})-g \text { in } \Omega \times \mathbb{T},
\end{equation}

\noindent
and for its $\partial \Omega_{q}$ according to Eq~(\ref{eq:mass_balance_N})

\begin{equation} \label{eq:mass_balance_N_pinn}
\Pi_{p_q} := -\mathcal{N}[\bm{\kappa}](\nabla p-\rho \mathbf{g}) \cdot \boldsymbol{n} - q_{D} \text { on } \partial \Omega_{q} \times \mathbb{T}.
\end{equation}



\noindent
\rev{Demanding the $\Pi$ terms above to be as close to zero as possible corresponds to fulfilling Eqs. (\ref{eq:linear_balance}), (\ref{eq:linear_balance_Omega_t}), (\ref{eq:mass_balance}), and (\ref{eq:mass_balance_N}).} 

\subsubsection*{Loss function definition}

The loss function applied with the PINN scheme is composed of two parts (here we use a mean squared error - $M S E$ as the metric). The error on the training data ($M S E_{tr}$) and the mean square value of the regularization term given by the physics-informed function ($M S E_{\Pi}$):

\begin{equation}\label{eq:loss_gen}
M S E=M S E_{tr}+M S E_{\Pi},
\end{equation}

\noindent
where

\begin{equation}\label{eq:loss_gen_pi}
M S E_{\Pi} = M S E_{\Pi_{\Omega}} + M S E_{\Pi_{\partial \Omega_N}},
\end{equation}

\noindent

\begin{equation}\label{eq:loss_gen_pi_omega}
M S E_{\Pi_{\Omega}} = M S E_{\Pi_{\bm{u}}} + M S E_{\Pi_{p}},
\end{equation}

\noindent
and
\begin{equation}\label{eq:loss_gen_pi_partial_omega}
M S E_{\Pi_{\partial \Omega_N}} = M S E_{\Pi_{\bm{u_{\sigma}}}} + M S E_{\Pi_{p_q}}.
\end{equation}

\noindent
where $M S E_{\Pi_{\bm{u}}}$, $ M S E_{\Pi_{\bm{u_{\sigma}}}}$, $M S E_{\Pi_{p}}$, and $M S E_{\Pi_{p_q}}$ correspond to the loss function of Eqs~(\ref{eq:linear_balance_pinn}), (\ref{eq:linear_balance_N_pinn}), (\ref{eq:mass_balance_pinn}), and (\ref{eq:mass_balance_N_pinn}), respectively. \rev{The Dirichlet boundary conditions given on $\partial \Omega_D$ with respect to the linear momentum Eq~(\ref{eq:linear_balance_D}) and mass balance Eq~(\ref{eq:mass_balance_D}) equations are automatically incorporated into the $M S E_{tr}$.} 


\rev{
A graphical presentation of how boundary points, initial points, as well as domain data,  are used for training the neural networks is provided in Fig~\ref{fig:loss_function_explain}. For the forward model, the set of points that constitutes $\partial \Omega_{D} \times \mathbb{T}$ and $\Omega \text { at } t = 0$ contributes to the $M S E_{tr}$ term of the loss function. We denote this set of training data as $N_b$. Besides, we have a set of collocation points that are sampled from $\partial \Omega_{N} \times \mathbb{T}$ and $\Omega \times \mathbb{T}$. These training data are used to minimize the $\Pi$ parts of the loss function, and we denote these data points as $N_{\Pi}$.
}

\rev{
For the inverse model, we will have a set of known/measured data points that can be sampled from the whole domain, i.e., $\Omega \times \mathbb{T}$. We refer to this set of training data as $N_{tr}$. All this data will contribute to both the $M S E_{tr}$ and the $M S E_{\Pi_{\Omega}}$ terms of the loss function. Also, we may include an extra set of collocation points (i.e., data where we would have no measured values), which would contribute only to the $M S E_{\Pi_{\Omega}}$ term.
}



\begin{figure}[!ht]
   \centering
        \includegraphics[width=9.0cm, height=7.0cm]{pictures/domain_test.pdf}
        \includegraphics[width=7.0cm, height=7.0cm]{pictures/domain_train.pdf}
   \caption{
{\bf Illustration of the parts of input space used for training the PINN:} (a) forward model and (b) inverse model. The collocation ($\Omega \times \mathbb{T}$), boundary and initial points ($\partial \Omega_{D} \times \mathbb{T}$ and $\Omega \text { at } t = 0$) are utilized in the forward model. However, only the training points ($\Omega \times \mathbb{T}$) are employed in the inverse model. $\mathrm{X}$ represents a set of spatial coordinates ($x$, $y$, and $z$), $t$ is a set of coordinates in the time domain, and $\mathrm{S}$ is a set of solution values ($\bm{u}$ and $p$) corresponding to $\mathrm{X}$ and $t$.}
   \label{fig:loss_function_explain}
\end{figure}


\subsubsection*{Training the PINN}

\rev{
The structure of a PINN architecture is shown in Fig \ref{fig:des_for_inv}. We assume a simple case of two inputs, $\mathrm{I_1}$ and $\mathrm{I_2}$, one output, $\mathrm{O}$, and one hidden layer with two neurons. Moreover, we here assume that $\Pi$ is a function of $\mathrm{O}$, the first derivative of $\mathrm{O}$ with respect to the set of inputs, and the set of physical parameters, $\theta$. With $\left(\cdot\right)_{\mathrm{tr}}$ we represent a training set, where the $\mathrm{O}$ values are known for given $\mathrm{I_{1,tr}}$ and $\mathrm{I_{2,tr}}$. This set is used to evaluate the $M S E_{tr}$ term of Eq~(\ref{eq:loss_gen}) at each training cycle. The $M S E_{\Pi}$ part of Eq~(\ref{eq:loss_gen}) is obtained by calculating $\Pi$ for the collocation points $\left(\cdot\right)_{\Pi}$; see bottom part of Fig \ref{fig:des_for_inv}. Calculating $\Pi$ involves knowing $\mathrm{O}$ and its derivatives with respect to $\mathrm{I_{1}}$ and $\mathrm{I_{2}}$. As these values are unknown for the collocation points, they are estimated using the current approximation of $\mathrm{O}(\mathrm{W},\mathrm{b},\mathrm{I_{1}},\mathrm{I_{2}})$ and its derivatives, which can be obtained using automatic differentiation \cite{griewank1989automatic,corliss2002automatic}. In this way, the regularization term $M S E_{\Pi}$ indirectly depends on the weights, $\mathrm{W}$, and biases, $\mathrm{b}$, of the neural architecture shown at the top of Fig \ref{fig:des_for_inv}. As a result, $\Pi$ can be written as a function of $\mathrm{W}$, $\mathrm{b}$, $\mathrm{I_{1,\Pi}}$, $\mathrm{I_{2,\Pi}}$, and $\theta$. If we denote this function, $\gamma$,  we have $\gamma\left( \mathrm{W},\mathrm{b},\mathrm{I_{1,\Pi}},\mathrm{I_{2,\Pi}}, \theta \right)$ $=$ $\Pi\left(O(\mathrm{W},\mathrm{b},\mathrm{I_{1,\Pi}},\mathrm{I_{2,\Pi}}),\frac{\partial\mathrm{O}}{\partial\mathrm{I_{1}}}, \frac{\partial\mathrm{O}}{\partial\mathrm{I_{2}}}, \theta \right)$.  Note that the specific mapping of $\mathrm{I_{1,\Pi}}$ and $\mathrm{I_{2,\Pi}}$ to $\Pi$ has been defined as the physical informed neural network in previous work \cite{raissi2019physics} as opposed to the neural network shown in the top of Fig \ref{fig:des_for_inv}.
}

\rev{
For both the forward and the inverse modeling cases, one trains the neural networks to establish a mapping from the input space given by $\mathrm{I_1}$ and $\mathrm{I_2}$ to the output space, $\mathrm{O}$, by minimizing $M S E_{tr}$ and $M S E_{\Pi}$. The essential differences between the two cases come down to the type of training examples being available to train the networks, as discussed in Fig  \ref{fig:loss_function_explain}, and whether the physical parameters ($\theta$) are known or not. For the forward modeling, we apply $\left(\cdot\right)_{\mathrm{tr}}$ to evaluate the $M S E_{tr}$ term and $\left(\cdot\right)_{\Pi}$ to evaluate the $M S E_{\Pi}$ term. The $\mathrm{W}$ and $\mathrm{b}$ are then adjusted to minimize the sum of these terms. One should emphasize that for the forward modeling, all the variables to be learned belong to the neural network depicted in the top of Fig \ref{fig:des_for_inv}.
}



\rev{
For the inverse modeling, the aim is to estimate $\theta$. We still, however, train a neural network to predict $\mathrm{O}$ (similar to the forward case). That is, we are not using a loss function that involves measuring a distance between estimated values of $\theta$ and their ground truth values. Instead, the reasoning behind solving the inverse problem is that we expect the unknown $\theta$ to converge towards their true values during training because we also allow $\theta$ to be adjusted along with $\mathrm{W}$ and $\mathrm{b}$. During a training phase, these variables are learned to minimize the combined sum of $M S E_{tr}$ and $M S E_{\Pi}$. Specifically, the variables are adjusted by backpropagating the errors as calculated by Eq. (\ref{eq:loss_gen}) \cite{hinton2006reducing,goodfellow2016deep}. Unlike the forward problem, the boundary and initial conditions are unknown and cannot be used to generate training examples. Instead, we provide training examples that, in real cases, would be obtained from measurements, which ideally correspond to solution points of the forward problem inside $\Omega$ (see Fig \ref{fig:loss_function_explain}b). 
}











\rev{
As discussed above, the solution of the inverse problem is based on training the neural network to establish a mapping from $\mathrm{I_1}$ and $\mathrm{I_2}$ to $\mathrm{O}$ as we do in the forward case. This means that the hyperparameters estimated from the forward modeling may act as qualified estimates for the hyperparameters in the inverse case, assuming the number of training examples in both cases are similar. By definition of the inverse problem, we do not know the true values of $\theta$, so we would have to assume that the exact values of $\theta$ have a limited influence on the hyperparameters. A more direct way of estimating the hyperparameters in the inverse case would be to divide the data into training and validation sets and then select the hyperparameters that minimize $MSE$ with respect to learning the mapping from $\mathrm{I_1}$ and $\mathrm{I_2}$ to $\mathrm{O}$. The downside of this is that we could then not spend all the measurement data on training the neural networks.
}






\begin{figure}[!ht]
   \centering
        \includegraphics[width=13.0cm, height=7.0cm]{pictures/des_for_inv_rev2.pdf}
   \caption{
{\bf Graphical illustration of how a traditional neural network is linked to a physics-informed neural network.} The set of $\theta$ represents unknown physical parameters that we want to estimate.}
   \label{fig:des_for_inv}
\end{figure}

\section*{Results and Discussion}

We first consider the results of forward and inverse models of the nonlinear diffusivity equation. Subsequently, we present the results of the nonlinear Biot's equations for both the forward and inverse models. 

\subsection*{Nonlinear diffusivity equation}
We assume $\alpha = 0$ to decouple the momentum and mass balance equations and focus on Eq~(\ref{eq:mass_balance}), which is reduced to

\begin{equation} \label{eq:mass_balance_decouple}
\phi c_{t} \frac{\partial p}{\partial t}-\nabla \cdot \mathcal{N}[\bm{\kappa}](\nabla p-\rho \mathbf{g})=g \text { in } \Omega \times \mathbb{T},
\end{equation}

\noindent
where $c_t$ is a total compressibility. The same boundary and initial conditions, Eqs~(\ref{eq:mass_balance_D}) to (\ref{eq:mass_balance_I}), are still valid. 


\rev{
We take $\Omega = \left[ 0,1\right]^1$, $\mathbb{T} = \left[0,1\right]$, and choose the exact solution in $\Omega$ as:
}
\begin{equation}\label{eq:darcy_exact_solution}
p(x,t) := \sin(x+t),
\end{equation}

\noindent
and $\mathcal{N}[\bm{\kappa}]$ as:

\begin{equation}  \label{eq:perm_darcy}
\mathcal{N}[\bm{\kappa}] := \bm{\kappa_0}p^2,
\end{equation}

\noindent
where $x$ and $t$ represent points in x-direction and time domain, respectively. The $\bm{\kappa_0}$ is assumed to be a scalar in this case, i.e., $\bm{\kappa_0}={\kappa_0}$. All the physical constants are set to $1.0$; and subsequently, $g$ is chosen as: 

\begin{equation}  \label{eq:darcy_source}
g(x,t) := \cos(x+t) + \sin(x+t)^3 - 2\cos(x+t)^2\sin(x+t),
\end{equation}

\noindent
to satisfy the exact solution. Furthermore, the homogeneous boundary conditions are applied to all boundaries and initial conditions using Eq~(\ref{eq:darcy_exact_solution}). Combining Eqs~(\ref{eq:mass_balance_decouple}) to (\ref{eq:darcy_source}), the physics-informed function ($\Pi$) for the nonlinear diffusivity equation is defined:

\begin{equation}\label{eq:pif_darcy_forward}
\Pi(x,t) :=\phi c_{t} \frac{\partial p}{\partial t}-{\kappa}_0\frac{\partial }{\partial x} p^2(\frac{\partial }{\partial x} p)-g.
\end{equation}

\rev{
We generate the solution based on an interval mesh having 2559 equidistant spatial intervals and 99 temporal ones; hence, in total, we have 256000 solution points, including the points on boundaries. Subsequently, we randomly draw $n$ training examples. Half of the remaining points are randomly selected as a validation set, and the remaining ones are used for testing. As an illustration, if we have 256000 solution points and use 100 examples to train the model. We then use 127950 examples for the validation and 127950 examples for the test set. 
}





\rev{For the nonlinear diffusivity equation, the forward modeling with a neural network should calculate the $p$ at any given $x$ and $t$ from provided values of $\phi$, $c_t$, and ${\kappa}_0$. For the inverse case, the aim is to infer $\phi$, $c_t$, and ${\kappa}_0$ from observed values of $x$, $t$, and $p$.} The architecture of the neural network corresponding to the top of Fig \ref{fig:des_for_inv} is illustrated in Fig~\ref{fig:darcy_nn}. We have two input nodes and one output node for this case. \rev{We use L-BFGS \cite{liu1989limited}; a quasi-Newton, fullbatch gradient-based optimization algorithm to minimize the loss function with stop criteria as $\frac{ \left|{M S E}^k - {M S E}^{k+1} \right| }{\mathrm{max} \left( \left|{M S E}^k\right|,\left|{M S E}^{k+1}\right|,1.0 \right )} <= 10^{-16}$ where $\left(\cdot\right)^{k}$ and $\left(\cdot\right)^{k+1}$ are previous and current iteration, respectively.} \rev{The L-BFGS algorithm has several advantages that are suitable for this study; for example, it is more stable compared to stochastic gradient descent (SGD) and can handle large batch sizes well \cite{liu1989limited, malouf2002comparison}.} This setting is used for all the simulations presented in this paper unless stated otherwise.


\begin{figure}[!ht]
   \centering
        \includegraphics[width=9.5cm, height=6.0cm]{pictures/darcy_nn.pdf}
   \caption{
{\bf Neural network architecture used for nonlinear diffusivity equation.} This network corresponds to the top part of Fig \ref{fig:des_for_inv}. There are two input nodes, $x$ and $t$, and one output node, $p$. The number of hidden layers, $N_{hl}$, and the number of neurons for each hidden layer, $N_n$, denote the hyperparameters.}
   \label{fig:darcy_nn}
\end{figure}



\rev{
To reiterate, the sampling strategies of the forward and inverse modelings are illustrated in Fig \ref{fig:loss_function_explain}. To be specific, in the forward case, we must determine the solution of the partial differential equation based on known boundary values for the time interval of $\mathbb{T} = \left(0,\tau\right]$ combined with the initial values at $t = 0$. These boundary and initial points are used to calculate the $M S E_{tr}$ term of Eq. \ref{eq:loss_gen}. The inner points act as collocation points and are used to calculate the $M S E_{\Pi}$ term of Eq. \ref{eq:loss_gen}. Note that for the forward case, we can pick as many points as we wish to calculate for both terms of the loss function. 
For the inverse model, we know a set of $x$ $\mathrm{and}$ $t$ and their corresponding values of $p$ in advance. These examples are used to calculate the $M S E_{tr}$ and $M S E_{\Pi}$ terms of Eq. \ref{eq:loss_gen}. In a real setting, the amount of training examples is limited by the available measurements. 
}



\subsubsection*{Verification of the forward model}

Following Eq (\ref{eq:loss_gen}) we obtain:

\begin{equation}\label{eq:loss_darcy_forward}
M S E=M S E_{b}+M S E_{\Pi},
\end{equation}

\noindent
where

\begin{equation}\label{eq:loss_p_darcy_forward}
M S E_{b}=\frac{1}{N_{b}} \sum_{i=1}^{N_{b}}\left|p\left(x_{b}^{i}, t_{b}^{i}\right)-p^{i}\right|^{2}
\end{equation}

\noindent
and

\begin{equation}\label{eq:loss_pif_darcy_forward}
M S E_{\Pi}=\frac{1}{N_{\Pi}} \sum_{i=1}^{N_{\Pi}}\left|\Pi\left(x_{\Pi}^{i}, t_{\Pi}^{i}\right)\right|^{2},
\end{equation}

\noindent
where $\left\{x_{b}^{i}, t_{b}^{i}, p^{i}\right\}_{i=1}^{N_{b}}$ refer to the initial and boundary training data on $p(x,t)$ while $\left\{x_{\Pi}^{i}, t_{\Pi}^{i}\right\}_{i=1}^{N_{\Pi}}$ denote the collocation points for $\Pi$, which are sampled using the Latin Hypercube method provided within the pyDOE package \cite{baudin2015pydoe}. The $N_b$ is the combined number of initial and boundary training data, and $N_{\Pi}$ is the number of collocation points.

\rev{
In Table \ref{table:darcy_forward_nu_nc}, we explore the accuracy of the PINN as a function of the number of training examples, i.e., different numbers of $N_b$ and $N_{\Pi}$ for a given size of the network; $N_{hl}$ and $N_n$ are fixed to four and ten, respectively. Due to the stochastic behavior of the training procedure of the neural networks, we calculate the results as an average over many realizations as it otherwise might be difficult to observe a clear trend of the relative $\mathcal{L}^2$ error as a function of the number of training examples (see Panel A of Table \ref{table:darcy_forward_nu_nc}). From Panels B and C of Table \ref{table:darcy_forward_nu_nc}, we observe that the average accuracy of the PINN is improved as $N_b$ and $N_{\Pi}$ are increased. Moreover, the combined number of $N_b$ and $N_{\Pi}$ to obtain an average value of the $\mathcal{L}^2$ error below $1.0 \times 10^{-2}$ is in the order of 100s. To illustrate the impact of $\Pi$ terms, we consider the case where we simply train a neural network to interpolate a feed-forward solution based on a known set of solution points. Without $\Pi$ regularization we obtain an average relative $\mathcal{L}^2$ error of $5.6 \times 10^{-2}$ based on 96 training examples while the average $\mathcal{L}^2$ error becomes less than $1.0 \times 10^{-2}$ when including $\Pi$. 
}


\rev{
\begin{table}[!ht]
\centering
\caption{
{\bf Diffusivity equation - forward modeling: PINN performance as function of training set size.} Relative ${\mathcal{L}^2}$ errors between the exact and predicted values of ${p}$ for the \rev{validation} set. The table shows the dependency on the number of the initial and boundary training data, $N_b$, and on the number of collocation points, $N_{\Pi}$. Here, the network architecture is fixed to 4 layers with 10 neurons per hidden layer.}
\subcaption*{Panel A: Average over 1 realization}
\begin{tabular}{|l|l|l|l|l|l|l|l|l|l|}
\hline
\backslashbox{$N_b$}{$N_{\Pi}$} 
          & 2     & 5     & 10    & 20    & 40    & 80    & 160   & 320   & 640 \\
    \hline
    2     & \cellcolor[rgb]{ .988,  .851,  .859}0.1010 & \cellcolor[rgb]{ .984,  .831,  .839}0.1150 & \cellcolor[rgb]{ .973,  .412,  .42}0.4060 & \cellcolor[rgb]{ .98,  .627,  .635}0.2570 & \cellcolor[rgb]{ .984,  .776,  .784}0.1530 & \cellcolor[rgb]{ .988,  .98,  .988}0.0108 & \cellcolor[rgb]{ .984,  .718,  .729}0.1930 & \cellcolor[rgb]{ .988,  .961,  .973}0.0236 & \cellcolor[rgb]{ .988,  .965,  .976}0.0216 \\
    \hline
    3     & \cellcolor[rgb]{ .988,  .945,  .957}0.0339 & \cellcolor[rgb]{ .988,  .898,  .91}0.0665 & \cellcolor[rgb]{ .984,  .843,  .855}0.1050 & \cellcolor[rgb]{ .988,  .976,  .988}0.0122 & \cellcolor[rgb]{ .839,  .882,  .945}0.0020 & \cellcolor[rgb]{ .988,  .855,  .867}0.0973 & \cellcolor[rgb]{ .957,  .965,  .988}0.0025 & \cellcolor[rgb]{ .988,  .949,  .961}0.0314 & \cellcolor[rgb]{ .988,  .984,  .996}0.0074 \\
    \hline
    6     & \cellcolor[rgb]{ .988,  .957,  .969}0.0248 & \cellcolor[rgb]{ .988,  .937,  .949}0.0401 & \cellcolor[rgb]{ .984,  .835,  .847}0.1110 & \cellcolor[rgb]{ .988,  .878,  .89}0.0807 & \cellcolor[rgb]{ .627,  .733,  .871}0.0012 & \cellcolor[rgb]{ .988,  .984,  .996}0.0057 & \cellcolor[rgb]{ .69,  .776,  .894}0.0014 & \cellcolor[rgb]{ .988,  .976,  .988}0.0116 & \cellcolor[rgb]{ .988,  .988,  1}0.0046 \\
    \hline
    12    & \cellcolor[rgb]{ .988,  .984,  .996}0.0072 & \cellcolor[rgb]{ .984,  .82,  .827}0.1230 & \cellcolor[rgb]{ .988,  .976,  .988}0.0126 & \cellcolor[rgb]{ .988,  .953,  .965}0.0296 & \cellcolor[rgb]{ .643,  .745,  .878}0.0012 & \cellcolor[rgb]{ .988,  .976,  .988}0.0119 & \cellcolor[rgb]{ .592,  .71,  .859}0.0010 & \cellcolor[rgb]{ .71,  .792,  .902}0.0015 & \cellcolor[rgb]{ .988,  .984,  .996}0.0076 \\
    \hline
    24    & \cellcolor[rgb]{ .988,  .973,  .984}0.0152 & \cellcolor[rgb]{ .776,  .839,  .925}0.0018 & \cellcolor[rgb]{ .533,  .667,  .839}0.0008 & \cellcolor[rgb]{ .529,  .663,  .835}0.0008 & \cellcolor[rgb]{ .49,  .635,  .824}0.0006 & \cellcolor[rgb]{ .725,  .804,  .906}0.0016 & \cellcolor[rgb]{ .537,  .671,  .839}0.0008 & \cellcolor[rgb]{ .478,  .627,  .82}0.0006 & \cellcolor[rgb]{ .569,  .694,  .851}0.0009 \\
    \hline
    48    & \cellcolor[rgb]{ .988,  .976,  .988}0.0124 & \cellcolor[rgb]{ .988,  .976,  .988}0.0127 & \cellcolor[rgb]{ .706,  .788,  .898}0.0015 & \cellcolor[rgb]{ .988,  .988,  1}0.0037 & \cellcolor[rgb]{ .988,  .988,  1}0.0044 & \cellcolor[rgb]{ .655,  .753,  .882}0.0013 & \cellcolor[rgb]{ .988,  .988,  1}0.0035 & \cellcolor[rgb]{ .38,  .561,  .784}0.0002 & \cellcolor[rgb]{ .988,  .988,  1}0.0026 \\
    \hline
    96    & \cellcolor[rgb]{ .988,  .984,  .996}0.0076 & \cellcolor[rgb]{ .584,  .706,  .859}0.0010 & \cellcolor[rgb]{ .925,  .945,  .976}0.0023 & \cellcolor[rgb]{ .702,  .788,  .898}0.0015 & \cellcolor[rgb]{ .988,  .988,  1}0.0026 & \cellcolor[rgb]{ .765,  .831,  .922}0.0017 & \cellcolor[rgb]{ .569,  .694,  .851}0.0009 & \cellcolor[rgb]{ .518,  .655,  .831}0.0007 & \cellcolor[rgb]{ .584,  .706,  .859}0.0010 \\
    \hline
    192   & \cellcolor[rgb]{ .988,  .988,  1}0.0036 & \cellcolor[rgb]{ .494,  .639,  .824}0.0006 & \cellcolor[rgb]{ .988,  .965,  .976}0.0205 & \cellcolor[rgb]{ .976,  .98,  .996}0.0025 & \cellcolor[rgb]{ .576,  .698,  .855}0.0010 & \cellcolor[rgb]{ .988,  .988,  1}0.0051 & \cellcolor[rgb]{ .467,  .62,  .816}0.0005 & \cellcolor[rgb]{ .522,  .659,  .835}0.0008 & \cellcolor[rgb]{ .38,  .561,  .784}0.0002 \\
    \hline
    384   & \cellcolor[rgb]{ .827,  .875,  .941}0.0020 & \cellcolor[rgb]{ .804,  .859,  .933}0.0019 & \cellcolor[rgb]{ .557,  .682,  .847}0.0009 & \cellcolor[rgb]{ .42,  .588,  .8}0.0004 & \cellcolor[rgb]{ .537,  .671,  .839}0.0008 & \cellcolor[rgb]{ .353,  .541,  .776}0.0001 & \cellcolor[rgb]{ .451,  .612,  .812}0.0005 & \cellcolor[rgb]{ .451,  .608,  .808}0.0005 & \cellcolor[rgb]{ .404,  .576,  .792}0.0003 \\
    \hline
\end{tabular}
\bigskip
\subcaption*{Panel B: Average over 3 realizations}
\begin{tabular}{|l|l|l|l|l|l|l|l|l|l|}
\hline
\backslashbox{$N_b$}{$N_{\Pi}$} 
          & 2     & 5     & 10    & 20    & 40    & 80    & 160   & 320   & 640 \\
    \hline
    2     & \cellcolor[rgb]{ .976,  .514,  .522}0.0972 & \cellcolor[rgb]{ .984,  .71,  .718}0.0596 & \cellcolor[rgb]{ .984,  .804,  .816}0.0417 & \cellcolor[rgb]{ .98,  .573,  .58}0.0859 & \cellcolor[rgb]{ .984,  .824,  .835}0.0380 & \cellcolor[rgb]{ .988,  .976,  .988}0.0086 & \cellcolor[rgb]{ .765,  .831,  .922}0.0041 & \cellcolor[rgb]{ .988,  .969,  .98}0.0105 & \cellcolor[rgb]{ .859,  .898,  .953}0.0050 \\
    \hline
    3     & \cellcolor[rgb]{ .98,  .631,  .643}0.0744 & \cellcolor[rgb]{ .98,  .678,  .69}0.0652 & \cellcolor[rgb]{ .976,  .514,  .522}0.0969 & \cellcolor[rgb]{ .984,  .804,  .816}0.0415 & \cellcolor[rgb]{ .988,  .933,  .941}0.0173 & \cellcolor[rgb]{ .988,  .882,  .894}0.0268 & \cellcolor[rgb]{ .988,  .886,  .898}0.0261 & \cellcolor[rgb]{ .988,  .976,  .988}0.0084 & \cellcolor[rgb]{ .988,  .98,  .992}0.0077 \\
    \hline
    6     & \cellcolor[rgb]{ .976,  .486,  .494}0.1020 & \cellcolor[rgb]{ .973,  .412,  .42}0.1160 & \cellcolor[rgb]{ .98,  .58,  .588}0.0841 & \cellcolor[rgb]{ .988,  .953,  .965}0.0133 & \cellcolor[rgb]{ .988,  .984,  .996}0.0071 & \cellcolor[rgb]{ .98,  .627,  .635}0.0756 & \cellcolor[rgb]{ .988,  .91,  .922}0.0217 & \cellcolor[rgb]{ .988,  .98,  .992}0.0078 & \cellcolor[rgb]{ .988,  .988,  1}0.0064 \\
    \hline
    12    & \cellcolor[rgb]{ .98,  .576,  .584}0.0853 & \cellcolor[rgb]{ .988,  .969,  .98}0.0099 & \cellcolor[rgb]{ .98,  .627,  .639}0.0751 & \cellcolor[rgb]{ .988,  .976,  .988}0.0088 & \cellcolor[rgb]{ .988,  .988,  1}0.0061 & \cellcolor[rgb]{ .984,  .784,  .792}0.0456 & \cellcolor[rgb]{ .702,  .784,  .898}0.0036 & \cellcolor[rgb]{ .965,  .973,  .992}0.0059 & \cellcolor[rgb]{ .675,  .769,  .89}0.0033 \\
    \hline
    24    & \cellcolor[rgb]{ .984,  .788,  .8}0.0447 & \cellcolor[rgb]{ .988,  .878,  .89}0.0272 & \cellcolor[rgb]{ .988,  .969,  .98}0.0104 & \cellcolor[rgb]{ .988,  .91,  .922}0.0214 & \cellcolor[rgb]{ .988,  .871,  .882}0.0290 & \cellcolor[rgb]{ .471,  .624,  .816}0.0015 & \cellcolor[rgb]{ .765,  .831,  .922}0.0041 & \cellcolor[rgb]{ .49,  .635,  .824}0.0017 & \cellcolor[rgb]{ .973,  .976,  .992}0.0060 \\
    \hline
    48    & \cellcolor[rgb]{ .988,  .918,  .929}0.0201 & \cellcolor[rgb]{ .714,  .796,  .902}0.0037 & \cellcolor[rgb]{ .435,  .6,  .804}0.0012 & \cellcolor[rgb]{ .988,  .902,  .914}0.0230 & \cellcolor[rgb]{ .745,  .816,  .914}0.0040 & \cellcolor[rgb]{ .988,  .984,  .996}0.0072 & \cellcolor[rgb]{ .769,  .831,  .922}0.0042 & \cellcolor[rgb]{ .969,  .973,  .992}0.0060 & \cellcolor[rgb]{ .576,  .698,  .855}0.0025 \\
    \hline
    96    & \cellcolor[rgb]{ .988,  .969,  .98}0.0104 & \cellcolor[rgb]{ .529,  .663,  .835}0.0020 & \cellcolor[rgb]{ .443,  .604,  .808}0.0013 & \cellcolor[rgb]{ .573,  .694,  .851}0.0024 & \cellcolor[rgb]{ .839,  .882,  .945}0.0048 & \cellcolor[rgb]{ .424,  .592,  .8}0.0011 & \cellcolor[rgb]{ .4,  .576,  .792}0.0009 & \cellcolor[rgb]{ .514,  .651,  .831}0.0019 & \cellcolor[rgb]{ .522,  .659,  .835}0.0020 \\
    \hline
    192   & \cellcolor[rgb]{ .988,  .925,  .937}0.0186 & \cellcolor[rgb]{ .455,  .612,  .812}0.0014 & \cellcolor[rgb]{ .545,  .675,  .843}0.0022 & \cellcolor[rgb]{ .373,  .553,  .78}0.0006 & \cellcolor[rgb]{ .38,  .561,  .784}0.0007 & \cellcolor[rgb]{ .424,  .588,  .8}0.0011 & \cellcolor[rgb]{ .373,  .553,  .78}0.0006 & \cellcolor[rgb]{ .353,  .541,  .776}0.0004 & \cellcolor[rgb]{ .388,  .565,  .788}0.0008 \\
    \hline
    384   & \cellcolor[rgb]{ .475,  .624,  .816}0.0015 & \cellcolor[rgb]{ .424,  .588,  .8}0.0011 & \cellcolor[rgb]{ .416,  .584,  .796}0.0010 & \cellcolor[rgb]{ .408,  .58,  .796}0.0009 & \cellcolor[rgb]{ .365,  .549,  .78}0.0005 & \cellcolor[rgb]{ .384,  .561,  .784}0.0007 & \cellcolor[rgb]{ .42,  .588,  .8}0.0010 & \cellcolor[rgb]{ .365,  .549,  .78}0.0005 & \cellcolor[rgb]{ .357,  .545,  .776}0.0005 \\
    \hline
\end{tabular}
\bigskip
\subcaption*{Panel C: Average over 24 realizations}
\begin{tabular}{|l|l|l|l|l|l|l|l|l|l|}
\hline
\backslashbox{$N_b$}{$N_{\Pi}$} 
          & 2     & 5     & 10    & 20    & 40    & 80    & 160   & 320   & 640 \\
    \hline
    2     & \cellcolor[rgb]{ .973,  .412,  .42}0.1002 & \cellcolor[rgb]{ .976,  .439,  .447}0.0965 & \cellcolor[rgb]{ .976,  .478,  .486}0.0904 & \cellcolor[rgb]{ .98,  .588,  .596}0.0732 & \cellcolor[rgb]{ .976,  .514,  .522}0.0849 & \cellcolor[rgb]{ .984,  .729,  .737}0.0514 & \cellcolor[rgb]{ .988,  .969,  .98}0.0139 & \cellcolor[rgb]{ .988,  .933,  .945}0.0194 & \cellcolor[rgb]{ .988,  .937,  .949}0.0188 \\
    \hline
    3     & \cellcolor[rgb]{ .976,  .494,  .502}0.0879 & \cellcolor[rgb]{ .98,  .561,  .569}0.0775 & \cellcolor[rgb]{ .976,  .506,  .514}0.0861 & \cellcolor[rgb]{ .984,  .729,  .737}0.0513 & \cellcolor[rgb]{ .984,  .831,  .843}0.0352 & \cellcolor[rgb]{ .984,  .831,  .843}0.0354 & \cellcolor[rgb]{ .988,  .89,  .902}0.0260 & \cellcolor[rgb]{ .988,  .988,  1}0.0111 & \cellcolor[rgb]{ .988,  .984,  .996}0.0117 \\
    \hline
    6     & \cellcolor[rgb]{ .98,  .6,  .612}0.0710 & \cellcolor[rgb]{ .984,  .714,  .722}0.0537 & \cellcolor[rgb]{ .984,  .718,  .729}0.0529 & \cellcolor[rgb]{ .988,  .902,  .914}0.0241 & \cellcolor[rgb]{ .988,  .965,  .976}0.0144 & \cellcolor[rgb]{ .988,  .89,  .898}0.0265 & \cellcolor[rgb]{ .988,  .961,  .973}0.0153 & \cellcolor[rgb]{ .671,  .765,  .886}0.0056 & \cellcolor[rgb]{ .988,  .988,  1}0.0107 \\
    \hline
    12    & \cellcolor[rgb]{ .984,  .71,  .722}0.0541 & \cellcolor[rgb]{ .984,  .831,  .843}0.0351 & \cellcolor[rgb]{ .984,  .816,  .827}0.0376 & \cellcolor[rgb]{ .988,  .863,  .875}0.0302 & \cellcolor[rgb]{ .969,  .973,  .992}0.0104 & \cellcolor[rgb]{ .988,  .906,  .918}0.0235 & \cellcolor[rgb]{ .988,  .988,  1}0.0111 & \cellcolor[rgb]{ .988,  .973,  .984}0.0135 & \cellcolor[rgb]{ .678,  .773,  .89}0.0058 \\
    \hline
    24    & \cellcolor[rgb]{ .988,  .855,  .863}0.0320 & \cellcolor[rgb]{ .988,  .91,  .922}0.0233 & \cellcolor[rgb]{ .988,  .969,  .98}0.0141 & \cellcolor[rgb]{ .988,  .859,  .871}0.0311 & \cellcolor[rgb]{ .988,  .957,  .969}0.0160 & \cellcolor[rgb]{ .643,  .745,  .878}0.0052 & \cellcolor[rgb]{ .698,  .78,  .894}0.0061 & \cellcolor[rgb]{ .588,  .706,  .859}0.0043 & \cellcolor[rgb]{ .561,  .686,  .847}0.0039 \\
    \hline
    48    & \cellcolor[rgb]{ .988,  .945,  .957}0.0178 & \cellcolor[rgb]{ .631,  .733,  .871}0.0050 & \cellcolor[rgb]{ .988,  .953,  .965}0.0162 & \cellcolor[rgb]{ .988,  .988,  1}0.0111 & \cellcolor[rgb]{ .565,  .69,  .851}0.0040 & \cellcolor[rgb]{ .467,  .624,  .816}0.0024 & \cellcolor[rgb]{ .596,  .71,  .859}0.0044 & \cellcolor[rgb]{ .569,  .694,  .851}0.0041 & \cellcolor[rgb]{ .478,  .627,  .82}0.0026 \\
    \hline
    96    & \cellcolor[rgb]{ .859,  .898,  .953}0.0087 & \cellcolor[rgb]{ .549,  .678,  .843}0.0037 & \cellcolor[rgb]{ .51,  .651,  .831}0.0031 & \cellcolor[rgb]{ .529,  .663,  .835}0.0034 & \cellcolor[rgb]{ .455,  .612,  .812}0.0022 & \cellcolor[rgb]{ .408,  .576,  .792}0.0015 & \cellcolor[rgb]{ .443,  .604,  .808}0.0021 & \cellcolor[rgb]{ .376,  .557,  .784}0.0010 & \cellcolor[rgb]{ .392,  .569,  .788}0.0012 \\
    \hline
    192   & \cellcolor[rgb]{ .718,  .8,  .906}0.0064 & \cellcolor[rgb]{ .62,  .729,  .871}0.0049 & \cellcolor[rgb]{ .424,  .588,  .8}0.0017 & \cellcolor[rgb]{ .38,  .561,  .784}0.0010 & \cellcolor[rgb]{ .38,  .561,  .784}0.0010 & \cellcolor[rgb]{ .384,  .565,  .788}0.0011 & \cellcolor[rgb]{ .388,  .565,  .788}0.0012 & \cellcolor[rgb]{ .373,  .553,  .78}0.0009 & \cellcolor[rgb]{ .369,  .553,  .78}0.0009 \\
    \hline
    384   & \cellcolor[rgb]{ .604,  .718,  .863}0.0046 & \cellcolor[rgb]{ .478,  .631,  .82}0.0026 & \cellcolor[rgb]{ .384,  .565,  .788}0.0011 & \cellcolor[rgb]{ .384,  .561,  .784}0.0011 & \cellcolor[rgb]{ .38,  .561,  .784}0.0010 & \cellcolor[rgb]{ .357,  .545,  .776}0.0007 & \cellcolor[rgb]{ .369,  .553,  .78}0.0008 & \cellcolor[rgb]{ .357,  .541,  .776}0.0007 & \cellcolor[rgb]{ .353,  .541,  .776}0.0006 \\
    \hline
\end{tabular}
\label{table:darcy_forward_nu_nc}
\end{table}
}


\rev{Next, we perform a sensitivity analysis to explore how the PINN performance depends on the hyperparameters, $N_{hl}$ and $N_{n}$. We do this using a fixed size of the training set with $N_b$ and $N_{\Pi}$ equal to 96 and 160, respectively. Ideally, the explored space of the hyperparameters should reveal a parameter set that produces a minimum of the loss function on a validation set. To deal with the stochastic nature of the training procedure, again, we calculate average values obtained over many realizations. The results are presented in Table \ref{table:darcy_forward_nl_nn}. From panel C, we observe that selecting $N_{hl}=6$ and $N_{n}=5$ corresponds to the best accuracy in the explored space. For a smaller and larger size of the architecture, the accuracy decreases, corresponding to underfitting and overfitting, respectively \cite{tetko1995neural}. We now apply the trained PINN model using $N_b = 96$, $N_{\Pi} = 160$, $N_{hl} = 6$, and $N_{n} = 5$ to the test set. The relative $\mathcal{L}^2$ error is 1.43 $\times$ $10^{-3}$, which is comparable to that of the validation set (see Table \ref{table:darcy_forward_nl_nn} - Panel C).}



\begin{table}[!ht]
\centering
\caption{
{\bf Diffusivity equation - forward modeling: PINN performance as function of hyperparameters.} Relative ${\mathcal{L}^2}$ errors between the exact and predicted values of ${p}$ for the \rev{validation} set. The table shows the dependency on the different number of hidden layers $N_{hl}$ and different number of neurons per layer $N_{n}$. Here, the total number of training and collocation points is fixed to $N_b$ = 96 and $N_{\Pi}$ = 160, respectively.}
\subcaption*{Panel A: Average over 1 realization}
\begin{tabular}{|l|l|l|l|l|l|l|}
\hline
\backslashbox{$N_{hl}$}{$N_n$}    & 2        & 5        & 10       & 20       & 40       & 80       \\ \hline
    2     & \cellcolor[rgb]{ .549,  .678,  .843}0.0009 & \cellcolor[rgb]{ .741,  .816,  .914}0.0015 & \cellcolor[rgb]{ .988,  .976,  .988}0.0062 & \cellcolor[rgb]{ .361,  .545,  .776}0.0003 & \cellcolor[rgb]{ .353,  .541,  .776}0.0003 & \cellcolor[rgb]{ .988,  .973,  .984}0.0088 \\
    \hline
    4     & \cellcolor[rgb]{ .988,  .988,  1}0.0033 & \cellcolor[rgb]{ .988,  .984,  .996}0.0036 & \cellcolor[rgb]{ .988,  .969,  .98}0.0098 & \cellcolor[rgb]{ .561,  .686,  .847}0.0009 & \cellcolor[rgb]{ .988,  .953,  .965}0.0144 & \cellcolor[rgb]{ .482,  .631,  .82}0.0007 \\
    \hline
    6     & \cellcolor[rgb]{ .373,  .553,  .78}0.0003 & \cellcolor[rgb]{ .988,  .984,  .996}0.0038 & \cellcolor[rgb]{ .8,  .855,  .933}0.0017 & \cellcolor[rgb]{ .988,  .945,  .957}0.0175 & \cellcolor[rgb]{ .988,  .918,  .929}0.0268 & \cellcolor[rgb]{ .553,  .682,  .847}0.0009 \\
    \hline
    8     & \cellcolor[rgb]{ .945,  .957,  .984}0.0021 & \cellcolor[rgb]{ .588,  .706,  .859}0.0010 & \cellcolor[rgb]{ .439,  .604,  .808}0.0005 & \cellcolor[rgb]{ .427,  .592,  .8}0.0005 & \cellcolor[rgb]{ .467,  .624,  .816}0.0006 & \cellcolor[rgb]{ .988,  .98,  .992}0.0058 \\
    \hline
    16    & \cellcolor[rgb]{ .988,  .875,  .886}0.0411 & \cellcolor[rgb]{ .38,  .561,  .784}0.0003 & \cellcolor[rgb]{ .988,  .965,  .976}0.0105 & \cellcolor[rgb]{ .498,  .643,  .827}0.0007 & \cellcolor[rgb]{ .988,  .945,  .957}0.0175 & \cellcolor[rgb]{ .863,  .898,  .953}0.0019 \\
    \hline
    32    & \cellcolor[rgb]{ .984,  .78,  .788}0.0735 & \cellcolor[rgb]{ .882,  .914,  .961}0.0019 & \cellcolor[rgb]{ .988,  .976,  .988}0.0073 & \cellcolor[rgb]{ .973,  .412,  .42}0.1970 & \cellcolor[rgb]{ .988,  .988,  1}0.0024 & \cellcolor[rgb]{ .988,  .973,  .984}0.0081 \\
    \hline
\end{tabular}
\bigskip
\subcaption*{Panel B: Average over 3 realizations}
\begin{tabular}{|l|l|l|l|l|l|l|}
\hline
\backslashbox{$N_{hl}$}{$N_n$}    & 2        & 5        & 10       & 20       & 40       & 80       \\ \hline
   2     & \cellcolor[rgb]{ .988,  .988,  1}0.0027 & \cellcolor[rgb]{ .678,  .773,  .89}0.0015 & \cellcolor[rgb]{ .867,  .902,  .957}0.0021 & \cellcolor[rgb]{ .988,  .988,  1}0.0042 & \cellcolor[rgb]{ .988,  .988,  1}0.0025 & \cellcolor[rgb]{ .988,  .988,  1}0.0026 \\
    \hline
    4     & \cellcolor[rgb]{ .353,  .541,  .776}0.0004 & \cellcolor[rgb]{ .988,  .988,  1}0.0025 & \cellcolor[rgb]{ .616,  .725,  .867}0.0013 & \cellcolor[rgb]{ .965,  .973,  .992}0.0024 & \cellcolor[rgb]{ .647,  .745,  .878}0.0014 & \cellcolor[rgb]{ .988,  .988,  1}0.0025 \\
    \hline
    6     & \cellcolor[rgb]{ .922,  .941,  .976}0.0023 & \cellcolor[rgb]{ .882,  .914,  .961}0.0021 & \cellcolor[rgb]{ .851,  .89,  .949}0.0020 & \cellcolor[rgb]{ .988,  .988,  1}0.0028 & \cellcolor[rgb]{ .988,  .988,  1}0.0026 & \cellcolor[rgb]{ .984,  .984,  .996}0.0025 \\
    \hline
    8     & \cellcolor[rgb]{ .988,  .988,  1}0.0049 & \cellcolor[rgb]{ .549,  .678,  .843}0.0010 & \cellcolor[rgb]{ .929,  .945,  .976}0.0023 & \cellcolor[rgb]{ .741,  .816,  .914}0.0017 & \cellcolor[rgb]{ .737,  .812,  .91}0.0016 & \cellcolor[rgb]{ .847,  .886,  .949}0.0020 \\
    \hline
    16    & \cellcolor[rgb]{ .973,  .412,  .42}0.5580 & \cellcolor[rgb]{ .988,  .988,  1}0.0029 & \cellcolor[rgb]{ .6,  .714,  .863}0.0012 & \cellcolor[rgb]{ .988,  .988,  1}0.0030 & \cellcolor[rgb]{ .827,  .875,  .941}0.0019 & \cellcolor[rgb]{ .988,  .988,  1}0.0027 \\
    \hline
    32    & \cellcolor[rgb]{ .988,  .875,  .886}0.1150 & \cellcolor[rgb]{ .988,  .957,  .969}0.0328 & \cellcolor[rgb]{ .988,  .988,  1}0.0036 & \cellcolor[rgb]{ .988,  .988,  1}0.0034 & \cellcolor[rgb]{ .988,  .98,  .992}0.0107 & \cellcolor[rgb]{ .91,  .933,  .973}0.0022 \\
    \hline
\end{tabular}
\bigskip
\subcaption*{Panel C: Average over 24 realizations}
\begin{tabular}{|l|l|l|l|l|l|l|}
\hline
\backslashbox{$N_{hl}$}{$N_n$}    & 2        & 5        & 10       & 20       & 40       & 80       \\ \hline
    2     & \cellcolor[rgb]{ .973,  .412,  .42}0.1660 & \cellcolor[rgb]{ .882,  .914,  .961}0.0022 & \cellcolor[rgb]{ .988,  .988,  1}0.0032 & \cellcolor[rgb]{ .988,  .988,  1}0.0025 & \cellcolor[rgb]{ .988,  .988,  1}0.0025 & \cellcolor[rgb]{ .741,  .812,  .91}0.0020 \\
    \hline
    4     & \cellcolor[rgb]{ .988,  .937,  .949}0.0171 & \cellcolor[rgb]{ .976,  .976,  .992}0.0024 & \cellcolor[rgb]{ .62,  .729,  .871}0.0018 & \cellcolor[rgb]{ .78,  .843,  .925}0.0021 & \cellcolor[rgb]{ .6,  .714,  .863}0.0018 & \cellcolor[rgb]{ .831,  .878,  .945}0.0022 \\
    \hline
    6     & \cellcolor[rgb]{ .988,  .988,  1}0.0032 & \cellcolor[rgb]{ .353,  .541,  .776}0.0014 & \cellcolor[rgb]{ .6,  .714,  .863}0.0018 & \cellcolor[rgb]{ .988,  .988,  1}0.0024 & \cellcolor[rgb]{ .745,  .816,  .914}0.0020 & \cellcolor[rgb]{ .933,  .949,  .98}0.0023 \\
    \hline
    8     & \cellcolor[rgb]{ .988,  .91,  .922}0.0253 & \cellcolor[rgb]{ .467,  .62,  .816}0.0016 & \cellcolor[rgb]{ .682,  .773,  .89}0.0019 & \cellcolor[rgb]{ .729,  .808,  .91}0.0020 & \cellcolor[rgb]{ .8,  .855,  .933}0.0021 & \cellcolor[rgb]{ .988,  .988,  1}0.0025 \\
    \hline
    16    & \cellcolor[rgb]{ .98,  .678,  .686}0.0912 & \cellcolor[rgb]{ .667,  .761,  .886}0.0019 & \cellcolor[rgb]{ .678,  .773,  .89}0.0019 & \cellcolor[rgb]{ .675,  .769,  .89}0.0019 & \cellcolor[rgb]{ .988,  .988,  1}0.0025 & \cellcolor[rgb]{ .988,  .988,  1}0.0033 \\
    \hline
    32    & \cellcolor[rgb]{ .984,  .835,  .847}0.0462 & \cellcolor[rgb]{ .988,  .969,  .98}0.0090 & \cellcolor[rgb]{ .988,  .984,  .996}0.0037 & \cellcolor[rgb]{ .988,  .988,  1}0.0027 & \cellcolor[rgb]{ .988,  .984,  .996}0.0036 & \cellcolor[rgb]{ .988,  .988,  1}0.0031 \\
    \hline
\end{tabular}
\label{table:darcy_forward_nl_nn}
\end{table}


\rev{
For the selected model ($N_{hl} = 6$, $N_{n} = 5$, $N_b$ = 96, and $N_{\Pi}$ = 160), we illustrate the behavior of the mean square training error of this model as function of the training iterations in Fig \ref{fig:darcy_loss_list}. One can observe that $M S E_{b}$ is higher than $M S E_{\Pi}$ in the beginning, but later it approaches zero faster than $M S E_{\Pi}$. The $MSE$ converges steadily without any oscillations.
}

\begin{figure}[!ht]
   \centering
        \includegraphics[width=8.0cm, height=8.0cm]{pictures/loss_list_darcy_forward.pdf}
   \caption{
{\bf Diffusivity equation - forward modeling: mean square training error plot.} This model uses $N_{hl} = 6$, $N_{n} = 5$, $N_b$ = 96, and $N_{\Pi}$ = 160. $M S E$, $M S E_{b}$, and $M S E_{\Pi}$ are calculated using Eqs. \ref{eq:loss_darcy_forward}, \ref{eq:loss_p_darcy_forward}, and \ref{eq:loss_pif_darcy_forward}, respectively.}
   \label{fig:darcy_loss_list}
\end{figure}


\subsubsection*{Inverse model of diffusivity equation}

We rewrite Eq~(\ref{eq:pif_darcy_forward}) in the parametrized form \cite{hesthaven2016certified} as presented below:

\begin{equation} \label{eq:darcy_para_inv}
\Pi(x,t) :=\theta_1 \frac{\partial p}{\partial t}-\theta_2\frac{\partial }{\partial x} p^2(\frac{\partial }{\partial x} p)-g,
\end{equation}

\noindent
where

\begin{equation} \label{eq:darcy_theta}
\begin{aligned}
\theta_1  = \phi c_{t}  \quad \text{and} \quad \theta_2  = \bm{\kappa}_0.
\end{aligned}
\end{equation}

\noindent
Unlike the forward mode, where the weights and biases are the unknown parameters to be learned, we now have two more unknowns, $\theta_1$ and $\theta_2$. Using Eq (\ref{eq:loss_gen}) we have:



\begin{equation}\label{eq:loss_p_darcy_inverse}
M S E_{{tr}}=\frac{1}{N_{{tr}}} \sum_{i=1}^{N_{{tr}}}\left|p\left(x_{{tr}}^{i} t_{{tr}}^{i}\right)-p^{i}\right|^{2}
\end{equation}

\noindent
and

\begin{equation}\label{eq:loss_pif_darcy_inverse}
M S E_{\Pi}=\frac{1}{N_{{tr}}} \sum_{i=1}^{N_{{tr}}}\left|\Pi\left(x_{{tr}}^{i}, t_{{tr}}^{i}\right)\right|^{2},
\end{equation}

\noindent
where $\left\{x_{{tr}}^{i}, t_{{tr}}^{i}, p^{i}\right\}_{i=1}^{N_{{tr}}}$ refers to the set of training data and $N_{tr}$ is the number of training data. \rev{ In contrast to the forward model, we do not need to specify specific collocation points to activate the $\Pi$ dependent loss term; here we can just use the training examples used to calculate $M S E_{tr}$}. All physical constants are set to one; hence, $\theta_1 = 1.0$ and $\theta_2 = 1.0$. \rev{Note that $\theta_1$ and $\theta_2$ are considered constant throughout the domain.}

\rev{
We use the hyperparameters obtained from the sensitivity analysis of the forward model (see Panel C of Table \ref{table:darcy_forward_nl_nn}), i.e., $N_{hl} = 6$ and $N_{n} = 5$. In Table \ref{table:darcy_inv_nn_nl}, we illustrate that this choice also yields the least $\mathcal{L}^2$ error of $p$ and percentage errors of $\theta_1$ and $\theta_2$ for the inverse problem, supporting the heuristic arguments presented above in the section "Training the PINN". Moreover, we find that the $\mathcal{L}^2$ error of $p$ and percentage error of $\theta_1$ and $\theta_2$ are not much different with different combinations of the hyperparameters. Note that the $\mathcal{L}^2$ error of $p$ and percentage errors of $\theta_1$ and $\theta_2$ presented in Table \ref{table:darcy_inv_nn_nl} represent average values over 24 realizations.
}

\begin{table}[htbp]
  \centering
  \caption{{\bf Diffusivity equation - inverse modeling: PINN performance as function of hyperparameters.} Relative ${\mathcal{L}^2}$ error of ${p}$ and percentage error of ${\theta_1}$ and ${\theta_2}$ for different number of hidden layers, ${N_{hl}}$, and different number of neurons per layer, ${N_{n}}$. The $N_{tr}$ is fixed at 250. Note that we pick the optimal hyperparameters, i.e., $N_{hl} = 6$ and $N_{n} = 5$ from the sensitivity analysis of the forward model. Results shown in this table are an average over 24 realizations.}
    \begin{tabular}{|l|l|l|l|l|l|l|}
\hline
&\backslashbox{$N_{hl}$}{$N_n$}    & 2        & 5        & 10\\
    \hline
    \multirow{3}[6]{*}{$p$} & 4     & \cellcolor[rgb]{ .988,  .922,  .933}0.0004 & \cellcolor[rgb]{ .824,  .871,  .941}0.0003 & \cellcolor[rgb]{ .988,  .988,  1}0.0003 \\
\cline{2-5}          & 6     & \cellcolor[rgb]{ .98,  .639,  .651}0.0007 & \cellcolor[rgb]{ .353,  .541,  .776}0.0002 & \cellcolor[rgb]{ .988,  .914,  .925}0.0004 \\
\cline{2-5}          & 8     & \cellcolor[rgb]{ .973,  .412,  .42}0.0010 & \cellcolor[rgb]{ .482,  .631,  .82}0.0002 & \cellcolor[rgb]{ .745,  .816,  .914}0.0003 \\
    \hline
    \multirow{3}[6]{*}{$\theta_1$} & 4     & \cellcolor[rgb]{ .988,  .988,  1}0.1888 & \cellcolor[rgb]{ .576,  .698,  .855}0.1357 & \cellcolor[rgb]{ .988,  .976,  .988}0.1969 \\
\cline{2-5}          & 6     & \cellcolor[rgb]{ .984,  .761,  .773}0.3200 & \cellcolor[rgb]{ .353,  .541,  .776}0.1065 & \cellcolor[rgb]{ .984,  .749,  .761}0.3261 \\
\cline{2-5}          & 8     & \cellcolor[rgb]{ .973,  .412,  .42}0.5195 & \cellcolor[rgb]{ .51,  .651,  .831}0.1272 & \cellcolor[rgb]{ .396,  .573,  .792}0.1125 \\
    \hline
    \multirow{3}[6]{*}{$\theta_2$} & 4     & \cellcolor[rgb]{ .988,  .988,  1}0.3443 & \cellcolor[rgb]{ .988,  .976,  .988}0.3562 & \cellcolor[rgb]{ .922,  .941,  .976}0.3289 \\
\cline{2-5}          & 6     & \cellcolor[rgb]{ .976,  .522,  .529}0.7121 & \cellcolor[rgb]{ .353,  .541,  .776}0.1912 & \cellcolor[rgb]{ .976,  .541,  .549}0.6946 \\
\cline{2-5}          & 8     & \cellcolor[rgb]{ .973,  .412,  .42}0.7949 & \cellcolor[rgb]{ .839,  .882,  .945}0.3088 & \cellcolor[rgb]{ .596,  .714,  .863}0.2504 \\
    \hline
    \end{tabular}%
  \label{table:darcy_inv_nn_nl}%
\end{table}%

\rev{
We depict the performance of the PINN model for solving the inverse problem as a function of $N_{tr}$ in Fig~\ref{fig:darcy_inv_find_N_f}. The reported error values of $\theta_1$ and $\theta_2$ are percentage errors, while the relative $\mathcal{L}^2$ error is shown for $p$. The error bars show the standard derivation ($\pm$1 SD) based on 24 realizations. We observe that a minimum of $N_{tr}=200$ is required by the PINN model to avoid substantial stochastic fluctuations in the estimated values of $\theta_1$ and $\theta_2$. Moreover, we observe that the PINN model with $N_{tr}>200$ provides average percentage errors of $\theta_1$ and $\theta_2$ less than 1 \%, but with a large standard deviation. Over 24 realizations of the trained PINN with $N_{tr}=200$, the minimum percentage errors are $2.31 \times 10^{-2}$ and $7.33 \times 10^{-2}$ and the maximum are $1.93 \times 10^{-1}$ and $4.92 \times 10^{-1}$, for $\theta_1$ and $\theta_2$, respectively.} \rev{To obtain 200 training examples in a real setting, i.e., lab experiments or field observations, is realistic; hence, the results illustrate the feasibility of PINN to solve the inverse problem based on a reasonably sized data set.} 



\begin{figure}[!ht]
   \centering
        \includegraphics[width=8.0cm, height=8.0cm]{pictures/darcy_inv_find_N_f.pdf}
   \caption{
{\bf Diffusivity equation - inverse modeling: PINN performance as function of training set size.} This figure shows the estimated error dependency on the amount of training data, ${N_{tr}}$. \rev{The error bars show mean and standard derivation ($\pm$ 1 SD) based on 24 realizations.} Note that there is no noise added in this investigation, and all physical constants are set to one, i.e., $\theta_1 =$ $\theta_2 = 1.0$. The reported error values of $\theta_1$ and $\theta_2$ are percentage errors while the relative $\mathcal{L}^2$ error is shown for $p$.}
   \label{fig:darcy_inv_find_N_f}
\end{figure}

Next, we perform a systematic study of the effect of \rev{additive noise} in data, which is created from the true data as follows \cite{raissi2019physics}:

\begin{equation}\label{eq:darcy_noise}
\bm{X}_{noise} = \bm{X}_{true} + \epsilon \mathcal{S}\left(\bm{X}_{true}\right) \mathcal{G}\left(0,1\right),
\end{equation}

\noindent
where $\bm{X}_{noise}$ and $\bm{X}_{true}$ is the vector of the data with and without noise, respectively. The $\epsilon$ determines the noise level, $\mathcal{S}\left(\cdot\right)$ represents a standard deviation operator, $\mathcal{G}\left(0,1\right)$ is a random value, which is sampled from the Gaussian distribution with mean and standard deviation of zero and one, respectively. The noise generated from this procedure is fully random and uncorrelated. The results are presented in Table \ref{table:darcy_inv_noise} \rev{as average values over 10 training realizations. We can observe that the error increases with the noise level ($\epsilon$) while the error decreases as the $N_{tr}$ is increased.} As expected, the PINN model requires more data to accurately approximate the unknown physical parameters when the noise level is high.


\begin{table}[!ht]
\centering
\caption{
{\bf Diffusivity equation - inverse modeling: PINN performance as function of noise.} This figure shows the average percentage errors of \rev{${\theta_1}$ and ${\theta_2}$ for different numbers of training data, ${N_{tr}}$, as function of the noise levels.} Here, the neural network architecture is kept fixed to 6 layers and 5 neurons per layer. The results are averages over 10 realizations.}
 \begin{tabular}{|l|l|l|l|l|l|l|}
\hline
\multicolumn{2}{|l|}{{\backslashbox{$N_{tr}$}{Noise ($\epsilon$)}} }                     & 0\% & 1\% & 5\% & 10\% \\ \hline
\multirow{4}[8]{*}{$\theta_1$} & 100   & \cellcolor[rgb]{ .506,  .651,  .831}0.17 & \cellcolor[rgb]{ .984,  .812,  .824}1.82 & \cellcolor[rgb]{ .976,  .451,  .459}4.40 & \cellcolor[rgb]{ .973,  .412,  .42}4.67 \\
\cline{2-6}          & 250   & \cellcolor[rgb]{ .486,  .635,  .824}0.15 & \cellcolor[rgb]{ .667,  .761,  .886}0.30 & \cellcolor[rgb]{ .988,  .929,  .941}1.00 & \cellcolor[rgb]{ .984,  .792,  .8}1.98 \\
\cline{2-6}          & 500   & \cellcolor[rgb]{ .443,  .604,  .808}0.12 & \cellcolor[rgb]{ .51,  .651,  .831}0.17 & \cellcolor[rgb]{ .988,  .961,  .973}0.77 & \cellcolor[rgb]{ .988,  .953,  .961}0.84 \\
\cline{2-6}          & 1000  & \cellcolor[rgb]{ .353,  .541,  .776}0.04 & \cellcolor[rgb]{ .416,  .584,  .796}0.09 & \cellcolor[rgb]{ .722,  .8,  .906}0.34 & \cellcolor[rgb]{ .988,  .937,  .949}0.94 \\
    \hline
    \multirow{4}[8]{*}{$\theta_2$} & 100   & \cellcolor[rgb]{ .463,  .62,  .816}0.22 & \cellcolor[rgb]{ .988,  .878,  .886}1.49 & \cellcolor[rgb]{ .976,  .49,  .498}4.35 & \cellcolor[rgb]{ .973,  .412,  .42}4.90 \\
\cline{2-6}          & 250   & \cellcolor[rgb]{ .482,  .631,  .82}0.24 & \cellcolor[rgb]{ .82,  .871,  .941}0.52 & \cellcolor[rgb]{ .984,  .835,  .843}1.80 & \cellcolor[rgb]{ .984,  .745,  .757}2.45 \\
\cline{2-6}          & 500   & \cellcolor[rgb]{ .467,  .62,  .816}0.23 & \cellcolor[rgb]{ .769,  .831,  .922}0.47 & \cellcolor[rgb]{ .988,  .945,  .957}0.99 & \cellcolor[rgb]{ .988,  .878,  .89}1.48 \\
\cline{2-6}          & 1000  & \cellcolor[rgb]{ .353,  .541,  .776}0.13 & \cellcolor[rgb]{ .537,  .671,  .839}0.28 & \cellcolor[rgb]{ .686,  .776,  .894}0.41 & \cellcolor[rgb]{ .988,  .973,  .984}0.79 \\
    \hline
\end{tabular}
\label{table:darcy_inv_noise}
\end{table}

\subsection*{Nonlinear Biot's equations}

\rev{From the nonlinear diffusivity equation section, we have shown that the PINN model can solve forward and inverse problems. We then progress to the multiphysics problem represented by the nonlinear Biot's equations. We take $\Omega = \left[ 0,1\right]^2$, $\mathbb{T} = \left[0,1\right]$, and choose the exact solution in $\Omega$ as:}

\begin{equation} \label{eq:biot_displacement_exact}
\bm{u}(x,y,t):=\left[ \begin{array}{cc}  { u} \\ { v}\end{array}\right] = \left[ \begin{array}{cc}  { \sin(x+y+t)} \\ { \cos(x+y+t)}\end{array}\right],
\end{equation}

\noindent
for the displacement variable where $u$ and $v$ are displacements in x- and y-direction, respectively. Note that as we focus on the 2-Dimensional domain; therefore, $\bm{u}(x,y,t)$ is composed of two spatial components. For the pressure variable, we choose 


\begin{equation}\label{eq:biot_pressure_exact}
p(x,y,t) := e^{(x+y+t)}.
\end{equation}

\noindent
Here $x$, $y$, and $t$ represent points in x-, y-direction, and time domain, respectively. The $\mathcal{N}[\bm{\kappa}]$ is chosen as:

\begin{equation}  \label{eq:perm_biot}
\mathcal{N}[\bm{\kappa}] := \bm{\kappa_0}e^{\varepsilon_{v}},
\end{equation}

\noindent
where $\bm{\kappa_0}$ represent initial rock matrix conductivity. Again, we assume $\bm{\kappa_0}$ to be a scalar in this case, i.e., $\bm{\kappa_0}={\kappa_0}$. The $\varepsilon_{v}$ is the total volumetric strain defined as:

\begin{equation} \label{eq:volumetric_str}
\varepsilon_{v} :=\operatorname{tr}({\bm{\varepsilon}})=\sum_{i=1}^{2} \varepsilon_{i i}.
\end{equation}

\noindent
The choice of $\mathcal{N}[\bm{\kappa}]$ function is selected to represent the change in a volumetric strain that affects the porous media conductivity, and it is adapted from \cite{Du2007, abou2013petroleum, kadeethum2019comparison}. All the physical constants are set to $1.0$; and subsequently, $\bm{f}$ is chosen as: 

\begin{equation} \label{eq:biot_displacement_source}
\bm{f}(x,y,t):=\left[ \begin{array}{cc}{f_u(x,y,t)} \\ {f_v(x,y,t)}\end{array}\right],
\end{equation}

\noindent
where

\begin{equation}\label{eq:biot_displacement_source1}
\begin{aligned}
f_u(x,y,t) & := -4.0\sin(x+y+t)-2.0\cos(x+y+t) -  e^{(x+y+t)},
\end{aligned}
\end{equation}

\noindent
and

\begin{equation}\label{eq:biot_displacement_source2}
\begin{aligned}
f_v(x,y,t) & := -4.0\cos(x+y+t)-2.0\sin(x+y+t) -  e^{(x+y+t)},
\end{aligned}
\end{equation}


\noindent
for the momentum balance equation, Eq~(\ref{eq:linear_balance}). The source term of the mass balance equation, Eq~(\ref{eq:mass_balance}), $g$ is chosen as: 


\begin{equation}  \label{eq:biot_mass_source}
\begin{aligned}
g(x,y,t) &:= \left( \cos \left( x+y+t \right) +\sin \left( x+y+t \right) -1
 \right) {{\rm e}^{\cos \left( x+y+t \right) -\sin \left( x+y+t
 \right) +x+y+t}} \\
 &-\cos \left( x+y+t \right) +{{\rm e}^{x+y+t}}-\sin
 \left( x+y+t \right), 
\end{aligned}
\end{equation}


\noindent
to satisfy the exact solution. Furthermore, the boundary conditions and initial conditions are applied using Eqs~(\ref{eq:biot_displacement_exact}) and (\ref{eq:biot_pressure_exact}). The $\Pi_{\bm{u}}(x,y,t)$ and $\Pi_{p}(x,y,t)$, Eqs (\ref{eq:linear_balance_pinn}) and (\ref{eq:mass_balance_pinn}), \rev{here act as the physics-informed function.} 

\rev{
We generate the exact solution points, Eqs \eqref{eq:biot_displacement_exact} and \eqref{eq:biot_pressure_exact}, based on a rectangular mesh ($\Omega = [0,1]^2$) with 99 equidistant intervals in both x- and y-direction, i.e. $\Delta x$ $=$ $\Delta y$. Using 49 equidistant temporal intervals, in total, we have 500000 examples. Similar to the diffusivity equation case, we draw $n$ training examples randomly. Subsequently, we split the remaining examples equally for validation and test sets. Again, assuming we have 500000 solution points for the sake of illustration, we use 100 examples to train the model; we then have 249950 examples for both the validation and the test sets.} 






\rev{
To recap, the forward modeling of Biot's system aims to predict the displacement ($\bm{u}$) and pressure ($p$) by specifying the initial and boundary conditions, collocation points, and as well as the physical parameters ($\mu_l$, $\lambda_l$, $\alpha$, $\phi$, $c_f$, $K_{s}$, and ${\kappa}_0$). The inverse modeling, however, aims to estimate the physical parameters from observed values of $\bm{u}$ and $p$ with their corresponding values of $x$, $y$, and $t$.
}

\rev{
}


\rev{In the case of the nonlinear Biot's equations, the architecture of the neural network corresponding to the top of Fig \ref{fig:des_for_inv} is presented in Fig~\ref{fig:biot_nn}. We have three input nodes and three output nodes for this case. For the forward modeling, the hyperparameters $N_{hl}$ and $N_n$ are found using a sensitivity analysis, and as argued in the above section, "Training the PINN," we apply the same hyperparameters for the inverse model. Again, we use L-BFGS; a quasi-Newton, fullbatch gradient-based optimization algorithm to minimize the loss function \cite{liu1989limited} for the forward model. For the inverse problem, we find that combining ADAM, stochastic gradient descent, and L-BFGS might provide faster convergence when training the neural network. Specifically, we use ADAM for the first 10000 iterations and then continue using L-BFGS until the stop criterion is met. Note that we apply the same stop criterion as described above for the diffusivity case. Since the ADAM is a first-order method compared to L-BFGS, which is a second-order model, ADAM is less computationally demanding. Initially, where the weights of the neural network are far from convergence, we speculate that the less computational effort of ADAM is an advantage. However, as we approach the minimum, L-BFGS is likely to provide a better estimate of the steepest descent. Whether these observations could be made when dealing with other types of partial differential equations is an open question, but the use of a similar combined optimization scheme has been reported in the literature for the case of the Navier-Stokes equations \cite{raissi2019physics}. }

\begin{figure}[!ht]
   \centering
        \includegraphics[width=9.5cm, height=6.0cm]{pictures/biot_nn.pdf}
   \caption{
{\bf Neural networks architecture used for nonlinear Biot's equations.} This figure corresponds to the top part of Fig \ref{fig:des_for_inv}. There are three inputs, $x$, $y$, and $t$, and three outputs, $u$, $v$, and $p$. The number of hidden layers, $N_{hl}$, and the number of neurons for each hidden layer, $N_n$, denote the hyperparameters.}
   \label{fig:biot_nn}
\end{figure}     




\subsubsection*{Verification of the forward model}

Again, we apply Eq. (\ref{eq:loss_gen}) and obtain:

\begin{equation}\label{eq:loss_biot_forward}
M S E=M S E_{b}+M S E_{\Pi_{\bm{u}}} + M S E_{\Pi_{p}}
\end{equation}

\noindent
where 

\begin{equation}\label{eq:loss_p_biot_forward}
\begin{aligned}
M S E_{b}&=\frac{1}{N_{b}} \sum_{i=1}^{N_{b}}(\left|u\left(x_{b}^{i}, y_{b}^{i}, t_{b}^{i} \right)-u^{i}\right|^{2} + \left|v\left(x_{b}^{i}, y_{b}^{i}, t_{b}^{i} \right)-v^{i}\right|^{2}\\ &+\left|p\left(x_{b}^{i}, y_{b}^{i}, t_{b}^{i} \right)-p^{i}\right|^{2}),
\end{aligned}
\end{equation}

\begin{equation}\label{eq:loss_pif_u_biot_forward}
M S E_{\Pi_{\bm{u}}}=\frac{1}{N_{\Pi_{\bm{u}}}} \sum_{i=1}^{N_{\Pi_{\bm{u}}}}\left|\Pi_{\bm{u}}\left(x_{\Pi_{\bm{u}}}^{i},y_{\Pi_{\bm{u}}}^{i}, t_{\Pi_{\bm{u}}}^{i}\right)\right|^{2},
\end{equation}

\noindent
and

\begin{equation}\label{eq:loss_pif_p_biot_forward}
M S E_{\Pi_{p}}=\frac{1}{N_{\Pi_{p}}} \sum_{i=1}^{N_{\Pi_{p}}}\left|\Pi_{p}\left(x_{\Pi_{p}}^{i},y_{\Pi_{p}}^{i}, t_{\Pi_{p}}^{i}\right)\right|^{2},
\end{equation}

\noindent
where $\left\{x_{b}^{i},y_{b}^{i},t_{b}^{i}, u^{i},v^{i},p^{i}\right\}_{i=1}^{N_{b}}$ refer to the initial and boundary training data. $\left\{x_{\Pi_{\bm{u}}}^{i}, y_{\Pi_{\bm{u}}}^{i},t_{\Pi_{\bm{u}}}^{i}\right\}_{i=1}^{N_{\Pi_{\bm{u}}}}$ and $\left\{x_{\Pi_p}^{i}, y_{\Pi_p}^{i},t_{\Pi_p}^{i}\right\}_{i=1}^{N_{\Pi_p}}$ \rev{ specify the collocation points for $\Pi_{\bm{u}}(x,y,t)$ and $\Pi_p(x,y,t)$, as defined in Eqs (\ref{eq:linear_balance_pinn}) and (\ref{eq:mass_balance_pinn}). Similar to the diffusivity equation case, these collocation points are sampled using the Latin Hypercube method \cite{baudin2015pydoe}.} $N_b$ denotes the number of initial and boundary training data, and $N_{\Pi_{\bm{u}}}$ and $N_{\Pi_p}$ are the number of collocation points for $\Pi_{\bm{u}}$ and $\Pi_p$, respectively. For the sake of simplification, in this investigation, we assume $N_{\Pi_{\bm{u}}}$ $=$ $N_{\Pi_p}$ $=$ $N_{\Pi}$ and $\left\{x_{\Pi_{\bm{u}}}^{i}, y_{\Pi_{\bm{u}}}^{i},t_{\Pi_{\bm{u}}}^{i}\right\}_{i=1}^{N_{\Pi_{\bm{u}}}}$ $=$ $\left\{x_{\Pi_p}^{i}, y_{\Pi_p}^{i},t_{\Pi_p}^{i}\right\}_{i=1}^{N_{\Pi_p}}$ $=$ $\left\{x_{\Pi}^{i}, y_{\Pi}^{i},t_{\Pi}^{i}\right\}_{i=1}^{N_{\Pi}}$.


\rev{In Fig \ref{fig:biot_exact_sol}, we illustrate an example of exact solutions of $u$, $v$, and $p$, and compare them to the prediction values obtained from our test set using the PINN trained with $N_b = 24$, $N_{\Pi} = 20$, $N_{hl} = 6$, and $N_{n} = 5$. This figure demonstrates that the PINN provides good approximations of the exact solutions. The dependency of the prediction accuracy on $N_b$ and $N_{\Pi}$ with the $N_{hl}$ and $N_n$ fixed to four and ten, respectively is illustrated in Table \ref{table:biot_forward_nu_nc}. Again, to deal with the stochastic behavior of the neural networks, we calculate an average of the relative $\mathcal{L}^2$ error over many realizations to obtain a clear pattern; see Panels B and C of Table \ref{table:biot_forward_nu_nc}. Similar to the diffusivity equation case, we observe the accuracy of the PINN is improved when $N_b$ and $N_{\Pi}$ are increased. We also note that the total amount of training examples required to achieve high accuracy, i.e., an average $\mathcal{L}^2$ error less than $1.0 \times 10^{-3}$, is in the order of 100s. Again, we can illustrate the impact of the $\Pi$ terms by considering the case where we train a neural network to interpolate a feed-forward solution based on a known set of solution points with and without making use of the regularization terms. Without using $\Pi$ we obtain an average relative $\mathcal{L}^2$ error of $1.1 \times 10^{-1}$ based on 96 solution points while the average $\mathcal{L}^2$ error becomes less than $1.0 \times 10^{-3}$ when including $\Pi$.}


\begin{figure}[!ht]
   \centering
        \includegraphics[width=5.0cm, height=4.0cm]{pictures/biot_u_sol_0.png}
        \includegraphics[width=5.0cm, height=4.0cm]{pictures/biot_u_sol_05.png}
        \includegraphics[width=5.0cm, height=4.0cm]{pictures/biot_u_sol_10.png}
        \includegraphics[width=5.0cm, height=4.0cm]{pictures/biot_v_sol_0.png}
        \includegraphics[width=5.0cm, height=4.0cm]{pictures/biot_v_sol_05.png}
        \includegraphics[width=5.0cm, height=4.0cm]{pictures/biot_v_sol_10.png}
        \includegraphics[width=5.0cm, height=4.0cm]{pictures/biot_p_sol_0.png}
        \includegraphics[width=5.0cm, height=4.0cm]{pictures/biot_p_sol_05.png}
        \includegraphics[width=5.0cm, height=4.0cm]{pictures/biot_p_sol_10.png}
   \caption{
{\bf Biot's equations - forward modeling: exact solutions (shown by surface plot) and 100 PINN predictions per time step using the test set (shown by black points).} The PINN was trained using $N_b = 24$, $N_{\Pi} = 20$, $N_{hl} = 6$, and $N_{n} = 5$. The top row illustrates the displacement in the x-direction, $u$, at (a) $t = 0.0$, (b) $t = 0.5$, and (c) $t = 1.0$. The middle row illustrates the displacement in the y-direction, $v$, at (d) $t = 0.0$, (e) $t = 0.5$, and (f) $t = 1.0$. The bottom row illustrates the pressure, $p$, at (g) $t = 0.0$, (h) $t = 0.5$, and (i) $t = 1.0$.}
   \label{fig:biot_exact_sol}
\end{figure}


\begin{table}[!ht]
\centering
\caption{
{\bf Biot's equations - forward modeling: PINN performance as function of training set size.} Sum of relative ${\mathcal{L}^2}$ errors between the exact and predicted values of ${u}$, ${v}$, and ${p}$ for the \rev{validation} set. The table shows the dependency on the number of the initial and boundary training data, $N_b$, and on the number of collocation points, $N_{\Pi}$. The hyperparemeters are fixed to 4 layers with 10 neurons per hidden layer.}
\subcaption*{Panel A: Average over 1 realization}
\begin{tabular}{|l|l|l|l|l|l|l|l|l|}
\hline
\backslashbox{$N_b$}{$N_{\Pi}$}             & 2      & 5      & 10     & 20     & 40     & 80     & 160    \\ \hline
   2     & \cellcolor[rgb]{ .976,  .482,  .49}0.4475 & \cellcolor[rgb]{ .976,  .51,  .518}0.4257 & \cellcolor[rgb]{ .973,  .412,  .42}0.5010 & \cellcolor[rgb]{ .98,  .682,  .694}0.2865 & \cellcolor[rgb]{ .98,  .686,  .698}0.2826 & \cellcolor[rgb]{ .988,  .863,  .871}0.1438 & \cellcolor[rgb]{ .984,  .792,  .804}0.1982 \\
    \hline
    3     & \cellcolor[rgb]{ .984,  .82,  .831}0.1770 & \cellcolor[rgb]{ .988,  .851,  .863}0.1508 & \cellcolor[rgb]{ .984,  .733,  .741}0.2464 & \cellcolor[rgb]{ .984,  .788,  .8}0.2009 & \cellcolor[rgb]{ .984,  .784,  .792}0.2059 & \cellcolor[rgb]{ .988,  .867,  .875}0.1410 & \cellcolor[rgb]{ .984,  .824,  .831}0.1751 \\
    \hline
    6     & \cellcolor[rgb]{ .976,  .533,  .541}0.4056 & \cellcolor[rgb]{ .976,  .463,  .471}0.4615 & \cellcolor[rgb]{ .988,  .937,  .949}0.0824 & \cellcolor[rgb]{ .416,  .584,  .796}0.0042 & \cellcolor[rgb]{ .988,  .886,  .898}0.1236 & \cellcolor[rgb]{ .988,  .875,  .886}0.1327 & \cellcolor[rgb]{ .984,  .835,  .843}0.1662 \\
    \hline
    12    & \cellcolor[rgb]{ .988,  .988,  1}0.0412 & \cellcolor[rgb]{ .988,  .933,  .945}0.0860 & \cellcolor[rgb]{ .988,  .953,  .965}0.0722 & \cellcolor[rgb]{ .404,  .576,  .792}0.0035 & \cellcolor[rgb]{ .357,  .545,  .776}0.0005 & \cellcolor[rgb]{ .369,  .553,  .78}0.0013 & \cellcolor[rgb]{ .761,  .827,  .918}0.0267 \\
    \hline
    24    & \cellcolor[rgb]{ .576,  .698,  .855}0.0148 & \cellcolor[rgb]{ .988,  .937,  .949}0.0824 & \cellcolor[rgb]{ .408,  .576,  .792}0.0037 & \cellcolor[rgb]{ .38,  .557,  .784}0.0019 & \cellcolor[rgb]{ .373,  .553,  .78}0.0015 & \cellcolor[rgb]{ .361,  .545,  .776}0.0007 & \cellcolor[rgb]{ .353,  .541,  .776}0.0003 \\
    \hline
    48    & \cellcolor[rgb]{ .988,  .976,  .988}0.0532 & \cellcolor[rgb]{ .42,  .588,  .8}0.0045 & \cellcolor[rgb]{ .392,  .569,  .788}0.0027 & \cellcolor[rgb]{ .361,  .545,  .776}0.0008 & \cellcolor[rgb]{ .38,  .557,  .784}0.0019 & \cellcolor[rgb]{ .353,  .541,  .776}0.0002 & \cellcolor[rgb]{ .502,  .647,  .827}0.0099 \\
    \hline
    96    & \cellcolor[rgb]{ .361,  .545,  .776}0.0007 & \cellcolor[rgb]{ .357,  .545,  .776}0.0005 & \cellcolor[rgb]{ .353,  .541,  .776}0.0003 & \cellcolor[rgb]{ .357,  .545,  .776}0.0006 & \cellcolor[rgb]{ .357,  .541,  .776}0.0004 & \cellcolor[rgb]{ .353,  .541,  .776}0.0001 & \cellcolor[rgb]{ .353,  .541,  .776}0.0003 \\
    \hline
\end{tabular}
\bigskip
\subcaption*{Panel B: Average over 3 realizations}
\begin{tabular}{|l|l|l|l|l|l|l|l|l|}
\hline
\backslashbox{$N_b$}{$N_{\Pi}$}             & 2      & 5      & 10     & 20     & 40     & 80     & 160    \\ \hline
2     & \cellcolor[rgb]{ .976,  .502,  .51}0.4738 & \cellcolor[rgb]{ .973,  .412,  .42}0.5502 & \cellcolor[rgb]{ .976,  .42,  .427}0.5458 & \cellcolor[rgb]{ .988,  .875,  .882}0.1493 & \cellcolor[rgb]{ .988,  .894,  .906}0.1310 & \cellcolor[rgb]{ .984,  .8,  .808}0.2138 & \cellcolor[rgb]{ .984,  .839,  .847}0.1804 \\
    \hline
    3     & \cellcolor[rgb]{ .98,  .651,  .663}0.3415 & \cellcolor[rgb]{ .988,  .902,  .914}0.1240 & \cellcolor[rgb]{ .98,  .635,  .643}0.3578 & \cellcolor[rgb]{ .98,  .694,  .706}0.3044 & \cellcolor[rgb]{ .988,  .863,  .875}0.1573 & \cellcolor[rgb]{ .988,  .906,  .918}0.1192 & \cellcolor[rgb]{ .988,  .867,  .875}0.1562 \\
    \hline
    6     & \cellcolor[rgb]{ .976,  .467,  .475}0.5052 & \cellcolor[rgb]{ .988,  .882,  .894}0.1423 & \cellcolor[rgb]{ .988,  .961,  .973}0.0733 & \cellcolor[rgb]{ .988,  .851,  .863}0.1688 & \cellcolor[rgb]{ .984,  .808,  .82}0.2066 & \cellcolor[rgb]{ .988,  .91,  .922}0.1158 & \cellcolor[rgb]{ .357,  .541,  .776}0.0005 \\
    \hline
    12    & \cellcolor[rgb]{ .988,  .98,  .992}0.0551 & \cellcolor[rgb]{ .988,  .98,  .992}0.0566 & \cellcolor[rgb]{ .408,  .58,  .796}0.0044 & \cellcolor[rgb]{ .467,  .62,  .816}0.0087 & \cellcolor[rgb]{ .988,  .988,  1}0.0472 & \cellcolor[rgb]{ .4,  .573,  .792}0.0037 & \cellcolor[rgb]{ .988,  .918,  .929}0.1108 \\
    \hline
    24    & \cellcolor[rgb]{ .988,  .988,  1}0.0498 & \cellcolor[rgb]{ .529,  .667,  .839}0.0135 & \cellcolor[rgb]{ .408,  .58,  .796}0.0045 & \cellcolor[rgb]{ .357,  .541,  .776}0.0005 & \cellcolor[rgb]{ .361,  .545,  .776}0.0009 & \cellcolor[rgb]{ .373,  .553,  .78}0.0016 & \cellcolor[rgb]{ .369,  .549,  .78}0.0013 \\
    \hline
    48    & \cellcolor[rgb]{ .396,  .573,  .792}0.0034 & \cellcolor[rgb]{ .373,  .553,  .78}0.0016 & \cellcolor[rgb]{ .682,  .773,  .89}0.0246 & \cellcolor[rgb]{ .365,  .549,  .78}0.0012 & \cellcolor[rgb]{ .384,  .565,  .788}0.0027 & \cellcolor[rgb]{ .353,  .541,  .776}0.0002 & \cellcolor[rgb]{ .369,  .553,  .78}0.0015 \\
    \hline
    96    & \cellcolor[rgb]{ .357,  .541,  .776}0.0004 & \cellcolor[rgb]{ .353,  .541,  .776}0.0002 & \cellcolor[rgb]{ .357,  .545,  .776}0.0006 & \cellcolor[rgb]{ .357,  .541,  .776}0.0005 & \cellcolor[rgb]{ .357,  .545,  .776}0.0006 & \cellcolor[rgb]{ .353,  .541,  .776}0.0001 & \cellcolor[rgb]{ .369,  .553,  .78}0.0014 \\
    \hline
\end{tabular}
\bigskip
\subcaption*{Panel C: Average over 27 realizations}
\begin{tabular}{|l|l|l|l|l|l|l|l|l|}
\hline
\backslashbox{$N_b$}{$N_{\Pi}$}             & 2      & 5      & 10     & 20     & 40     & 80     & 160    \\ \hline
2     & \cellcolor[rgb]{ .973,  .412,  .42}0.5320 & \cellcolor[rgb]{ .976,  .471,  .482}0.4828 & \cellcolor[rgb]{ .976,  .514,  .522}0.4473 & \cellcolor[rgb]{ .984,  .753,  .765}0.2458 & \cellcolor[rgb]{ .984,  .776,  .788}0.2253 & \cellcolor[rgb]{ .984,  .729,  .741}0.2660 & \cellcolor[rgb]{ .98,  .678,  .686}0.3100 \\
    \hline
    3     & \cellcolor[rgb]{ .976,  .494,  .502}0.4660 & \cellcolor[rgb]{ .976,  .545,  .553}0.4211 & \cellcolor[rgb]{ .98,  .631,  .639}0.3503 & \cellcolor[rgb]{ .984,  .835,  .847}0.1781 & \cellcolor[rgb]{ .988,  .851,  .863}0.1644 & \cellcolor[rgb]{ .988,  .847,  .859}0.1661 & \cellcolor[rgb]{ .984,  .824,  .835}0.1870 \\
    \hline
    6     & \cellcolor[rgb]{ .98,  .569,  .576}0.4031 & \cellcolor[rgb]{ .984,  .835,  .847}0.1765 & \cellcolor[rgb]{ .988,  .882,  .894}0.1380 & \cellcolor[rgb]{ .988,  .976,  .988}0.0583 & \cellcolor[rgb]{ .988,  .941,  .953}0.0871 & \cellcolor[rgb]{ .988,  .945,  .957}0.0852 & \cellcolor[rgb]{ .988,  .886,  .898}0.1340 \\
    \hline
    12    & \cellcolor[rgb]{ .988,  .918,  .929}0.1086 & \cellcolor[rgb]{ .988,  .969,  .98}0.0641 & \cellcolor[rgb]{ .988,  .988,  1}0.0471 & \cellcolor[rgb]{ .427,  .592,  .8}0.0059 & \cellcolor[rgb]{ .576,  .698,  .855}0.0169 & \cellcolor[rgb]{ .569,  .694,  .851}0.0164 & \cellcolor[rgb]{ .612,  .722,  .867}0.0195 \\
    \hline
    24    & \cellcolor[rgb]{ .988,  .984,  .996}0.0525 & \cellcolor[rgb]{ .647,  .745,  .878}0.0220 & \cellcolor[rgb]{ .482,  .631,  .82}0.0101 & \cellcolor[rgb]{ .541,  .671,  .839}0.0142 & \cellcolor[rgb]{ .416,  .584,  .796}0.0052 & \cellcolor[rgb]{ .369,  .553,  .78}0.0016 & \cellcolor[rgb]{ .365,  .549,  .78}0.0012 \\
    \hline
    48    & \cellcolor[rgb]{ .431,  .596,  .804}0.0061 & \cellcolor[rgb]{ .388,  .565,  .788}0.0030 & \cellcolor[rgb]{ .365,  .549,  .78}0.0013 & \cellcolor[rgb]{ .365,  .549,  .78}0.0013 & \cellcolor[rgb]{ .365,  .549,  .78}0.0012 & \cellcolor[rgb]{ .353,  .541,  .776}0.0003 & \cellcolor[rgb]{ .361,  .545,  .776}0.0009 \\
    \hline
    96    & \cellcolor[rgb]{ .376,  .557,  .784}0.0022 & \cellcolor[rgb]{ .353,  .541,  .776}0.0005 & \cellcolor[rgb]{ .357,  .545,  .776}0.0008 & \cellcolor[rgb]{ .357,  .541,  .776}0.0006 & \cellcolor[rgb]{ .353,  .541,  .776}0.0004 & \cellcolor[rgb]{ .357,  .541,  .776}0.0007 & \cellcolor[rgb]{ .353,  .541,  .776}0.0005 \\
    \hline
\end{tabular}
\label{table:biot_forward_nu_nc}
\end{table}



\rev{In Table \ref{table:biot_forward_nl_nn}, we present a sensitivity analysis of $N_{hl}$ and $N_{n}$ with a fixed size of the training set; $N_b=96$ and $N_{\Pi}=160$. Once again, the observed trend becomes more apparent by averaging over many training realizations. We can now identify an extremum in the explored space of the hyperparameters corresponding to have $N_{hl} = 6$ and $N_{n} = 20$; see Panel C of Table \ref{table:biot_forward_nl_nn}. We then apply all of the 27 PINN networks trained with that choice of hyperparameters to the test set. We obtain the average $\mathcal{L}^2$ error of the test set to be 9.05 $\times$ $10^{-5}$, which is comparable to that of the validation set, 9.08 $\times$ $10^{-5}$.}



\begin{table}[!ht]
\centering
\caption{
{\bf Biot's equations - forward modeling: PINN performance as function of hyperparameters.} Sum of relative ${\mathcal{L}^2}$ errors between the exact and predicted values of ${u}$, ${v}$, and ${p}$ for the \rev{validation} set. The table shows the dependency on the different number of hidden layers, $N_{hl}$, and different number of neurons per layer, $N_{n}$. Here, the total number of training and collocation points is fixed to $N_b$ = 96 and $N_{\Pi}$ = 160, respectively.}
\subcaption*{Panel A: Average over 1 realization}
\begin{tabular}{|l|l|l|l|l|l|l|}
\hline


\backslashbox{$N_{hl}$}{$N_n$}   & 2      & 5      & 10     & 20     & 40     & 80     \\ \hline
 2     & \cellcolor[rgb]{ .98,  .643,  .651}0.56506 & \cellcolor[rgb]{ .988,  .98,  .992}0.01617 & \cellcolor[rgb]{ .988,  .984,  .996}0.00782 & \cellcolor[rgb]{ .831,  .878,  .945}0.00045 & \cellcolor[rgb]{ .502,  .647,  .827}0.00021 & \cellcolor[rgb]{ .988,  .988,  1}0.00115 \\
    \hline
    4     & \cellcolor[rgb]{ .988,  .918,  .929}0.12080 & \cellcolor[rgb]{ .988,  .988,  1}0.00316 & \cellcolor[rgb]{ .392,  .569,  .788}0.00013 & \cellcolor[rgb]{ .482,  .631,  .82}0.00019 & \cellcolor[rgb]{ .49,  .639,  .824}0.00020 & \cellcolor[rgb]{ .431,  .596,  .804}0.00016 \\
    \hline
    6     & \cellcolor[rgb]{ .984,  .71,  .718}0.45949 & \cellcolor[rgb]{ .988,  .976,  .988}0.02069 & \cellcolor[rgb]{ .929,  .945,  .976}0.00052 & \cellcolor[rgb]{ .988,  .988,  1}0.00060 & \cellcolor[rgb]{ .353,  .541,  .776}0.00010 & \cellcolor[rgb]{ .494,  .639,  .824}0.00020 \\
    \hline
    8     & \cellcolor[rgb]{ .988,  .902,  .914}0.14333 & \cellcolor[rgb]{ .988,  .984,  .996}0.00971 & \cellcolor[rgb]{ .973,  .412,  .42}0.93757 & \cellcolor[rgb]{ .871,  .906,  .957}0.00048 & \cellcolor[rgb]{ .682,  .773,  .89}0.00034 & \cellcolor[rgb]{ .42,  .588,  .8}0.00015 \\
    \hline
    16    & \cellcolor[rgb]{ .988,  .906,  .914}0.14052 & \cellcolor[rgb]{ .988,  .902,  .914}0.14110 & \cellcolor[rgb]{ .788,  .847,  .929}0.00041 & \cellcolor[rgb]{ .494,  .639,  .824}0.00020 & \cellcolor[rgb]{ .494,  .639,  .824}0.00020 & \cellcolor[rgb]{ .482,  .631,  .82}0.00019 \\
    \hline
    32    & \cellcolor[rgb]{ .98,  .62,  .627}0.60485 & \cellcolor[rgb]{ .988,  .973,  .984}0.03188 & \cellcolor[rgb]{ .988,  .984,  .996}0.00866 & \cellcolor[rgb]{ .988,  .984,  .996}0.00972 & \cellcolor[rgb]{ .6,  .714,  .863}0.00028 & \cellcolor[rgb]{ .757,  .824,  .918}0.00039 \\
    \hline
\end{tabular}
\bigskip
\subcaption*{Panel B: Average over 3 realizations}
\begin{tabular}{|l|l|l|l|l|l|l|}
\hline
\backslashbox{$N_{hl}$}{$N_n$}   & 2      & 5      & 10     & 20     & 40     & 80     \\ \hline
2     & \cellcolor[rgb]{ .98,  .702,  .71}0.30761 & \cellcolor[rgb]{ .988,  .98,  .992}0.01074 & \cellcolor[rgb]{ .988,  .988,  1}0.00346 & \cellcolor[rgb]{ .635,  .741,  .875}0.00019 & \cellcolor[rgb]{ .765,  .831,  .922}0.00026 & \cellcolor[rgb]{ .643,  .745,  .878}0.00020 \\
    \hline
    4     & \cellcolor[rgb]{ .988,  .863,  .875}0.13684 & \cellcolor[rgb]{ .988,  .98,  .992}0.01101 & \cellcolor[rgb]{ .894,  .922,  .965}0.00032 & \cellcolor[rgb]{ .353,  .541,  .776}0.00006 & \cellcolor[rgb]{ .38,  .561,  .784}0.00007 & \cellcolor[rgb]{ .396,  .573,  .792}0.00008 \\
    \hline
    6     & \cellcolor[rgb]{ .988,  .855,  .867}0.14338 & \cellcolor[rgb]{ .988,  .949,  .961}0.04415 & \cellcolor[rgb]{ .612,  .725,  .867}0.00018 & \cellcolor[rgb]{ .471,  .624,  .816}0.00012 & \cellcolor[rgb]{ .42,  .588,  .8}0.00009 & \cellcolor[rgb]{ .376,  .557,  .784}0.00007 \\
    \hline
    8     & \cellcolor[rgb]{ .976,  .541,  .553}0.47702 & \cellcolor[rgb]{ .988,  .976,  .988}0.01634 & \cellcolor[rgb]{ .651,  .749,  .878}0.00020 & \cellcolor[rgb]{ .431,  .596,  .804}0.00010 & \cellcolor[rgb]{ .396,  .573,  .792}0.00008 & \cellcolor[rgb]{ .376,  .557,  .784}0.00007 \\
    \hline
    16    & \cellcolor[rgb]{ .98,  .678,  .69}0.33020 & \cellcolor[rgb]{ .988,  .863,  .875}0.13525 & \cellcolor[rgb]{ .988,  .988,  1}0.00081 & \cellcolor[rgb]{ .988,  .988,  1}0.00041 & \cellcolor[rgb]{ .49,  .635,  .824}0.00013 & \cellcolor[rgb]{ .502,  .647,  .827}0.00013 \\
    \hline
    32    & \cellcolor[rgb]{ .976,  .549,  .561}0.46854 & \cellcolor[rgb]{ .973,  .412,  .42}0.61381 & \cellcolor[rgb]{ .98,  .627,  .639}0.38437 & \cellcolor[rgb]{ .988,  .922,  .929}0.07503 & \cellcolor[rgb]{ .988,  .988,  1}0.00087 & \cellcolor[rgb]{ .443,  .604,  .808}0.00010 \\
    \hline
\end{tabular}
\bigskip
\subcaption*{Panel C: Average over 27 realizations}
\begin{tabular}{|l|l|l|l|l|l|l|}
\hline
\backslashbox{$N_{hl}$}{$N_n$}   & 2      & 5      & 10     & 20     & 40     & 80     \\ \hline
2     & \cellcolor[rgb]{ .98,  .694,  .702}0.24203 & \cellcolor[rgb]{ .988,  .961,  .973}0.02952 & \cellcolor[rgb]{ .506,  .647,  .827}0.00180 & \cellcolor[rgb]{ .4,  .573,  .792}0.00062 & \cellcolor[rgb]{ .361,  .545,  .776}0.00021 & \cellcolor[rgb]{ .365,  .549,  .78}0.00023 \\
    \hline
    4     & \cellcolor[rgb]{ .976,  .498,  .506}0.39596 & \cellcolor[rgb]{ .988,  .922,  .933}0.06005 & \cellcolor[rgb]{ .675,  .769,  .89}0.00368 & \cellcolor[rgb]{ .376,  .557,  .784}0.00039 & \cellcolor[rgb]{ .353,  .541,  .776}0.00010 & \cellcolor[rgb]{ .357,  .541,  .776}0.00014 \\
    \hline
    6     & \cellcolor[rgb]{ .98,  .588,  .596}0.32610 & \cellcolor[rgb]{ .988,  .902,  .91}0.07829 & \cellcolor[rgb]{ .361,  .545,  .776}0.00018 & \cellcolor[rgb]{ .353,  .541,  .776}0.00009 & \cellcolor[rgb]{ .353,  .541,  .776}0.00010 & \cellcolor[rgb]{ .357,  .541,  .776}0.00014 \\
    \hline
    8     & \cellcolor[rgb]{ .976,  .467,  .475}0.42070 & \cellcolor[rgb]{ .984,  .843,  .855}0.12314 & \cellcolor[rgb]{ .376,  .557,  .784}0.00037 & \cellcolor[rgb]{ .353,  .541,  .776}0.00013 & \cellcolor[rgb]{ .353,  .541,  .776}0.00010 & \cellcolor[rgb]{ .353,  .541,  .776}0.00012 \\
    \hline
    16    & \cellcolor[rgb]{ .973,  .412,  .42}0.46364 & \cellcolor[rgb]{ .984,  .843,  .855}0.12222 & \cellcolor[rgb]{ .988,  .878,  .89}0.09613 & \cellcolor[rgb]{ .988,  .933,  .945}0.05190 & \cellcolor[rgb]{ .353,  .541,  .776}0.00012 & \cellcolor[rgb]{ .353,  .541,  .776}0.00011 \\
    \hline
    32    & \cellcolor[rgb]{ .976,  .478,  .486}0.41372 & \cellcolor[rgb]{ .976,  .514,  .522}0.38439 & \cellcolor[rgb]{ .976,  .514,  .522}0.38473 & \cellcolor[rgb]{ .988,  .949,  .961}0.03998 & \cellcolor[rgb]{ .988,  .933,  .945}0.05080 & \cellcolor[rgb]{ .988,  .984,  .996}0.01052 \\
    \hline
\end{tabular}
\label{table:biot_forward_nl_nn}
\end{table}

\begin{comment}
\rev{
\begin{table}[!ht]
\centering
\caption{
{\bf Biot's equations - forward mdoeling: Relative $\mathcal{L}^2$ error of the \rev{test} set between the predicted and the exact solution, $u(x, y, t)$, $v(x, y, t)$, and $p(x, y, t)$ for different number of hidden layers $N_{hl}$ and different number of neurons per layer $N_{n}$. \rev{Here, the total number of training and collocation points is fixed to $N_b$ = 96 and $N_{\Pi}$ = 160, respectively, and the result is an average over 27 realizations}}}
\begin{tabular}{|l|l|l|l|l|l|l|}
\hline
\backslashbox{$N_{hl}$}{$N_n$}   & 2      & 5      & 10     & 20     & 40     & 80     \\ \hline
    2     & \cellcolor[rgb]{ .98,  .694,  .702}2.42E-01 & \cellcolor[rgb]{ .988,  .961,  .973}2.95E-02 & \cellcolor[rgb]{ .506,  .647,  .827}1.80E-03 & \cellcolor[rgb]{ .4,  .573,  .792}6.19E-04 & \cellcolor[rgb]{ .361,  .545,  .776}2.07E-04 & \cellcolor[rgb]{ .365,  .549,  .78}2.29E-04 \\
    \hline
    4     & \cellcolor[rgb]{ .976,  .498,  .51}3.96E-01 & \cellcolor[rgb]{ .988,  .922,  .933}6.01E-02 & \cellcolor[rgb]{ .675,  .769,  .89}3.68E-03 & \cellcolor[rgb]{ .376,  .557,  .784}3.93E-04 & \cellcolor[rgb]{ .353,  .541,  .776}9.82E-05 & \cellcolor[rgb]{ .357,  .541,  .776}1.36E-04 \\
    \hline
    6     & \cellcolor[rgb]{ .98,  .588,  .596}3.26E-01 & \cellcolor[rgb]{ .988,  .902,  .91}7.82E-02 & \cellcolor[rgb]{ .361,  .545,  .776}1.78E-04 & \cellcolor[rgb]{ .353,  .541,  .776}9.05E-05 & \cellcolor[rgb]{ .353,  .541,  .776}9.51E-05 & \cellcolor[rgb]{ .357,  .541,  .776}1.38E-04 \\
    \hline
    8     & \cellcolor[rgb]{ .976,  .471,  .478}4.20E-01 & \cellcolor[rgb]{ .984,  .843,  .855}1.23E-01 & \cellcolor[rgb]{ .376,  .557,  .784}3.71E-04 & \cellcolor[rgb]{ .353,  .541,  .776}1.28E-04 & \cellcolor[rgb]{ .353,  .541,  .776}1.02E-04 & \cellcolor[rgb]{ .353,  .541,  .776}1.24E-04 \\
    \hline
    16    & \cellcolor[rgb]{ .973,  .412,  .42}4.64E-01 & \cellcolor[rgb]{ .984,  .843,  .855}1.22E-01 & \cellcolor[rgb]{ .988,  .878,  .89}9.62E-02 & \cellcolor[rgb]{ .988,  .933,  .945}5.19E-02 & \cellcolor[rgb]{ .353,  .541,  .776}1.16E-04 & \cellcolor[rgb]{ .353,  .541,  .776}1.15E-04 \\
    \hline
    32    & \cellcolor[rgb]{ .976,  .478,  .486}4.14E-01 & \cellcolor[rgb]{ .976,  .514,  .522}3.84E-01 & \cellcolor[rgb]{ .976,  .514,  .522}3.85E-01 & \cellcolor[rgb]{ .988,  .949,  .961}4.00E-02 & \cellcolor[rgb]{ .988,  .933,  .945}5.07E-02 & \cellcolor[rgb]{ .988,  .984,  .996}1.05E-02 \\
    \hline
\end{tabular}
\label{table:biot_forward_nl_nn_test}
\end{table}
}
\end{comment}

\rev{
In Fig \ref{fig:biot_loss_list}, we show the behavior of the different loss terms as function of the training iterations, when using $N_{hl} = 6$, $N_{n} = 20$, $N_b = 96$, and $N_{\Pi} = 160$. Similar to the diffusivity equation, we observe that $M S E_{b}$ is generally higher than $M S E_{\Pi_{\bm{u}}}$ and $M S E_{\Pi_{p}}$. Furthermore, $M S E_{\Pi_{\bm{u}}}$ and $M S E_{\Pi_{p}}$ are comparable. Unlike the diffusivity equation case (see Fig \ref{fig:darcy_loss_list}), we observe minor oscillations of $M S E$, $M S E_{b}$, $M S E_{\Pi_{\bm{u}}}$, and $M S E_{\Pi_{p}}$ during convergence.
}

\begin{figure}[!ht]
   \centering
        \includegraphics[width=8.0cm, height=8.0cm]{pictures/loss_list_biot.pdf}
   \caption{
{\bf Biot's equations - forward modeling: mean square training error plot.} This model uses $N_{hl} = 6$, $N_{n} = 20$, $N_b = 96$, and $N_{\Pi} = 160$. $M S E$, $M S E_{b}$, $M S E_{\Pi_{\bm{u}}}$, and $M S E_{\Pi_{p}}$ are calculated using Eqs. \ref{eq:loss_biot_forward}, \ref{eq:loss_p_biot_forward}, \ref{eq:loss_pif_u_biot_forward}, and \ref{eq:loss_pif_p_biot_forward}, respectively.}
   \label{fig:biot_loss_list}
\end{figure}

\subsubsection*{Inverse model of Biot's equations}

We rewrite Eq~(\ref{eq:linear_balance_pinn}) to the
parametrized form as:

\begin{equation}\label{eq:biot_mom_inv}
\Pi_{\bm{u}} = \nabla \cdot \left [ 2 \theta_1 \bm{\varepsilon}(\bm{u})+\theta_2  \boldsymbol{u} \bm{I} \right ] -\theta_3 \nabla \cdot p \bm{I} - \bm{f} \text { in } \Omega \times \mathbb{T},
\end{equation}

\noindent
and Eq~(\ref{eq:mass_balance_pinn}) as:

\begin{equation} \label{eq:biot_mass_inv}
\Pi_p = \theta_4 \frac{\partial p}{\partial t}+\theta_3 \frac{\partial \nabla \cdot \boldsymbol{u}}{\partial t}- \theta_5 \nabla \cdot e^{\varepsilon_{v}}(\nabla p-\rho \mathbf{g}) - g \text { in } \Omega \times \mathbb{T},
\end{equation}

\noindent
where

\begin{equation} \label{eq:biot_theta}
\begin{aligned}
\theta_1 = \mu_l \text{,} \quad \theta_2 = \lambda_l \text{,} \quad
\theta_3 = \alpha \text{,} \quad
\theta_4 = \phi c_{f}+\frac{\alpha-\phi}{K_{s}}  \text{, } \text{and} \quad 
\theta_5 = {\kappa}_0.
\end{aligned}
\end{equation}

\noindent
\rev{The $\Pi$ terms now have five additional unknown parameters, $\theta_1$, $\theta_2$, $\theta_3$, $\theta_4$, and $\theta_5$ that along with the weights and biases of the neural network are adjusted during the training of the network.} Once again we apply Eq. (\ref{eq:loss_gen}) and obtain:

\begin{equation}\label{eq:loss_biot_inv}
M S E=M S E_{tr}+M S E_{\Pi_{\bm{u}}} + +M S E_{\Pi_{p}}
\end{equation}

\noindent
where 

\begin{equation}\label{eq:loss_p_biot_inv}
\begin{aligned}
M S E_{tr}&=\frac{1}{N_{tr}} \sum_{i=1}^{N_{tr}}(\left| u\left(x_{tr}^{i}, y_{tr}^{i}, t_{tr}^{i} \right)-u^{i}\right|^{2} + \left|v\left(x_{tr}^{i}, y_{tr}^{i}, t_{tr}^{i} \right)-v^{i}\right|^{2} \\ &+\left|p\left(x_{tr}^{i}, y_{tr}^{i}, t_{tr}^{i} \right)-p^{i}\right|^{2}),
\end{aligned}
\end{equation}

\noindent

\begin{equation}\label{eq:loss_pif_u_biot_inv}
M S E_{\Pi_{\bm{u}}}=\frac{1}{N_{\Pi_{\bm{u}}}} \sum_{i=1}^{N_{\Pi_{\bm{u}}}}\left|\Pi_{\bm{u}}\left(x_{\Pi_{\bm{u}}}^{i},y_{\Pi_{\bm{u}}}^{i}, t_{\Pi_{\bm{u}}}^{i}\right)\right|^{2},
\end{equation}

\noindent
and

\begin{equation}\label{eq:loss_pif_p_biot_inv}
M S E_{\Pi_{p}}=\frac{1}{N_{\Pi_{p}}} \sum_{i=1}^{N_{\Pi_{p}}}\left|\Pi_{p}\left(x_{\Pi_{p}}^{i},y_{\Pi_{p}}^{i}, t_{\Pi_{p}}^{i}\right)\right|^{2},
\end{equation}

\noindent
where $\left\{x_{tr}^{i},y_{tr}^{i},t_{tr}^{i}, u^{i},v^{i},p^{i}\right\}_{i=1}^{N_{tr}}$ refer to the set of training data. In contrast to the forward model, we \rev{can apply the training points as collocation points when calculating the terms given by Eqs (\ref{eq:biot_mom_inv}) and (\ref{eq:biot_mass_inv})}. To recap, the $N_{tr}$ denotes the number of training data. Similar to the forward model, all physical constants are set to one, i.e., $\theta_1$, $\theta_2$, $\theta_3$, $\theta_4$, and $\theta_5$ are equal to one. \rev{Note that these $\theta$ are constant throughout the domain.}


\rev{
We use the optimal hyperparameters, i.e., $N_{hl} = 6$ and $N_{n} = 20$ from the sensitivity analysis of the forward model (see Panel C of Table \ref{table:biot_forward_nl_nn}). In Table \ref{table:biot_inv_nn_nl}, we can observe that this choice also yields the least  $\mathcal{L}^2$ error with respect to both $p$, $u$, and $v$ and the lowest percentage error for the unknown physical parameters ($\theta_1$, $\theta_2$, $\theta_3$, $\theta_4$, and $\theta_5$). Moreover, we observe that the $\mathcal{L}^2$ error of the output space and percentage error of the unknown physical parameters are not much different with different combinations of the hyperparameters. 
}

\begin{table}[htbp]
  \centering
  \caption{{\bf Biot's equations - inverse modeling: PINN performance as function of hyperparameters.} Relative ${\mathcal{L}^2}$ error of ${p}$, ${u}$, and ${v}$ and percentage error of ${\theta_1}$, ${\theta_2}$, ${\theta_3}$, ${\theta_4}$, and ${\theta_5}$ for different number of hidden layers, ${N_{hl}}$, and different number of neurons per layer, ${N_{n}}$. The $N_{tr}$ is fixed at 250. Note that we pick the optimal hyperparameters, i.e., $N_{hl} = 6$ and $N_{n} = 20$ from the sensitivity analysis of the forward model. Results shown in this table are an average over 27 realizations.}
    \begin{tabular}{|l|l|l|l|l|l|l|}
\hline
&\backslashbox{$N_{hl}$}{$N_n$}    & 2        & 5        & 10\\
    \hline
    \multirow{3}[6]{*}{$p$} & 4     & \cellcolor[rgb]{ .906,  .929,  .969}0.00006 & \cellcolor[rgb]{ .427,  .592,  .8}0.00004 & \cellcolor[rgb]{ .98,  .584,  .592}0.00008 \\
\cline{2-5}          & 6     & \cellcolor[rgb]{ .973,  .412,  .42}0.00009 & \cellcolor[rgb]{ .353,  .541,  .776}0.00004 & \cellcolor[rgb]{ .988,  .988,  1}0.00006 \\
\cline{2-5}          & 8     & \cellcolor[rgb]{ .976,  .506,  .514}0.00008 & \cellcolor[rgb]{ .988,  .988,  1}0.00006 & \cellcolor[rgb]{ .871,  .902,  .957}0.00005 \\
    \hline
    \multirow{3}[6]{*}{$u$} & 4     & \cellcolor[rgb]{ .522,  .659,  .835}0.00019 & \cellcolor[rgb]{ .455,  .612,  .812}0.00016 & \cellcolor[rgb]{ .976,  .533,  .541}0.00047 \\
\cline{2-5}          & 6     & \cellcolor[rgb]{ .973,  .412,  .42}0.00051 & \cellcolor[rgb]{ .353,  .541,  .776}0.00013 & \cellcolor[rgb]{ .988,  .922,  .929}0.00036 \\
\cline{2-5}          & 8     & \cellcolor[rgb]{ .761,  .827,  .918}0.00027 & \cellcolor[rgb]{ .988,  .925,  .937}0.00036 & \cellcolor[rgb]{ .988,  .988,  1}0.00034 \\
    \hline
    \multirow{3}[6]{*}{$v$} & 4     & \cellcolor[rgb]{ .424,  .592,  .8}0.00042 & \cellcolor[rgb]{ .475,  .624,  .816}0.00045 & \cellcolor[rgb]{ .98,  .627,  .635}0.00098 \\
\cline{2-5}          & 6     & \cellcolor[rgb]{ .973,  .412,  .42}0.00114 & \cellcolor[rgb]{ .353,  .541,  .776}0.00038 & \cellcolor[rgb]{ .988,  .945,  .957}0.00074 \\
\cline{2-5}          & 8     & \cellcolor[rgb]{ .78,  .843,  .925}0.00060 & \cellcolor[rgb]{ .988,  .878,  .89}0.00079 & \cellcolor[rgb]{ .988,  .988,  1}0.00071 \\
    \hline
    \multirow{3}[6]{*}{$\theta_1$} & 4     & \cellcolor[rgb]{ .392,  .569,  .788}0.20180 & \cellcolor[rgb]{ .49,  .639,  .824}0.28665 & \cellcolor[rgb]{ .973,  .412,  .42}1.06223 \\
\cline{2-5}          & 6     & \cellcolor[rgb]{ .988,  .961,  .973}0.72512 & \cellcolor[rgb]{ .353,  .541,  .776}0.16844 & \cellcolor[rgb]{ .976,  .549,  .557}0.97916 \\
\cline{2-5}          & 8     & \cellcolor[rgb]{ .514,  .655,  .831}0.30534 & \cellcolor[rgb]{ .988,  .988,  1}0.70613 & \cellcolor[rgb]{ .976,  .525,  .537}0.99229 \\
    \hline
    \multirow{3}[6]{*}{$\theta_2$} & 4     & \cellcolor[rgb]{ .569,  .69,  .851}0.65507 & \cellcolor[rgb]{ .573,  .694,  .851}0.66560 & \cellcolor[rgb]{ .973,  .412,  .42}3.56789 \\
\cline{2-5}          & 6     & \cellcolor[rgb]{ .984,  .78,  .792}2.23747 & \cellcolor[rgb]{ .353,  .541,  .776}0.23053 & \cellcolor[rgb]{ .98,  .624,  .635}2.80360 \\
\cline{2-5}          & 8     & \cellcolor[rgb]{ .851,  .89,  .949}1.21152 & \cellcolor[rgb]{ .988,  .988,  1}1.47513 & \cellcolor[rgb]{ .984,  .722,  .729}2.45137 \\
    \hline
    \multirow{3}[6]{*}{$\theta_3$} & 4     & \cellcolor[rgb]{ .353,  .541,  .776}0.02968 & \cellcolor[rgb]{ .537,  .671,  .839}0.05016 & \cellcolor[rgb]{ .98,  .647,  .659}0.23851 \\
\cline{2-5}          & 6     & \cellcolor[rgb]{ .973,  .412,  .42}0.33292 & \cellcolor[rgb]{ .353,  .541,  .776}0.02935 & \cellcolor[rgb]{ .988,  .988,  1}0.10100 \\
\cline{2-5}          & 8     & \cellcolor[rgb]{ .98,  .678,  .69}0.22593 & \cellcolor[rgb]{ .988,  .949,  .961}0.11796 & \cellcolor[rgb]{ .745,  .816,  .914}0.07371 \\
    \hline
    \multirow{3}[6]{*}{$\theta_4$} & 4     & \cellcolor[rgb]{ .365,  .549,  .78}0.04834 & \cellcolor[rgb]{ .478,  .627,  .82}0.06162 & \cellcolor[rgb]{ .973,  .412,  .42}0.19135 \\
\cline{2-5}          & 6     & \cellcolor[rgb]{ .976,  .498,  .506}0.18104 & \cellcolor[rgb]{ .353,  .541,  .776}0.04666 & \cellcolor[rgb]{ .988,  .855,  .863}0.13785 \\
\cline{2-5}          & 8     & \cellcolor[rgb]{ .867,  .902,  .957}0.10711 & \cellcolor[rgb]{ .98,  .663,  .671}0.16124 & \cellcolor[rgb]{ .988,  .988,  1}0.12127 \\
    \hline
    \multirow{3}[6]{*}{$\theta_5$} & 4     & \cellcolor[rgb]{ .392,  .569,  .788}0.19710 & \cellcolor[rgb]{ .427,  .592,  .8}0.21622 & \cellcolor[rgb]{ .976,  .518,  .525}1.10700 \\
\cline{2-5}          & 6     & \cellcolor[rgb]{ .973,  .412,  .42}1.24043 & \cellcolor[rgb]{ .353,  .541,  .776}0.17633 & \cellcolor[rgb]{ .988,  .91,  .922}0.60595 \\
\cline{2-5}          & 8     & \cellcolor[rgb]{ .984,  .749,  .761}0.81199 & \cellcolor[rgb]{ .988,  .988,  1}0.50331 & \cellcolor[rgb]{ .776,  .839,  .925}0.39558 \\
    \hline
    \end{tabular}%
  \label{table:biot_inv_nn_nl}%
\end{table}%

\rev{
The performance of the PINN model as a function of the number of training examples $N_{tr}$ is depicted in Fig~\ref{fig:biot_inv_find_N_f}. We can observe that the stochastic variations in the estimated values of $\theta_1$, $\theta_2$, $\theta_3$, $\theta_4$ and $\theta_5$, in general, are reduced the more training examples we apply. Moreover, the relative $\mathcal{L}^2$ errors of $u$, $v$, and $p$ are always less than 0.01 \%. \rev{Using in the order of 1000 examples, the average estimation error of the physical parameters is in the order of 1 \%, but as with the diffusivity case, there is a large variation between the trained PINN models. Within the 27 realizations of the trained PINN models with $N_{tr} = 1000$ the percentage errors of $\theta_1$, $\theta_2$, $\theta_3$, $\theta_4$ and $\theta_5$ varied between $0.05$, $0.02$, $0.03$, $0.04$, and $0.03$ and  $2.51$, $6.14$, $4.75$, $8.12$, and $3.60$, respectively. Having 1000 training examples in actual cases, i.e., lab experiments or field observations, is realistic. Hence, also for the Biot's equations, we observe the feasibility of the PINN model to handle the inverse problem by estimating the unknown physical parameters using a reasonably sized data set. }
}





\begin{figure}[!ht]
   \centering
        \includegraphics[width=8.0cm, height=8.0cm]{pictures/biot_inv_find_N_f_lam.pdf}
        \includegraphics[width=8.0cm, height=8.0cm]{pictures/biot_inv_find_N_f_pri.pdf}
   \caption{
{\bf Biot's equations - inverse modeling: PINN performance as function of training set size.} This figure shows the estimated error dependency on the amount of training data, ${N_{tr}}$. (a) $\theta_1$, $\theta_2$, $\theta_3$, $\theta_4$, and $\theta_5$ and (b) $u$, $v$, and $p$. \rev{The error bars show mean and standard derivation ($\pm$ 1 SD) based on 27 realizations.} The reported error values of $\theta_1$, $\theta_2$, $\theta_3$, $\theta_4$, and $\theta_5$ are percentage errors while the relative $\mathcal{L}^2$ error is shown for $u$, $v$, and $p$. Note that there is no noise in this investigation.}
   \label{fig:biot_inv_find_N_f}
\end{figure}


\rev{Next, we perform a systematic study of the effect of noise in data, which is created utilizing Eq~(\ref{eq:darcy_noise}). The results are presented in Table \ref{table:biot_inv_noise} for $\theta_1$, $\theta_2$, $\theta_3$, $\theta_4$, and $\theta_5$. \rev{Note that these results are an average over ten realizations.} We can observe that the percentage error of $\theta_1$, $\theta_2$, $\theta_3$, $\theta_4$, and $\theta_5$ increase along with the noise level ($\epsilon$), but as $N_{tr}$ is increased, the impact of the noise is reduced.}  





\begin{table}[!ht]
\centering
\caption{
{\bf Biot's equations - inverse modeling: PINN performance as function of noise.} This figure shows the average percentage errors of ${\theta_1}$, ${\theta_2}$, ${\theta_3}$, ${\theta_4}$, and ${\theta_5}$ for different numbers of training data ${N_{tr}}$ corrupted by different noise levels (${\epsilon}$). Here, the neural network architecture is kept fixed to 6 layers and 20 neurons per layer. These results are an average over 10 realizations.}
\begin{tabular}{|l|l|l|l|l|l|}

\hline
\multicolumn{2}{|l|}{{\backslashbox{$N_{tr}$}{Noise ($\epsilon$)}} }                     & 0\% & 1\% & 5\% & 10\% \\ \hline
    \multirow{5}[10]{*}{$\theta_1$} & 1000  & \cellcolor[rgb]{ .435,  .596,  .804}0.91 & \cellcolor[rgb]{ .988,  .957,  .969}6.40 & \cellcolor[rgb]{ .988,  .937,  .949}7.15 & \cellcolor[rgb]{ .973,  .412,  .42}29.55 \\
\cline{2-6}          & 1500  & \cellcolor[rgb]{ .443,  .604,  .808}0.98 & \cellcolor[rgb]{ .945,  .957,  .984}4.64 & \cellcolor[rgb]{ .988,  .973,  .984}5.70 & \cellcolor[rgb]{ .984,  .757,  .769}14.86 \\
\cline{2-6}          & 2000  & \cellcolor[rgb]{ .4,  .576,  .792}0.67 & \cellcolor[rgb]{ .608,  .718,  .863}2.17 & \cellcolor[rgb]{ .988,  .984,  .996}5.26 & \cellcolor[rgb]{ .984,  .733,  .741}15.95 \\
\cline{2-6}          & 2500  & \cellcolor[rgb]{ .38,  .561,  .784}0.52 & \cellcolor[rgb]{ .839,  .882,  .945}3.87 & \cellcolor[rgb]{ .988,  .98,  .992}5.33 & \cellcolor[rgb]{ .984,  .82,  .831}12.24 \\
\cline{2-6}          & 5000  & \cellcolor[rgb]{ .353,  .541,  .776}0.30 & \cellcolor[rgb]{ .498,  .643,  .827}1.39 & \cellcolor[rgb]{ .765,  .831,  .922}3.32 & \cellcolor[rgb]{ .988,  .929,  .941}7.60 \\
    \hline
    \multirow{5}[10]{*}{$\theta_2$} & 1000  & \cellcolor[rgb]{ .518,  .659,  .835}1.48 & \cellcolor[rgb]{ .984,  .725,  .737}10.26 & \cellcolor[rgb]{ .98,  .639,  .647}12.30 & \cellcolor[rgb]{ .973,  .412,  .42}17.50 \\
\cline{2-6}          & 1500  & \cellcolor[rgb]{ .408,  .58,  .796}0.84 & \cellcolor[rgb]{ .988,  .871,  .882}6.95 & \cellcolor[rgb]{ .988,  .933,  .945}5.44 & \cellcolor[rgb]{ .976,  .541,  .549}14.53 \\
\cline{2-6}          & 2000  & \cellcolor[rgb]{ .498,  .643,  .827}1.37 & \cellcolor[rgb]{ .988,  .98,  .992}4.38 & \cellcolor[rgb]{ .988,  .878,  .89}6.73 & \cellcolor[rgb]{ .98,  .588,  .596}13.48 \\
\cline{2-6}          & 2500  & \cellcolor[rgb]{ .4,  .573,  .792}0.80 & \cellcolor[rgb]{ .949,  .961,  .984}3.95 & \cellcolor[rgb]{ .769,  .831,  .922}2.91 & \cellcolor[rgb]{ .988,  .898,  .91}6.31 \\
\cline{2-6}          & 5000  & \cellcolor[rgb]{ .353,  .541,  .776}0.51 & \cellcolor[rgb]{ .604,  .718,  .863}1.97 & \cellcolor[rgb]{ .765,  .831,  .922}2.90 & \cellcolor[rgb]{ .808,  .859,  .933}3.14 \\
    \hline
    \multirow{5}[10]{*}{$\theta_3$} & 1000  & \cellcolor[rgb]{ .804,  .859,  .933}1.19 & \cellcolor[rgb]{ .988,  .882,  .894}2.85 & \cellcolor[rgb]{ .98,  .686,  .694}5.06 & \cellcolor[rgb]{ .973,  .412,  .42}8.10 \\
\cline{2-6}          & 1500  & \cellcolor[rgb]{ .671,  .765,  .886}0.87 & \cellcolor[rgb]{ .859,  .898,  .953}1.32 & \cellcolor[rgb]{ .988,  .91,  .922}2.52 & \cellcolor[rgb]{ .984,  .71,  .718}4.80 \\
\cline{2-6}          & 2000  & \cellcolor[rgb]{ .553,  .682,  .847}0.58 & \cellcolor[rgb]{ .871,  .906,  .957}1.35 & \cellcolor[rgb]{ .988,  .914,  .925}2.49 & \cellcolor[rgb]{ .984,  .769,  .78}4.11 \\
\cline{2-6}          & 2500  & \cellcolor[rgb]{ .357,  .545,  .776}0.11 & \cellcolor[rgb]{ .467,  .62,  .816}0.38 & \cellcolor[rgb]{ .949,  .961,  .984}1.54 & \cellcolor[rgb]{ .984,  .769,  .78}4.12 \\
\cline{2-6}          & 5000  & \cellcolor[rgb]{ .353,  .541,  .776}0.10 & \cellcolor[rgb]{ .412,  .584,  .796}0.24 & \cellcolor[rgb]{ .988,  .98,  .992}1.72 & \cellcolor[rgb]{ .988,  .925,  .937}2.36 \\
    \hline
    \multirow{5}[10]{*}{$\theta_4$} & 1000  & \cellcolor[rgb]{ .4,  .573,  .792}0.46 & \cellcolor[rgb]{ .988,  .925,  .937}4.70 & \cellcolor[rgb]{ .984,  .718,  .729}8.35 & \cellcolor[rgb]{ .973,  .412,  .42}13.75 \\
\cline{2-6}          & 1500  & \cellcolor[rgb]{ .396,  .569,  .788}0.43 & \cellcolor[rgb]{ .882,  .914,  .961}3.01 & \cellcolor[rgb]{ .988,  .898,  .906}5.22 & \cellcolor[rgb]{ .984,  .71,  .718}8.53 \\
\cline{2-6}          & 2000  & \cellcolor[rgb]{ .357,  .541,  .776}0.23 & \cellcolor[rgb]{ .878,  .91,  .961}2.99 & \cellcolor[rgb]{ .988,  .961,  .969}4.12 & \cellcolor[rgb]{ .98,  .62,  .631}10.08 \\
\cline{2-6}          & 2500  & \cellcolor[rgb]{ .357,  .545,  .776}0.24 & \cellcolor[rgb]{ .435,  .6,  .804}0.64 & \cellcolor[rgb]{ .988,  .859,  .871}5.88 & \cellcolor[rgb]{ .98,  .663,  .675}9.34 \\
\cline{2-6}          & 5000  & \cellcolor[rgb]{ .353,  .541,  .776}0.20 & \cellcolor[rgb]{ .4,  .573,  .792}0.45 & \cellcolor[rgb]{ .62,  .729,  .871}1.62 & \cellcolor[rgb]{ .984,  .824,  .831}6.53 \\
    \hline
    \multirow{5}[10]{*}{$\theta_5$} & 1000  & \cellcolor[rgb]{ .722,  .8,  .906}3.46 & \cellcolor[rgb]{ .988,  .882,  .894}9.22 & \cellcolor[rgb]{ .98,  .671,  .682}16.15 & \cellcolor[rgb]{ .973,  .412,  .42}24.52 \\
\cline{2-6}          & 1500  & \cellcolor[rgb]{ .631,  .737,  .875}2.70 & \cellcolor[rgb]{ .973,  .976,  .992}5.64 & \cellcolor[rgb]{ .984,  .733,  .741}14.16 & \cellcolor[rgb]{ .98,  .561,  .569}19.75 \\
\cline{2-6}          & 2000  & \cellcolor[rgb]{ .373,  .557,  .784}0.44 & \cellcolor[rgb]{ .757,  .824,  .918}3.76 & \cellcolor[rgb]{ .984,  .812,  .824}11.56 & \cellcolor[rgb]{ .98,  .627,  .639}17.55 \\
\cline{2-6}          & 2500  & \cellcolor[rgb]{ .4,  .573,  .792}0.66 & \cellcolor[rgb]{ .58,  .698,  .855}2.23 & \cellcolor[rgb]{ .988,  .957,  .969}6.82 & \cellcolor[rgb]{ .988,  .949,  .961}7.09 \\
\cline{2-6}          & 5000  & \cellcolor[rgb]{ .353,  .541,  .776}0.24 & \cellcolor[rgb]{ .471,  .624,  .816}1.28 & \cellcolor[rgb]{ .812,  .863,  .937}4.26 & \cellcolor[rgb]{ .988,  .984,  .996}5.91 \\
    \hline
\end{tabular}
\begin{flushleft} As the noise ($\epsilon$) is increased, the error increases as expected.
\end{flushleft}
\label{table:biot_inv_noise}
\end{table}

\rev{All calculations were carried out using a XeonE5\_2660v3 processor with a single thread. As an example of the Biot's equations, the CPU time for training the neural networks using $N_{tr} = 10000$ and $15000$ with no noise are 128037 seconds and 186154 seconds, respectively. Note that the reported values are obtained from the model trained using the combined ADAM and L-BFGS. Using L-BFGS alone, the CPU times of model with $N_{tr} = 10000$ and $15000$ are 222681 and 294768 seconds, respectively.} 




\section*{Conclusion}

\rev{
This paper studies the application of physics-informed neural networks (PINN) for solving the nonlinear diffusivity and Biot's equations in the context of forward and inverse modelings. The following conclusions are drawn:
}

\begin{itemize}

\item \rev{PINN can be used to solve the forward modeling problem for the nonlinear diffusivity and Biot's equations, at least for the type of geometries considered in this paper. The displacement and pressure variables of our test sets could be predicted with an average $\mathcal{L}^2$ error of 9.05 $\times$ $10^{-5} \pm 3.1 \times 10^{-4}$ based on 27 realizations.  
}


\item \rev{For the inverse modeling cases, PINN can predict all of the unknown physical parameters with an average percentage error of around 1\%; however, the stochastic variations from one PINN implementation to the next is quite large. Using, for instance, 1000 training examples in the Biot's equation case, the percentage error of the estimated physical parameters over 27 PINN models could vary from $0.02$ to $8.12$. Increasing the number of training examples reduces this problem. Still, our results indicate it would be essential to do an average over PINN models with different random initialization of the weights and biases. This process may lead to the requirement of more processing power. This challenge might be even higher when applying PINN to more complex geometries and heterogeneous materials.
}




\item \rev{For the inverse modeling, PINN is tolerant to a noise level up to 5 \% (the estimation error of physical parameters is approximately less than 15 \%.). Again, this requires that one does an average over several PINN realizations. As expected, the result improves when the number of training examples is increased.
}

\item \rev{We have presented arguments on why the hyperparameters selection process for the forward case is likely to be applicable to the inverse case. For the cases considered here, this was confirmed experimentally. However, this should be explored in more detail by investigating the use of PINN for other types of nonlinear partial differential equations. 
}

\end{itemize}

\rev{
Finally, in terms of future work, the capability of the physics-informed neural networks should be tested in the case where the input data is incomplete, i.e., $\bm{u}$ and $p$ are not available at the same spatial and temporal coordinates. Besides, one could investigate the potential benefits of training networks using mini-batches. Moreover, smarter initialization of the weights and biases (based on transfer learning principles) could potentially be employed to increase the speed and accuracy of the training procedure \cite{yosinski2014transferable}.
}



\section*{Acknowledgments}
The research leading to these results has received funding from the Danish Hydrocarbon Research and Technology Centre under the Advanced Water Flooding program.

\section*{Author Contributions}
\begin{flushleft}
\textbf{Conceptualization:} Teeratorn Kadeethum, Thomas M Jørgensen, Hamidreza M Nick \par

\textbf{Formal Analysis:} Teeratorn Kadeethum \par

\textbf{Funding Acquisition:} Hamidreza M Nick \par

\textbf{Software:} Teeratorn Kadeethum \par

\textbf{Supervision:} Thomas M Jørgensen, Hamidreza M Nick \par

\textbf{Validation:} Thomas M Jørgensen, Hamidreza M Nick \par

\textbf{Writing – Original Draft Preparation:} Teeratorn Kadeethum, Thomas M Jørgensen \par

\textbf{Writing – Review and Editing:} Teeratorn Kadeethum, Thomas M Jørgensen, Hamidreza M Nick  \par
\end{flushleft}
\nolinenumbers


\bibliographystyle{unsrt}  
