2.3 Coupling GIPL to GSSHA

thermal and permafrost model that is used to compute a one-dimensional (vertical)

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

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

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

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

<sup>∂</sup>t xl ð Þ ; <sup>τ</sup>

temperature distribution, an appropriate ground temperature profile based on the

�

where g is the geothermal gradient, a small constant (km�<sup>1</sup>

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

θð Þ¼ x; t ηð Þ x

<sup>∂</sup><sup>τ</sup> <sup>¼</sup> <sup>∇</sup> � ð Þ k xð Þ ; <sup>t</sup> <sup>∇</sup>t xð Þ ; <sup>τ</sup> (4)

C xð Þ ; s ds þ Lθð Þ x; t (5)

), and θ(x,t) is the volumetric unfrozen

t xð Þ¼ <sup>u</sup>; τ tair (6)

<sup>∂</sup><sup>x</sup> <sup>¼</sup> <sup>g</sup> (7)

t xð Þ¼ ; τ t0ð Þ x (8)

1, t≥t�

, t < t�

a tj j � � <sup>t</sup> �<sup>b</sup>

η(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,

). H(x,t) is an

), C(x,s) is the

). For the initial

(9)

soil temperature profile over time using static values of soil moisture.

Equation [28] is the basis of the GIPL numerical model:

enthalpy function defined as

Hydrology - The Science of Water

condition:

geothermal gradient as

point location is used:

and snow cover.

106

volumetric heat capacity (MJ m�<sup>3</sup> K�<sup>1</sup>

<sup>∂</sup>H xð Þ ; <sup>t</sup>

t is the temperature. k(x,t) is a thermal conductivity (W m�<sup>1</sup> K�<sup>1</sup>

ðt 0

H xð Þ¼ ; t

where L is the volumetric latent heat of freeze/thaw (MJ m�<sup>3</sup>

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 included in either.

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 shown in Figure 2.

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

#### Figure 1.

GIPL as a frozen ground/permafrost component in GSSHA.

In GSSHA, only the surficial aquifer is simulated, so the saturated groundwater domain is down to the first confining layer in the subsurface. This is typically on the order of a few to hundreds of meters deep. The unsaturated zone domain is any soil above the saturated zone. The unsaturated domain is dynamic in both space and time and can vary from no domain (groundwater table is at or above the soil surface) to the depth of the surficial aquifer, depending on groundwater conditions. The unsaturated zone is further divided into four regions, corresponding to the

Because of the differences in domains, and requirements for solution, in the coupled framework, the GIPL domain and discretization are independent of the saturated and unsaturated soil water domains and discretizations. The linkage of computational nodal discretized information from GIPL to GSSHA and vice-versa is

The primary effect of soil temperature on the soil water domain is that freezing soils result in much lower hydraulic conductivities of the soil. In the unsaturated zone, the vertical soil hydraulic conductivity is computed from the relative satura-

The relative fraction of liquid water of the total soil moisture, SE, can be com-

where n, m, and α are the van Genuchten parameters; t is the soil temperature in °C. SE is always 1 for temperatures above 0°C. For temperatures below �10°C, the

As soil freezes, pathways through the soil, pores, close reducing the ability of the soil to transmit water. This results in a reduction of the soil hydraulic conductivity. In the unfrozen portion of the soil, an exponential response in effective hydraulic conductivity has been measured for freezing/thawing mineral and organic soils [33]. The temperature-adjusted relative saturation, SE, can be used to compute the soil temperature-adjusted hydraulic conductivity as the sum of the hydraulic conductivity of the unfrozen pores and the hydraulic conductivity of the frozen pores:

for T ≤0°

SEln Kð Þþ<sup>t</sup> ð Þ <sup>1</sup>�SE Kf (11)

C (10)

, Kt is the hydraulic

.

<sup>1</sup> <sup>þ</sup> ð Þ <sup>∝</sup>j j <sup>1</sup>:22:<sup>t</sup> <sup>n</sup> <sup>m</sup>

Ksoil ðÞ¼ t e

conductivity for SE = 1, and Kf is the frozen hydraulic conductivity (SE = 0). The frozen soil has a very little capacity to transmit water. In GSSHA Kf is set to 10�<sup>6</sup>

The effect of soil temperature on the saturated groundwater solution is considered by adjusting the transmissivity of the groundwater based on the frozen water content in each groundwater cell. Because of soil structure, lateral hydraulic conductivities are essentially 10 to 100 times larger than vertical hydraulic conductivities. In practical terms, frozen saturated media has little to no ability to laterally

In the soil thermal calculations, the full profile of soil temperatures allows the

determination of discrete soil thermal cells that are below freezing and thus

where Ksoil(t) is the effective hydraulic conductivity in m s�<sup>1</sup>

2.3.3 Soil heat transfer effect on effective groundwater transmissivity

transport water compared to the unfrozen media.

109

A, B, and C soil horizons, as well as the deeper groundwater media.

tion (SE), a nondimensional parameter that varies between 0 and 1.

SE <sup>¼</sup> <sup>1</sup>

2.3.2 Soil temperature effect on hydraulic conductivity

Process Modeling of Soil Thermal and Hydrological Dynamics

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

puted from the soil temperature as [32]

value of SE is assumed to be 0.

shown in Figure 2.

Figure 2. GSSHA and GIPL coupling architecture.

state. Instead, GIPL performs soil temperature computations utilizing the time and space varying soil moistures computed using the Richards equation. The time step of GIPL is set to the overall GSSHA model time step, which is on the order of minutes, as opposed to days. The GSSHA input formats are used to specify 3D distributed soil parameter values for the GIPL solution. Several thermohydrodynamic formulations and modeling concepts are implemented to link and exchange the information in GIPL and GSSHA.

#### 2.3.1 Linking soil temperature and soil water computational nodes

The solution domain of the GIPL soil temperature model overlaps in a somewhat complex manner with both the saturated and unsaturated soil water movement domains in GSSHA. In no-flux lower GIPL boundary condition, the GIPL domain must extend very deep into the soil/permafrost, as much as 1000 m or more.

Process Modeling of Soil Thermal and Hydrological Dynamics DOI: http://dx.doi.org/10.5772/intechopen.84414

In GSSHA, only the surficial aquifer is simulated, so the saturated groundwater domain is down to the first confining layer in the subsurface. This is typically on the order of a few to hundreds of meters deep. The unsaturated zone domain is any soil above the saturated zone. The unsaturated domain is dynamic in both space and time and can vary from no domain (groundwater table is at or above the soil surface) to the depth of the surficial aquifer, depending on groundwater conditions. The unsaturated zone is further divided into four regions, corresponding to the A, B, and C soil horizons, as well as the deeper groundwater media.

Because of the differences in domains, and requirements for solution, in the coupled framework, the GIPL domain and discretization are independent of the saturated and unsaturated soil water domains and discretizations. The linkage of computational nodal discretized information from GIPL to GSSHA and vice-versa is shown in Figure 2.

#### 2.3.2 Soil temperature effect on hydraulic conductivity

The primary effect of soil temperature on the soil water domain is that freezing soils result in much lower hydraulic conductivities of the soil. In the unsaturated zone, the vertical soil hydraulic conductivity is computed from the relative saturation (SE), a nondimensional parameter that varies between 0 and 1.

The relative fraction of liquid water of the total soil moisture, SE, can be computed from the soil temperature as [32]

$$\mathcal{S}\_E = \left(\frac{1}{1 + (\infty|1.22.t|)^n}\right)^m \text{ for } T \le \mathbf{0}^\circ \text{C} \tag{10}$$

where n, m, and α are the van Genuchten parameters; t is the soil temperature in °C. SE is always 1 for temperatures above 0°C. For temperatures below �10°C, the value of SE is assumed to be 0.

As soil freezes, pathways through the soil, pores, close reducing the ability of the soil to transmit water. This results in a reduction of the soil hydraulic conductivity. In the unfrozen portion of the soil, an exponential response in effective hydraulic conductivity has been measured for freezing/thawing mineral and organic soils [33]. The temperature-adjusted relative saturation, SE, can be used to compute the soil temperature-adjusted hydraulic conductivity as the sum of the hydraulic conductivity of the unfrozen pores and the hydraulic conductivity of the frozen pores:

$$K\_{soil}\left(t\right) = e^{S\_E \ln(K\_t) + (1 - S\_E)K\_f} \tag{11}$$

where Ksoil(t) is the effective hydraulic conductivity in m s�<sup>1</sup> , Kt is the hydraulic conductivity for SE = 1, and Kf is the frozen hydraulic conductivity (SE = 0). The frozen soil has a very little capacity to transmit water. In GSSHA Kf is set to 10�<sup>6</sup> .

#### 2.3.3 Soil heat transfer effect on effective groundwater transmissivity

The effect of soil temperature on the saturated groundwater solution is considered by adjusting the transmissivity of the groundwater based on the frozen water content in each groundwater cell. Because of soil structure, lateral hydraulic conductivities are essentially 10 to 100 times larger than vertical hydraulic conductivities. In practical terms, frozen saturated media has little to no ability to laterally transport water compared to the unfrozen media.

In the soil thermal calculations, the full profile of soil temperatures allows the determination of discrete soil thermal cells that are below freezing and thus

state. Instead, GIPL performs soil temperature computations utilizing the time and space varying soil moistures computed using the Richards equation. The time step of GIPL is set to the overall GSSHA model time step, which is on the order of minutes, as opposed to days. The GSSHA input formats are used to specify 3D distributed soil parameter values for the GIPL solution. Several thermo-

hydrodynamic formulations and modeling concepts are implemented to link and

The solution domain of the GIPL soil temperature model overlaps in a somewhat complex manner with both the saturated and unsaturated soil water movement domains in GSSHA. In no-flux lower GIPL boundary condition, the GIPL domain must extend very deep into the soil/permafrost, as much as 1000 m or more.

exchange the information in GIPL and GSSHA.

Figure 2.

108

GSSHA and GIPL coupling architecture.

Hydrology - The Science of Water

2.3.1 Linking soil temperature and soil water computational nodes

incapable of transmitting saturated soil water laterally. However, in the 2D saturated flow calculations, there is no vertical discretization of the domain. To account for the effect of the frozen sections in the groundwater media, the depth of flow (b) [Eq. (3)] is adjusted by subtracting the total frozen portions of the unsaturated zone, contiguous or not, to compute an effective flow depth (beffective). The effect of this reduced the transmissivity in Eq. (2), which is the primary control on flow in the equation.

The depth of the unfrozen saturated media in GSSHA is identified by determining which soil thermal computational nodes correspond to the saturated media depth and then adding up all unfrozen,T > 0, cells to determine the unfrozen saturated flow depth (beffective). If the groundwater surface is between frozen and unfrozen thermal cells, then the portion unfrozen is determined by interpolation. This avoids the overestimation of the effective saturated depth.

Once the effective saturated depth is calculated, local�/grid-based groundwater transmissivity is defined as.

$$T = K\_{soil} \, ^\* \, b\_{\text{effective}} \tag{12}$$

tionalmap.gov/viewer.html, was employed to develop the study area model shown

Process Modeling of Soil Thermal and Hydrological Dynamics

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

Figure 4(a) shows the soils of the study area, and the soil description, as per Rieger et al. [35], is shown in Table 1. The study area vegetation-type distribution, based on a unified statewide system for classifying vegetation in Alaska [36], is

CPCRW lies in the interior climatic division of Alaska characterized by low annual precipitation, low cloudiness, low humidity, and large diurnal and annual temperature variations. In the study region, the 30-year normals [37] show that January is typically the coldest month, with a mean temperature of 24.4°C, and July, the warmest month, with a mean temperature of 17.1°C. The average annual precipitation in this region is 285 mm, and the wettest months are June, July, and

in Figure 3.

111

Figure 3. Study area.

shown in Figure 4(b).

3.2.2 Climate of the study area

#### 3. Study area and data

#### 3.1 Description of the study site

Figure 3(c) shows the location of the study area with elevation and gridded model domain. The study area is in the Caribou-Poker Creeks Research Watershed (CPCRW), 48 km north of Fairbanks with latitude of 65°10<sup>0</sup> N and longitude of 147° 30<sup>0</sup> W Alaska. The CPCRW, which encompasses an area of 101.5 km<sup>2</sup> within the boreal forest, is part of the Long-Term Ecological Research (LTER) network. Parts of this watershed are underlain by discontinuous permafrost [34].

There is a weather station on the summit of Caribou Peak as shown in Figure 3 which is called CPEAK and is at an elevation of 773 m. This station has a 10-meter tower with various atmospheric sensors and ground temperature thermistors at several depths (http://www.lter.uaf.edu/data/metadata/id/442/inline/1).

As shown in Figure 3, A 10 by 10 GSSHA/GIPL gridded model was developed for simulations in the study area. This study area model domain and location is selected as to include the CPEAK station so that the observations from the station could be used in model development and validation.

#### 3.2 Data types and sources

The Caribou-Poker Creeks Research Watershed (CPCRW) is one of the study sites for the Long-Term Ecological Research (LTER) project where long-term observed data are maintained by the Institute of Arctic Biology at the University of Alaska, Fairbanks, and made available online (http://www.lter.uaf.edu/data). The website has a search filter page where the hydrology and climate data for the study sites is available. The data from CPEAK was employed as initial condition of soil profile temperature and soil moisture and hydrometeorological forcings/input to run the model.

#### 3.2.1 Topography, soil, and land use

A 50-meter digital elevation model (DEM), obtained from the National Elevation Dataset (NED) and downloaded through the National Map Viewer http://na

Process Modeling of Soil Thermal and Hydrological Dynamics DOI: http://dx.doi.org/10.5772/intechopen.84414

incapable of transmitting saturated soil water laterally. However, in the 2D saturated flow calculations, there is no vertical discretization of the domain. To account for the effect of the frozen sections in the groundwater media, the depth of flow (b) [Eq. (3)] is adjusted by subtracting the total frozen portions of the unsaturated zone, contiguous or not, to compute an effective flow depth (beffective). The effect of this reduced the transmissivity in Eq. (2), which is the primary control on flow in

The depth of the unfrozen saturated media in GSSHA is identified by determin-

Once the effective saturated depth is calculated, local�/grid-based groundwater

Figure 3(c) shows the location of the study area with elevation and gridded model domain. The study area is in the Caribou-Poker Creeks Research Watershed (CPCRW), 48 km north of Fairbanks with latitude of 65°10<sup>0</sup> N and longitude of 147° 30<sup>0</sup> W Alaska. The CPCRW, which encompasses an area of 101.5 km<sup>2</sup> within the boreal forest, is part of the Long-Term Ecological Research (LTER) network. Parts

There is a weather station on the summit of Caribou Peak as shown in Figure 3 which is called CPEAK and is at an elevation of 773 m. This station has a 10-meter tower with various atmospheric sensors and ground temperature thermistors at several depths (http://www.lter.uaf.edu/data/metadata/id/442/inline/1).

As shown in Figure 3, A 10 by 10 GSSHA/GIPL gridded model was developed for simulations in the study area. This study area model domain and location is selected as to include the CPEAK station so that the observations from the station

The Caribou-Poker Creeks Research Watershed (CPCRW) is one of the study sites for the Long-Term Ecological Research (LTER) project where long-term observed data are maintained by the Institute of Arctic Biology at the University of Alaska, Fairbanks, and made available online (http://www.lter.uaf.edu/data). The website has a search filter page where the hydrology and climate data for the study sites is available. The data from CPEAK was employed as initial condition of soil profile temperature and soil moisture and hydrometeorological forcings/input to

A 50-meter digital elevation model (DEM), obtained from the National Elevation Dataset (NED) and downloaded through the National Map Viewer http://na

<sup>T</sup> <sup>¼</sup> Ksoil <sup>∗</sup> beffective (12)

ing which soil thermal computational nodes correspond to the saturated media depth and then adding up all unfrozen,T > 0, cells to determine the unfrozen saturated flow depth (beffective). If the groundwater surface is between frozen and unfrozen thermal cells, then the portion unfrozen is determined by interpolation.

This avoids the overestimation of the effective saturated depth.

of this watershed are underlain by discontinuous permafrost [34].

could be used in model development and validation.

the equation.

transmissivity is defined as.

Hydrology - The Science of Water

3. Study area and data

3.2 Data types and sources

3.2.1 Topography, soil, and land use

run the model.

110

3.1 Description of the study site

tionalmap.gov/viewer.html, was employed to develop the study area model shown in Figure 3.

Figure 4(a) shows the soils of the study area, and the soil description, as per Rieger et al. [35], is shown in Table 1. The study area vegetation-type distribution, based on a unified statewide system for classifying vegetation in Alaska [36], is shown in Figure 4(b).

#### 3.2.2 Climate of the study area

CPCRW lies in the interior climatic division of Alaska characterized by low annual precipitation, low cloudiness, low humidity, and large diurnal and annual temperature variations. In the study region, the 30-year normals [37] show that January is typically the coldest month, with a mean temperature of 24.4°C, and July, the warmest month, with a mean temperature of 17.1°C. The average annual precipitation in this region is 285 mm, and the wettest months are June, July, and

The initial soil moisture condition obtained from Caribou-Poker Creeks Research Watershed for 2002-5-1, 1 AM [41], is shown in Figure 6. This initial soil moisture condition of the soil profile is employed by the numerical solution of the soil moisture, Richards equation, in the unsaturated vadose zone (Eq. (1)).

Process Modeling of Soil Thermal and Hydrological Dynamics

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

The model parameter values, distributed on grids horizontally and on soil layers vertically, for processes, such as overland flow, infiltration, evapotranspiration, and

Soil heat conductivity (W/m K) Soil volumetric heat capacity (J/m3 K)

3.0 2,800,000

4.2 Model parameter values

Figure 5.

Figure 6.

Table 2.

113

Initial soil moisture condition.

Soil thermal parameter values.

Initial soil temperature condition.
