Numerical Simulation of Heat Transport Problems in Porous Media Coupled with Water Flow Using the Network Method

: In the present work, a network model for the numerical resolution of the heat transport problem in porous media coupled with a water ﬂow is presented. Starting from the governing equations, both for 1D and 2D geometries, an equivalent electrical circuit is obtained after their spatial discretization, so that each term or addend of the differential equation is represented by an electrical device: voltage source, capacitor, resistor or voltage-controlled current source. To make this possible, it is necessary to establish an analogy between the real physical variables of the problem and the electrical ones, that is: temperature of the medium and voltage at the nodes of the network model. The resolution of the electrical circuit, by means of the different circuit resolution codes available today, provides, in a fast, simple and precise way, the exact solution of the temperature ﬁeld in the medium, which is usually represented by abaci with temperature-depth proﬁles. At the end of the article, a series of applications allow, on the one hand, to verify the precision of the numerical tool by comparison with existing analytical solutions and, on the other, to show the power of calculation and representation of solutions of the network models presented, both for problems in 1D domains, typical of scenarios with vertical ﬂows, and for 2D scenarios with regional ﬂow.


Introduction
The topic of heat transport coupled to water flow in porous media has aroused great interest among researchers in recent decades [1,2], finding in the literature a great variety of temperature-depth profiles, both for regional [3,4] and exclusively vertical flows [5,6] and for a wide variety of cases in the soil surface temperature condition [7][8][9][10][11]. The first advances in this discipline, whose processes are governed by coupled, non-linear PDEs derived from momentum, mass and energy conservation laws, date from the middle of the 20th century [12][13][14]. The patterns of solutions to this complex phenomenon are fundamental in many fields of engineering, such as the exploitation of geothermal energy extraction systems [15,16], saltwater intrusion [17], underground spread of pollutants [18] and petroleum engineering [19], among others.
In this research, the study of the heat transport problem partially coupled to a constant velocity field (either solved or previously known) corresponds to a linear problem, except for the surface boundary condition when it is a time-dependent function. The harmonic temperature conditions on the soil surface (which are the usual ones) give rise to a heat transport problem that is coupled with the movement of the fluid. These continuously changing conditions, although seasonally repetitive (day-night cycle; winter-summer cycle) determine equally changing temperature-depth profiles, with a certain inertia that is determined by the thermal parameters of the ground (conductivity and specific heat). Despite the changing conditions on the soil surface, it is possible to speak of stationary temperature-depth profiles, both on a daily and annual scale.
Within the great variety of numerical techniques that exist for solving these types of problems, it is worth highlighting the network simulation method [20], a finite vol-modifying the permeability of the soil, which would lead to changes in the fluid velocity in those scenarios governed by hydraulic potentials (in the scenarios studied here, it is the fluid velocity itself that governs the problem).
In this article, starting from the governing equations of the heat transport problem partially coupled to a known and constant velocity field (with vertical or horizontal direction), the design of a network model that will allow the resolution of this type of problem, both in 1D and 2D rectangular geometries, is presented. Subsequently, and after verifying the results obtained with our tool with the existing solutions in the reference literature [7][8][9][10], a series of selected applications (both for regional flow and exclusively vertical upflow) will be simulated in which different functions (dependent on time) will be assumed for the boundary condition of the soil surface temperature, obtaining the temperature-depth profiles once the stationary situation has been reached.

Mathematical Model
In this section, the equations and boundary conditions that govern different scenarios that address the problem of heat transport coupled to a water flow will be presented.
The physical model of the problem is that of a 2D rectangular geometry scenario that represents a fully saturated porous medium in which fluid flow and heat processes take place in such a way that the latter is coupled to the water flow (for vertical flow it will be enough to consider a one-dimensional domain). That is, the solution to the thermal problem depends on and is coupled to the solution of the mechanical problem.
The mechanical problem concerns the velocity field of water, which, for the purposes of this paper, will be uniform. That is, the water velocity will always be constant both in modulus, direction and sense, Equation (1). This consideration, although at first it seems quite strong, in certain groundwater flow scenarios it can be considered without committing major errors, as occurs in different real cases such as the discharge of an aquifer to the bottom of the sea or a lake (where the flow is exclusively vertical and it is usual for the velocity to be constant for long periods of time) or in regional (horizontal) flows [36], where it is also common to find large areas of aquifers with constant average velocities. On the other hand, the fact of knowing the velocity field (since it is imposed) makes the influence of viscosity (which is affected by temperature) negligible (fluid viscosity would affect velocities in those scenarios that were governed by hydraulic potential, but this is not the case with the problems discussed in this paper).
On the other hand, the thermal problem refers to the temperature field at each point of the porous medium, and is defined by the set of Equations (2)- (6) [1][2][3].
The mathematical model, according to the nomenclature of Figure 1, is as follows: v(x, y, t) = q 0 (1) (ρ e c e ) ∂T ∂t − k m ∇ 2 (T) + ρ w c w q 0 ∇(T) = 0 (2) T (x=0,y,t) = T left and T (x=L,y,t) = T right (4) T T (x, y=0, t) = T m + T sin sin 2π P t T (x,y,t=0) = T s (8) where Equation (2) represents the balance of heat transport in the medium (governing equation particularized to the case of constant water velocity). Equation (3) represents a constant temperature condition (or Dirichlet condition) at the bottom of the domain, also called the first-class condition, while Equation (4) represents the same type of condition but in the lateral boundaries of the domain. Equations (5)-(7) define the temperature condition on the ground surface in the form of a constant value, a step function or a sinusoidal, respectively. L and H are the dimensions of the rectangular domain. Finally, Equation (8) represents the initial condition defined by the initial temperature in the entire soil domain. Regarding the temperature conditions on the soil surface, Equations (5)- (7), it is important to note that, in principle, the only feasible boundary condition on this contour is the sinusoidal temperature distribution, Equation (7), equivalent to 24 h temperature variation for day and night. However, in this work, the conditions of constant (Equation (5)) and stepped (Equation (6)) temperature have also been included, for various reasons. The Regarding the temperature conditions on the soil surface, Equations (5)- (7), it is important to note that, in principle, the only feasible boundary condition on this contour is the sinusoidal temperature distribution, Equation (7), equivalent to 24 h temperature variation for day and night. However, in this work, the conditions of constant (Equation (5)) and stepped (Equation (6)) temperature have also been included, for various reasons. The first one is to illustrate that with the network method it is very easy to implement any boundary condition. In fact, in future research we plan to work with other surface temperature functions (even with tabulated data) whose implementation in the network model is just as simple. On the other hand, the constant temperature condition can sometimes serve as a simplification to consider an average daily temperature, although it could also refer to a monthly or annual average. For its part, the stepped temperature condition could be assimilated to a simplification similar to the previous one, but instead of having a single daily average temperature, there are two (an average maximum temperature for the day and an average minimum temperature for the night).

Mathematical Model for Vertical Flow
The exclusively vertical flow of water in porous aquifers is a fairly common situation to find in practice [37]. Thus, for example, this type of flow occurs in situations such as the discharge of a river into an aquifer [38], two aquifers separated by an aquitard [39] or in those cases in which an aquifer discharges directly onto the bottom of the sea or a lake [2].
The governing Equation (2) becomes, for one-dimensional vertical flow of an incompressible fluid through homogeneous porous media, as follows: On the other hand, depending on the boundary and initial conditions adopted, three variants of the problem with exclusive vertical flow are defined, governed by Equation (9).

Constant Surface Temperature
It is a scenario where the temperature of the ground surface (T g ) remains constant throughout the time, being different from the initial temperature of the soil mass (T s ). It is also necessary to set a temperature at the bottom of the domain as a boundary condition (T o ). Table 1 shows the values that these boundary and initial thermal conditions, Equations (3), (5) and (8), will take, for this case of vertical water flow and constant temperature on the ground surface, in the later applications section. Table 1. Boundary and initial thermal conditions for constant surface temperature and vertical flow scenario.

Parameter
Value Units In this scenario, the surface temperature varies alternately between two values (simulating a period with a colder surface temperature, and a warmer one). To do this, it is necessary to define a stepped temperature function, Figure 2, which will be imposed as a boundary condition on the ground surface.

Stepped Surface Temperature
In this scenario, the surface temperature varies alternately between two values ( ulating a period with a colder surface temperature, and a warmer one). To do this, necessary to define a stepped temperature function, Figure 2, which will be imposed boundary condition on the ground surface.  Table 2 summarizes the boundary and initial thermal conditions, Equations (3 and (8), for the case of vertical flow of water with stepped variation of the temperatur the ground surface that will be shown in the applications section.   Table 2 summarizes the boundary and initial thermal conditions, Equations (3), (6) and (8), for the case of vertical flow of water with stepped variation of the temperature on the ground surface that will be shown in the applications section. Table 2. Boundary and initial thermal conditions for stepped surface temperature and vertical flow scenario.

Parameter
Value Units This time, the surface temperature varies following a sinusoidal function (simulating what can be, in certain situations, the temperature fluctuation throughout a day). On this occasion, a sinusoidal function will be imposed as a boundary condition on the ground surface, Figure 3. Table 3 summarizes the boundary and initial thermal conditions, Equations (3), (7) and (8), for the case of vertical flow of water with sinusoidal temperature on the ground surface that will be numerically solved in the applications section. occasion, a sinusoidal function will be imposed as a boundary condition on the ground surface, Figure 3.  Table 3 summarizes the boundary and initial thermal conditions, Equations (3), (7) and (8), for the case of vertical flow of water with sinusoidal temperature on the ground surface that will be numerically solved in the applications section. Table 3. Boundary and initial thermal conditions for sinusoidal surface temperature and vertical flow scenario.

Parameter
Value

Mathematical Model for Regional Flow
In hydrogeology, a regional flow is commonly called that whose movement of water is exclusively horizontal [3,4]. This type of water movement is typically found in underground hydrology, especially in aquifer extensions that are relatively far from the discharge area [40].
As described in [41], the general equation that describes the simultaneous flow of fluid and heat in the ground [37] is: while for isotropic thermal media and two-dimensional rectangular geometry, Equation (10) reduces to: However, the vertical flow is zero, so Equation (11) is further simplified (note that, given the hypothesis of uniform water flow, vx has been replaced by qo,x): In the applications section, the sinusoidal surface temperature variant of the problem governed by Equation (12) will be presented.  Table 3. Boundary and initial thermal conditions for sinusoidal surface temperature and vertical flow scenario.

Mathematical Model for Regional Flow
In hydrogeology, a regional flow is commonly called that whose movement of water is exclusively horizontal [3,4]. This type of water movement is typically found in underground hydrology, especially in aquifer extensions that are relatively far from the discharge area [40].
As described in [41], the general equation that describes the simultaneous flow of fluid and heat in the ground [37] is: while for isotropic thermal media and two-dimensional rectangular geometry, Equation (10) reduces to: However, the vertical flow is zero, so Equation (11) is further simplified (note that, given the hypothesis of uniform water flow, v x has been replaced by q o,x ): In the applications section, the sinusoidal surface temperature variant of the problem governed by Equation (12) will be presented.

Sinusoidal Surface Temperature
Again, the surface temperature is described by a sinusoidal function, which will be imposed as a boundary condition on the ground surface. Table 4 summarizes the boundary and initial thermal conditions, Equations (3), (4), (7) and (8), for the application case of regional flow of water with sinusoidal temperature on the ground surface. Table 4. Boundary and initial thermal conditions for sinusoidal surface temperature and regional flow scenario.

Parameter
Value Units Free condition

Network Model Design
In order to design a network model (or circuit) equivalent to a certain process to be studied, it is necessary, first of all, to establish an analogy between the physical variables of the problem and the different electrical quantities. For the purposes of this research, we chose to establish an equivalence between the physical variable (real) temperature of the fluid (T) and the variable (analogous) electric potential (V). In this way, the spatially discretized equations of the mathematical model and the equations of the network model (corresponding to the analogous variables) for a volume element coincide.
The network model development procedure consists of reticulating the space into volume elements or elementary cells to later define the set of differential equations in finite differences (leaving time as a continuous variable) to be applied in these elementary cells. Once the corresponding equivalence between the real variables of the problem and the electrical variables has been established, each term of the differential equations in finite spatial differences can be identified with an electrical device (resistor, current source, voltage source, capacitor, etc.), which is arranged between the different nodes of the elementary cell and is balanced at the central node. It is important to note that, since the thermal (and where appropriate, hydraulic) characteristics of the material are included in the expressions of the electrical devices that make up each elementary cell (which represents a certain portion of soil), with the network method it would be possible to model and simulate non-homogenous porous media, since this technique allows designing circuits with cells of different properties, which would be implemented through simple programming routines. These types of media could include (i) layered stratified soils, (ii) changing properties with depth, (iii) soils with local anomalies and (iv) thermal and hydraulic properties imported from a database, among others.
Once the electrical circuit is designed, its resolution is carried out, which is quick and easy to obtain thanks to the different circuit resolution software available today [21], some of them free to access [22]. The precision of the results will be solely at the expense of the fineness of the chosen mesh: the higher the number of cells, the more accurate the model will be. In general, the network method requires very few cells to achieve highly accurate [33] results, with negligible relative errors, below 1%.

Elementary Cell for Vertical Flow
For the case of vertical flow, Equation (9) can be expressed in spatial finite differences (retaining time as a continuous variable) and, following the nomenclature of Figure 4, remains as follows: From the point of view of the network method, each addend of Equation (13) can be considered an electric current: which balance at the central node of the elementary cell, being: In the network method, those linear terms such as J R+∆ and J R−∆ can be implemented by resistive-type devices. For this reason, and following Ohm's law, I = V/R, the values of the resistors can be expressed as: These elements, as shown in Figure 4, are located between the central node and one of the ends, which is where the potential fall occurs (in the physical analogy, the temperature). For its part, the linear term J c,e , which includes a time derivative, can be implemented as a capacitor (C i ) of ρ e c e capacitance.
From the point of view of the network method, each addend of Equation (13) can be considered an electric current: which balance at the central node of the elementary cell, being: In the network method, those linear terms such as J ∆ and J ∆ can be implemented by resistive-type devices. For this reason, and following Ohm's law, I = V/R, the values of the resistors can be expressed as: These elements, as shown in Figure 4, are located between the central node and one of the ends, which is where the potential fall occurs (in the physical analogy, the temperature). For its part, the linear term J , , which includes a time derivative, can be implemented as a capacitor (C ) of ρ c capacitance.
On the other hand, J , represents a non-linear term, which must be implemented in the circuit as voltage-controlled current source. In this element, the current is specified directly, by means of the following expressions: where T ∆ and T ∆ are directly read at the extreme nodes of the corresponding cell i. This device, in the same way as the capacitor, must be implemented between the central node of the cell and the common ground node.  On the other hand, J G,w represents a non-linear term, which must be implemented in the circuit as voltage-controlled current source. In this element, the current is specified directly, by means of the following expressions: where T i+∆ and T i−∆ are directly read at the extreme nodes of the corresponding cell i. This device, in the same way as the capacitor, must be implemented between the central node of the cell and the common ground node.

Elementary Cell for Regional Flow
For the case of regional flow, Equation (12) is expressed in spatial finite differences (retaining time as a continuous variable and following the nomenclature of Figure 5) as follows: Again, each addend of Equation (18) can be considered an electric current: which balance at the central node of the elementary cell, being: As seen before, the linear terms J R,i+∆ and J R,i−∆ can be implemented by resistors, whose values will be specified by: and will be placed between the central node and the right and left extreme nodes ( Figure 5), respectively. For their part, the terms J R,j+∆ and J R,j−∆ are also modeled with both resistors, located between the central node and the upper and lower extreme nodes, respectively. Their values will come to be specified by: As for the case of vertical flow, the linear storage term J c,e is implemented by means of a capacitor (C i,j ), whose capacitance value is ρ e c e . Finally, the non-linear term J G,w will be implemented by the voltage-controlled current source G w,i,j , defined by: where T i+∆,j and T i−∆,j are directly read at the right and left extreme nodes of the corresponding cell i,j. This device is implemented between the central node of the cell and the common ground node, in the same way as the capacitor.

Boundary and Initial Conditions
Both for the cases of regional flow and exclusively vertical flow, the only initial condition to implement is that of the initial ground temperature (T ), Equation (8). In the network method, this condition is very easy to establish, since it is only necessary to indicate what is the initial voltage at which the capacitors (C or C , ) are charged.
Regarding the boundary conditions, in those cases where the temperature of the contour remains invariant, Equations (3), (4) and (5), it is enough to place a voltage source (V or V , ) of constant value, connected from the corresponding contour node to the common ground node, as can be seen in Figure 6. For those cases in which it is desired that the temperature of the soil surface varies in a stepwise manner, the network method allows the implementation of a stepped function in a simple way, by means of a function commonly called PULSE [42,43]. In this case, the value of the voltage source (V or V , ) will vary with time, as shown in Figure 2, between T and T , depending on the values assigned to the parameters (Table 2) necessary to define the PULSE function: T , T , P , T and P.

Boundary and Initial Conditions
Both for the cases of regional flow and exclusively vertical flow, the only initial condition to implement is that of the initial ground temperature (T s ), Equation (8). In the network method, this condition is very easy to establish, since it is only necessary to indicate what is the initial voltage at which the capacitors (C i or C i,j ) are charged.
Regarding the boundary conditions, in those cases where the temperature of the contour remains invariant, Equations (3)-(5), it is enough to place a voltage source (V i or V i,j ) of constant value, connected from the corresponding contour node to the common ground node, as can be seen in Figure 6.

Boundary and Initial Conditions
Both for the cases of regional flow and exclusively vertical flow, the only initial condition to implement is that of the initial ground temperature (T ), Equation (8). In the network method, this condition is very easy to establish, since it is only necessary to indicate what is the initial voltage at which the capacitors (C or C , ) are charged.
Regarding the boundary conditions, in those cases where the temperature of the contour remains invariant, Equations (3), (4) and (5), it is enough to place a voltage source (V or V , ) of constant value, connected from the corresponding contour node to the common ground node, as can be seen in Figure 6. For those cases in which it is desired that the temperature of the soil surface varies in a stepwise manner, the network method allows the implementation of a stepped function in a simple way, by means of a function commonly called PULSE [42,43]. In this case, the value of the voltage source (V or V , ) will vary with time, as shown in Figure 2, between T and T , depending on the values assigned to the parameters (Table 2) necessary to define the PULSE function: T , T , P , T and P. For those cases in which it is desired that the temperature of the soil surface varies in a stepwise manner, the network method allows the implementation of a stepped function in a simple way, by means of a function commonly called PULSE [42,43]. In this case, the value of the voltage source (V i or V i,j ) will vary with time, as shown in Figure 2, between T 1 and T 2 , depending on the values assigned to the parameters ( Table 2) necessary to define the PULSE function: T D , T R , P W , T F and P.

Contour cell
Finally, for the case of sinusoidal variation of the ground surface temperature (Figure 3), the specific sinusoidal function [31,32] has been used, which defines the instantaneous value of the potential (that is, the temperature) provided by the voltage source. In this case, the parameters to be specified are: T m and T sin (Tables 3 and 4).

Verification of the Model and Applications
In this section, a total of four applications will be presented that will serve both to validate the precision of our models and to illustrate a series of scenarios that can be found in real practical situations of heat transport in saturated porous media in which exists, at the same time, a flow of water.

Verification of the Model. Constant Surface Temperature and Different Vertical Flow Rates
This first application aims to show the high precision that is achieved with the network method in obtaining the numerical solution to this type of problem. The model addressed is that of exclusively vertical flow with constant temperature on the ground surface, as described in Section 2.1.1. For this scenario, Table 5 summarizes the geometric parameters of the problem, as well as the thermal properties of the porous medium, while Table 1 does it with the conditions of initial temperature of the soil and in the lower and upper contours.  Figure 7, the values of the temperature T are represented as a function of the depth Y, once the steady state has been reached and for different values and directions of the vertical velocity v y : upward flow (positive v y ), downward flow (negative v y ) and nonexistent flow (v y equal to 0). As can be seen, and due to the fact that the temperature of the bottom edge is lower than that of the ground surface, for the cases of upward flow there is a cooling of the medium (taking the case of null flow as a reference), the more pronounced, the higher the velocity v y , while for the cases of downward flow, what is observed is a warming of the medium.
The results presented here coincide with the analytical solutions presented by Bredehoeft and Papaopulos [7] and satisfactorily approximate the data collected in the in situ study carried out by Duque et al. [2] in the coastal lagoon of Ringkøbing Fjord (Denmark), as shown in Table 6. As can be seen, the maximum relative error is 7.53% when compared with the real data of Duque et al. [2], although this is considerably reduced to 4.77% when compared with their fitted curve. On the other hand, when we compare with the analytical solution proposed by Bredehoeft and Papaopulos [7], the maximum relative error does not exceed 0.47%. Without a doubt, these are very low errors for the engineering field, which shows that the solutions obtained with our tool are highly accurate and valid for the simulation of real scenarios.
It should be noted that the number of cells into which the medium has been divided is 100 (Table 5), a mesh for which values of the temperature variable have been obtained with maximum relative errors below 0.5%, in comparison with the temperatures provided by Bredehoeft and Papaopulos [7], Table 6. A mesh sensitivity analysis is presented in Table 7, which shows the maximum relative error made by our approximation (taking as reference the analytical solutions [7]) as a function of the number of cells, as well as calculation times. As can be seen, the network method achieves very precise solutions with undemanding grids, since with meshes above 50 cells, the maximum relative error does not exceed 1%. On the other hand, above 100 cells, more demanding meshes hardly reduce the error, while computing times are significantly affected. Therefore, we consider a 100 cell mesh to be suitable. The results presented here coincide with the analytical solutions presented by Bredehoeft and Papaopulos [7] and satisfactorily approximate the data collected in the in situ study carried out by Duque et al. [2] in the coastal lagoon of Ringkøbing Fjord (Denmark), as shown in Table 6. As can be seen, the maximum relative error is 7.53% when compared with the real data of Duque et al. [2], although this is considerably reduced to 4.77% when compared with their fitted curve. On the other hand, when we compare with the analytical solution proposed by Bredehoeft and Papaopulos [7], the maximum relative error does not exceed 0.47%. Without a doubt, these are very low errors for the engineering field, which shows that the solutions obtained with our tool are highly accurate and valid for the simulation of real scenarios.  It should be noted that the number of cells into which the medium has been divided is 100 (Table 5), a mesh for which values of the temperature variable have been obtained

Stepped Surface Temperature and Different Vertical Flow Rates
In this case, the model addressed is the one described in Section 2.1.2, with exclusively vertical flow and stepped temperature variation on the soil surface. Table 8 shows the geometric parameters of the problem and the thermal properties of the ground for this case. On the other hand, Table 2 summarizes the temperature conditions in the lower and upper contours (with the different parameters necessary to define the stepped temperature function), as well as the initial temperature of the medium.  Figure 8 schematizes how the steady state is reached, at a certain depth y i , after the start of the problem, for the case of stepped surface temperature. Once this state is reached, it is observed that the temperature is not constant, but fluctuates between certain values in wave form, the variable ∆T y,max being the amplitude of that wave. Stepped surface temperature and vertical flow. Figure 9 represents the amplitude ( ∆T , ) that the temperature steady wave reaches, as a function of the depth Y, for different values and directions of the vertical velocity v : upward flow (positive v ), downward flow (negative v ) and non-existent flow (v equal to 0). As can be seen (for the specific values of the initial and boundary conditions of temperature, Table 2), as in the previous application, a cooling of the medium occurs when the water flow is upward (taking the v = 0 case as a reference), of greater magnitude the greater the velocity v , while a warming of the medium is manifested when the flow is downward. The mesh used for this application was also 100 cells, which has proven sufficient to achieve high precision with the network method in this type of problem. Stepped surface temperature and vertical flow. Figure 9 represents the amplitude (∆T y,max ) that the temperature steady wave reaches, as a function of the depth Y, for different values and directions of the vertical velocity v y : upward flow (positive v y ), downward flow (negative v y ) and non-existent flow (v y equal to 0). As can be seen (for the specific values of the initial and boundary conditions of temperature, Table 2), as in the previous application, a cooling of the medium occurs when the water flow is upward (taking the v y = 0 case as a reference), of greater magnitude the greater the velocity v y , while a warming of the medium is manifested when the flow is downward. The mesh used for this application was also 100 cells, which has proven sufficient to achieve high precision with the network method in this type of problem. flow (v equal to 0). As can be seen (for the specific values of the initial and boundary conditions of temperature, Table 2), as in the previous application, a cooling of the medium occurs when the water flow is upward (taking the v = 0 case as a reference), of greater magnitude the greater the velocity v , while a warming of the medium is manifested when the flow is downward. The mesh used for this application was also 100 cells, which has proven sufficient to achieve high precision with the network method in this type of problem. Figure 9. Temperature steady wave amplitude as a function of depth for the stepped surface temperature case with different vertical flow rates. Figure 9. Temperature steady wave amplitude as a function of depth for the stepped surface temperature case with different vertical flow rates.

Sinusoidal Surface Temperature and Different Vertical Flow Rates
Now, the scenario addressed corresponds to a case of exclusively vertical flow with sinusoidal temperature variation on the ground surface (as described in Section 2.1.3), whose geometric parameters and thermal properties of the porous medium are summarized in Table 9. The initial and boundary temperature conditions are summarized in Table 3. Table 9. Parameters of the scenario with sinusoidal surface temperature. Vertical flow.

Parameter
Value Units k m 0.84 J/(sm • C) ρ e c e 3.55 × 10 6 J/(m 3 • C) ρ w c w 4.18 × 10 6 J/(m 3 • C) H 1 m n 100 Cell height (∆y) 0.01 m Figure 10 shows schematically how the steady state is reached (at a certain depth y i ) after the start of the problem, for the case of sinusoidal surface temperature. In the stationary situation, it is observed how the temperature is not constant, but oscillates in the form of a wave between certain values, the variable ∆T y,max being the amplitude of such a wave. Figure 11 represents the amplitude (∆T y,max ) that the temperature steady wave reaches, as a function of the depth Y, for different values and directions of the vertical velocity v y . As in the previous applications (and for the specific values, reflected in Table 3, of the initial and boundary temperature conditions), upflows induce a cooling of the medium (taking the v y = 0 case as a reference), of greater magnitude as we increase v y . For downflows, the effect is the opposite (increased temperatures). The simulations were performed again with a 100 cell grid.
Cell height (∆y) 0.01 m Figure 10 shows schematically how the steady state is reached (at a certain depth y ) after the start of the problem, for the case of sinusoidal surface temperature. In the stationary situation, it is observed how the temperature is not constant, but oscillates in the form of a wave between certain values, the variable ∆T , being the amplitude of such a wave.  ) that the temperature steady wave reaches, as a function of the depth Y, for different values and directions of the vertical velocity v . As in the previous applications (and for the specific values, reflected in Table  3, of the initial and boundary temperature conditions), upflows induce a cooling of the medium (taking the v = 0 case as a reference), of greater magnitude as we increase v . For downflows, the effect is the opposite (increased temperatures). The simulations were performed again with a 100 cell grid.

Sinusoidal Surface Temperature and Regional Flow
Finally, we address a scenario where the temperature at the ground surface varies sinusoidally but where the flow is, on this occasion, horizontal (as described in Section 2.2.1). The boundary conditions and the initial soil temperature are summarized in Table  4, while the medium thermal properties and the geometry of the domain (2D) are collected in Table 10. In this sense, it should be noted that a scenario with horizontal water flow from left to right has been chosen, so that the heat existing in the left lateral contour of the problem is transported by the water flow to the porous medium. Regarding its dimensions, a length large enough has been chosen so that there are areas of the porous medium that are not affected by this boundary condition. In this way, we can analyze the effect that the variation in surface temperature has on the porous medium both in the vicinity of the lateral heat source and in distant locations. The chosen grid has been 100 cells in the flow direction (horizontal) by 20 cells vertically (as there is no flow component in that direction, the mesh can be much less demanding). Table 10. Parameters of the scenario with sinusoidal surface temperature. Regional flow.

Sinusoidal Surface Temperature and Regional Flow
Finally, we address a scenario where the temperature at the ground surface varies sinusoidally but where the flow is, on this occasion, horizontal (as described in Section 2.2). The boundary conditions and the initial soil temperature are summarized in Table 4, while the medium thermal properties and the geometry of the domain (2D) are collected in Table 10. In this sense, it should be noted that a scenario with horizontal water flow from left to right has been chosen, so that the heat existing in the left lateral contour of the problem is transported by the water flow to the porous medium. Regarding its dimensions, a length large enough has been chosen so that there are areas of the porous medium that are not affected by this boundary condition. In this way, we can analyze the effect that the variation in surface temperature has on the porous medium both in the vicinity of the lateral heat source and in distant locations. The chosen grid has been 100 cells in the flow direction (horizontal) by 20 cells vertically (as there is no flow component in that direction, the mesh can be much less demanding). Table 10. Parameters of the scenario with sinusoidal surface temperature. Regional flow.

Parameter
Value Units For this case, with a daily period (P) of the sinusoidal wave of surface temperature, it is observed as the amplitude (∆T y,max ) that reaches the steady wave of temperature in the ground; it is not very significant (less than 8% of the amplitude of the surface temperature sinusoidal wave at a distance of H/4 from the surface), being even lower at those depths far from the surface. temperature, it is observed as the amplitude (∆T , ) that reaches the steady wave of temperature in the ground; it is not very significant (less than 8% of the amplitude of the surface temperature sinusoidal wave at a distance of H/4 from the surface), being even lower at those depths far from the surface. This is due to the fact that the heat drag component has a great importance compared to the diffusive one (the regional flow velocity, v , is quite high while the heat conductivity of the rock-fluid matrix, k , presents a relatively low value), so that, for the values of lateral temperature (entry to the domain) and sinusoidal surface function, Table 4, the variation of the porous medium temperature due to the effect of the surface temperature is minimized by the drag effects.
However, when the period (P) of the surface temperature sinusoidal wave is annual, the amplitude (∆T , ) reached by the steady temperature wave in the medium does become significant (about 75% of the amplitude of the surface temperature sinusoidal wave, at a distance of H/4 from the surface), descending, logically, as we move away from the surface, Figure 15. On this occasion, the increase in the period of the surface temperature wave contributes to the diffusive process, so that the drag and diffusion components are more balanced in this case.     This is due to the fact that the heat drag component has a great importance compared to the diffusive one (the regional flow velocity, v x , is quite high while the heat conductivity of the rock-fluid matrix, k m , presents a relatively low value), so that, for the values of lateral temperature (entry to the domain) and sinusoidal surface function, Table 4, the variation of the porous medium temperature due to the effect of the surface temperature is minimized by the drag effects.
However, when the period (P) of the surface temperature sinusoidal wave is annual, the amplitude (∆T y,max ) reached by the steady temperature wave in the medium does become significant (about 75% of the amplitude of the surface temperature sinusoidal wave, at a distance of H/4 from the surface), descending, logically, as we move away from the surface, Figure 15. On this occasion, the increase in the period of the surface temperature wave contributes to the diffusive process, so that the drag and diffusion components are more balanced in this case.

Conclusions
The network method has been revealed as an optimal tool for the numerical modeling and simulation of heat transport problems with water flow coupling in porous media. Its high versatility has made it possible to approach, in a simple way, both one-dimensional models with exclusively vertical heat transport and water flow and more complex regional flow models where it is necessary to resort to two-dimensional geometries.
From the analogy established between the temperature of the porous medium and the electrical potential of the equivalent circuit, the main physical phenomena that occur in the problem, in the form of terms or addends of the differential governing equation, are perfectly reproduced in the network model through the different electrical components that make up the equivalent circuit. Thus, the capacitors have served to implement both the initial temperature of the porous medium and the heat storage process of the ground throughout the simulation. For their part, electrical resistors govern the heat diffusion phenomenon, based on the heat conductivity of the rock-fluid matrix and the size of the elementary cell. The heat drag effects caused by the coupled water flow have been modeled using current sources, which collect properties such as density and specific heat, as well as the velocity of the water flow. Finally, the temperature boundary conditions (constant, stepped or sinusoidal) are easily implemented by voltage sources, duly placed in the corresponding contour nodes.
The network model proposed here has been successfully verified through a 1D application in which the water flow is exclusively vertical. The temperature-depth profiles obtained have been compared with the analytical solutions available in the scientific literature, verifying that the solutions provided by our tool coincide with these, leaving the relative errors dependent on the size of the chosen mesh: relative errors below 1% for meshes above 50 cells, decreasing below 0.5% when the number of cells increases to 100. For all this, the network method is revealed as a powerful numerical tool, in addition to being fast and simple, in the simulation of this type of problem.

Conclusions
The network method has been revealed as an optimal tool for the numerical modeling and simulation of heat transport problems with water flow coupling in porous media. Its high versatility has made it possible to approach, in a simple way, both one-dimensional models with exclusively vertical heat transport and water flow and more complex regional flow models where it is necessary to resort to two-dimensional geometries.
From the analogy established between the temperature of the porous medium and the electrical potential of the equivalent circuit, the main physical phenomena that occur in the problem, in the form of terms or addends of the differential governing equation, are perfectly reproduced in the network model through the different electrical components that make up the equivalent circuit. Thus, the capacitors have served to implement both the initial temperature of the porous medium and the heat storage process of the ground throughout the simulation. For their part, electrical resistors govern the heat diffusion phenomenon, based on the heat conductivity of the rock-fluid matrix and the size of the elementary cell. The heat drag effects caused by the coupled water flow have been modeled using current sources, which collect properties such as density and specific heat, as well as the velocity of the water flow. Finally, the temperature boundary conditions (constant, stepped or sinusoidal) are easily implemented by voltage sources, duly placed in the corresponding contour nodes.
The network model proposed here has been successfully verified through a 1D application in which the water flow is exclusively vertical. The temperature-depth profiles obtained have been compared with the analytical solutions available in the scientific liter-ature, verifying that the solutions provided by our tool coincide with these, leaving the relative errors dependent on the size of the chosen mesh: relative errors below 1% for meshes above 50 cells, decreasing below 0.5% when the number of cells increases to 100. For all this, the network method is revealed as a powerful numerical tool, in addition to being fast and simple, in the simulation of this type of problem.
Finally, and once the precision of our model has been demonstrated, a series of scenarios, both in 1D (vertical flow) and 2D (horizontal flow) domains and with different temperature conditions of the ground surface, has been addressed in order to illustrate the versatility of our tool, obtaining the ranges of values (in the form of maximum amplitudes) between which the temperature fluctuates at a certain depth and/or position of the porous medium.
In future research, the model presented here could be extended to problems where the properties of the porous medium change over time (swelling and erosion processes), soils with non-homogenous thermal and hydraulic characteristics or other surface temperature functions, different from those presented here, including tabulated data from actual temperature distributions.

Acknowledgments:
We would like to thank Fundación Séneca for the scholarship 21271/FPI/19. Fundación Séneca. Región de Murcia (Spain) was awarded to José Antonio Jiménez-Valera, which has allowed us to carry out this investigation.