R/swe.deltasnow.R

In swemod: Snow Water Equivalent Modeling with Parametric and Statistical Models

swe.deltasnow <- function(data, rho.max=440, rho.null=107, c.ov=0.004838623, k.ov=0.1996423, k.exp=0.02634508, tau=0.02, timestep = 24, verbose = FALSE) {
 
  if(!inherits(data,"data.frame"))
    stop("swe.deltasnow: data must be of class data.frame")
  
  if(!any((is.element(colnames(data), c("hs","date")))))
    stop("swe.deltasnow: data must contain at least two columns named 'hs' and 'date'")
  
  Hobs <- data$hs
  if(any(is.na(Hobs)))
    stop("swe.deltasnow: snow depth data must not be NA")
  if(!all(Hobs >= 0))
    stop("swe.deltasnow: snow depth data must not be negative")
  if(!is.numeric(Hobs))
    stop("swe.deltasnow: snow depth data must be numeric")
  if(Hobs[1] != 0)
    stop("swe.deltasnow: snow depth observations must start with 0")
  
  if(!inherits(data$date,"character"))
    stop("swe.deltasnow: date column must be of class character")
  
  z <- zoo(Hobs,as.Date(data$date))
  if(!is.regular(z, strict = TRUE))
    stop("swe.deltasnow: date column must be strictly regular")
  
  #comptype <- match.arg(compactiontype)
  #if(length(comptype) == 0) 
  comptype <- "exp"
  
  snowpack.dd <- NULL
  
  # check coefficients, e.g. tau > 0, etc.
  
  H        <- c()                              # modeled total height of snow at any day [m]
  SWE      <- c()                              # modeled total SWE at any day [kg/m2]
  ly       <- 1                                # layer number [-]
  ts       <- timestep * 3600                  # timestep between observations [s]
  
  #rho.max  <- 405                              # maximum density of a snow layer [kg/m3]
  #rho.null <- 90                               # fresh snow density [kg/m3]
 
  # compaction
  #k.exp    <- 0.03                            # compaction factor for exponential compaction [m3/kg]
  
  # overburden
  #c.ov     <- 0.01                            # overburden factor due to fresh snow or rain
  #k.ov     <- 0.14                            # 
  
  # tau     <- 0.02                           # uncertainty bound [m]
  
  # preallocate matrix as days X layers   
  ly.tot   <- length(which(Hobs>0))          # maximum number of layers [-]
  day.tot  <- length(Hobs)                    # total days from first to last snowfall [-]
  h        <- matrix(0,ly.tot,day.tot)         # modeled height of snow in all layers [m]
  swe      <- matrix(0,ly.tot,day.tot)         # modeled swe in all layers [kg/m2]
  age      <- matrix(0,ly.tot,day.tot)         # age of modeled layers [days]

  # helper
  m        <- rep("",day.tot)                  # vector of (verbose) messages
  prec     <- 10^-10                           # precision for arithmetic comparisons
  
 
  
  
  
  #-------------------------------------------------------------------------
  # compaction of snowpack
  #-------------------------------------------------------------------------
  compactH <- function(h.d, swe.d, age.d, ly.tot, ly, comptype, k.exp, rho.max, ts, prec){
    # h.d=h[,t];swe.d=swe[,t];age.d=age[,t]
    
    #snowpack.dd <- NULL
      
    # .d  -> today
    # .dd -> tomorrow
    
    # compute overburden for each layer
    swe.hat.d <- c()
    for(i in 1:(ly.tot-1)){
      swe.hat.d[i] <- sum(swe.d[(i+1):ly.tot])
    }
    swe.hat.d[ly.tot] <- 0
    
    snowpack.d   <- data.frame(h  = h.d, swe = swe.d, swe.hat = swe.hat.d, age = age.d)
    H.d          <- sum(snowpack.d$h)
    
    #snowpack.dd  <<- data.frame(t(apply(snowpack.d, 1, metamorphic.compaction, H.d, comptype))) # compress today's snowpack layerwise
    #colnames(snowpack.dd) <<- c("h","swe","age","rho")
    #rownames(snowpack.dd) <<- paste0("dd.layer",1:nrow(snowpack.dd))
    
    a <- data.frame(t(apply(snowpack.d[(1:ly),], 1, metamorphic.compaction, H.d, comptype, k.exp, rho.max, ts, prec)))
    b <- data.frame(rep(0,ly.tot-ly),rep(0,ly.tot-ly),rep(0,ly.tot-ly),rep(0,ly.tot-ly))
    colnames(a) <- colnames(b) <- c("h","swe","age","rho")
    snowpack.dd <<- rbind(a,b)
    rownames(snowpack.dd) <<- paste0("dd.layer",1:nrow(snowpack.dd))
    
    
    return(snowpack.dd)
  }
  
  
  
  #-------------------------------------------------------------------------
  # compaction of individual snow layers without additional load 
  # today's values are compacted to tomorrow's values
  #-------------------------------------------------------------------------
  metamorphic.compaction <- function(x, H.d, comptype, k.exp, rho.max, ts, prec){
    
    g        <- 9.81  # gravity [ms-2]
    eta.null <- 8.5e6 # effective compactive viscosity of snow for "zero-density" [Pa s]
    
    # .d  -> today
    # .dd -> tomorrow
    age.d <- ifelse(x[1] == 0, 0, x[4])
    switch(comptype,
           # "power" = {
           #   h.dd <- x[1] * (age.d + 1)^-k.pow
           #   h.dd <- ifelse(x[2]/h.dd > rho.max, x[2]/rho.max, h.dd)
           # },
           "exp" = {
             h.dd <- x[1]/(1 + (x[3] * g * ts)/eta.null * exp(-k.exp * x[2]/x[1]))
             h.dd <- ifelse(x[2]/h.dd > rho.max, x[2]/rho.max, h.dd)
             
           },
           # "maxexp" = { 
           #   h.dd <- x[1]/(1 + (x[3] * g * ts)/eta.null * exp(-k.maxexp * x[2]/x[1])/(rho.max - x[2]/x[1]))
           #   
           # },
           {
             stop(paste0("'",comptype,"'"," is no valid type for snow layer compaction"))
           }
    )
    h.dd <- ifelse(x[1] == 0, 0, h.dd)
    swe.dd  <- x[2]
    age.dd  <- ifelse(x[1] == 0, 0, age.d + 1)
    rho.dd  <- ifelse(x[1] == 0, 0, swe.dd/h.dd)
    rho.dd  <- ifelse(rho.max - rho.dd < prec, rho.max, rho.dd)
    return(cbind(h = h.dd, swe = swe.dd, age = age.dd, rho = rho.dd))
  } 
  
  
  
  #-------------------------------------------------------------------------
  # assigns snowpack properties of tomorrow
  #-------------------------------------------------------------------------
  assignH <- function(sp.dd, h, swe, age, H, SWE, t, day.tot){
    if(t < day.tot){
      h[,t+1]   <- sp.dd$h
      swe[,t+1] <- sp.dd$swe
      age[,t+1] <- sp.dd$age
      
      H[t+1]   <- sum(h[,t+1])
      SWE[t+1] <- sum(swe[,t+1])
    }
    return(list(h=h, swe=swe, age=age, H=H, SWE=SWE))
  }
  
  
  
  drenchH <- function(t, ly, ly.tot, day.tot, Hobs, h, swe, age, H, SWE, rho.max, c.ov, k.ov, k.exp, ts, prec, comptype, m){
    Hobs.d=Hobs[t]; h.d=h[,t]; swe.d=swe[,t]; age.d=age[,t]
    
    # if rain is available: event <- "rain"
    event <- "melt"
    
    switch(event,
           melt = {
             
             msg(m,t,paste("melt "))
             
             runoff <- 0
             
             for(i in ly:1){
               
               if( sum(h.d[-i]) + swe.d[i]/rho.max - Hobs.d >= prec ){
                 h.d[i] <- swe.d[i]/rho.max 
               } else {
                 h.d[i] <- swe.d[i]/rho.max + abs(sum(h.d[-i]) + swe.d[i]/rho.max - Hobs.d)
                 break
               }
               
             }
             
             # all layers have rho.max
             if( all(rho.max - swe.d[1:ly]/h.d[1:ly] <= prec) ){
               msg(m,t,paste("no further compaction "))
               
               # produce runoff if sum(h.d) - Hobs.d is still > 0
               #if ( sum(h.d) - Hobs.d > prec ){
               msg(m,t,paste("runoff "))
               # decrease swe from all layers?
               # or beginning with lowest?
               #swe.d[1:ly] <- swe.d[1:ly] - (sum(h.d) - Hobs.d) * rho.max
               scale <- Hobs.d/sum(h.d)
               runoff <- (sum(h.d) - Hobs.d) * rho.max  # excess is converted to runoff [kg/m2]
               h.d <- h.d * scale                       # all layers are compressed (and have rho.max) [m]
               swe.d <- swe.d * scale
               
               # }
               
               
             } else {
               msg(m,t,paste("compaction "))
             }
             
             h[,t]   <- h.d
             swe[,t] <- swe.d
             age[,t] <- age.d
             H[t]    <- sum(h[,t])
             SWE[t]  <- sum(swe[,t])
             
             # no further compaction possible
             #snowpack.tomorrow <- cbind(h = h.d, swe = swe.d, age = age.d, rho = swe.d/h.d)
             #colnames(snowpack.tomorrow) <- c("h","swe","age","rho")
             snowpack.tomorrow <- compactH(h[,t], swe[,t], age[,t], ly.tot, ly, comptype, k.exp, rho.max, ts, prec)
             
             # set values for next day
             rl <- assignH(snowpack.tomorrow, h, swe, age, H, SWE, t, day.tot)
             h <- rl$h
             swe <- rl$swe
             age <- rl$age
             H <- rl$H
             SWE <- rl$SWE
             
             
           },
           # rain = {
           #   # daily rain amount is given as R.d
           #   # or
           #   # daily climatological probability of rain is given
           #   #  # rain amount after Bertle, 1966
           #   #  R.d <- 2.11 * SWE.d * (1 - Hobs.d/H.d)
           #   # ...
           #   #
           #   #----------------------------
           #   # model a rain on snow event 
           #   #----------------------------
           #   
           #   
           #   # compaction of all layers by rain overburden: sigma.r [0:1]
           #   sigma.r <- 1 - c.ov * R.d * exp(-k.ov * snowpack.dd$rho/(rho.max - snowpack.dd$rho))
           #   h.d     <- sigma.r * h.d
           #   
           #   # layerwise distribution of R.d
           #   # start with youngest layer ly, add as much R.d to that layer as it needs to reach rho.max
           #   # then go one layer deeper into the snowpack
           #   # repeat until at least one layer can not reach rho.max, then stop
           #   # OR 
           #   # all layers have reached rho.max, then stop and produce runoff with remaining R.d
           #   i <- ly
           #   while(R.d > 0){
           #     #if(R.d < 0)
           #     #  swe.d[i]  <- swe.d[i] + R.d
           #     #print(paste(i,R.d))
           #     delta.swe <- h.d[i] * rho.max - swe.d[i] # layer tolerates this swe amount to reach rho.max
           #     if(delta.swe > R.d){
           #       delta.swe <- R.d
           #       
           #     } 
           #     swe.d[i]  <- swe.d[i] + delta.swe
           #     R.d  <- R.d - delta.swe
           #     #if (i > 1) 
           #     i <- i - 1
           #     #msg(t,paste(" ",round(R.d,2)))
           #     if(i<=0 & R.d > 0){
           #       msg(m,t,paste("produce runoff from melted snow"))
           #       break
           #       
           #     }
           #     
           #     #stopifnot(i>0) 
           #   }
           #},
           {
             stop(paste0("event '",event,"' is not defined."))
           }
           
           
    )
    
    return(list(h=h, swe=swe, age=age, H=H, SWE=SWE))
  }
  
  
  
  scaleH <- function(t, ly, ly.tot, day.tot, deltaH, Hobs, h, swe, age, H, SWE, rho.max, k.exp, ts, prec, comptype, m){
    Hobs.d=Hobs[t]; h.d=h[,t]; swe.d=swe[,t]; age.d=age[,t]; deltaH.d = deltaH
    
    if( deltaH.d >= 0 ){
      
      msg(m,t,paste("stretching within tau "))
      scale <- Hobs.d / sum(h.d)
      h.d <- h.d * scale
      
    } else {
      
      msg(m,t,paste("compaction within tau "))
      
      runoff <- 0
      
      # scale.min << 1 means that layer can be strongly compressed
      # scale.min -> 1 means that layer can be compressed only a little
      #scale.min <- swe.d/h.d / rho.max # density ratio [-]
      scale.min <- ifelse(rho.max - swe.d/h.d < prec, 1, swe.d/h.d / rho.max)
      
      # added-up height of all layers that already have maximum density 
      # and cannot be scaled (in this case: compressed) anymore [m]
      non.scalable.h <- sum(h.d[which(scale.min >= 1)])
      
      # scaling factor for all layers with densities lower rho.max, 
      # i.d. those that can still be scaled (compressed or streched) [-]
      # scale < 1 as long as modelled H is larger than observed Hobs 
      #scale <- (Hobs.d - non.scalable.h) / (sum(h.d) - non.scalable.h)
      #scale <- ifelse(sum(h.d) == non.scalable.h, 1, (Hobs.d - non.scalable.h) / (sum(h.d) - non.scalable.h))
      if(sum(h.d) == non.scalable.h){
        scale <- ifelse(Hobs.d > sum(h.d), Hobs.d / sum(h.d), sign(Hobs.d - non.scalable.h) * 1)
      } else {
        scale <- (Hobs.d - non.scalable.h) / (sum(h.d) - non.scalable.h)
      }
      
      
      # Check, if there is any layer that would have to be compressed more than it can be compressed.
      # (They are close to rho.max, but have not quite reached rho.max.)
      # This condition can only be fulfilled if compression must take place,
      # because streching (i.d. Hobs.d > sum(h.d)) necessitates "scale > 1" 
      # and - as "scale.min" always is between 0 and 1 -   # the condition is not fulfilled.)
      # If all layers can be compressed sufficiently OR if the layer must be streched: see "else"...  
      if ( any(na.omit(scale.min) > scale) ){
        
        # Loop as long as compressed model snow heights add up 
        # to more than observed snow height.
        # Exit loop when modelled heights all together equal observed heights.
        nr.comp <- 0
        #while ( sum(h.d) > Hobs.d ){
        while ( sum(h.d) - Hobs.d > prec ){
          
          #print(paste(nr.comp,": ",sum(h.d)," - ",Hobs.d))
          
          if(nr.comp>20)
            stop("error in while")
          
          # scale becomes zero or negative, when all layers, that cannot be scaled since they are too dense already,
          # add up to more than the observed snow height. (Although if there are some other layers that can still be compressed.)
          # In this case runoff muss be created.
          if (scale <= 0){
            idx <- which(scale.min < 1)              # index for those layers that still can be compressed []
            h.d[idx]  <- h.d[idx] * scale.min[idx]   # Layers that still can be compressed are compressed to maximum density.
            # Now we end up with modelled layers all at maximum density, and and their heights add up to more than the observed.
            
            scale <- Hobs.d/sum(h.d)                 # scaling factor (between 0 and 1) to force modelled sum(h) to observed Hobs
            msg(m,t,paste("runoff "))
            runoff <- (sum(h.d) - Hobs.d) * rho.max  # excess is converted to runoff [kg/m2]
            h.d <- h.d * scale                       # all layers are compressed (and have rho.max) [m]
            swe.d <- swe.d * scale                   # all layer's SWEs are reduced proportionally; "new sum of SWE" equals "old sum of SWE + runoff" [kg/m2]
            
            
            # The following condition means scale > 0, and as we are in the "if-loop" scale < 1 (see above). 
            # Therefore, scale is [0-1] for this "else"-condition.
            # This means there are still layers that can be compressed 
            # AND the incompressible do not add up to more than observed.
            # So let's compress and see...
          } else {
            msg(m,t,paste("stepwise compaction "))
            idx <- which(scale.min >= scale)         # index for those layers that cannot be compressed so much anymore []
            h.d[idx]  <- h.d[idx]  * scale.min[idx]  # Layers that hardly can be compressed anymore are compressed to maximum density.
            h.d[-idx] <- h.d[-idx] * scale           # Layers that can be compressed as much as they should are compressed.
            # Now we end up with modelled layers at maximum density and others that aren't so dense.
            # Still, there might be some more compression necessary, so let's calculate the new references.
            
            # new minimum values by which a layer can now be compressed. (More and more get a 1 here.)
            #scale.min <- swe.d/h.d / rho.max
            scale.min <- ifelse(rho.max - swe.d/h.d < prec, 1, swe.d/h.d / rho.max)
            
            # added-up height of all layers that already have maximum density and 
            # cannot be compressed anymore [m]. (This value increases with every loop-iteration)
            non.scalable.h <- sum(h.d[which(scale.min >= 1)]) 
            
            # scaling factor to do the loop again until modelled snow height equals observed snow height
            scale <- (Hobs.d - non.scalable.h) / (sum(h.d) - non.scalable.h)
            
            # if(sum(h.d) == non.scalable.h){
            #   scale <- sign(Hobs.d - non.scalable.h) * 1
            # } else {
            #   scale <- (Hobs.d - non.scalable.h) / (sum(h.d) - non.scalable.h)
            # }
            
            
          }
          nr.comp <- nr.comp + 1
          
          
        }
        #msg(t,paste0(" (",nr.comp," steps)"))
        
        # If all layers can be compressed sufficiently 
        # OR have to be streched since Hobs.d > sum(h.d)
      } else {
        
        # scale layers to observed snow depth
        # i.e. compress or stretch depending on scale
        if(verbose){
          txt <- ifelse(scale > 1, "stretching within tau", "compaction within tau")
          msg(m,t,paste(txt))
        }
        
        h.d <- h.d * scale # "new" layer's heights; afterwards sum(h.d) = Hobs.d!
        #idx <- which(scale.min < 1) 
        #h.d[idx] <- h.d[idx] * scale  
        
      } 
      
    }
    
    h[,t]   <- h.d
    swe[,t] <- swe.d
    age[,t] <- age.d
    H[t]    <- sum(h[,t])
    SWE[t]  <- sum(swe[,t])
    
    # compact actual day 
    # if all layers already have maximum density rho.max
    # the snowpack will not be changed by the following step
    snowpack.tomorrow <- compactH(h[,t], swe[,t], age[,t], ly.tot, ly, comptype, k.exp, rho.max, ts, prec)
    
    # set values for next day
    rl <- assignH(snowpack.tomorrow, h, swe, age, H, SWE, t, day.tot)
    h <- rl$h
    swe <- rl$swe
    age <- rl$age
    H <- rl$H
    SWE <- rl$SWE
    
    return(list(h=h, swe=swe, age=age, H=H, SWE=SWE))
    
  }
  
  
  
  #-------------------------------------------------------------------------
  # keep track of messages in a vector
  #-------------------------------------------------------------------------
  msg <- function(m,t,strg){
    if(verbose){
      cat(paste(strg))
      if(is.null(m[t])){
        m[t] <- strg  
      } else {
        m[t] <- paste(m[t],strg)
      }
    }  
  }
  
  
  
  
  
  
  for (t in 1:day.tot) {
    
    msg(m,t,paste("day",t,": "))
    
   
    # snowdepth = 0, no snow cover
    if( Hobs[t] == 0 ){                                          
      H[t]   <- 0
      SWE[t] <- 0
      #RHO[t] <- 0
      
      # there is snow
    } else if( Hobs[t] > 0 ){
      
      # first snow in/during season
      if( Hobs[t-1] == 0 ){
        ly <- 1
        msg(m,t,paste("produce layer",ly))
        age[ly,t]     <- 1
        h[ly,t]       <- Hobs[t]   
        H[t]          <- Hobs[t]
        swe[ly,t]     <- rho.null * Hobs[t]
        SWE[t]        <- swe[ly,t]
        #RHO[t]        <- SWE[t]/H[t]
        
        # compact actual day 
        snowpack.tomorrow <- compactH(h[,t], swe[,t], age[,t], ly.tot, ly, comptype, k.exp, rho.max, ts, prec)
        
        # set values for next day
        rl <- assignH(snowpack.tomorrow, h, swe, age, H, SWE, t, day.tot)
        h <- rl$h
        swe <- rl$swe
        age <- rl$age
        H <- rl$H
        SWE <- rl$SWE
        
        
        # non-first day of snow covered period
      } else if ( Hobs[t-1] > 0 ){
        
        deltaH <- Hobs[t] - H[t]
        
        if( deltaH > tau ){
          msg(m,t,paste("create new layer",ly+1))

          sigma.s <- 1 - c.ov * (deltaH) * rho.null * exp(-k.ov * snowpack.dd$rho/(rho.max - snowpack.dd$rho))
          h[,t]     <- sigma.s * h[,t]
          swe[,t]   <- swe[,t-1]
          age[(1:ly),t]   <- age[(1:ly),t-1] + 1
          H[t]      <- sum(h[,t])
          SWE[t]    <- sum(swe[,t])
          #RHO[t]    <- SWE[t]/H[t]
          
          # only for new layer
          ly            <- ly + 1
          h[ly,t]       <- Hobs[t] - H[t]
          swe[ly,t]     <- rho.null * h[ly,t]
          age[ly,t]     <- 1
          
          # recompute
          H[t]   <- sum(h[,t])
          SWE[t] <- sum(swe[,t])
          #RHO[t] <- SWE[t]/H[t]
          
          # compact actual day 
          #print(paste(t,H[t]))
          snowpack.tomorrow <- compactH(h[,t], swe[,t], age[,t], ly.tot, ly, comptype, k.exp, rho.max, ts, prec)
          
          # set values for next day
          rl <- assignH(snowpack.tomorrow, h, swe, age, H, SWE, t, day.tot)
          h <- rl$h
          swe <- rl$swe
          age <- rl$age
          H <- rl$H
          SWE <- rl$SWE
          
          
          # no mass gain or loss, but scaling
        } else if( deltaH >= -tau & deltaH <= tau ) {
          msg(m,t,paste("scaling: "))
          #scaleH(Hobs[t], h[,t], swe[,t], age[,t], deltaH)
          rl <- scaleH(t, ly, ly.tot, day.tot, deltaH, Hobs, h, swe, age, H, SWE, rho.max, k.exp, ts, prec, comptype, m)
          h <- rl$h
          swe <- rl$swe
          age <- rl$age
          H <- rl$H
          SWE <- rl$SWE
  
          
          
        } else if ( deltaH < -tau ){
          
          msg(m,t,paste("drenching: "))
          #drenchH(Hobs[t], h[,t], swe[,t], age[,t])
          rl <- drenchH(t, ly, ly.tot, day.tot, Hobs, h, swe, age, H, SWE, rho.max, c.ov, k.ov, k.exp, ts, prec, comptype, m)
          h <- rl$h
          swe <- rl$swe
          age <- rl$age
          H <- rl$H
          SWE <- rl$SWE
          
          
        } else {
          msg(m,t,"??")
        }
        
      }
      
      
    } 
    msg(m,t,"\n")
  }
  
  return(SWE)
  
}