A Deep Learning Galerkin Method for the Closed-Loop Geothermal System

There has been an arising trend of adopting deep learning methods to study partial differential equations (PDEs). This article is to propose a Deep Learning Galerkin Method (DGM) for the closed-loop geothermal system, which is a new coupled multi-physics PDEs and mainly consists of a framework of underground heat exchange pipelines to extract the geothermal heat from the geothermal reservoir. This method is a natural combination of Galerkin Method and machine learning with the solution approximated by a neural network instead of a linear combination of basis functions. We train the neural network by randomly sampling the spatiotemporal points and minimize loss function to satisfy the differential operators, initial condition, boundary and interface conditions. Moreover, the approximate ability of the neural network is proved by the convergence of the loss function and the convergence of the neural network to the exact solution in L^2 norm under certain conditions. Finally, some numerical examples are carried out to demonstrate the approximation ability of the neural networks intuitively.


Introduction
The PDEs are used to describe a variety of physical phenomena such as fluid dynamics, quantum mechanics, and elasticity. Many subsurface energy and environmental applications are controlled by multiphase fluid flow in porous media, including hydrocarbon recovery from oil and gas reservoirs [3,8], geological carbon dioxide storage [4,6,7], geothermal energy extraction [14,18,19]. The management of such processes has been increasingly relying on high fidelity physics-based numerical simulations. As one main form of renewable green energy, geothermal energy is environmental, independent of the weather or climate condition and widely distributed, hence it can significantly contribute to the sustainable energy application [13,5,9]. The closed-loop heat exchanger, one major technique to extract the geothermal energy, circulates working fluid down to the rock mass with high temperatures and then back to the surface through a continuous and closed-loop pipe [10,12,17,20]. There has been some studies shows that working fluid in the closed-loop geothermal system can be water, superheated steam, carbon dioxide, and so on [11,16,17].
In this article, we consider the closed-loop system for the geothermal reservoir with porous media flows [15,21,22,23]. This model combines the governing equations of the fluid flow in the pipe, the porous media flow in the geothermal reservoir, and the heat flow into a coupled system through the heat exchanging conditions and the no-fluid-communication conditions on the interface between the pipe and porous media. The proposed model is a multi-physics interface problem. Traditionally, the numerical simulators discretize the governing equations of multi-physics interface problem by finite difference/volume/element methods [26,27,28,29,30,31,32] resulting in large scale algebraic system of equations that may require significant computational resources and advanced numerical solvers and preconditioners [24,25].
During the past decade significant advances have been made in deep learning approaches [2] in both theory and application. With recent advances in CPU and GPU-based parallel processing, deep learning has achieved great success in numerous applications including image recognition, natural language processing, and automated driving to name a few. As a result, many mathematicians have introduced neural networks into the PDEs precisely because multilayer feedforward networks are proven to be universal approximators for the PDEs [33,34]. More specifically, once the network structure is determined, any order derivatives of the neural network can be obtained analytically. Coupled with the automatic differentiation technique, some studies have shown that neural networks can be applied to solve the PDEs [35][36][37][38][39][40][41][42][43][44][45][46][47][48][49][50][51]. Besides, many effective algorithms are proposed to solve some high-dimensional PDEs in [58][59][60][61][62][63][64][65]. Sirignano and Spiliopoulos [36] proposed the DGM to address the curse of dimensionality problem when solving high-dimensional PDEs. There are also some significant studies, Karpatne et al. [52] introduced Physics-guided Neural Networks, and successfully applied it in lake temperature modeling. Wang et al. [53] developed the Theory-guided Neural Network to reliably predict subsurface single-phase flow in porous media. Yan et al. [54] described a gradient-based deep neural network constrained by the physics related to multiphase flow in porous media.
Inspired by the framework of [36], we have successfully applied the DGM for solving the secondorder linear elliptic equations and stokes problems in [55][56][57]. In these studies, there is few about coupled multi-physics model and multi-physics interface problem. Facing the challenge, the focus of this work is to explore the use of DNN as a universal approximator to learn the porous media flow, fluid flow and heat flow in the closed-loop system. Our algorithm is meshfree and parallel for the equations and variables, the neural network is trained on batches of randomly sampled time and space points to satisfy the differential operators, initial condition, boundary and interface conditions. Moreover, we obtain the convergence of the loss function and the neural network, in addition, the performance of the method is demonstrated by some numerical experiments.
The remaining of this paper is organized as follows. In the next section, we discuss the new coupled multi-physics model and the DGM algorithm. Section 3 is devoted to illustrate the convergence of the loss function, while we present the proof of another theorem to guarantee the convergence of neural network's solution in Section 4. Finally, some numerical results are provided in Section 5.

Equations and preliminaries
As showed in Figure 1, the global domain Ω consists of two subdomains Ω f and Ω p , where Ω f , Ω p ∈ R 2 are open, regular, simply connected, and bounded by Lipschitz continuous boundaries of ∂Ω f \ I and ∂Ω p \ I.
Here Ω f Ω p = ∅,Ω f Ω p = I andΩ f Ω p =Ω.n f andn p are the unit normal vectors which satisfy the condition ofn f = −n p on the interface I. They point outward from the free flow region Ω f and porous media flow region Ω p respectively. The time frame is considered in [0, T ]. We assume that the fluid flow with heat transfer is governed by the combination of the Navier-Stokes equations with the heat equations in Ω f : 2) 3) Here u f , p f , θ f denotes the free fluid flow region velocity vector field, pressure, and temperature respectively. The unit vector ξ denote the direction of the gravitational acceleration. f f and Υ f and the external force terms. P r is the Prandtl number, Ra represents the Rayleigh number, k f > 0 is the thermal conductivity. Γ M , Γ E denotes the Dirichlet and Neumann boundary condition respectively in the pipe region boundaries where Γ M Γ E = ∂Ω f \ I.
We assume the porous media region Ω p is homogeneous and isotropic. The porous media flow with heat transfer can be governed by the following Darcy's law coupled with heat equations: 14) Here u p , p p and θ p denote the porous medium fluid flow region velocity vector field, pressure, and temperature respectively. The source term is denoted by Υ p , Da is the Darcy number which represents the relative effect of the permeability of the medium versus the cross-sectional area, k p is the thermal conductivity. Γ N , Γ Z denote the Dirichlet and Neumann boundary condition respectively in the porous media region boundaries where Γ N Γ Z = ∂Ω p \ I. • Continuity of temperature across the interface: • Continuity of heat flux across the interface:

18)
• No-communication and no-slip condition for the free flow on the interface: where τ denote the unit tangential vector along the interface.
• No-porous media flow passing through the interface: In order to derive the variational formulation for the model problem (2.1)-(2.16), we define the following spaces as follows: : Define a product space We also define the solenoidal spaces The trilinear form is defined as follows: In a similar manner, we define another two trilinear forms for any (t f , t p ) ∈ W I and (θ f , θ p ) ∈ W T : For simplicity of notations, let us denote where the (u f , p f , θ f ) is the analytic solution of the system (2.1)-(2.8) and the (u p , p p , θ p ) is the analytic solution of the system (2.9)-(2.16).

The DGM algorithm
We apply the multilayer feed forward networks, which consists of an input layer, one or more hidden layers and one output layer, to approximate solution of the equations (2.1)-(2.16). The approximation capabilities of neural network architectures have recently been investigated by many authors. Especially, in [60], they have ascertained that the standard multilayer feedforward networks with activation function ψ can approximate any continuous function well on arbitrary compact subsets X of Ω, whenever activation function ψ is continuous, bounded and nonconstant. For convenience, we will specifically contrive our results only for the instance where there is only one hidden layer, each hidden unit has the same activation function ψ i and one output unit. The simplified diagram of this neural network is shown in Figure 2. The set of all the functions are described by such a network with n hidden units for velocity, pressure and temperature of Ω f and Ω p , respectively. Especially, is a product space.
• C n * : the class of neural networks consisting of a single hidden layer with n hidden units and d random points x j , j = 1, 2, · · · , d and one output unit, for convenience, we still use this notation to represent multilayer neural networks; • β k i , σ k j , σ k ji are weights, c k,i are the thresholds of the network.
) are the parameters of the neural network.
Assuming that (U f , P f , Θ f ) and (U p , P p , Θ p ) are neural network solutions to the system (2.1)-(2.8) and (2.9)-(2.16), respectively. Then the loss function of the two systems are defined by , (2.23) and There are some notes for the loss function: where ν is a positive probability density on y ∈ Y; } are the parameters of the neural network, ν i (i = 1, 2, 3) are the positive probability densities of (Ω f , Ω p ), ∂Ω \ I and I, respectively; • J i (i = f, p, I) is defined to measure the difference between the output of the network with the known input. The whole deep learning process is to find appropriate parameters Λ such that J is as close to 0 as possible.
Step 2. At the random points s n = {(x n , t n ), (z n , τ n )}, we compute the square error G(Λ n ; s n ) by Step 3. Apply the stochastic gradient descent step at s n by Here, the learning rate η n decreases as n increases.
Step 4. Repeat step 3 until G(Λ n ; s n ) meets the convergence condition

Algorithm 1 Deep Learning Galerkin Method.
Input: training data, s n = {(x n , t n ), (z n , τ n )}, max number of iterations M; Output: Neural network predicted solutions, 1: Specify the training set s n = {(x n , t n ), (z n , τ n )} for the initial/boundary conditions, domain and equations; 3: Specify a loss function by summing the initial/boundary conditions, the residual of the PDE; 4: Train the neural network to find the best parameters Λ by minimizing the loss function J in (2.27); 5: Get the predicted solution on the domain; 6: Return: which has at most polynomial growth on u, uniformly with respect to x. This means that For the (U p , u p ), (Θ f , θ f ), (Θ p , θ p ), (P f , p f ), (P p , p p ), the above properties are utilized. Lemma 1 . If u f ∈ Y f , then there exists a positive constant C > 0 such that The weak formulation of the closed-loop geothermal system (2.1)-(2.16) has a unique weak solution (u f , u p ; Lemma 3 . By Theorem 3 of [33] we know that there exist functions ] d are uniformly 2-dense on compact sets of C ∞ (0, T ; L 2 (Ω f )) and C ∞ ((0, T ; L 2 (Ω p )), respectively. This means that for u f ∈ C ∞ (0, T ; L 2 (Ω f )), u p ∈ C ∞ (0, T ; L 2 (Ω p ) and 1 > 0, 2 > 0, there exist • Θ f and Θ p ∈ C n θ I (ω I ) are uniformly 2-dense on compact sets of C ∞ (0, T ; L 2 (Ω f ) × L 2 (Ω p )), the meaning is that for θ f , θ p ∈ C ∞ (0, T ; L 2 (Ω f ) × L 2 (Ω p )) and 3 > 0, 4 > 0, there exist • P f ∈ C n p f (ψ f ), P p ∈ C n pp (ψ p ) are uniformly 2-dense on compact sets of C 2 (0, T ; L 2 (Ω f )) and C 2 (0, T ; L 2 (Ω p )), respectively. This means that for p f ∈ C 2 (0, T ; L 2 (Ω f )), p p ∈ C 2 (0, T ; L 2 (Ω p )) and 5 > 0, 6 . (3.1) In the following we will estimate the terms on the right hand side of (3.1) one by one as follows. First, we have by applying the hypothesis (H1) and Lemma 3. Similarly, we can obtain

Add (3.3) and (3.2) to (3.1), we can get that
It is similar to J f , .

(3.5)
Then we estimate the right hand side of (3.5) using the same method, and Plugging the above estimates into (3.5), we obtain that There is J f + J p + J I ≤ K , i.e., J ≤ K by (3.4), (3.8) and (3.9).

Convergence of the neural network to the unique solution
This subsection contains convergence results of neural networks to the unique solution of the closed-loop geothermal system as n → ∞. During the whole parallel computing process, each variable can be calculated independently. For convenience of explanation, we take the same network structure for each variable. Inspired by the thought of the Galerkin Method, each neural network (U n f , U n p , P n f , P n p , Θ n f , Θ n p ) satisfies the system  16). Moreover, when the sequence (U n f , U n p , Θ n f , Θ n p ) is uniformly bounded and equicontinuous, the neural networks (U n f , U n p , Θ n f , Θ n p ) converges strongly to (u f , u p , θ f , θ p ) in L 2 (Ω f ) and L 2 (Ω p ), respectively. Furthermore, (U n f , U n p , Θ n f , Θ n p ) uniformly converges to (u f , u p , θ f , θ p ) in Ω f and Ω p , respectively. The weak formulations of the coupled system (2.1)-(2.16) is to find (U n f , U n p ; P n f , P n Using the Poincaré inequality and Young's inequality can lead to Let Θ n f = Θ n f (0, x), W n f = ∂Θ n f (0,x) ∂t in (3.21), we obtain that Combine the Poincaré, Cauchy-Schwarz and Young's inequalities, there is and integrating the equation (3.27) in time, we get that This means that {Θ n f }, {Θ n p } are uniformly bounded in L 2 (0, T ; H 1 (Ω)). Plugging the equations (3.17) and (3.18), Integrating the equation (3.29) in time and using the Gronwall inequality can lead to and 2P rRaξ . (3.36) Plugging the above estimates into (3.31), we obtain that .
)ds, (3.38) where This shows that  In the following we will estimate the terms on the right hand side of (3.39) one by one as follows. .

Numerical results
In this section, we present two numerical experiments to validate the proposed DGM for the closedloop geothermal system. In section 4.1, an analytic solution is provided to show the convergence of the neural network solutions and the impact of the hidden layers on the convergence. Section 4.2 is conducted to investigate a benchmark problem for thermal convection in a squared cavity.

Convergence and stability tests.
In this examples, we consider the problem (2.1)-(2.16) on a 2D spatial domain consisted of two subdomains. One is the free fluid flow region Ω f = [0, 1] × [1,2], and the other one is the porous media region Ω p = [0, 1] × [0, 1]. The selected exact solution for the model is given by (4.1) The initial condition, boundary condition, and forcing term of the model problem (2.1) to (2.16) are determined by the exact solutions of the model problem. Choose the parameters value a = 1.0, P r = 1.0, Ra = 1.0, k f = k p = 1.0, Da = 1.0 and γ = 10 5 . For convenience, we denote the L 2 error (errL 2 ) and relative L 2 error (errRL 2 ) as follows where U i and Θ i (i = 1, 2, · · · , m) are the neural network solutions of velocity and temperature, respectively. And u i and θ i (i = 1, 2, · · · , m) are the corresponding exact solutions of velocity and temperature.
We train five architectures of the neural networks on the training set with 900 sample points, these neural networks include one to five hidden layers respectively and each hidden layer has 16 units. Table 1 and Figure 3 demonstrate the approximation ability of the DGM for the closed-loop geothermal system (2.1)-(2.16) on 2D spatial domain. Moreover, the value of the loss function J(U ) decreases gradually while the number of hidden layers increases, which is consistent with the result of theorem 1. Besides, combining the L 2 error and relative L 2 error, we find that when the loss function goes to zero, the solution of the neural network converges to the exact solution, which is the content expressed in theorem 2. Obviously, Figure 4 and Figure 5 show the good agreement between the exact solution and neural networks results with the increasing of the hidden layers. These results suggest that the neural network is close to the exact solution with the increasing of hidden layers, which is consistent with the theory.   Left boundary is considered as heated wall with θ f = θ p = 1 and right boundary is cold with θ f = θ p = 0. On the interface between the two regions, the interface conditions 2.17-2.20, which are proposed for the model, are utilized.