Next Article in Journal
Impact of Main Pipe Flow Velocity on Leakage and Intrusion Flow: An Experimental Study
Next Article in Special Issue
Bayesian Hierarchical Model Uncertainty Quantification for Future Hydroclimate Projections in Southern Hills-Gulf Region, USA
Previous Article in Journal
Sediment Grain-Size Characteristics and its Sources of Ten Wind-Water Coupled Erosion Tributaries (the Ten Kongduis) in the Upper Yellow River
Previous Article in Special Issue
Flash Flood Simulation for Ungauged Catchments Based on the Distributed Hydrological Model
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

Catchment Hydrological Modeling with Soil Thermal Dynamics during Seasonal Freeze-Thaw Cycles

1
Coastal and Hydraulics Laboratory, U.S. Army Engineer Research and Development Center, 3909 Halls Ferry Road, Vicksburg, MS 39180-6199, USA
2
Geophysical Institute, University of Alaska Fairbanks, Fairbanks, AK 99775, USA
*
Author to whom correspondence should be addressed.
Water 2019, 11(1), 116; https://doi.org/10.3390/w11010116
Submission received: 16 November 2018 / Revised: 29 December 2018 / Accepted: 2 January 2019 / Published: 10 January 2019
(This article belongs to the Special Issue Catchment Modelling)

Abstract

:
To account for the seasonal changes in the soil thermal and hydrological dynamics, the soil moisture state physical process defined by the Richards Equation is integrated with the soil thermal state defined by the numerical model of phase change based on the quasi-linear heat conductive equation. The numerical model of phase change is used to compute a vertical soil temperature profile using the soil moisture information from the Richards solver; the soil moisture numerical model, in turn, uses this temperature and phase, information to update hydraulic conductivities in the vertical soil moisture profile. Long-term simulation results from the test case, a head water sub-catchment at the peak of the Caribou Poker Creek Research Watershed, representing the Alaskan permafrost active region, indicated that freezing temperatures decreases infiltration, increases overland flow and peak discharges by increasing the soil ice content and decaying the soil hydraulic conductivity exponentially. Available observed and the simulated soil temperature comparison analysis showed that the root mean square error for the daily maximum soil temperature at 10-cm depth was 4.7 °C, and that for the hourly soil temperature at 90-cm and 300-cm was 0.17 °C and 0.14 °C, respectively.

1. Introduction

Both the Intergovernmental Panel on Climate Change latest full assessment [1] and the U.S. National Climate Assessment [2] unequivocally state that climate change is occurring due to anthropogenic factors, that significant warming has already occurred, and that both data trends and climate modeling indicate that climate change and the problems associated with climate change will only worsen in the future. Of the many issues related to climate change, climate change directly affects soil temperature and trends in soil temperature. After analyzing 50 years of soil temperature profile data at 5, 10, 20, 50, 100, and 150 cm depths across Canada, Qian et al. [3] found a soil temperature warming trend of +0.30 °C/decade. Warming was especially prevalent during spring. Qian et al. [3] also found that the increasing trend of soil temperature was directly associated with the increasing trends in air temperature. A change in soil temperature has a determinative role on the response of biochemical processes that control the soil biological cycle of plant nutrient and carbon use and production [4]. A positive trend in the soil temperature has a determinative role on the emission of greenhouse gases in Alaska and other arctic and sub-arctic ecosystems [5,6]. Melvin et al. [7] showed widespread damages to Alaska public infrastructure from increased soil temperature due to climate change.
Seasonal changes in air temperature lead to analogous changes in the soil thermal regime, thereby affecting the hydrological response [8,9]. Plot-scale studies by Dunne and Black [9], and Stähli et al. [10] showed that storm runoff generation processes of a watershed are altered as the soil water phase transitions from a non-frozen state to a freezing state with reduced infiltration and enhanced runoff. Without physical processes that explicitly simulate such seasonal changes in hydrological regime, calibrated/effective parameter values alone may not be sufficient to simulate the rainfall-runoff response, especially in a long-term numerical simulation of projected future global warming scenarios in higher latitude regions [8]. It is also significant to note that the freezing condition is highly variable in space and time [11], which a simple conceptual model, or even a physics-based numerical model lacking the required process, may fail to account for. In order to fully analyze the effects of climate change on soil thermal dynamics, and associated hydrologic response for the future climate change scenarios on the response of watersheds with permafrost and/or a seasonally frozen soil regime, understanding and application of soil thermal dynamics interaction with soil hydrological dynamics are required.
The purpose of this paper is to describe the development of a distributed, numerical model that accounts for both physics of soil thermodynamics and hydrodynamics including the physical linkage between thermodynamics and hydrodynamics of soil horizons/layers and the atmosphere. This numerical model could then be employed to explicitly simulate the soil thermal state and its effects on hydrological dynamics during changing seasons. The architecture of the numerical models’ linkage development is achieved by two way coupling of the soil moisture physical state accounting of Richards mathematical formulation in the US Army’s physics-based distributed hydrology model, Gridded Surface Subsurface Hydrologic Analysis (GSSHA) model [12,13] and soil thermal regime simulation capability of the Geophysical Institute Permafrost Laboratory (GIPL) numerical model [14,15] of the 1D quasi-linear heat equation with phase change between liquid water and ice. The coupled numerical schemes provide a clear insight of the temporal and spatial variability of the soil thermal state during changing seasons and the effects on hydrological processes that are significant for long-term simulation, such as infiltration and evapotranspiration.
The coupled model was deployed in the headwater sub-catchment, 0.2 km2, of the Caribou Poker Creek Research Watershed which included the soil properties and hydro-meteorological measurement site at the peak of the Caribou Poker Creek Research Watershed called CPEAK. The temperature calibration and verification and the coupled model induced thermodynamic and hydrodynamic interaction results’ analysis were done in this 0.2 km2 headwater sub-catchment. Finally, the 0.2 km2 headwater sub-catchment thermodynamic and hydrodynamic interaction results and analysis were employed to explain contrasting spatial and temporal variation in the observed discharge in bigger watershed scales of CPCRW during 2002 summer and fall seasons.

2. Study Area and Data

The test case selected in this study illustrates modeling permafrost active area with GIPL coupled in GSSHA.

2.1. Description of the Study Site

Figure 1 shows the test case, study area, location encompassing watershed and elevations, and model domain and grid.
The study area is contained in the in the Caribou Poker Creek Research (CPCRW), located 48 km north of Fairbanks 65°10′ N, 147°30′ W Alaska. The CPCRW site is part of the LTER (Long-Term Ecological Research) network. Parts of this watershed are underlain by permafrost, where the maximum seasonal thaw depth thickness is about 0.52 m at a low elevation point near the confluence of Poker and Caribou Creeks [16]. The CPCRW watershed encompasses an area of 101.5 km2 and is located within the boreal forest area.
As shown in Figure 1, the site includes a weather station on the summit of Caribou Peak, CPEAK, which is at an elevation of 773 m. The CPEAK station consists of a 10-m tower with various atmospheric sensors and ground temperature thermistors at several depths (http://www.lter.uaf.edu/data/metadata/id/442/inline/1). The site is located on the broad summit area of Caribou Peak, right at the tree-line elevation. There are no higher peaks nearby to limit the exposure of this site.
A 10 × 10 GSSHA/GIPL model grid was developed for simulations in the test case. The study area model domain, Figure 1, includes the CPEAK station so that the observations from the station could be used in model development and validation.

2.2. Data Acquisition

The Institute of Arctic Biology at the University of Alaska, Fairbanks maintains a Long-Term Ecological Research (LTER) site (funded by NSF) called the Bonanza Creek LTER. The Caribou-Poker Creeks Research Watershed (CPCRW) is one of the study sites for this LTER project where long-term monitoring data are collected and made available online (http://www.lter.uaf.edu/data). Clicking on Access Data, Study Sites Catalog in the Data menu of this website brings up a search filter page for the Study Sites Catalog. In this study the data from CPEAK were employed as the initial condition, and forcings used to run and validate the test case model.

2.2.1. Topography

The study area model shown in Figure 1 was developed from a 50 m digital elevation model (DEM) obtained from the National Elevation Dataset (NED), the primary elevation data product of the USGS. The data were downloaded through the National Map Viewer http://nationalmap.gov/viewer.html.

2.2.2. Soil and Land Use

Figure 2a shows the soils of the study area. A description of the soil properties, according to Rieger et al. [17], is shown in Table 1. In Table 1, literature-defined soil thermal properties’ value [18,19,20] and the soil hydraulic conductivity value [21] are also listed. Figure 2b shows the vegetation in the study area based on a unified statewide system for classifying vegetation in Alaska [22].

2.2.3. Climate of the Study Area

The study area lies in the interior climatic division of Alaska, a region that is characterized by large diurnal and annual temperature variations, low annual precipitation, low cloudiness and low humidity. In the study region, the 30-year normals [23] show that January is typically the coldest month, with a mean temperature of −24.4 °C, while July, the warmest month, has a mean of 17.1 °C. The annual precipitation in this region averages 285 mm, over half of which occurs in the months of June, July and August [24] and with 30% falling as snow from October until mid-April [23]. Figure 3 shows the monthly accumulated precipitation where Figure 3a is the monthly accumulated precipitation from the year 2001 to 2010 and Figure 3b is the 10-year average of monthly accumulated precipitation shown in Figure 3a. The monthly accumulated precipitation is derived from the available hourly precipitation data from a tipping bucket sensor at the CPEAK station. Figure 3c is the monthly average temperature from the year 2001 to 2010 and Figure 3d is the 10-year average of monthly average temperature shown in Figure 3c. The monthly average temperature is derived from the available hourly temperature data at the CPEAK station. Figure 3 shows that January is typically the coldest month and July is the warmest month with the maximum amount of precipitation.
The model employed CPEAK hydro meteorological data (http://www.lter.uaf.edu/data) for the month of May of the year 2002 including: relative humidity (%) from 1993 to 2016, air temperature (°C) from 1993 to 2016, soil temperature (°C) from 1998 to 2016, solar radiation (W h m−2) from 1999 to 2016, wind Speed (kts) from 1998 to 2005, barometric pressure (in Hg) from 2000 to 2005, and tipping bucket rainfall rates (mm/hour) from 1993 to 2016. The contact information for the data is at http://www.lter.uaf.edu/contact.

3. Methodology

3.1. Process Model Development of the Seasonal Change of the Soil Moisture Phase

To take account of the seasonal change of the soil moisture phase, water and ice phase, the GIPL, soil thermal regime simulation model, is employed to calculate a soil temperature profile in every lateral two-dimensional GSSHA computational element employing the soil moisture information from GSSHA hydrologic simulations; GSSHA, in turn, employs the temperature and water/ice phase change information from GIPL to adjust hydraulic conductivities for both the unsaturated soil flow and saturated groundwater flow [13]. This two-way coupling increases computational accuracy in both models by providing additional information and required processes which were not previously included in either model.
GSSHA takes account of important hydrological/processes such as evapotranspiration, infiltration, snow accumulation and melting, overland flow, saturated groundwater flow etc. based on physics and distributed in two-dimensional grids [25,26,27,28,29]. GIPL solves the one-dimensional (1D) non-linear heat equation with phase change employing an implicit, finite-difference, numerical scheme. GIPL treats the process of soil freezing and thawing in accordance with relationships between the soil unfrozen water content and temperature. The soil thermal state provided by GIPL is a point process, and is solved for each cell in the 2D grid.

3.2. Numerical Model of Soil Heat Transfer

The Stefan problem [30,31] with phase change is the problem of thawing or freezing via conduction of heat. GIPL solves the Stefan problem employing the enthalpy formulation. The 1D, vertical, quasi-linear heat conductive equation [32] is the basis of the GIPL numerical model:
H ( x , t ) τ = ( k ( x , t ) t ( x , τ ) )
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 time and t is temperature. k(x,t) is a thermal conductivity (Wm−1K−1). H(x,t) is an enthalpy function defined as:
H ( x , t ) = 0 t C ( x , s ) d s + L θ ( x , t )
where L is the volumetric latent heat of freeze/thaw (MJm−3), C(x,s) is volumetric heat capacity (MJm−3K−1) where s is the temperature integral term and θ (x,t) is volumetric unfrozen water content (fraction of the total water content). Boundary and initial conditions are required in solving Equation (1). The boundary condition on the upper extent of the domain corresponds to the near land surface air layer. Embedding a seasonal snow layer into the air layer is allowed by the fictitious domain formulation [33]. The upper boundary condition is defined as the Dirichlet type boundary condition.
t ( x u , τ ) = t a i r
where tair is a daily averaged air temperature. The lower boundary is set as the geothermal gradient as:
t ( x l , τ ) x = g
where g is geothermal gradient, a small constant (Km−1) usually within the range of 0.01–0.04 (Km−1) across areas with no geothermal anomalies [34]. For the initial temperature distribution, an appropriate ground temperature profile based on the point location is used.
t ( x , τ ) = t 0 ( x )
The unfrozen water content θ (x,t) is defined as:
θ ( x , t ) = η ( x ) { 1 , t t * a | c t | b , t < t *
η (x) is a volumetric soil moisture content. a and b are dimensionless positive constants [35]. The constant t * is a freezing point depression, the temperature at which ice begins to form in the soil. For simplification, we assumed the freezing point as zero degrees Celsius. The depth and time variation of the unfrozen water content θ ( x , t ) depend on hydrologic forcing and soil type. A detailed mathematical description of the model and numerical solution methods of Equation (1) can be found in Marchenko et al. [15], Nicolsky et al. [20], and Sergueev et al. [32].
GIPL input data include soil thermal properties, lithological data and vegetative cover, climate data and snow cover. The boundary conditions of the climate data are defined by Equations (3) and (4).

3.3. Numerical Model of Soil Moisture

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. The Richards Equation [36] 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 [37,38]:
C m ( ψ ) ψ τ z [ K s o i l ( ψ ) ( ψ z 1 ) ] W = 0
where Cm, the specific moisture capacity, is the slope of the soil water characteristic curve defined as rate of change of volumetric moisture content with respect to the matric head, ψ is the soil capillary head (cm), z is the vertical coordinate (cm) (downward positive), τ is time (h), K s o i l ( ψ ) (cm) is the effective hydraulic conductivity and W is a flux term added for sources and sinks (cm h−1), such as ET and infiltration. The head-based form is valid for both saturated and unsaturated conditions [39].
Heads for each cell are first computed using an implicit central-difference in space and forward difference solution and then flux updating is performed to ensure mass balance. The variables Ksoil and Cm from Equation (7) are non-linear parameters and depend on the water content of each cell. Ksoil and Cm are calculated employing Brooks and Corey [40] equations, as extended by Huston and Cass [41]. As infiltration, runoff and ET are dependent on the soil moisture state of the soil column, they are computed as part of the overall solution of Richards Equation, as described in Downer [37] and Downer and Ogden [38]. The proper application of the Richards Equation in GSSHA was described in Downer and Ogden [42]. Actual evapotranspiration is computed from the potential evapotranspiration by adjusting the potential evapotranspiration for the soil moisture in each cell (GSSHAWIKI, 2011, http://www.gsshawiki.com, accessed on 30 November 2018). Actual evapotranspiration depends on the soil moisture, hydraulic properties of the soil and plant characteristics.

3.4. Linking Soil Temperature and Soil Water Computational Nodes

The solution domain of GIPL soil temperature model overlaps in a somewhat complex manner with both the saturated and unsaturated soil water movement domains in GSSHA. If there is 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 soil horizons, as well as the deeper groundwater media as shown in Figure 4. A soil horizon is a layer having different physical, chemical or biological characteristics from the layers above and beneath. In Figure 4, soil horizons for thermodynamic and hydrodynamic process may be the same or different as per the significance of the physical, chemical and biological characteristics for individual process.
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 4. In the soil thermal profile of Figure 4, we set 56.7 °C as the highest land surface temperature as per: https://wmo.asu.edu/content/world-highest-temperature, retrieved 27 December 2018. In the soil thermal profile of Figure 4, we set −93.2 °C as the lowest land surface temperature as per: https://www.nasa.gov/content/goddard/nasa-usgs-landsat-8-satellite-pinpoints-coldest-spots-on-earth, retrieved 27 December 2018.

3.5. Linking Soil Thermodynamics with Soil Moisture Hydrodynamics

The linkage of the soil thermal and water movement solutions in GSSHA facilitates the temperature domain solution adjusting the soil thermal properties based on the varying soil moisture and the soil water movement domain solutions adjusting for the varying soil temperature condition.

3.6. Linking Soil Temperature and 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 non-dimensional 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 [43]:
S E = ( 1 1 + ( | 1.22 t | ) n ) m f o r   T 0   ° C
where n, m and α are the Van-Genuchten-Parameters (derived from literature [44]); t is 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 [43]. 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 s o i l ( t ) = e S E l n ( K t ) + ( 1 S E ) K f
where Ksoil(t) is the effective hydraulic conductivity in m s−1; Kt is the hydraulic conductivity for SE = 1 and Kf is the frozen hydraulic conductivity for SE = 0. In practice, the contribution from the frozen portion of the soil, (1 − SE)Kf is quite small and is often neglected.

4. Model Development of the Study Area

A 10 × 10 grid was developed to cover the study area, including the CPEAK station. The model includes precipitation, overland flow, unsaturated zone computations using the 1D Richards Equation, ET using the Penman-Monteith equation [45,46], and 1D soil thermal computations using GIPL. The grid with elevations, soils, and vegetation are shown in Figure 1 and Figure 2a,b respectively.

4.1. Initial Condition

The coupled model included the numerical solution of the soil moisture using Richards Equation in the unsaturated vadose zone, Equation (7), and the numerical solution of the soil thermal state using quasi-linear heat conductive equation, Equation (1). Both of those numerical solutions, Equations (1) and (7), require the initial condition. For this study, initial conditions for soil temperature were obtained from observations from CPEAK ground temperature thermistors at several depths on 2002-5-1, 1 AM, Figure 5a. Initial soil moistures were derived from observations at the Caribou-Poker Creeks Research Watershed for 2002-5-1, 1 AM [47], Figure 5b. In lieu of soil temperature and moisture measurements, the model itself can be used to establish the proper initial conditions given a sufficient prior period of observed hydro-meteorological forcing input. Downer and Ogden [38] and Senarath et al. [48] showed how the GSSHA model could not only estimate soil moistures with depth, but could actually be used to nullify the effects of assuming an incorrect initial soil moisture state by allowing the model to run through at least one saturating precipitation event and some subsequent drying period before using model results. Without either type of observation, the effect of the initial condition on model results would have to be assessed utilizing sensitivity analysis. To ensure that the initial temperature conditions did not influence the results, the maximum difference of soil temperatures at all depth levels between two successive annual cycles needs to be less than 0.01 °C [14]. In this study, we employed observed soil temperature as the initial condition and did not conduct a long-term equilibrium temperature profile assessment.

4.2. Overall Model Processes and Parameter Values

4.2.1. Thermodynamic Process

The soil thermal parametric values represent Alaskan woodland and tundra ecosystem sites in the permafrost active region. Table 2 shows the parameters and values employed in this test case model where the dominant soil type is silt loam as shown in Figure 2a. These parameter values are as per the literature [17,18,19].
As Equation (8) is deployed in the thermodynamic process, the Van-Genuchten-Parameters’ values were employed as per [44] were n, m and α are 2, 0.5 and 1, respectively, for the silt dominant soil texture.

4.2.2. Routing Process

Diffusive wave approach was deployed to route the flow. The Manning’s roughness parameter values for the routing model as shown in the Table 3 were employed from the literature [48,49]. The land use types defined in Table 3 are from Figure 2b.

4.2.3. Infiltration Process

The soil moisture and infiltration process is defined by the Richards Equation defined by Equation (7). The soil parameter values for the infiltration model as shown in the Table 4 were employed from published values for silt loam soil [21]. Brooks and Corey [40] equations, as extended by Huston and Cass [41] were employed as the soil water retention model in the infiltration process.

4.2.4. Evapotranspirational Process for the Long-Term Simulation

Hydrological modeling of seasonal change requires the model to be set in the long-term simulation mode. To take account of the depletion of the soil moisture after/before rainfall events in the long-term simulation, the evapotranspirational process defined by the Penman-Monteith equation was deployed. ET parameters were estimated based on land cover. For the Penman-Monteith method, values of land surface albedo, vegetation height (for aerodynamic resistance term), vegetation canopy resistance (for stomatal control of the loss of water), and vegetation transmission coefficient (light penetration through canopy) are needed. Table 5 shows those parameter values.
All the parameter values defined in Table 5 were taken from the GSSHA wiki https://www.gsshawiki.com. Table 5 albedo and vegetation height values are defined in Eagleson [50]; canopy resistance values are defined in Monteith [46]; and the transmission coefficient values are defined in Sutton [51].

5. Result and Discussion

5.1. Observed and Simulated Soil Profile Temperature

Figure 6 shows the comparison of the observed and simulated soil temperature at various depths over the simulation period, 1–31 May 2002. May is a transition between spring and summer in the study area when the temperatures are increasing above freezing. As shown in Figure 6, the air temperature has the most significant influence in the near surface soil layer as shown by the simulated and the observed soil temperature at 10-cm depth time series trend line, in Figure 6a. This is because the upper boundary condition for the surface of the soil is the air temperature, as shown by Equation (3), not the actual surface of the soil. This results in the model being overly sensitive to the value of air temperature very near the surface boundary. The effect of this error dissipates as the simulation computation moves deeper into the soil. As the soil layer depth increased, air temperature influence on soil thermo-dynamics diminished along with the increase in the time lag influence as shown in Figure 6b. It was also found that the soil thermal conductivity, the dominant thermo-dynamic, parameter value was most sensitive and effective for the near surface simulated temperature to capture the frequency of the observed soil temperature and the air temperature. In Figure 6a, the top 10-cm depth temperature time series shows that the maximum amplitude of the daily temperature was better simulated than the minimum amplitude. The observed and simulated maximum daily soil temperature at 10-cm depth is compared in Figure 6c. The root mean errors in Figure 6b are 0.17 °C and 0.14 °C for 90-cm and 300-cm, respectively, and in Figure 6c, 4.6 °C. While the model tends to capture the long-term evolution of soil temperature, the model underestimated daily temperature fluctuations in the top soil layer.

5.2. Catchment Hydrological Dynamics under Freezing and Thawing

5.2.1. Effective Hydraulic Conductivity with and without Thermodynamics

Figure 7 shows the comparison, with and without thermodynamics, simulation evolution of the effective hydraulic conductivity at 10 cm depth starting on 1 May 2002. The initial condition for the soil profile, shown in Figure 5a, is below freezing for the entire profile. So, for the case of soil hydraulic conductivity with thermal adjustments, the hydraulic conductivity is reduced from the beginning, corresponding to the initial condition, which was below freezing.
Figure 7 also show the corresponding air temperature with is continuously below zero degrees C for the first 107 h. That is the reason the simulated effective hydraulic conductivity with thermodynamic process in that continuous negative air temperature period is almost impermeable. More specifically, this impermeability is the result of reduced effective saturation in the soil due to the soil being completely frozen as shown in Figure 6. Figure 7 shows that the effective hydraulic conductivity without thermodynamic process maintained its value well above zero.
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 consistent with other observations of freezing and thawing of mineral and organic soils [32]. The simulated effective hydraulic conductivity in Figure 7 is from Equation (9), where the effective soil hydraulic conductivity Ksoil at a given temperature (t) is a function of hydraulic conductivity of the unfrozen soil and the effective saturation SE.

5.2.2. Hydrologic Runoff Response

Figure 8 shows the comparison of GSSHA simulated infiltration and discharge, at the outlet of the case study watershed of 0.2 km2 (headwater sub-catchment) shown in Figure 1, with and without taking account of 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 6, 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 7, which results in decreased infiltration and increased runoff in response to the precipitation event, as shown in Figure 8 graph, with thermodynamics. On the other hand, loss of frozen soil and permafrost or without taking account of thermodynamics will lead to enhanced connectivity between the surface and ground water storage regimes and decreased overland flow as shown by the graph without thermodynamics in Figure 8.

5.3. Discussion

In this thermo-hydrodynamic numerical modelling effort, the analysis was limited to the immediate effects of the soil temperature on the soil hydraulic conductivity around the measured soil location. We did not attempt to analyze the watershed wide effects due to the many complications related to doing so. The CPCRW watershed has very complex hydrology with many hydrologic processes, related to moss covered peatland and wetlands surface/water groundwater interaction, etc. Even for a summer event when and where thermodynamic effects are null, CPCRW hydrology is observed to be a complex process due to significant presence of moss covered peatland and wetland [52]. Moss-covered peatlands and wetlands significantly influence the hydraulic conductivity, infiltration, evapotranspiration and discharge [53]. In an active permafrost area there is also a seasonal patterns of freezing and thawing of the peat [53]. Teasing out the effects of frozen soil on the hydrology with these complicating factors is a very difficult process and this study has not met this level of validation, which would be a follow-on study. To isolate the immediate effects of the soil temperature on soil hydraulic conductivity, in this work we limited our analysis to a small sub-watershed of the scale of 0.2 km2, located at the CPCRW peak as shown in Figure 1. Moss, peatland and wetland are absent in this 0.2 km2 sub-watershed at the peak. The study site includes a temperature and hydro meteorological data measurement station, called CPEAK. There is no runoff information at CPCRW at this scale. However, we can use the runoff from our small watershed to assess the impact that the change in soil temperature and soil water state have on the hydrology. Simulation results indicate that the numerical runoff and infiltration results at this 0.2 km2 scale directly correspond to potential seasonal variation of discharge at larger watershed scales of CPCRW that have measured discharge data. The following paragraphs in this section demonstrate how this 0.2 km2 numerical runoff result help explain the measured variation, temporal and spatial, of discharge at larger scales within CPCRW.
Figure 9a shows the hourly air temperature at CPEAK from July 2002 through October 2002.
Figure 9b shows the corresponding hourly observed discharge (available at http://www.lter.uaf.edu/data) for C2 and C3 watersheds which are the sub-watersheds of CPCRW [54,55]. Two distinct differences in discharges from C2 and C3 can be observed:
  • Although C2 and C3 are the watersheds of similar size, 5.2 km2 and 5.7 km2, respectively, C3 discharge is significantly higher than C2 discharge.
  • From 19 September 2002, C2 discharge rapidly increases to a value roughly equal to C3 discharge.
While the C2 catchment has negligible permafrost, 3.5%, the C3 basin has 53.2% permafrost coverage. Permafrost is a permanently frozen soil. The soil thermal/hydrodynamic relations explained here could help explain the response in both basins. In an attempt to explain the above distinct behavior of C2 and C3 basin discharge, we simulated the CPEAK, 0.2 km2 watershed, for September 2002. Figure 10 shows the simulated hydrological and thermal variables for the event of 19–20 September 2002. Figure 10a shows the verification of the simulated soil temperature at 10 cm depth with the observed temperature at 10 cm depth with the root mean error of 0.2. Figure 9a and Figure 10a show a significant drop in soil temperature reaching to −3.8 °C during 19–20 September 2002, freezing the soil water content. Because of the frozen soil water content, the thermo-hydrodynamic model simulated effective hydraulic conductivity dropped to a value of 0.00003 cm h−1 at the peak of the precipitation event in Figure 10b. The simulated effective hydraulic conductivity value without taking into account of thermodynamics, was 0.02 cm h−1 at the peak of the precipitation event in Figure 10b, a difference of 4 orders of magnitude. Because of this drop in hydraulic conductivity value the thermo-hydrodynamic model simulation resulted in decreased infiltration and increased runoff, as shown in Figure 10c,d respectively. This numerical model’s increased runoff response in freezing soil conditions is the representation of the physical process that influence discharge at the C3 basin under freezing condition.

6. Conclusions

Seasonal changes lead to changes in the soil thermal regime. The seasonal change in the thermal regime, freezing and thawing, of soil has a direct impact on the hydrological response of a watershed due to changes in the void space in the soil pore-space. Therefore, to represent the watershed rainfall runoff partitioning mechanism in this freezing and thawing scenario, a hydrological model should capture the seasonal rise of soil temperatures and thaw of frozen soils; thereby, intuitively representing the transitioning of effective soil hydraulic conductivity.
This study demonstrates a coupled architecture for simulating the interaction between soil thermodynamics and hydrodynamics. The coupling architecture linked the thermo-dynamic GIPL model into GSSHA’s hydro-dynamic modeling framework. The water and ice phase change condition and its subsequent influence in the vedose zone hydrology is numerically described and experimentally demonstrated/deployed and verified in the headwater region at the peak of the Caribou Poker Creek Research Watershed.
In this study, the months of May and September of the year 2002 were employed in the long-term simulation as May is a transition between spring and summer when the temperatures are increasing above freezing and September is a transition between autumn and winter when the temperatures are decreasing below freezing. The coupled numerical schemes provided a clear insight of the temporal and spatial variability of soil thermal state during these changing seasons and the effects on hydrological processes that are significant for long-term simulation, such as infiltration and runoff. A comparative analysis, under the freezing condition, showed decreased effective hydraulic conductivity, decreased infiltration and increased runoff. The simulated results, agreeing the observed ones, showed that the air temperature has the most significant influence in the near surface soil layer and the effect diminished along with the increase in the soil profile depth. It was also found that the soil thermal conductivity, the dominant thermo-dynamic parameter value, was most sensitive and effective for the near-surface simulated temperature to capture the frequency of the observed soil temperature and the air temperature.

Author Contributions

Software development, research conceptualization, methodology, analysis, writing, reviewing and editing: N.R.P., C.W.D. and S.M.

Funding

Article Processing Charges (APCs) for the publication were supported by the Flood and Coastal Systems Research Program, US Army Engineer Research and Development Center.

Acknowledgments

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 Caribou-Poker Creeks Research Watershed (CPCRW) in Alaska. Constructive comments from two anonymous reviewers are greatly appreciated.

Conflicts of Interest

The authors declare no conflict of interest.

References

  1. IPCC. Climate Change 2014: Synthesis Report. Contribution of Working Groups I, II and III to the Fifth Assessment Report of the Intergovernmental Panel on Climate Change; Core Writing Team, Pachauri, R.K., Meyer, L.A., Eds.; IPCC: Geneva, Switzerland, 2014; 151p. [Google Scholar]
  2. U.S. Global Change Research Program. Climate Science Special Report: Fourth National Climate Assessment, Volume I; Wuebbles, D.J., Fahey, D.W., Hibbard, K.A., Dokken, D.J., Stewart, B.C., Maycock, T.K., Eds.; USGCRP: Washington, DC, USA, 2017; 470p. [Google Scholar] [CrossRef]
  3. Qian, B.; Gregorich, E.G.; Gameda, S.; Hopkins, D.W.; Wang, X.L. Observed soil temperature trends associated with climate change in Canada. J. Geophys. Res. 2011, 116. [Google Scholar] [CrossRef] [Green Version]
  4. Schipper, L.A.; Hobbs, J.K.; Rutledge, S.; Arcus, V.L. Thermodynamic theory explains the temperature optima of soil microbial processes and high Q10 values at low temperatures. Glob. Chang. Biol. 2014, 20, 3578–3586. [Google Scholar] [CrossRef] [PubMed]
  5. Zhuang, Q.; Melillo, J.M.; McGuire, A.D.; Kicklighter, D.W.; Prinn, R.G.; Steudler, P.A.; Felzer, B.S.; Hu, S. Net emissions of CH4 and CO2 in Alaska: Implications for the region’s greenhouse gas budget. Ecol. Appl. 2007, 17. [Google Scholar] [CrossRef]
  6. Davidson, E.A.; Janssens, I.A. Temperature sensitivity of soil carbon decomposition and feedbacks to climate change. Nature 2006, 440, 65–173. [Google Scholar] [CrossRef] [PubMed]
  7. Melvina, A.M.; Larsenb, P.; Boehlertc, B.; Neumannc, J.E.; Chinowskye, P.; Espinete, X.; Martinichf, J.; Baumannc, M.S.; Rennelsc, L.; Bothnerc, A.; et al. Climate change damages to Alaska public infrastructure and the economics of proactive adaptation. Proc. Natl. Acad. Sci. USA 2017, 114, E122–E131. [Google Scholar] [CrossRef] [PubMed]
  8. Interagency Climate Adaptation Team. Adapting to Climate Change in Minnesota, Minnesota Pollution Control Agency; Document Number: p-gen4 -0; Interagency Climate Adaptation Team: MN, USA, 2017. [Google Scholar]
  9. Dunne, T.; Black, R.D. Runoff processes during snowmelt. Water Resour. Res. 1971, 7, 1160–1172. [Google Scholar] [CrossRef]
  10. Stähli, M.; Jansson, P.E.; Lundin, L.C. Soil moisture redistribution and infiltration in frozen sandy soils. Water Resour. Res. 1999, 35, 95–103. [Google Scholar] [CrossRef] [Green Version]
  11. Stähli, M. Hydrological significance of soil frost for pre-alpine areas. J. Hydrol. 2017, 546, 90–102. [Google Scholar] [CrossRef]
  12. Downer, C.W.; Ogden, F.L. Gridded Surface Subsurface Hydrologic Analysis (GSSHA) User’s Manual, Version 1.43 for Watershed Modeling System 6.1; ERDC/CHL SR-06-1; System Wide Water Resources Program, Coastal and Hydraulics Laboratory, U.S. Army Corps of Engineers, Engineer Research and Development Center: Vicksburg, MS, USA, 2006; 207p. [Google Scholar]
  13. Pradhan, N.R.; Downer, C.W.; Marchenko, S.; Lijedahl, A.; Douglas, T.A.; Byrd, A. Development of a Coupled Framework for Simulating Interactive Effects of Frozen Soil Hydrological Dynamics in Permafrost Regions; ERDC TR-13-15; U.S. Army Engineer Research and Development Center: Vicksburg, MS, USA, 2013. [Google Scholar]
  14. Jafarov, E.E.; Marchenko, S.S.; Romanovsky, V.E. Numerical modeling of permafrost dynamics in Alaska using a high spatial resolution dataset. Cryosphere Discuss. 2012, 6, 89–124. [Google Scholar] [CrossRef]
  15. Marchenko, S.; Romanovsky, V.; Tipenko, G. Numerical modeling of spatial permafrost dynamics in Alaska. In Proceedings of the Nineth International Conference on Permafrost, Willey, Institute of Northern Engineering, University of Alaska, Fairbanks, AK, USA, 29 June–3 July 2008. [Google Scholar]
  16. Bolton, W.R.; Hinzman, L.D.; Yoshikawa, K. Stream flow studies in a watershed underlain by discontinuous peramafrost. In Water Resources in Extreme Environments; Kane, D.L., Ed.; American Water Resources Association: Anchorage, AK, USA, 2000; pp. 31–36. [Google Scholar]
  17. Rieger, S.; Furbush, C.E.; Schoephorster, D.; Sumnmerfield, H.; Geiger, L.C. Soils of the Caribou-Poker Creeks Research Watershed, Alaska; CRREL Technical Report 236, AD 744451; Cold Regions Research and Engineering Lab: Hanover, NH, USA, 1972. [Google Scholar]
  18. Farouki, O.T. The thermal properties of soils in cold regions. Cold Reg. Sci. Technol. 1981, 5, 67–75. [Google Scholar] [CrossRef]
  19. Konovalov, A.A.; Roman, L.T. The thermophysical properties of peat soils. Soil Mech. Found. Eng. 1973, 10, 179–181. [Google Scholar] [CrossRef]
  20. Nicolsky, D.J.; Romanovsky, V.E.; Tipenko, G.S. Using in-situ temperature measurements to estimate saturated soil thermal properties by solving a sequence of optimization problems. Cryosphere 2007, 1, 41–58. [Google Scholar] [CrossRef] [Green Version]
  21. Rawls, W.J.; Brakensiek, D. Prediction of soil water properties for hydrologic modelling. In Watershed Management in the Eighties; Ward, T.J., Ed.; American Society of Civil Engineers: Reston, VI, USA, 1985; pp. 293–299. [Google Scholar]
  22. Viereck, L.A.; Dyrness, C.T. A Preliminary Classification Systems for Vegetation of Alaska; USDA Pacific Northwest Forest and Range Experiment Station, General Technical Report PNW-106; USDA: Washington, DC, USA, 1980.
  23. U.S. Department of Commerce. Monthly Normals of Temperature, Precipitation, and Heating and Cooling Degree-Days, 1941–1970; National Oceanic and Atmospheric Adminstration: Asheville, NC, USA, 1973.
  24. U.S. Department of Commerce. Local Climatological Data Annual Summary with Comparative Data, Fairbanks, Alaska; National Oceanic and Atmospheric Adminstration: Asheville, NC, USA, 1980.
  25. Downer, C.W.; Skahill, B.E.; Graulau-Santiago, J.A.; Weston, D.; Pradhan, N.R.; Byrd, A.R. Gridded Surface Subsurface Hydrologic Analysis Modeling for Analysis of Flood Design Features at the Picayune Strand Restoration Project; ERDC/CHL TR-15-X; U.S. Army Engineer Research and Development Center: Vicksburg, MS, USA, 2015. [Google Scholar]
  26. Massey, T.C.; Pradhan, N.R.; Byrd, A.R.; Cresitello, D.E. USACE-ERDC coastal storm modelling systems in support of hurricane sandy operations. Flood Risk Manag. Newslett. 2013, 6, 2–3. [Google Scholar]
  27. Ogden, F.L.; Pradhan, N.R.; Downer, C.W.; Zahner, J.A. Relative importance of impervous area, drainage density, width function, and subsurface storm drainage on flood runoff from an urbanized catchment. Water Resour. Res. 2011, 47, W12503. [Google Scholar] [CrossRef]
  28. Pradhan, N.R.; Downer, C.W.; Johnson, B.E. A physics based hydrologic modeling approach to simulate non-point source pollution for the purposes of calculating TMDLs and designing abatement measures. In Practical Aspects of Computational Chemistry III; Leszczynski, J., Shukla, M.K., Eds.; Springer: New York, NY, USA, 2014; Chapter 9; pp. 249–282. [Google Scholar]
  29. Downer, C.W.; Pradhan, N.R.; Byrd, A. Modeling Subsurface Storm and Tile Drain Systems in GSSHA with Superlink; ERDC TR-14-11; U.S. Army Engineer Research and Development Center: Vicksburg, MS, USA, 2014. [Google Scholar]
  30. Alexiades, V.; Solomon, A.D. Mathematical Modeling of Melting and Freezing Processes; Hemisphere: Washington, DC, USA, 1993; 325p. [Google Scholar]
  31. Verdi, C. Numerical Aspects of Parabolic Free Boundary and Hysteresis Problems; Lecture Notes in Mathematics; Springer: New York, NY, USA, 1994; pp. 213–284. [Google Scholar]
  32. Sergueev, D.; Tipenko, G.; Romanovsky, V.; Romanovskii, N. Mountain permafrost thickness evolution under the influence of long-term climate fluctuations (results from numerical simulation). In Proceedings of the 8th International Conference on Permafrost, Zürich, Switzerland, 21–25 July 2003; Taylor and Francis: Philadelphia, PA, USA, 2003; pp. 1017–1021. [Google Scholar]
  33. Marchuk, G.I.; Kuznetsov, Y.A.; Matsokin, A.M. Fictitious domain and domain decomposition methods. Sov. J. Num. Anal. Math. Modell. 1986, 1, 1–86. [Google Scholar] [CrossRef]
  34. Pollack, H.N.; Hurter, S.J.; Johnson, J.R. Heat flow from the earth’s interior: Analysis of the global data set. Rev. Geophys. 1993, 31, 267–280. [Google Scholar] [CrossRef]
  35. Lovell, C.W. Temperature effects on phase composition and strength of partially frozen soil. Highw. Res. Board Bull. 1957, 168, 74–95. [Google Scholar]
  36. Richards, L.A. Capillary conduction of liquids in porous mediums. Physics 1931, 1, 318–333. [Google Scholar] [CrossRef]
  37. Downer, C.W. Identification and Modeling of Important Stream Flow Producing Processes in Watersheds. Ph.D. Thesis, University of Connecticut, Storrs, CT, USA, 2002. [Google Scholar]
  38. Downer, C.W.; Ogden, F.L. Appropriate vertical discretization of Richards’ equation for two-dimensional watershed-scale modeling. Hydrol. Process. 2004, 18, 1–22. [Google Scholar] [CrossRef]
  39. Haverkamp, R.; Vauclin, M.; Touma, J.; Wierenga, P.J.; Vachaud, G. A Comparison of Numerical Simulation Models for One-Dimensional Infiltration. Soil Sci. Soc. Am. J. 1977, 41, 2. [Google Scholar] [CrossRef]
  40. Brooks, R.J.; Corey, A.T. Hydraulic Properties of Porous Media; Hydrology Paper 3; Colorade State University: Fort Collins, CO, USA, 1964. [Google Scholar]
  41. Huston, J.L.; Cass, A. A retentivity function for use in soil-water simulation models. J. Soil Sci. 1987, 38, 105–113. [Google Scholar]
  42. Downer, C.W.; Ogden, F.L. Prediction of runoff and soil moisture at the watershed scale: Effects of model complexity and parameter assignment. Water Resour. Res. 2003, 39, 1045. [Google Scholar] [CrossRef]
  43. Zhang, Y.; Carey, S.K.; Quinton, W.L.; Janowicz, J.R.; Pomeroy, J.W.; Flerhinger, G.N. Comparison of algorithms and parameterisations for infiltration into organic-covered permafrost soils. Hydrol. Earth Syst. Sci. 2010, 14, 729–750. [Google Scholar] [CrossRef] [Green Version]
  44. Schulla, J. Model Description WaSiM. 2017. Available online: http://www.wasim.ch/downloads/doku/wasim/wasim_2017_en.pdf (accessed on 2 January 2019).
  45. Monteith, J. Evaporation and environment. Symp. Soc. Exp. Biol. 1965, 19, 205–234. [Google Scholar] [PubMed]
  46. Monteith, J. Evaporation and surface temperature. R. J. Q. Meteorol. Soc. 1981, 107, 1–27. [Google Scholar] [CrossRef]
  47. Chapin, F.S.; Ruess, R.W. Caribou-Poker Creeks Research Watershed: Hourly Soil Moisture at Varying Depths from 2000 to Present. 2018. Available online: https://doi.org/10.6073/pasta/5207a6643bae4f9792aa21180bc1c08d (accessed on 9 December 2018).
  48. Senarath, S.U.S.; Ogden, F.L.; Downer, C.W.; Sharif, H.O. On the calibration and verification of two-dimensional, distributed, Hortonian, continuous watershed models. Water Resour. Res. 2000, 36, 1495–1510. [Google Scholar] [CrossRef] [Green Version]
  49. Engman, E.T. Roughness coefficients for routing surface runoff. J. Irrig. Drain. Eng. 1986, 112, 39–53. [Google Scholar] [CrossRef]
  50. Eagleson, P.S. Dynamic Hydrology; McGraw-Hill: New York, NY, USA, 1970; 461p. [Google Scholar]
  51. Sutton, O.G. Micrometeorology; McGraw Hill: New York, NY, USA, 1953; p. 68. [Google Scholar]
  52. Downer, C.W.; Wahl, M. EGU2017-10350 Conceptualizing Peatlands in a Physically-Based Spatially Distributed Hydrologic Model; European Geophysical Union: Munich, Germany, 2017. [Google Scholar]
  53. Whitfield, P.H.; Kamp, G.V.D.; St-Hilaire, A. Introduction to peatlands special issue: Improving hydrological prediction in Canadian peatlands. Can. Water Resour. J. 2009, 34, 303–310. [Google Scholar] [CrossRef]
  54. Bolton, W.R.; Hinzman, L.; Yoshikawa, K. Water balance dynamics of three small catchments in a Sub-Arctic boreal forest. In Northern Research Basins Water Balance (Proceedings of a Workshop Held at Victoria, Canada, March 2004); Kane, D.L., Yang, D., Eds.; IAHS: London, UK, 2004; pp. 213–223. [Google Scholar]
  55. Haugen, R.K.; Slaughter, C.W.; KEHowe, K.E.; Dingman, S.L. Hydrology and Climatology of the Caribou-Poker Creeks Research Watershed, Alaska; CRREL Report 82-26; U.S. Army Cold Regions Research and Engineering Laboratory: Hanover, NH, USA, 1982. [Google Scholar]
Figure 1. The test case, study area, location encompassing watershed and elevations, and model domain and grid.
Figure 1. The test case, study area, location encompassing watershed and elevations, and model domain and grid.
Water 11 00116 g001
Figure 2. Land use and soil map of the study area: (a) soil type and (b) land use type.
Figure 2. Land use and soil map of the study area: (a) soil type and (b) land use type.
Water 11 00116 g002
Figure 3. Graphical information of the monthly precipitation and temperature.
Figure 3. Graphical information of the monthly precipitation and temperature.
Water 11 00116 g003
Figure 4. Thermodynamics coupling into the hydrodynamics.
Figure 4. Thermodynamics coupling into the hydrodynamics.
Water 11 00116 g004
Figure 5. Initial conditions: (a) Initial soil temperature profile for the thermodynamic numerical simulation and (b) Initial soil moisture content profile for the unsaturated zone numerical simulation.
Figure 5. Initial conditions: (a) Initial soil temperature profile for the thermodynamic numerical simulation and (b) Initial soil moisture content profile for the unsaturated zone numerical simulation.
Water 11 00116 g005
Figure 6. Time-series of observed and simulated temperature at: (a) 10-cm depth, (b) 90-cm and 300-cm depths and (c) maximum values at 10-cm depth, (d) corresponding hourly precipitation rate.
Figure 6. Time-series of observed and simulated temperature at: (a) 10-cm depth, (b) 90-cm and 300-cm depths and (c) maximum values at 10-cm depth, (d) corresponding hourly precipitation rate.
Water 11 00116 g006
Figure 7. Hydraulic conductivity under freezing and thawing soil active layer.
Figure 7. Hydraulic conductivity under freezing and thawing soil active layer.
Water 11 00116 g007
Figure 8. Hydrologic response, with and without thermodynamics, to precipitation events: (a) precipitation events (b) infiltration comparison (c) discharge comparison.
Figure 8. Hydrologic response, with and without thermodynamics, to precipitation events: (a) precipitation events (b) infiltration comparison (c) discharge comparison.
Water 11 00116 g008
Figure 9. Caribou Poker Creek Research Watershed (CPCRW) (a) Observed discharge Observed discharge of C2 and C3 catchments, sub-catchments of CPCRW shown in Figure 1 (b) air temperature from Cpeak station shown in Figure 1.
Figure 9. Caribou Poker Creek Research Watershed (CPCRW) (a) Observed discharge Observed discharge of C2 and C3 catchments, sub-catchments of CPCRW shown in Figure 1 (b) air temperature from Cpeak station shown in Figure 1.
Water 11 00116 g009
Figure 10. Catchment hydrological modeling with and without soil thermal dynamics during September 2002; (a) simulated and observed soil temperature (b) precipitation (c) simulated infiltration with and without soil thermodynamics (d) simulated discharge with and without thermodynamics.
Figure 10. Catchment hydrological modeling with and without soil thermal dynamics during September 2002; (a) simulated and observed soil temperature (b) precipitation (c) simulated infiltration with and without soil thermodynamics (d) simulated discharge with and without thermodynamics.
Water 11 00116 g010aWater 11 00116 g010b
Table 1. Soil textural property of the study area.
Table 1. Soil textural property of the study area.
Soil SeriesUSDA TextureDrainageSoil Heat Conductivity (W/mK)Soil Volumetric Heat Capacity (J/m3K)Saturated Hydraulic Conductivity (cm/hour)
FairplaySilt loam and gravelly silt loamModerately well drained0.1–2.81,600,000–2,800,0000.6
OlnesSilt loam and very gravelly silt loamWell drained0.1–2.81,600,000–2,800,0000.6
Table 2. Soil thermal parameter values.
Table 2. Soil thermal parameter values.
Soil Heat Conductivity (W/mK)Soil Volumetric Heat Capacity (J/m3K)
3.02,800,000
Table 3. Manning’s roughness values.
Table 3. Manning’s roughness values.
Land Use TypeManning’s Roughness Values
Coniferous open0.17
Coniferous woodland0.19
Deciduous closed0.2
Shrub tall0.25
Table 4. Soil parameter values for the Richards infiltration scheme.
Table 4. Soil parameter values for the Richards infiltration scheme.
Infiltration ParameterSoil Layer
TopMiddleLower
Soil thickness (m)0.200.50-
Saturated hydraulic conductivity (m/s)0.000080.000080.00008
Pore size distribution index0.600.6940.694
Bubbling Pressure head (m)0.080.060.06
Porosity (m3/m3)0.420.400.40
Residual soil moisture content (m3/m3)0.040.0450.045
Table 5. Evapotranspiration parameter values.
Table 5. Evapotranspiration parameter values.
Land Use TypeAlbedoVegetation Height (m)Canopy Resistance (s/m)Transmission Coefficient
Coniferous open0.15101200.18
Coniferous woodland0.15101200.18
Deciduous closed0.2122000.18
Shrub tall0.21.31500.5

Share and Cite

MDPI and ACS Style

Pradhan, N.R.; Downer, C.W.; Marchenko, S. Catchment Hydrological Modeling with Soil Thermal Dynamics during Seasonal Freeze-Thaw Cycles. Water 2019, 11, 116. https://doi.org/10.3390/w11010116

AMA Style

Pradhan NR, Downer CW, Marchenko S. Catchment Hydrological Modeling with Soil Thermal Dynamics during Seasonal Freeze-Thaw Cycles. Water. 2019; 11(1):116. https://doi.org/10.3390/w11010116

Chicago/Turabian Style

Pradhan, Nawa Raj, Charles W. Downer, and Sergei Marchenko. 2019. "Catchment Hydrological Modeling with Soil Thermal Dynamics during Seasonal Freeze-Thaw Cycles" Water 11, no. 1: 116. https://doi.org/10.3390/w11010116

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop