62. Infiltrating Heat Front
An analytical solution used in chemical engineering for modeling concentration in packed-bed reactors , commonly referred to as a “Danckwerts” (or “third-type”) boundary condition, is adapted here for confirming the accuracy of heat transport as solved by the UZE package in MODFLOW 6. The analytical solution has the following characteristics: (1) it solves for total energy flux (advection and conduction) along a 1-dimensional (1D) profile [i.e., there is no conduction (“thermal bleeding”) with the surrounding materials], (2) heat flux can only enter the active domain with the inflow (in this case, the infiltration) and it does not require the temperature at the boundary to be equal to the temperature of the infiltration, and (3) it solves for total heat flux (instead of temperature) throughout the 1D domain.
The Danckwerts analytical solution takes the following form:
\(q_{T_z} = q_{T_0} + \dfrac{1}{2} \left( q_{T_{infil}} - q_{T_0} \right) \left( \textit{erfc} \biggl\{ \dfrac{z-\nu t}{2 \sqrt{Dt}} \biggr\} + \textit{exp} \Bigl\{ \dfrac{\nu z}{D}\Bigr\} \cdot \textit{erfc} \biggl\{ \dfrac{z+\nu t}{2 \sqrt{Dt}} \biggr\} \right)\)
where \(q_{T_z}\) is the heat flux and depth \(z\), \(q_{T_0}\)is the infiltrating heat flux at \(t=0\), \(q_{T_{infil}}\) is the amount of infiltrating heat flux at \(t > 0\), \(\textit{erfc}\) is the complementary error function, \(z\) is the distance from infiltrating heat flux boundary, which, in this example is the depth below land surface, \(t\) is time (in days), \(\nu\) represents the “thermal convection velocity” determined from,
\(\nu=q \cdot \dfrac{\rho_w C_{p_w}}{S_{w_z} \theta \rho_w C_{p_w} + \left( 1-\theta \right) \rho_s C_{p_w}}\)
with \(q\) equal to the volumetric infiltration rate, \(\rho_w\) is the density of water, \(C_{p_w}\) is the heat capacity of water, \(S_{w_z}\) is the saturation at depth \(z\), \(\theta\) is the water content, \(rho_s\) is the density of a aquifer solids, and \(C_{p_w}\) is the heat capacity of the aquifer solids. Additionally, \(D\) represents the bulk thermal diffusivity,
\(D=\dfrac{k_{T_{bulk}}}{S_{w_z} \theta \rho_w C_{p_w} + \left(1 - \theta \right) \rho_s C_{p_w}}\)
and \(k_{T_{bulk}}\) is the bulk thermal conductivity represented by,
\(k_{T_{bulk}} = S_{w_z} \theta k_{T_w} + \left(1-\theta \right) k_{T_s}\)
and \(k_{T_w}\) and \(k_{T_s}\) are the thermal conductivities of the water and aquifer material, respectively.
62.1. Example description
A 1D model grid of the unsaturated zone is used to simulate the downward migration of an infiltrating heat front (Figure 62.1). Steady flow conditions are simulated with the GWF model. The UZF package simulates flow through the unsaturated zone with a constant infiltration rate of 0.01 \(\frac{m}{d}\). An initial temperature of 10\(^{\circ}C\) is specified for the entire model domain. After steady flow and transport conditions are established with a quasi-steady-state stress period, the temperature of the infiltration is increased to 20\(^{\circ}C\) resulting in an energy source loading of 9.68\(\frac{J}{s}\) (calculated from \(q_{infil} \cdot T_{infil} \cdot \rho_w \cdot C_{p_w}\)). Other pertinent parameter values are provided in Table 62.1.
A constant head boundary is placed at the bottom of the model to remove water that recharges the water table in order to prevent the water table from rising into the unsaturated column (Figure 62.1).
Figure 62.1 View of 1-dimensional model setup. Total thickness of the unsaturated zone is 10 \(m\) and is discretized with 100 cells that are each 10 \(cm\) thick.
Parameter |
Value |
|---|---|
Number of layers (\(-\)) |
101 |
Number of rows (\(-\)) |
1 |
Number of columns (\(-\)) |
1 |
Number of simulated periods (\(-\)) |
2 |
Cell width (\(m\)) |
1.0 |
Cell length (\(m\)) |
1.0 |
Cell thickness (\(m\)) |
0.1 |
Initial head (\(m\)) |
0.05 |
Top of the model grid (\(m\)) |
10.00 |
Initial temperature (\(^{\circ}C\)) |
10.0 |
Advection scheme (\(-\)) |
Upstream |
Longitudinal mechanical dispersion term (\(m\)) |
0.0 |
Porosity (\(-\)) |
0.2 |
Density of dry solid aquifer material (\(\frac{kg}{m^3}\)) |
1500.0 |
Heat capacity of dry solid aquifer material (\(\frac{J}{kg \cdot ^{\circ} C}\)) |
760.0 |
Density of water (\(\frac{kg}{m^3}\)) |
1000.0 |
Heat capacity of water (\(\frac{J}{kg \cdot ^{\circ} C}\)) |
4183.0 |
Thermal conductivity of water (\(\frac{W}{m \cdot ^{\circ} C}\)) |
0.5918 |
Thermal conductivity of solid aquifer material (\(\frac{W}{m \cdot ^{\circ} C}\)) |
0.27 |
Infiltration rate (\(\frac{m}{d}\)) |
0.01 |
Contant head at the model outlet (\(m\)) |
0.05 |
Residual water content of the unsaturated zone (\(-\)) |
0.0001 |
Saturated water content of the unsaturated zone (\(-\)) |
0.20 |
Initial water content of the unsaturated zone (\(-\)) |
0.055 |
Brooks-Corey epsilon parameter (\(-\)) |
4.0 |
Vertical hydraulic conductivity of the unsaturated zone (\(\frac{m}{d}\)) |
1.0 |
Initial temperature in simulation domain (\(^{\circ} C\)) |
10.0 |
Temperature of infiltrating water (\(^{\circ} C\)) |
20.0 |
62.2. Example Results
The heat front that migrates downward through the unsaturated zone as a result of energy loading associated with the infiltration is shown in Figure 62.2 at 10, 50, and 100 days. Figure 62.3 shows the same comparison broken out into 3 subplots, but further parses the downward heat migration into its advective (dark green) and conductive (light green) components. Where the temperature gradients are steepest, the conductive flux of energy plays a more prominent role in the downward migration of heat.
Figure 62.2 Comparison of simulated migration of infiltrating heat front to an analytical solution
Figure 62.3 The same infiltrating heat fronts as shown in Figure 62.2, but highlights the advective and conductive heat fluxes separately
62.3. Jupyter Notebook
The Jupyter notebook used to create the MODFLOW 6 input files for this example and post-process the results is: