Thornthwaite-Mather water balance

bh.R

# balanco hídrico
# Baseado em: 
# Pereira, Antonio Roberto. (2005). Simplificado o balanço hídrico de Thornthwaite-Mather.
# Bragantia, 64(2), 311-313. https://dx.doi.org/10.1590/S0006-87052005000200019

# Posse, GO
# https://edisciplinas.usp.br/pluginfile.php/3433083/mod_resource/content/1/Meteoro_aula11-Balanco_Hidrico.pdf
dados <- data.frame(p = c(271, 215, 230, 119, 20, 9, 5, 12, 30, 123, 223, 280),
                    etp = c(116, 97, 104, 88, 78, 63, 62, 90, 94, 109, 106, 106))
# campina Grande
# Pereira (2005)
dados2 <- data.frame(p = c(41, 55, 100, 129, 95, 107, 124, 58, 38, 17, 19, 21),
                     etp = c(108, 109, 115, 107, 95, 80, 62, 78, 77, 102, 108, 117))


# BH Climatológico
# Roda até estabilizar, com erro menor que 1 mm
# ou seja, último passo de tempo + 1 ≃ primeiro passo de tempo
# pode não ser adequado para rodadas diárias
# daí é melhor usar o sequencial
# espera um data frame com duas colunas
# p: precipitacao (mm)
# etp: evapotranspiração potencial (mm)
# cad: capacidade de água disponível
calc_bh_clim <- function(dados, cad = 100) {
  obs <- length(dados$p)
  dados$p_etp <- dados$p - dados$etp
  
  dados$arm[1] <- cad
  dados$alt <- -99
  dados$etr <- -99
  dados$def <- -99
  dados$exd <- -99
  
  # rodada inicial, para preencher arm
  for (m in 2:obs) {
    if (dados$p_etp[m] < 0) {
      dados$arm[m] <- min(cad, dados$arm[m-1] * exp(dados$p_etp[m]/cad))
    } else {
      dados$arm[m] <- min(cad, dados$arm[m-1] + dados$p_etp[m])
    }
  }
  
  # calc 13o mes ~ semelhante ao início
  arm13 <- ifelse(dados$p_etp[1] < 0,
                  min(cad, dados$arm[obs] * exp(dados$p_etp[1]/cad)),
                  min(cad, dados$arm[obs] + dados$p_etp[1])
                  )
  
  # iteracao até estabilizar
  count <- 1
  while (abs(dados$arm[1] - arm13) > 1) {
    count <- count + 1
    dados$arm[1] <- arm13
    for (m in 2:obs) {
      if (dados$p_etp[m] < 0) {
        dados$arm[m] <- min(cad, dados$arm[m-1] * exp(dados$p_etp[m]/cad))
      } else {
        dados$arm[m] <- min(cad, dados$arm[m-1] + dados$p_etp[m])
      }
    }
    arm13 <- ifelse(dados$p_etp[1] < 0,
                    min(cad, dados$arm[obs] * exp(dados$p_etp[1]/cad)),
                    min(cad, dados$arm[obs] + dados$p_etp[1])
    )
  }
  
  cat(paste("Armazenamento estabilizado com", count, "rodadas"))
  
  # finalizando o BH
  for (m in 1:obs) {
    if (m == 1) {
      dados$alt[m] <- dados$arm[m] - dados$arm[obs]
    } else{
      dados$alt[m] <- dados$arm[m] - dados$arm[m-1]
    }
  }
  
  dados$etr <- dados$etp
  dados$etr[dados$p_etp < 0] <- dados$p[dados$p_etp < 0] + abs(dados$alt[dados$p_etp < 0])
  
  dados$def <- dados$etp - dados$etr
  
  dados$exd <- 0
  dados$exd[dados$arm == cad] <- dados$p_etp[dados$arm == cad] - dados$alt[dados$arm == cad]
  
  return(dados)
}


# BH Sequencial
# Inicia rodada com CAD a um determinado percentual
# melhor é descartar primeiro ano para dar tempo de estabilizar o armazenamento
# espera um data frame com duas colunas
#     p: precipitacao (mm)
#     etp: evapotranspiração potencial (mm)
# cad: capacidade de água disponível
# pcad: porcentual da cad completa no primeiro passo de tempo
calc_bh_seq <- function(dados, cad = 100, pcad = 0.7) {
  obs <- length(dados$p)
  dados$p_etp <- dados$p - dados$etp
  
  dados$arm[1] <- cad * pcad
  dados$alt <- -99
  dados$etr <- -99
  dados$def <- -99
  dados$exd <- -99
  
  for (m in 2:obs) {
    if (dados$p_etp[m] < 0) {
      dados$arm[m] <- min(cad, dados$arm[m-1] * exp(dados$p_etp[m]/cad))
    } else {
      dados$arm[m] <- min(cad, dados$arm[m-1] + dados$p_etp[m])
    }
  }
  
  # finalizando o BH
  for (m in 1:obs) {
    if (m == 1) {
      dados$alt[m] <- dados$arm[m] - dados$arm[obs]
    } else{
      dados$alt[m] <- dados$arm[m] - dados$arm[m-1]
    }
  }
  
  dados$etr <- dados$etp
  dados$etr[dados$p_etp < 0] <- dados$p[dados$p_etp < 0] + abs(dados$alt[dados$p_etp < 0])
  
  dados$def <- dados$etp - dados$etr
  
  dados$exd <- 0
  dados$exd[dados$arm == cad] <- dados$p_etp[dados$arm == cad] - dados$alt[dados$arm == cad]
  
  return(dados)
}

calc_bh_clim(dados)

calc_bh_clim(dados2, cad = 125)