The Effect of Groundwater Level on Forced Convection Heat Transfer from a Buried Cylinder

The effect of groundwater level on forced convection heat transfer from a cylinder of 1 m in diameter buried at a depth of 10 m is studied numerically. The porosity and permeability of the soil are fixed at 0.5 and 1.0 × 10–11 m2, respectively, throughout the present study, while the inlet velocity U0 and the position of the phreatic surface (h), measured from the center of the cylinder, are varied in the ranges from 10–9 to 10–3 m s−1 and from − 4 to 4 m, respectively. Under these conditions, the heat fluxes on the cylinder are computed numerically. Various approximate solutions are also proposed and compared with those obtained numerically. The heat transfer results are also presented in nondimensional form for generality. The paper is divided into three parts, depending on the relative positions of the phreatic surface to the heated cylinder. First, the heat transfer rates are analyzed when the phreatic surface is above the cylinder, so that the cylinder is completely submerged in the water. In the second case, the phreatic surface is assumed to lie below the cylinder, so that the cylinder is left in the water-free space. In the third case, the phreatic surface is in contact with the cylinder surface, so that the cylinder is partially submerged in groundwater. For the above respective cases, three different approximate analyses are presented, and their validities and limitations are evaluated in comparison with the numerical calculations.


Gravitational acceleration vector (m s −2 ) h
Position of phreatic surface measured from the center of the cylinder (m) h Average heat transfer coefficient (W m −2 K −1 ) k Thermal conductivity (W m −1 K −1 ) K Permeability of porous medium (m 2 ) Nu Local Nusselt number: Dynamic pressure (N m −2 ) Pe Péclet number U 0 d∕ (-) q ′′ Local heat flux on the cylinder (W m −2 ) q ′′ Average heat flux on the cylinder (W m −2 ) Heat flux vector (W m −2 ) s Arc length of circular cylinder (m) T Temperature (K) u x-Component of u (m s −1 ) v y-Component of u (m s −1 ) Groundwater velocity vector (m s −1 ) U 0 Incoming groundwater velocity (m s −1 )

Introduction
Flow in porous media is a very important research subject because of its wide-ranging applications, such as geophysical processes to advanced technologies, e.g., fuel cell development as reported by Rosli et al. (2009). In geophysics, the dynamics of fluid flows in porous media was first studied in connection with groundwater hydrology, as summarized by such authors as Bear (1979) and Hillel (1982). Geothermal energy development is another closely related example, where the idea of transport in porous media plays an important role. The geothermal energy is stored in the deep underground permeable layer, called geothermal reservoir, in the form of hot brine, and it must be tapped to the Earth's surface for electricity generation. The geothermal reservoir's dynamical behavior necessary for sustainable steam production is often analyzed with a model of saturated porous medium. Early development of geothermal energy and its conversion to electricity was well documented by Di Pippo (1980). More recent information on this subject is also available in Di Pippo (2015). Early work on heat transfer and convection in geothermal reservoirs was conducted by Cheng (1978), and Garg and Kassoy (1981). Recent understanding on deep underground reservoirs, however, is that they are more properly modeled as a fracture network system rather than a porous medium, as demonstrated by Zimmerman and Bodvarsson (1996), and more recently Aghajannezhad et al. (2021). Because of high temperature in the reservoir, steam-liquid two-phase convection is possible, as shown by Straus (1977, 1979), and a recent review on this subject is available in Valavanides (2018). Low-enthalpy geothermal energy can be also tapped with ground-coupled heat pumps (GCHP). This direct use of thermal energy for building heating and agricultural greenhouses is a fast-growing industry today, as reviewed by Sarbu and Sebarchievici (2016), and amply documented in the handbook from ASHRAE (American Society of Heating and Air-Conditioning Engineers) (2021).
Heat transfer from a heated object placed in the ground is often encountered in both high-and low-enthalpy geothermal energy developments. For example, Cheng and Minkowycz (1977) investigated natural convection from a hot intrusion in a geothermal system. Merkin (1979) studied natural convection from a general two-dimensional heated body. Since then, there have appeared numerous research reports on natural convection around heated bodies in porous media with various different flow models, such as: non-Darcy flow models; non-Newtonian fluid-filled porous media; non-uniform permeability distributions; conduction and convection conjugate problems. The aforementioned problems are compiled in reference books such as those by Bejan and Kraus (2003), Nield and Bejan (2013) and Vafai (2015).
Although it is less prevalent, forced convection heat transfer from a heated circular cylinder embedded in a porous medium is an important subject in connection with borehole heat exchangers (BHE) employed in GCHP, as described by Ngo and Lai (2009). Groundwater flow influence on BHE performance was examined by such authors as Diao et al. (2004), Katsura et al. (2006), Lee and Lam (2007), Wang et al. (2009), Holmberg et al. (2016, Li et al. (2020) and Casasso and Sethi (2014). According to these authors, the presence of groundwater flow noticeably enhances BHE performance. In order to understand the mechanism of this enhancement, forced convection heat transfer from an isothermal cylinder may serve as a basis for the analysis.
Early work on forced convection heat transfer from a heated cylinder dates back to Cheng (1982), who employed the Darcy model for fluid flow through a porous matrix and the local thermal equilibrium condition between the porous matrix and saturated fluid. However, he did not report the average heat transfer coefficients, which are often more important from the engineering point of view. Integrating the local values of heat flux over the cylinder surface, the average heat transfer coefficients were presented in the book by Nield and Bejan (2013). Forced convection heat transfer from a heated elliptic cylinder was analyzed by Kimura (1988a), which plays a role to bridge heat transfer from a flat plate and that from a circular cylinder. Kimura (1989) also considered transient heat transfer from an isothermal circular cylinder. When a constant flux is specified over a circular cylinder, Kimura (1988b) carried out a boundary-layer analysis and devised groundwater velocimetry. A groundwater velocimeter based on this idea was constructed and reported by Kimura (2015).
Experimental work on the heat transfer in a packed bed was performed by Nasr et al. (1994). Their main interest was the enhancement of heat transfer from an isothermal cylinder by high-conductivity packed-bed materials. The validity of the above-mentioned averaged heat transfer characteristics was marginally verified in their experiments. The forced convection in a packed bed is mainly focused on a high fluid velocity regime. Therefore, it is natural that the research has been directed toward heat transfer analyses using various non-Darcy flow models, and the eventual relaxation of the local thermal equilibrium (LTE) condition between the fluid and the porous matrix. These non-Darcy flow models and local thermally non-equilibrium (LTNE) models require the introduction of new parameters. Therefore, the mathematical system that describes the conservation laws becomes increasingly more complex, and we are still far from a complete understanding of heat transfer behavior in those models. Nonetheless, research reports on packed beds have been published by many authors. For example, heat transfer from an isothermal cylinder based on various extended Darcy models was reported on by such authors as Thevenin and Sadaoui (1995) and Layeghi and Nouri-Borujerdi (2004). It is expected that the LTE condition between the fluid and solid phases will become less valid as the heat transport speeds in the fluid and solid become significantly different, which is conceivable when the flow speed increases. This local thermally non-equilibrium problem was investigated with the Darcy flow model by Rees et al. (2003) and Wong et al. (2004). Al-Sumaily et al. (2012) presented a numerical analysis with the Brinkman-extended non-Darcy flow model and LTNE condition.
The above brief review on forced convection from a heated circular cylinder reveals the basic assumption common to all: the porous medium is completely saturated with flowing fluid. On the other hand, McKibbin and Kimura (2014) were the first to propose a problem in which the heated cylinder was partially submerged in groundwater, i.e., the cylinder was positioned near the interface between a water-saturated lower layer and an air-saturated upper layer, as shown in Fig. 1. The problem was inspired by the need to evaluate the heat transfer from a buried circular tank filled with water, which can serve as a possible heat reservoir for heat pump systems in the heating and cooling of residential or office buildings. In such a case, the tank is frequently either fully or partially submerged in the groundwater. In order to define the heat transfer problem from the water-filled tank, it is necessary to introduce two additional parameters: the relative positions of the groundwater level and ground surface to the buried cylinder. This makes the proposed problem rather intractable. Although the paper demonstrated the mathematical power to compute the heat transfer for several sets of the above two geometric parameters, it suffers from the lack of a general picture of heat transfer characteristics as the parameters involved are varied.
In the present paper, we relax one of the geometric conditions, i.e., the relative position of the ground surface to the tank, which means that the tank is positioned far from the ground surface. This simplification makes the problem much more tractable. Therefore, the present paper is only concerned with the groundwater level relative to the tank. We further model the tank as an isothermally heated circular cylinder in the two-dimensional plane. The present paper is divided into the following three parts: (a) the cylinder is completely submerged in shallow groundwater; (b) the cylinder is positioned above the groundwater surface, also called the phreatic surface, but only by a small distance; (c) the cylinder is partially submerged in groundwater. We show in this paper a general picture of heat transfer characteristics for the respective categories. Case (a) is illustrated schematically in Fig. 1.
We start with the mathematical formulation and the definition of the domain to be analyzed, along with the boundary conditions. Then, we conduct numerical computations to generate heat transfer results for a given set of parameters in the respective cases. We also provide simple analyses to delineate physical reasons for the observed heat transfer characteristics. Lastly, the conclusions of the present work are drawn.

The Governing Equations
In this section, we describe a mathematical system derived from the mass, momentum and energy conservation laws in a saturated porous medium. We are concerned with the heat transfer in the ground, where the groundwater velocity is usually small, typically 0.001 to 1 mm s −1 in Darcy velocity, and the average particle size composing the porous matrix is typically less than 1 mm. The permeability is usually of the order of 10 −11 m 2 . Due to the above conditions, it is generally believed that the Darcy law and the LTE condition are guaranteed. The conservation equations for mass, momentum and energy in vectorial form are: Here, we employ the Brinkman-extended momentum equation. In the above equations, is a velocity vector, is the vector of gravitational acceleration, is the identity matrix, and P is the dynamic pressure. , and c p are the fluid density, viscosity and specific heat, respectively. The effective thermal conductivity of the porous medium is defined as where the subscripts s and f indicate the porous matrix and the saturating fluid, respectively, with as the porosity and K as the permeability.

The Physical Properties and Parameters Employed in the Present Paper
The physical properties of air and water at T 0 = 293.15 K (inlet water temperature) are used for the present numerical computations. The cylinder temperature is fixed at T h = 303.15 K . The physical properties of the porous matrix are taken from the data generally accepted in groundwater hydrology, e.g., Bear (1979), and they are listed in Table 1. In a typical sand and gravel layer in an aquifer, it is known that the porosity varies from 0.4 to 0.5, and the mean grain size is between 0.1 and 1 mm. The permeability is in the order of 10 −11 m 2 . It is worth noting that the Carman-Kozeny formula, example, see Nield and Bejan (2013), gives 2.78 × 10 −11 m 2 for spherical beads for which d p = 0.1 mm, and the porosity of = 0.5 . Assuming a typical permeable aquifer, the parameters involved in the present work are summarized in Table 2.

The Computational Domain and the Boundary Conditions
The groundwater level is always positioned close to the cylinder. As shown in Fig. 2, the domain of analysis is 20 m × 30 m, and the center of the circular cylinder (1 m in diameter) is positioned at the mid-height and 10 m from the left boundary. The thermal boundary conditions on the top, bottom and left boundaries are fixed at a constant temperature, and the right boundary is adiabatic. The groundwater enters into the lower layer of the domain through the left boundary at a specified velocity, while the top (phreatic surface) of the water-saturated layer is stress-free and the bottom is non-slip. The right boundary is set to the usual outlet condition with a constant dynamic pressure. The groundwater surface is horizontal, and the capillary fringe is ignored; on the   other hand, the upper layer is saturated with dry air, which is assumed to be completely stagnant, and the velocity vector there is set to zero in the momentum and energy equations. Because of the local thermal equilibrium assumption, we solve a single energy equation, with two different fluid properties, i.e., the lower layer has the thermal properties of a water-sand mixture, while the upper layer has the properties of an air-sand mixture. Therefore, the air-saturated upper layer is essentially modeled as a rectangular solid region. Referring to Fig. 2, the domain of analysis and the associated boundary conditions are mathematically expressed as:

Numerical Methodology and Validation
The set of the governing Eqs.
(1)-(3) were solved together with the above boundary conditions. The COMSOL Multiphysics software Ver.5.4, a finite-element-based numerical solver, was employed to solve the problem. Non-uniform elements were generated over the region external to the cylinder, with mesh refinement being employed near the surface of the cylinder. The total number of elements was 140,540. The local and average Nusselt numbers were compared with those obtained by Thevenin and Sadaoui (1995) for the Darcy numbers Da = 10 −2 , 10 −4 , the Péclet numbers Pe = 1 − 10 3 and the porosity = 0.9 . The local Nusselt number distribution on the cylinder surface at Pe = 7.0 is shown in Fig. 3a. Good agreement is observed. The slight difference that is apparent may be due to the rather coarse grid used by Thevenin and Sadaoui (1995). As regards the average Nusselt numbers, the two results are almost identical, as shown in Fig. 3b. It is also worth noting that for Da = 10 −4 the average Nusselt numbers are very close to the boundary-layer solution with the Darcy flow model, which reads as Nu = 1.015Pe 1∕2 (Nield and Bejan 2013).

Numerical Results
The temperature fields near the cylinder for four different water velocities are shown in Fig. 4. In these cases, the cylinder is completely submerged, and the phreatic surface is Fig. 3 The present numerical results for local and average Nusselt numbers compared with those by Thevenin and Sadaoui (1995); a local Nusselt numbers at Pe = 7, b average Nusselt numbers as a function of the Péclet number for two different Darcy numbers positioned at 1.0 m from the cylinder center. In the figure, the vertical position of the cylinder center is shown with a dashed-dotted line, and the phreatic surface is indicated with a dashed line. When the water velocity is less than U 0 = 10 −7 m s −1 , it is evident that the flow effects are not very visible, and that the heat transfer is dominated by conduction. On the other hand, when the velocity exceeds U 0 = 10 −5 m s −1 , a thin thermal boundary along the cylinder is clearly visible. It should be also noted that only the temperature in the water-saturated zone is strongly influenced by the heated cylinder. This situation is further pronounced at U 0 = 10 −3 m s −1 . The average Nusselt numbers for all of the combinations of velocity and phreatic surface position are given in Table 3.
In parallel with the temperature fields, we also show the local heat fluxes over the cylinder in Fig. 5. Figure 5a, b clearly indicates that the heat transfer rates are functions of the phreatic surface levels at small groundwater velocities. On the other hand, Fig. 5c, d, for large groundwater velocities, shows no dependence on the phreatic surface level. This further enforces our previous argument.

Effects of the Phreatic Surface
The series of temperature fields at different groundwater velocities in Fig. 4 prompts us to look in more depth at the relative thickness between the thermal boundary layer and the  submerged depth. When the submerged depth is given as a condition, there is a critical velocity, beyond which the thermal boundary layer is confined to the water-saturated zone, and consequently the presence of the phreatic surface does not influence the heat transfer rate over the cylinder. In this case, the heat transfer can be obtained in the same way as a case in which the region is completely saturated with groundwater.
Here, we try to determine the critical velocity when the depth of the submerged cylinder is given. Since this is only a rough argument, we are only concerned with the average value of the boundary-layer thickness. According to the boundary-layer theory, the average thickness over the cylinder is inversely proportional to the square root of the Péclet number (Nield and Bejan 2013). Therefore, when the depth is equal to the thermal boundarylayer thickness , the corresponding velocity, which is expressed as a function of the Péclet number, is the critical velocity; in this case, where a is an arbitrary O(1) constant. The value of a will be determined by comparing the average Nusselt number Nu ∞ with h = 10 m and the numerical results for Nu with an arbitrary submerged depth. We compute the relative differences of the two defined by Eq. (6) and show them in Fig. 6: In the figure, it should be noted that the Péclet number appearing in the abscissa is divided by 2.74 in order to avoid fractional numbers in the abscissa scales. In this way, we can easily relate the value of the Peclet number with the actual incoming water velocity, i.e., Pe/2.74 = 1 corresponds to U 0 = 10 −6 m s −1 . It seems that a = 2 in Eq. (5) fits well in order to divide the range into small and large ΔNu regions. For example, when the submerged depth is equal to the radius of the cylinder, i.e., h/d = 1 in Fig. 1, the heat transfer can be estimated by the boundary-layer theory as long as the Péclet number is greater than 16.

Numerical Results
The numerical computations are performed when the phreatic surface is positioned below the cylinder. In this case, h takes a negative value and it must satisfy the condition h < −d∕2 . In Fig. 7, the temperature fields are shown for two different groundwater velocities, U 0 = 10 −6 m s −1 and 10 −3 m s −1 , when the phreatic surface is set to h = −1 m , where the dashed-dotted and dashed lines indicate the cylinder center and the phreatic surface, respectively, as before. When the velocity is U 0 = 10 −6 m s −1 , a temperature rise near the cylinder can be observed both in the water-and the air-saturated layers. On the other hand, no temperature rise is observed in the water-saturated layer when the velocity is 10 −3 m s −1 . This is because the flowing groundwater quickly sweeps away the conducting heat from above, when the velocity is large. Therefore, the phreatic surface is kept at a nearly constant temperature equal to the incoming water temperature.

Approximate Evaluation of Heat Transfer Rate Due to Mirror-Image Method
The above observation of the temperature fields suggests the possibility of obtaining an analytical solution for large groundwater velocities. In that case, the phreatic surface can be assumed to be at a constant temperature equal to that of the incoming water. This is the problem of heat conduction in a semi-infinite space, where a constant temperature cylinder is placed above the reflection boundary at a different constant temperature. The problem can be solved with the mirror-image method, as described in the book by Eckert and Drake (1972). The average heat flux over the cylinder is given in nondimensional form by In Fig. 8, we show local heat flux variations obtained by numerical computations. When h = − 4 m, the heat fluxes are nearly independent of the groundwater velocity. This is understandable, because the magnitude of the temperature gradient becomes smaller as one moves away from the heated object. Therefore, the thermal property there is less influential on the heat flux over the cylinder. On the other hand, when h = − 1 m, the groundwater velocity greatly influences heat transfer rates over the cylinder. In this case, contrary to the former example, the temperatures very close to the cylinder are significantly altered by the presence of groundwater. The heat flux becomes large as the velocity becomes large. This is because the nearly constant temperature on the phreatic surface at a large velocity forms a steep temperature gradient under the cylinder. In Fig. 9, we plot the average heat flux as a function of the position of the phreatic surface, in both dimensional and nondimensional forms. The dashed line is obtained by using Eq. (7), whereas the other lines are obtained numerically for the different velocities and phreatic surface levels. It is clearly seen that Eq. (7) and the numerical computations give nearly identical results when the water velocity is greater than U 0 = 10 −4 m s −1 . At U 0 = 10 −5 m s −1 , the values predicted by Eq. (7) are still useful if errors of about 20% are accepted. However, at U 0 = 10 −6 m s −1 , the numerical solutions start deviating greatly from those given by Eq. (7). As the groundwater velocity decreases, the position of the phreatic surface becomes less influential on the resulting heat transfer. This, as discussed before, is due to the fact that the heat transfer is essentially determined by the temperature gradient in the region close to the cylinder, and the far-field temperatures become less significant. It is also instructive to compute the conductive heat flux on an isothermal cylinder placed in an air-saturated infinite porous medium with the physical properties in Table 1 and the same temperature difference at a very large time. According to Carslaw and Jaeger (1946), the average heat flux on the

Numerical Results
The third case is a situation in which the cylinder is partially submerged in the groundwater. In the present paper, this is equivalent to saying that the parameter h is in the range −0.5 m < h < 0.5 m . In Fig. 10, the temperature fields near the cylinder for U 0 = 1.0 × 10 −3 (m s −1 ) for three different positions of phreatic surface level are shown. A large rise in temperature in the dry upper layer and thin thermal boundary layer in the lower water-saturated zone is common to all three figures. This implies that the thermal resistance over the cylinder facing to the dry upper layer is much greater than that on the portion facing to the water-saturated lower layer. The figures indicate that this is true irrespective of the phreatic surface levels. This point is further demonstrated in Fig. 11, where the local heat fluxes are plotted as a function of the angle θ. In the graphs, the heat fluxes on the portion wetted with groundwater are generally much greater than those on the dry portion. This feature is more pronounced as the groundwater velocity increases. However, even at the velocity of U 0 = 1.0 × 10 −6 (m s −1 ) , a substantial difference between the two is still visible.
An additional interesting observation is made when h < 0 . As seen in Fig. 11c, a spikelike heat flux appears at the circumferential position at which the groundwater first strikes the cylinder. There, the incoming cold groundwater meets the hot cylinder wall with almost no deceleration in the velocity. This situation forms a singularity in the temperature gradient at that position; on the other hand, for h > 0 , the incoming groundwater must reverse  stagnant region is much smaller, and the heat flux starts rising sharply, yet no spike-like heat flux is produced, as seen in Fig. 11b. In summary, the numerically obtained average heat fluxes over the cylinder are listed with various combinations of the groundwater velocity and the phreatic surface level in Table 4.

Rough Estimate of Average Heat Transfer Rate
The above observations prompt us to formulate a method that can approximately evaluate the average heat flux on the cylinder, or the total heat flow from the cylinder. Suppose that we have a given condition for the incoming groundwater velocity, and that the cylinder is completely submerged. If this is the case, the average Nusselt number is given by Nu = 1.015Pe 1∕2 . If the cylinder is partially submerged, the first approximation is to neglect the heat flow on the cylinder wall contacting with the upper dry layer, compared with the heat flow on the lower wet layer, as seen in Fig. 11. That is to say, the heat flow takes place only on the portion where the cylinder wall is contacting with the water-saturated lower layer. In doing so, we first define the fraction of the arc length, s, over the total peripheral length, In order to compute the average heat flux, we simply multiply the average heat flux from the completely submerged case with the fraction defined by Eq. (8), i.e., The results obtained from Eq. (9) are listed in Table 5. The above results compare well with the numerical ones listed in Table 4, especially when both the groundwater velocity and the wetted fraction are large. This is quite understandable, because the basis of the approximate method is Nu = 1.015Pe 1∕2 , which is derived from the boundary-layer assumptions. Also, since the heat flow from the portion that is in contact with the dry layer is completely neglected in Eq. (9), the approximation must be more accurate as the fraction is large and close to unity. This translates into that the average heat flux values in the upper left region of the two tables compare well, whereas those in the lower right region of the two compare rather poorly. Nonetheless, the simplicity of Eq. (9) should be appreciated. The above observations can be summarized graphically in Fig. 12.

Concluding Remarks
Forced convection heat transfer from a circular cylinder placed in soil in the presence of the phreatic surface has been investigated. The two-dimensional space subject to analysis is composed of two horizontal layers, where the upper one is saturated with stagnant air, and the lower is saturated with flowing groundwater. The ground surface is fixed 10 m above the cylinder, and the problem is addressed for the three different phreatic surface levels relative to the cylinder: (a) the cylinder is completely submerged in shallow groundwater; (b) the cylinder is positioned above the groundwater surface, but only by a small distance; (c) the cylinder is partially submerged in groundwater. The physical properties of the porous matrix are chosen from those in a typical permeable aquifer.  When the heated cylinder is completely submerged, the critical groundwater velocity is determined by Eq. (5), beyond which the average heat transfer can be computed as that in an infinite groundwater-saturated medium. 2. When the heated cylinder is situated above the phreatic surface level, a simple formula to compute the average heat transfer rate is proposed. Its validity has been extensively tested via numerical simulations. It has been proved that the formula in Eq. (7) gives good predictions when the groundwater velocity is such that Pe > 27.4. 3. When the heated cylinder is partially submerged, the heat transfer is essentially controlled by the arc length contacting with the water-saturated layer. This observation allows us to introduce the simple formula of Eq. (9), i.e., the heat transfer rate for the completely submerged case multiplied by a wetted arc length ratio to the circular perimeter. Extensive numerical verifications have been conducted, and it is found that the formula gives good predictions, particularly when the groundwater velocity and the wetted arc length ratio are large.