HcGroundImplicit: Numerical solution of ground heat conduction
In geocryology/PermafrostTools: PermafrostTools

Description Usage Arguments Details Value Author(s) Examples

Implicitly solves ground heat conduction based on finite- difference Crank-Nicolson scheme.

1
2
3

HcGroundImplicit(dt, mat, bc.lower, bc.upper, unfrozen.type = "INTERVAL",
  unfrozen.par = c(1), bcutype = "NEUMANN", bcltype = "NEUMANN",
  layers.sno = c(0), type.sno = "CALONNE", type.gnd = "COSENZA")

`dt`	Time step [s].
`mat`	Data frame with relevant ground properties.
`bc.lower`	Lower boundary condition as either [W m-2] (Neumann) or [C or K] (Dirichlet). Array with dimension (N).
`bc.upper`	Upper boundary condition as either [W m-2] (Neumann) or [C or K] (Dirichlet). Array with dimension (N).
`unfrozen.type`	Character string indicating the type of function to be used; the standard is "INTERVAL". The full options are: "DALLAMICO": Dall’Amico, M., Endrizzi, S., Gruber, S., & Rigon, R. (2011). A robust and energy-conserving model of freezing variably-saturated soil. The Cryosphere, 5(2), 469–484. doi:10.5194/tc-5-469-2011 Typical values of van Genuchten parameters are found in Table 2 of Gubler, S., Endrizzi, S., Gruber, S., & Purves, R. S. (2013). Sensitivities and uncertainties of modeled ground temperatures in mountain environments. Geoscientific Model Development, 6(4), 1319–1336. doi:10.5194/gmd-6-1319-2013 "MOTTAGHY": Mottaghy, D., & Rath, V. (2006). Latent heat effects in subsurface heat transport modelling and their impact on palaeotemperature reconstructions. Geophysical Journal International, 164(1), 236-245. "INTERVAL": Phase change takes place in an interval of specified with (unfrozen.par) below 0C and at a constant rate.
`unfrozen.par`	Parameter set for the chosen unfrozen water function. "DALLAMICO": unfrozen.par[1]: van Genuchten alpha [mm-1], unfrozen.par[2]: van Genuchten n [-], unfrozen.par[3]: residual water content [m3/m3]. The saturated water content is given by the input (mat$wat) and can thus vary with depth. "MOTTAGHY": unfrozen.par[1]: width of freezing interval [K], unfrozen.par[2]: omega. "INTERVAL" unfrozen.par[1]: width of freezing interval [K].
`bcutype`	Type of upper boundary condition used. Can be "NEUMANN" (standard, fixed heat flux) or "DIRICHLET" (fixed temperature).
`bcltype`	Type of upper boundary condition used. Can be "NEUMANN" (standard, fixed heat flux) or "DIRICHLET" (fixed temperature).
`layers.sno`	Index to identify snow layers within 'mat'.
`type.sno`	Character string indicating the type of parameterization to be used for snow; the standard is CALONNE. Available are: CALONNE: Calonne, N., Flin, F., Morin, S., Lesaffre, B., du Roscoat, S. R., & Geindreau, C. (2011). Numerical and experimental investigations of the effective thermal conductivity of snow. Geophysical Research Letters, 38(23), doi:10.1029/2011GL049234 "YEN": Yen, Y.-C. (1981). Review of thermal properties of snow, ice and sea ice (34 pages). Hanover, NH, USA. COSENZA: Cosenza, P., Guerin, R., & Tabbagh, A. (2003). Relationship between thermal conductivity and water content of soils using numerical modelling. European Journal of Soil Science, 54(3), 581–588. doi:10.1046/j.1365-2389.2003.00539.x STURM: Sturm, M., J. Holmgren, M. König, and K. Morris (1997), The thermal conductivity of seasonal snow, Journal of Glaciology, 43(143), 26–41. JORDAN: Jordan, R. E., Andreas, E. L., & Makshtas, A. P. (1999). Heat budget of snow-covered sea ice at North Pole 4. Journal of Geophysical Research, 104(C4), 7785. doi:10.1029/1999JC900011 WILLIAMS: Figure 4.11 in "The Frozen Earth: fundamentals of geocryology" by P. Williams and M. Smith, 1989.
`type.gnd`	Character string indicating the type of parameterization to be used for ground; the standard is "COSENZA". Available are: COSENZA: Cosenza, P., Guerin, R., & Tabbagh, A. (2003). Relationship between thermal conductivity and water content of soils using numerical modelling. European Journal of Soil Science, 54(3), 581–588. doi:10.1046/j.1365-2389.2003.00539.x GEOMETRIC: Geometric mean, intermediate mixed conductivity model, approximation of randomly oriented consituent elements. ARITHMETIC: Arithmetic mean, high mixed conductivity model, approximation of consituent elements layered parallel to temperature gradient. HARMONIC: Harmonic mean, low mixed conductivity model, approximation of consituent elements layered normal to temperature gradient.

Before the solution, ground and snow properties are updated. This means that the proportions of ice and liquid water are updated and that the thermal conductivity and heat capacity of the soil are changed accordingly. Phase change is parameterised as apparent heat capacity based on the derivative of the freezing characteristic curve. The example below shows how ground properties are described.

Returns input array with Tj solved for next time step and temperature- dependent properties updated.

Stephan Gruber <stephan.gruber@carleton.ca>

#== Settings ==================
dt <- 3600    #time step [s]
ns <-  240    # how many time steps to compute before saving? 
st <- 4 * 365 # length of run [days]
dzmin <- 0.02 # minimal z-spacing ("how fine is the grid") [m]
zmax  <- 25   # depth of lowermost node center ("how large is the domain") [m]
base  <- 1.5  # resolution reduction ("how strong is the grid corsened with depth")
              # 1: equal spacing, >1 growing 
Ti <- -1      # initial temperature [C]

type.gnd      <- "COSENZA"  # parameterization for soil thermal conductivity 
unfrozen.type <- "DALLAMICO" # paramterization for unfrozen water content
unfrozen.par  <- c(0.001, 1.4, 0.05)  # parameter(s) for unfrozen water content parameterization
bcutype       <- "DIRICHLET"# upper boundary condition type
bcltype       <- "NEUMANN"  # lower boundary condition type
layers.sno    <- c(0)       # index of nodes that are snow
type.sno      <- "CALONNE"  # parameterization for snow thermal conductivity

#== Preparation ================
st    <- st * 3600 * 24 # convert to [s]
t     <- seq(from = 0, to = st, by = dt) #vector of time [s]
nt    <- length(t)
soild <- SoilDiscretize(dzmin, zmax, base)
nz    <- length(soild$z) #number of soil discretizations

#make boundary condition
bc.upper <- seq(from = -1, to = 3, length.out = nt) #linear warming
bc.lower <- rep(0.05, nt)

#make material properties and discretization
mat <- data.frame( zj = -soild$z, # depth of node [m]
                  dz = soild$dz, # with of node [m]
                  Tj = rep( Ti, nz), # initial temperature [C]
                 soi = rep(0.5, nz), # volumetric proportion of soil matrix [0-1]
                 wat = rep(0.0, nz), # volumetric proportion of total water [0-1]
                 liq = rep(  0, nz), # calculate later, volumetric proportion of liquid water [0-1]
                 ice = rep(  0, nz), # calculate later, volumetric proportion of ice [0-1]                
               k.soi = rep(2.5, nz), # thermal conductivity of soil matrix [W m-1 K-1]
               c.soi = rep(2e6, nz)  # volumetric heat capacity of soil matrix [J m-3 K-1]
)

#== Target arrays to hold results =========
out.Tj  <- rbind(NULL,t(mat$Tj))
out.liq <- rbind(NULL,t(mat$liq))
out.t   <- t[0]
out.dz  <- mat$dz

#== Loop over time and solve =======
for (i in 1:nt) {
  mat <- HcGroundImplicit(dt, mat, bc.lower[i], bc.upper[i], unfrozen.type = unfrozen.type, 
                          unfrozen.par = unfrozen.par, bcutype = bcutype, bcltype = bcltype,
                          layers.sno = layers.sno, type.sno = type.sno, type.gnd = type.gnd)
#output only in specified interval
check <- i/ns - floor(i/ns)
  if (check == 0) {
    out.Tj  <- rbind( out.Tj,t(mat$Tj))
    out.liq <- rbind(out.liq,t(mat$liq))
    out.t   <- t[i]
    print(paste(i,nt))   
  }		
}

#== Plotting ================     
for (n in 1:147) {
 plot(out.Tj[n,],mat$zj, type="l",xlim=c(-1.2 , 3))     
 abline(v=0)
}