Process Modeling of Soil Thermal and Hydrological Dynamics

Nawa Raj Pradhan, Charles W. Downer and Sergei Marchinko

## Abstract

To explicitly simulate the soil thermal state effects on hydrological response, the soil thermal regime, frozen soil, and permafrost simulation capability in the Geophysical Institute Permafrost Laboratory (GIPL) model have been included into the physically based, distributed watershed model Gridded Surface Subsurface Hydrologic Analysis (GSSHA). The GIPL model is used to compute a vertical soil temperature profile in every lateral two-dimensional GSSHA computational element using the soil moisture information from hydrologic simulations in GSSHA; GSSHA, in turn, uses this temperature and phase, ice content, and information to adjust hydraulic conductivities for both the vertical unsaturated soil flow and lateral saturated groundwater flow. This two-way coupling increases computational accuracy in both models by providing additional information and processes not previously included in either. The soil moisture physical state is defined by the Richards Equation, and the soil thermal state is defined by the numerical model of phase change based on quasi-linear heat conduction equation. Results from the demonstration site, a head water sub-catchment at the peak of the Caribou-Poker Creeks Research Watershed, representing Alaskan woodland and tundra ecosystem in permafrost-active region, indicated that freezing temperatures reduce soil thermal conductivity and soil storage capacity, thereby increasing overland flow and peak discharges.

Keywords: soil thermodynamics, soil hydrodynamics, soil temperature, soil moisture, soil hydraulic conductivity, soil infiltration

#### 1. Introduction

Frozen soil has the potential to significantly affect hydrological response from hydrologic units of any scale, from the plot size to globally. In the United States, land surface freezing is an important hydrologic factor during winter in low temperature prevailing areas [1]. Approximately 60% of the Northern Hemisphere land surface is frozen in winter [2]. Frozen ground, heavy rains, and rapid snowmelt provide favorable conditions for extreme flooding [3]. Given the potential for significant consequences, frozen soil should be considered an important hydrologic factor in analysis of the potential for flooding, especially in the early spring when and where frozen ground, heavy rains, and rapid snowmelt prevail. Frozen soil hydrologic effects at larger scales, regional and global, include contribution of net freshwater to the Arctic Ocean [4], affecting ocean salinity and global thermohaline circulation [5]. While global and regional hydrologic and climate studies require a realistic representation of the hydraulic and thermal properties [6], a physical process-based representation of these properties starts from a small, point, or field scale. Plot-scale studies by Dunne and Black [7] and Stähli et al. [8], among others, show that a catchment's rainfall runoff generation process is altered when infiltration is reduced during the transition of the soil water phase toward the freezing state. As the distribution of frozen soils is highly variable in both space and time [9], to fully understand both the occurrence of frozen soil and the subsequent impact on hydrology requires a physics-based, distributed analytical tool that provides a two-way coupling of the soil thermodynamic and hydrologic state to provide insight into the temporal and spatial variability of soil thermal state and its effects on the time and space distributions of the hydrological response.

2.1.1 Unsaturated zone model

(cm h�<sup>1</sup>

∂ ∂x

105

Txx ∂h ∂x 

of the Richards Equation is solved [20, 21]:

Cmð Þ ψ

rated and unsaturated conditions [22].

2.1.2 Saturated groundwater model

þ ∂ ∂x

problem can be described as [20]

∂ ∂x

2.2 Soil thermodynamic model

Txy ∂h ∂y 

Kxxb <sup>∂</sup>Ews ∂x 

where T is the transmissivity (m<sup>2</sup> s

∂ψ <sup>∂</sup><sup>τ</sup> � <sup>∂</sup> ∂z

Process Modeling of Soil Thermal and Hydrological Dynamics

DOI: http://dx.doi.org/10.5772/intechopen.84414

[23] equations, as extended by Huston and Cass [24].

flow equations as described by Pinder and Bredehoeft [25]:

þ ∂ ∂y

> þ ∂ ∂y

Details of the solution can be found in Downer [20].

The unsaturated, or vadose, zone controls the flux of water from the land surface to the saturated groundwater zone and determines the partitioning of water into runoff, infiltration, ET, and groundwater recharge. GSSHA provides many ways to analyze these processes including an integrated solution of soil water movement and state described by the Richards Equation [19]. The Richards Equation [19] is considered the most physically correct mathematical representation of the vadose zone. While flow in the vadose zone is in three dimensions, the predominant direction of flow is vertical. In GSSHA, the 1D, vertical, head-based form

Ksoilð Þ ψ

∂ψ <sup>∂</sup><sup>z</sup> � <sup>1</sup>

� W ¼ 0 (1)

where Cm is the specific moisture capacity, ψ is the soil capillary head (cm), z is the vertical coordinate (cm) (downward positive), τ is time (h), Ksoilð Þ ψ (cm) is the effective hydraulic conductivity, and W is a flux term added for sources and sinks

), such as ET and infiltration. The head-based form is valid in both satu-

Heads for each cell as first computed using an implicit central difference in space and forward difference solution and then flux updating are performed to ensure mass balance. The variables Ksoil and Cm in Eq. (1) are nonlinear and depend on the water content of each cell. Ksoil and Cm are calculated employing Brooks and Corey

While groundwater flow can be free surface or confined, and may be threedimensional, GSSHA solves a 2D lateral solution of the free surface groundwater

> þ ∂ ∂y

Tyy ∂h ∂y 

¼ S

∂Ews

¼ S ∂h

), h is the hydraulic head (m), S is the

<sup>∂</sup><sup>τ</sup> <sup>þ</sup> W xð Þ ; <sup>y</sup>; <sup>τ</sup> (2)

<sup>∂</sup><sup>τ</sup> <sup>þ</sup> W xð Þ ; <sup>y</sup>; <sup>τ</sup> (3)

).

Tyx ∂h ∂x 

�1

storage term (dimensionless), and W is the flux term for sources and sinks (m s�<sup>1</sup>

Downer [20] simplifies Eq. (2) by assuming no directional heterogeneity in the transmissivity terms, expressing transmissivity as the product of the hydraulic conductivity of the media (Ksoil) and the depth of the saturated media (b).

Substituting surface water elevation (Ews = h + datum) for the head, the free surface

Kyyb <sup>∂</sup>Ews ∂y 

GIPL was developed at the Geophysical Institute Permafrost Laboratory for the explicit purpose of determining soil temperature profiles and locating the areas of frozen soil in the profile, including the permafrost. GIPL is a stand-alone 1D soil

In this chapter we describe the development of a coupled soil thermodynamic/ hydrologic model that accounts for the physical processes of thermodynamics, the interaction of the thermodynamics with hydrodynamics, and the variability of freezing condition and runoff in space and time utilizing two well-known and widely applied physics-based models for soil thermodynamics and watershed hydrology, GIPL [10, 11], and GSSHA [12, 13]. This coupled model is demonstrated on a contrived watershed, around a previous GIPL test site [13]. In this demonstration study, the coupled model was deployed in the headwater sub-catchment of the Caribou-Poker Creeks Research Watershed.

## 2. Development of thermo-hydrodynamic model

The coupled model was developed from two widely known, documented, and applied based models, with applications for permafrost and watershed analysis, the Geophysical Institute Permafrost Laboratory (GIPL) model and the Gridded Surface Subsurface Hydrologic Analysis (GSSHA) model. A brief description of the models and the relevant processes are described in the following sections.

#### 2.1 Watershed hydrological model

The GSSHA model is the basis for the overall model development. The GSSHA model was chosen because it is a published, fully distributed, physics-based, watershed model with wide applicability for computing watershed system states, such as soil moisture, snow cover, overland flow depth, stream discharge, and groundwater heads. The GSSHA works on a uniform grid. Any number of point processes, such as distributed precipitation, snow accumulation and melting, evapotranspiration, infiltration, groundwater recharge, etc., can be computed [14–18]. Point processes can be integrated to produce the system response, such as discharge, groundwater flow, etc. GSSHA is an option-driven model, and many processes can be simulated with multiple methods. The focus here is on the subsurface processes.

GSSHA simulates both saturated and unsaturated water movements. In reality, subsurface flow is three-dimensional with complex layering, varying free surface, confined aquifers, and other complications. In GSSHA, the subsurface is represented in a physically explicit manner consistent with the main purpose of GSSHA, which is to simulate surface water flows and the interaction of the surface and subsurface systems. Saturated and unsaturated flows are solved in separate domains. The two domains are linked with head and flux boundaries.

#### 2.1.1 Unsaturated zone model

global thermohaline circulation [5]. While global and regional hydrologic and climate studies require a realistic representation of the hydraulic and thermal properties [6], a physical process-based representation of these properties starts from a small, point, or field scale. Plot-scale studies by Dunne and Black [7] and Stähli et al. [8], among others, show that a catchment's rainfall runoff generation process is altered when infiltration is reduced during the transition of the soil water phase toward the freezing state. As the distribution of frozen soils is highly variable in both space and time [9], to fully understand both the occurrence of frozen soil and the subsequent impact on hydrology requires a physics-based, distributed analytical tool that provides a two-way coupling of the soil thermodynamic and hydrologic state to provide insight into the temporal and spatial variability of soil thermal state and its effects on the time and space distributions of the

In this chapter we describe the development of a coupled soil thermodynamic/ hydrologic model that accounts for the physical processes of thermodynamics, the interaction of the thermodynamics with hydrodynamics, and the variability of freezing condition and runoff in space and time utilizing two well-known and widely applied physics-based models for soil thermodynamics and watershed hydrology, GIPL [10, 11], and GSSHA [12, 13]. This coupled model is demonstrated on a contrived watershed, around a previous GIPL test site [13]. In this demonstration study, the coupled model was deployed in the headwater sub-catchment of the

The coupled model was developed from two widely known, documented, and applied based models, with applications for permafrost and watershed analysis, the Geophysical Institute Permafrost Laboratory (GIPL) model and the Gridded Surface Subsurface Hydrologic Analysis (GSSHA) model. A brief description of the models and the relevant processes are described in the following sections.

The GSSHA model is the basis for the overall model development. The GSSHA

model was chosen because it is a published, fully distributed, physics-based, watershed model with wide applicability for computing watershed system states, such as soil moisture, snow cover, overland flow depth, stream discharge, and groundwater heads. The GSSHA works on a uniform grid. Any number of point processes, such as distributed precipitation, snow accumulation and melting, evapotranspiration, infiltration, groundwater recharge, etc., can be computed [14–18]. Point processes can be integrated to produce the system response, such as discharge, groundwater flow, etc. GSSHA is an option-driven model, and many processes can be simulated with multiple methods. The focus here is on the

GSSHA simulates both saturated and unsaturated water movements. In reality, subsurface flow is three-dimensional with complex layering, varying free surface, confined aquifers, and other complications. In GSSHA, the subsurface is represented in a physically explicit manner consistent with the main purpose of GSSHA, which is to simulate surface water flows and the interaction of the surface and subsurface systems. Saturated and unsaturated flows are solved in separate

domains. The two domains are linked with head and flux boundaries.

hydrological response.

Hydrology - The Science of Water

Caribou-Poker Creeks Research Watershed.

2.1 Watershed hydrological model

subsurface processes.

104

2. Development of thermo-hydrodynamic model

The unsaturated, or vadose, zone controls the flux of water from the land surface to the saturated groundwater zone and determines the partitioning of water into runoff, infiltration, ET, and groundwater recharge. GSSHA provides many ways to analyze these processes including an integrated solution of soil water movement and state described by the Richards Equation [19]. The Richards Equation [19] is considered the most physically correct mathematical representation of the vadose zone. While flow in the vadose zone is in three dimensions, the predominant direction of flow is vertical. In GSSHA, the 1D, vertical, head-based form of the Richards Equation is solved [20, 21]:

$$\mathbf{C}\_{m}(\boldsymbol{\psi})\frac{\partial \boldsymbol{\psi}}{\partial \boldsymbol{\tau}} - \frac{\partial}{\partial \boldsymbol{z}} \left[ K\_{sol}(\boldsymbol{\psi}) \left( \frac{\partial \boldsymbol{\psi}}{\partial \boldsymbol{z}} - \mathbf{1} \right) \right] - \boldsymbol{W} = \mathbf{0} \tag{1}$$

where Cm is the specific moisture capacity, ψ is the soil capillary head (cm), z is the vertical coordinate (cm) (downward positive), τ is time (h), Ksoilð Þ ψ (cm) is the effective hydraulic conductivity, and W is a flux term added for sources and sinks (cm h�<sup>1</sup> ), such as ET and infiltration. The head-based form is valid in both saturated and unsaturated conditions [22].

Heads for each cell as first computed using an implicit central difference in space and forward difference solution and then flux updating are performed to ensure mass balance. The variables Ksoil and Cm in Eq. (1) are nonlinear and depend on the water content of each cell. Ksoil and Cm are calculated employing Brooks and Corey [23] equations, as extended by Huston and Cass [24].

#### 2.1.2 Saturated groundwater model

While groundwater flow can be free surface or confined, and may be threedimensional, GSSHA solves a 2D lateral solution of the free surface groundwater flow equations as described by Pinder and Bredehoeft [25]:

$$\frac{\partial}{\partial \mathbf{x}} \left( T\_{\text{xx}} \frac{\partial \mathbf{h}}{\partial \mathbf{x}} \right) + \frac{\partial}{\partial \mathbf{x}} \left( T\_{\text{xy}} \frac{\partial \mathbf{h}}{\partial \mathbf{y}} \right) + \frac{\partial}{\partial \mathbf{y}} \left( T\_{\text{yx}} \frac{\partial \mathbf{h}}{\partial \mathbf{x}} \right) + \frac{\partial}{\partial \mathbf{y}} \left( T\_{\text{yy}} \frac{\partial \mathbf{h}}{\partial \mathbf{y}} \right) = \mathbf{S} \frac{\partial \mathbf{h}}{\partial \mathbf{r}} + \mathbf{W}(\mathbf{x}, \mathbf{y}, \mathbf{r}) \tag{2}$$

where T is the transmissivity (m<sup>2</sup> s �1 ), h is the hydraulic head (m), S is the storage term (dimensionless), and W is the flux term for sources and sinks (m s�<sup>1</sup> ).

Downer [20] simplifies Eq. (2) by assuming no directional heterogeneity in the transmissivity terms, expressing transmissivity as the product of the hydraulic conductivity of the media (Ksoil) and the depth of the saturated media (b). Substituting surface water elevation (Ews = h + datum) for the head, the free surface problem can be described as [20]

$$\frac{\partial}{\partial \mathbf{x}} \left( K\_{\mathbf{x}\mathbf{x}} b \frac{\partial E\_{w\mathbf{s}}}{\partial \mathbf{x}} \right) + \frac{\partial}{\partial \mathbf{y}} \left( K\_{\mathbf{y}\mathbf{y}} b \frac{\partial E\_{w\mathbf{s}}}{\partial \mathbf{y}} \right) = \mathbf{S} \frac{\partial E\_{w\mathbf{s}}}{\partial \mathbf{z}} + \mathbf{W}(\mathbf{x}, \mathbf{y}, \mathbf{z}) \tag{3}$$

Details of the solution can be found in Downer [20].

#### 2.2 Soil thermodynamic model

GIPL was developed at the Geophysical Institute Permafrost Laboratory for the explicit purpose of determining soil temperature profiles and locating the areas of frozen soil in the profile, including the permafrost. GIPL is a stand-alone 1D soil

thermal and permafrost model that is used to compute a one-dimensional (vertical) soil temperature profile over time using static values of soil moisture.

The Stefan problem [26, 27] with phase change is the problem of thawing or freezing via conduction of heat. GIPL solves the Stefan problem employing the enthalpy formulation. One-dimensional, vertical, quasi-linear heat conduction Equation [28] is the basis of the GIPL numerical model:

$$\frac{\partial H(\mathbf{x},t)}{\partial \boldsymbol{\pi}} = \nabla \cdot \left( k(\mathbf{x},t) \nabla t(\mathbf{x},\boldsymbol{\pi}) \right) \tag{4}$$

2.3 Coupling GIPL to GSSHA

DOI: http://dx.doi.org/10.5772/intechopen.84414

Process Modeling of Soil Thermal and Hydrological Dynamics

included in either.

shown in Figure 2.

Figure 1.

107

GIPL as a frozen ground/permafrost component in GSSHA.

The linkage of the soil thermal and water movement solutions in GSSHA facilitates the temperature domain solution based on the varying soil moisture and the soil water movement domain solutions adjusted for the varying soil temperature condition. In linking the GIPL model to GSSHA, GIPL is essentially included in GSSHA as another point process. The GIPL model is used to compute a vertical soil temperature profile in every lateral two-dimensional GSSHA computational element using the soil moisture information from hydrologic simulations in GSSHA; GSSHA, in turn, uses this temperature and water phase information to adjust hydraulic conductivities for both the vertical unsaturated soil flow and lateral saturated groundwater flow. This two-way coupling increases computational accuracy in both models by providing additional information and processes not previously

As shown in Figure 1, the GSSHA model provides the spatial variability of land surface and hydrodynamic parameters such as air temperature, soil properties, subsurface soil moisture state, etc. to GIPL. GIPL employs the information provided by GSSHA to update the soil thermal state and pass it back to GSSHA. This thermal state information is employed to determine whether the soils are frozen or not. This frozen or unfrozen information is employed by GSSHA to adjust hydraulic conductivities, hydraulic transmissivity, and soil moisture state and saturation levels used in water flow computations. These updated soil saturation and soil moisture state are then employed by the GIPL model to produce new thermal state profiles. This exchange of information continues throughout simulation duration, as

The implementation of GIPL into GSSHA provides a robust soil temperature/ moisture tool with many improvements compared to the stand-alone version of GIPL. In GSSHA, the 1D GIPL model is used to compute a soil temperature profile in every 2D GSSHA grid cell, providing a 3D map of soil temperature/water state. As implemented in GSSHA, GIPL no longer utilizes a static time step nor soil water

where x is a vertical spatial variable which ranges between xu, upper depth of the computational unit, and xL, lower depth of the computational unit. τ is the time and t is the temperature. k(x,t) is a thermal conductivity (W m�<sup>1</sup> K�<sup>1</sup> ). H(x,t) is an enthalpy function defined as

$$H(\mathbf{x}, t) = \int\_0^t \mathbf{C}(\mathbf{x}, \mathbf{s}) d\mathbf{s} + L\theta(\mathbf{x}, t) \tag{5}$$

where L is the volumetric latent heat of freeze/thaw (MJ m�<sup>3</sup> ), C(x,s) is the volumetric heat capacity (MJ m�<sup>3</sup> K�<sup>1</sup> ), and θ(x,t) is the volumetric unfrozen water content (fraction of the total water content). Boundary and initial conditions are required in solving Eq. (4). The boundary condition on the upper extent of the domain corresponds to the near land surface air layer. Embedding seasonal snow layer into the air layer is allowed by the fictitious domain formulation [29]. The upper boundary condition is defined as the Dirichlet-type boundary condition:

$$t(\mathfrak{x}\_{\mathfrak{u}}, \mathfrak{r}) = t\_{\mathfrak{a}\mathfrak{r}} \tag{6}$$

where tair is a daily averaged air temperature. The lower boundary is set as the geothermal gradient as

$$\frac{\partial \mathbf{t}(\mathbf{x}\_l, \mathbf{r})}{\partial \mathbf{x}} = \mathbf{g} \tag{7}$$

where g is the geothermal gradient, a small constant (km�<sup>1</sup> ). For the initial temperature distribution, an appropriate ground temperature profile based on the point location is used:

$$t(\mathfrak{x}, \mathfrak{r}) = t\_0(\mathfrak{x}) \tag{8}$$

The unfrozen water content θ(x,t) is defined as:

$$\theta(\mathbf{x},t) = \eta(\mathbf{x}) \begin{cases} \mathbf{1}, & t \ge t\_\* \\ a|t\_\*-t|^{-b}, & t < t\_\* \end{cases} \tag{9}$$

η(x) is a volumetric soil moisture content. a and b are dimensionless positive constants [30]. The constant t� is a freezing point depression, the temperature at which ice begins to form in the soil. The depth and time variation of the unfrozen water content θð Þ x; t depends on hydrologic forcing and soil type. The numerical solution in GIPL is an implicit, finite difference scheme. A detailed mathematical description of the model and numerical solution methods of Eq. (4) can be found in Marchenko et al. [11], Sergueev et al. [28], and Nicolsky et al. [31]. GIPL input data include soil thermal properties, lithological data and vegetative cover, climate data, and snow cover.
