R/hc_ground_implicit.R

In geocryology/PermafrostTools: PermafrostTools

# =============================================================================
#'
#' @title Numerical solution of ground heat conduction 
#'
#' @description Implicitly solves ground heat conduction based on finite-
#'              difference Crank-Nicolson scheme.   
#'
#' @details  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.
#'
#' @param dt Time step [s].
#' @param mat Data frame with relevant ground properties.  
#'
#' @param bc.lower Lower boundary condition as either [W m-2] (Neumann) or [C or K] 
#'            (Dirichlet). Array with dimension (N).  
#' @param bc.upper Upper boundary condition as either [W m-2] (Neumann) or [C or K] 
#'            (Dirichlet). Array with dimension (N). 
#' @param bcutype Type of upper boundary condition used. Can be "NEUMANN" 
#'                (standard, fixed heat flux) or "DIRICHLET" (fixed temperature).  
#' @param bcltype Type of upper boundary condition used. Can be "NEUMANN"
#'                (standard, fixed heat flux) or "DIRICHLET" (fixed temperature). 
#' 
#' @param 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.  	
#'
#' @param 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].  
#' @param 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.
#' 
#' @param 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.
#'
#' @param layers.sno Index to identify snow layers within 'mat'.
#
#' @return Returns input array with Tj solved for next time step and temperature-
#'         dependent properties updated.
#' 
#' @export
#' @examples
#' #== 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)
#' } 
#' 
#' @author Stephan Gruber <stephan.gruber@@carleton.ca>
#'
# =============================================================================

HcGroundImplicit <- function(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") {
                             	
  #make zero volumetric heat source term
  gj <- mat$Tj * 0
 
  #TIME LOOP
  for (j in 1:length(bc.upper)) {
 	#compute ground properties
    mat <- Unfrozen(mat, unfrozen.type, unfrozen.par)
    mat <- TermCond(mat, layers.sno = layers.sno, type.sno = type.sno, type.gnd = type.gnd)
    mat <- HeatCapacity(mat)
    #mat <- SnowProp(mat, snow.type, snow.par)  # not ready yet
 	  
    #compute new solution (transpose to get from array to data frame)  	
  	mat$Tj <- Cranknicolson(mat$Tj, mat$kj, mat$cj, gj, mat$zj, dt, 
  	                          bc.lower[j], bc.upper[j], 
  	                          bcutype = bcutype, bcltype = bcltype)[]                                                                                              
  } 	  	  	  	
  return(mat)
}