Process Modeling of Soil Thermal and Hydrological Dynamics

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.

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.

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, 15, 16, 17, 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

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]:

$$C_m\left(\psi \right)\frac{\partial \psi }{\partial \tau }-\frac{\partial }{\partial z}\left[K_{\mathit{\text{soil}}}\left(\psi \right)\left(\frac{\partial \psi }{\partial z}-1\right)\right]-W=0$$ E1

where C_m is the specific moisture capacity, $\psi$ is the soil capillary head (cm), z is the vertical coordinate (cm) (downward positive), $\tau$ is time (h), $K_{\mathit{\text{soil}}}\left(\psi \right)$ (cm) is the effective hydraulic conductivity, and W is a flux term added for sources and sinks (cm h⁻¹), 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 K_soil and C_m in Eq. (1) are nonlinear and depend on the water content of each cell. K_soil and C_m 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 three-dimensional, GSSHA solves a 2D lateral solution of the free surface groundwater flow equations as described by Pinder and Bredehoeft [25]:

$$\frac{\partial }{\partial x}\left(T_{\mathit{xx}}\frac{\partial h}{\partial x}\right)+\frac{\partial }{\partial x}\left(T_{\mathit{xy}}\frac{\partial h}{\partial y}\right)+\frac{\partial }{\partial y}\left(T_{\mathit{yx}}\frac{\partial h}{\partial x}\right)+\frac{\partial }{\partial y}\left(T_{\mathit{yy}}\frac{\partial h}{\partial y}\right)=S\frac{\partial h}{\partial \tau }+W\left(x,y,\tau \right)$$ E2

where T is the transmissivity (m² s⁻¹), h is the hydraulic head (m), S is the storage term (dimensionless), and W is the flux term for sources and sinks (m s⁻¹).

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 (K_soil) and the depth of the saturated media (b). Substituting surface water elevation (E_ws = h + datum) for the head, the free surface problem can be described as [20]

$$\frac{\partial }{\partial x}\left(K_{\mathit{xx}}b\frac{\partial E_{\mathit{ws}}}{\partial x}\right)+\frac{\partial }{\partial y}\left(K_{\mathit{yy}}b\frac{\partial E_{\mathit{ws}}}{\partial y}\right)=S\frac{\partial E_{\mathit{ws}}}{\partial \tau }+W\left(x,y,\tau \right)$$ E3

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\left(x,t\right)}{\partial \tau }=\nabla ·\left(k\left(x,t\right)\nabla t\left(x,\tau \right)\right)$$ E4

where x is a vertical spatial variable which ranges between x_u, upper depth of the computational unit, and x_L, lower depth of the computational unit. τ is the time and t is the temperature. k(x,t) is a thermal conductivity (W m⁻¹ K⁻¹). H(x,t) is an enthalpy function defined as

$$H\left(x,t\right)={\int }_0^{t}C\left(x,s\right)\mathit{ds}+\mathit{L\theta }\left(x,t\right)$$ E5

where L is the volumetric latent heat of freeze/thaw (MJ m⁻³), C(x,s) is the volumetric heat capacity (MJ m⁻³ K⁻¹), and $\theta$ (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\left(x_u,\tau \right)=t_{\mathit{air}}$$ E6

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

$$\frac{\partial t\left(x_l,\tau \right)}{\partial x}=g$$ E7

where g is the geothermal gradient, a small constant (km⁻¹). For the initial temperature distribution, an appropriate ground temperature profile based on the point location is used:

$$t\left(x,\tau \right)=t_0\left(x\right)$$ E8

The unfrozen water content $\theta$ (x,t) is defined as:

$$\theta \left(x,t\right)=\eta \left(x\right)\left\{\begin{array}{cc}1,& t\ge t_{*}\\ a{\left|t_{*}-t\right|}^{-b},& t<t_{*}\end{array}\right.$$ E9

$\eta$ (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 $\theta \left(x,t\right)$ 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.

2.3. Coupling GIPL to 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 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.

Figure 1.

GIPL as a frozen ground/permafrost component in GSSHA.

Figure 2.

GSSHA and GIPL coupling architecture.

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 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 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.

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 (S_E), a nondimensional parameter that varies between 0 and 1.

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

$$S_E={\left(\frac{1}{1+{\left(\propto \left|1.22.t\right|\right)}^n}\right)}^m\mathit{\text{for}}T\le 0^{°}C$$ E10

where n, m, and α are the van Genuchten parameters; t is the soil temperature in °C. S_E is always 1 for temperatures above 0°C. For temperatures below −10°C, the value of S_E 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, S_E, 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_{\mathit{\text{soil}}}\left(t\right)=e^{S_E}$$ E11

where K_soil(t) is the effective hydraulic conductivity in m s⁻¹, K_t is the hydraulic conductivity for S_E = 1, and K_f is the frozen hydraulic conductivity (S_E = 0). The frozen soil has a very little capacity to transmit water. In GSSHA K_f is set to 10⁻⁶.

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 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 (b_effective). 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 (b_effective). 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.

$$\mathit{T}={{\mathit{K}}_{\text{soil}}}^{\ast }{\mathit{b}}_{\text{effective}}$$ E12

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′ N and longitude of 147°30′ W Alaska. The CPCRW, which encompasses an area of 101.5 km² within the boreal forest, is part of the Long-Term Ecological Research (LTER) network. Parts of this watershed are underlain by discontinuous permafrost [34].

Figure 3.

Study area.

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://nationalmap.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).

Figure 4.

Soil types and vegetation types of the study area.

Soil series	USDA texture	Drainage
Fairplay	Silt loam and gravelly silt loam	Moderately well drained
Olnes	Silt loam and very gravelly silt loam	Well drained

Table 1.

Soil properties of the study area.

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 August [38]. Annual snowfall averages 1692 mm and commonly covers the ground from October to April [37]. The model employed CPEAK hydrometeorological data of the year 2002 included relative humidity (%), wind speed (kts), air temperature (°F), barometric pressure (in Hg), solar radiation (W h m⁻²), and tipping bucket rainfall rates and are available at http://www.lter.uaf.edu/data.

4. Model development of the study area

A 10 by 10 gridded model of the study area, shown in Figure 3, includes precipitation, overland flow, unsaturated zone computations using the 1D Richards equation, ET using the Penman-Monteith Equation [39, 40], and 1D soil thermal computations using GIPL.

4.1. Initial soil moisture and soil temperature conditions

The initial temperature condition which is from CPEAK ground temperature thermistors at several depths on 2002-5-1, 1 AM, is shown in Figure 5. The numerical solution of the soil thermal state using quasi-linear heat conduction equation (Eq. (4)) employs this initial temperature condition of the soil profile.

Figure 5.

Initial soil temperature condition.

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)).

Figure 6.

Initial soil moisture condition.

4.2. Model parameter values

The model parameter values, distributed on grids horizontally and on soil layers vertically, for processes, such as overland flow, infiltration, evapotranspiration, and thermodynamics, are assigned based on the soil texture, shown in Figure 4a, and land use, shown in Figure 4b.

Table 2 shows the soil thermal parameters with values, representing the silt loam [31, 42, 43] of Alaskan woodland and tundra ecosystem sites in permafrost-active region, employed in the thermodynamics process.

Soil heat conductivity (W/m K)	Soil volumetric heat capacity (J/m^3 K)
3.0	2,800,000

Table 2.

Soil thermal parameter values.

Table 3 shows Manning’s roughness parameter values, representing the vegetation type of Figure 4(b) [44, 45], employed in the routing process.

Land use type	Manning’s roughness values
Coniferous open	0.17
Coniferous woodland	0.19
Deciduous closed	0.2
Shrub tall	0.25

Table 3.

Manning’s roughness values.

Table 4 shows the soil hydraulic parameters with values, representing the silt loam ([46], https://www.gsshawiki.com) of Figure 4(a), employed in the infiltration and soil water retention process.

Infiltration parameter	Soil layer
	Top	Middle	Lower
Soil thickness (cm)	20	50	—
Saturated hydraulic conductivity (cm/h)	0.5	0.5	0.5
Pore size distribution index	0.60	0.694	0.694
Wetting front suction head (cm)	8.0	6.0	6.0
Porosity (m^3/m^3)	0.42	0.40	0.40
Residual soil moisture content (m^3/m^3)	0.04	0.045	0.045

Table 4.

Soil parameter values for the Richards infiltration scheme.

In this study, Penman-Monteith method was employed as evapotranspiration process. Table 5 shows the Penman-Monteith parameters that includes vegetation transmission coefficient (light penetration through canopy), values of land surface albedo, vegetation height (for aerodynamic resistance term), and vegetation canopy resistance (for stomatal control of the loss of water). Table 5 shows the ET parameter values, representing the vegetation type of Figure 4(b), employed in the evapotranspiration process.

Land use type	Albedo	Vegetation height (cm)	Canopy resistance (s/m)	Transmission coefficient
Coniferous open	0.15	10	120	0.18
Coniferous woodland	0.15	10	120	0.18
Deciduous closed	0.2	12	200	0.18
Shrub tall	0.2	1.3	150	0.5

Table 5.

Evapotranspiration parameter values.

All the parameter values defined in Table 5 were taken from the literature, which are also defined in the GSSHA wiki https://www.gsshawiki.com. The literature values for albedo and vegetation height are defined in Eagleson [47]. The literature values for canopy resistance are defined in Monteith [39]. The literature values for the transmission coefficient are defined in Sutton [48].

5. Result and discussion

5.1. Soil thermodynamics

The simulation period was from May 1st to May 31st of 2002, a period during which air temperatures are beginning to rise above freezing. The daily maximum soil temperature obtained from the simulation is compared with the observed one in Figure 7. Both the observed and simulated temperatures in Figure 7 is effective at a depth of 10 cm in the soil profile. The root-mean-square error for this daily maximum soil temperature was 4.7°C. It was found that the soil thermal conductivity parameter value was most sensitive and effective for the near-surface simulated temperature. It was also found that the temperature of the near-surface soil layer is directly influenced by the air temperature.

Figure 7.

Comparison of the time series of observed and simulated temperature.

5.2. Soil’s effective hydraulic conductivity

Figure 8 shows the simulation evolution of the effective hydraulic conductivity starting on May 1st of 2002. There are no observed hydraulic conductivity values to compare to, but the model intuitively represents the condition. Initially, the hydraulic conductivity in Figure 8 is very low, almost zero in the first 113 hours, the result of reduced effective saturation in the soil due to the soil being completely frozen (Figure 7). As the air and soil temperatures rise, the soil begins to thaw, with a resulting increase in hydraulic conductivity. The exponential rise of hydraulic conductivity is a consistent observation from freezing/thawing mineral and organic soils [33]. The simulated effective hydraulic conductivity in Figure 8 is from Eq. (11), where the effective soil hydraulic conductivity K_soil at a given temperature (t) is a function of hydraulic conductivity of the unfrozen soil and the effective saturation S_E.

Figure 8.

Hydraulic conductivity under freezing and thawing soil active layer.

5.3. Hydrologic runoff response

Figure 9 shows the comparison of GSSHA simulated discharge with and without taking into account freezing and thawing soil properties in the study area. While there are no measurements of runoff from the site, the results are consistent with the results presented for air and soil temperature and resulting hydraulic conductivity shown above. The freezing soil temperature (Figure 7) leads to increased coverage of frozen soil which in turn leads to less soil pore water storage. The reduced soil pore water storage capacity leads to decrease hydraulic conductivity as shown in Figure 8, which results in a flashier response to the precipitation event as shown in Figure 9 graph with thermodynamics. On the other hand, loss of frozen soil and permafrost or without taking into account thermodynamics will lead to enhanced connectivity between the surface and groundwater storage regimes and decreased overland flow as shown by the graph without thermodynamics in Figure 9.

Figure 9.

Hydrograph with and without thermodynamics.

5.4. Discussion

Unless the simulation discharge represents the distributed physical runoff process in a realistic way, even a well-calibrated simulated discharge at a catchment outlet may only be a right answer for wrong reasons. The modeling and the simulations in this study have explicitly taken into account frozen soil as an important hydrologic factor. The simulation results showed the variability of freezing condition in space and time. How this variability of freezing condition in space and time would affect the distribution of overland runoff is particularly important from the concern of climate change and land use change effect in local hydrology [7, 49]. Figure 10 shows the model simulated distribution of spatially distributed overland runoff, at the peak of the precipitation event of May 7th and May 8th of 2002 in CPEAK station, by taking into account thermodynamics and without taking into account thermodynamics. It is clear from Figure 10 that the simulation significantly underestimated the runoff when thermodynamic process was absent.

Figure 10.

Comparison of distribution of overland runoff with and with thermodynamic process in the model simulation.

6. Recommendations and conclusions

A coupled framework for simulating the interaction between soil temperature, including permafrost, and hydrology was developed by incorporating the soil temperature and permafrost model GIPL into the distributed, physics-based hydrologic model GSSHA. This chapter describes the architecture for numerically linking the GIPL thermodynamic model into GSSHA’s hydrodynamic modeling framework. Deploying this enhanced capability showed that GSSHA hydrodynamics include soil moisture saturation feedback in the vadose zone as well as the corresponding soil ice content effects on hydraulic conductivity and transmissivity. In this study the coupled model was deployed in the headwater region at the peak of the Caribou-Poker Creeks Research Watershed that included the hydrometeorological and soil property measurement station called CPEAK.

The model captured the seasonal rise of soil temperatures and thaw of frozen soils. The model showed intuitively correct representations of soil hydraulic conductivity and runoff, consistent with the observed rise of soil temperatures. Numerical simulations showed the hydrologic importance of frozen soils with the implication that climate change could have large effects on hydrology as air and soil temperatures rise near the poles resulting in loss of permafrost and increase seasonal thawing of soils.

The main use of the “process modeling of soil thermal and hydrological dynamics” has been to generate spatial and temporal dataset of permafrost distribution and ground temperature dynamics as well as the active layer thickness, the depth of the seasonal thaw. Such dataset would be useful in the assessments of a wide range of thermo-hydrodynamic related fields, including ecology, climatology, and socio-economy, in the cold regions.

Acknowledgments

Pradhan et al. [13] was the initial hydrodynamic and thermodynamic model integration effort which was then funded by Strategic Environmental Research and Development Program (SERDP) under Project Number 11 RC-2110. Theoretical and technical contributions of Dr. Anna Liljedahl, at the University of Alaska Fairbanks are significant for that initial hydrodynamic and thermodynamic model integration effort. Article Processing Charges (APCs) for the publication were supported by the U.S. Army Corps of Engineers Flood and Coastal Storm Damage Reduction Research and Development Program under the technical direction of Dr. Julie Rosati and Mr. Ian E. Floyd. The Bonanza Creek Long-Term Ecological Research program (BNZ LTER) at the University of Alaska Fairbanks is thanked for making and sharing the hydrology and climate data set of the CPCRW in Alaska. Constructive comments from anonymous reviewers are greatly appreciated.