1 Introdução e Contexto

Este Documento registra a evolução metodológica da pesquisa de Iniciação Científica (PIBIC 2025/2026) sobre estimativa de vazões na Bacia do Alto Paraguai usando sensoriamento remoto.

O Pantanal, maior área úmida tropical do mundo, integra a Bacia do Alto Paraguai e apresenta dinâmica hidrológica marcada pela sazonalidade da precipitação e por pulsos de inundação com atraso de dois a seis meses em relação ao pico chuvoso (Hamilton, Sippel, e Melack 1996). A rede fluviométrica é esparsa, tornando atraente o uso de satélites para estimar vazão.

Este trabalho propõe e avalia modelos alternativos de estimativa de vazão baseados em produtos de sensoriamento remoto, utilizando a relação entre as Anomalias de Armazenamento Terrestre de Água (TWSA) do satélite GRACE e a descarga medida em Porto Murtinho.

Pergunta da Pesquisa: É possível estimar a vazão mensal em Porto Murtinho usando produtos gerados de satélite?

2 Linha do Tempo

Etapa 1
Balanço Hídrico Tradicional — \(Q = P - ET - \Delta S\) com GRACE, CHIRPS/ERA5 e MOD16/SSEBop via GEE.
Etapa 2
Fechamento impossível — Sentido inverso funcionou, mas usar \(\Delta S\) do GRACE para estimar Q não funcionou.
Etapa 3
Relação exponencial — \(Q = \alpha \cdot e^{\beta \cdot TWSA}\) (Duvvuri e Beighley 2023). R² ≈ 0,48 — promissor mas limitado para a região umida do pantanal (fatores P e ET não são considerados no modelo).
Etapa 4
Modelos estendidos — Incorporação de P (ERA5) e ET (SSEBop). Melhora moderada no modelo de (Duvvuri e Beighley 2023).
Etapa 5 - Atual
Análise de montante — Separação da vazão incremental do Pantanal descontando sub-bacias do Planalto com lag por CCF.

3 Dados Utilizados

3.1 Produtos de Sensoriamento Remoto

Todos foram adquiridos via Google Earth Engine (GEE), extraídos como médias espaciais da bacia, exportados como séries temporais mensais em .csv.

3.1.1 Precipitação

3.1.1.1 CHIRPS v2.0 — Climate Hazards Group InfraRed Precipitation with Station data

Combina estimativas de temperatura de topo de nuvens (satélites geoestacionários) com dados de pluviômetros terrestres. Alta resolução espaço-temporal, muito utilizado em bacias tropicais sul-americanas. Os valores diários são somados para totais mensais em mm/mês. Desenvolvido pelo grupo CHG da UC Santa Barbara.

3.1.1.2 ERA5-Land — ECMWF Reanalysis v5 Land

Reanálise atmosférica do Centro Europeu de Previsão do Tempo (ECMWF). Combina dados observados com modelos físicos de circulação atmosférica para reconstruir o histórico atmosférico global. Os dados de precipitação vêm em metros e são multiplicados por 1000 para converter em milímetros.

3.1.2 Evapotranspiração

3.1.2.1 MOD16 — MODIS Global Evapotranspiration

Baseado na equação de Penman-Monteith, que estima a ET a partir de energia disponível na superfície, resistência estomática da vegetação e déficit de pressão de vapor. Utiliza dados MODIS de albedo, índice de área foliar (LAI) e temperatura de superfície. Os dados têm fator de escala de 0,1 — multiplicado para obter mm/mês.

3.1.2.2 SSEBop — Simplified Surface Energy Balance (Operational)

Modelo operacional da USGS baseado no balanço simplificado de energia na superfície. Usa temperatura de superfície terrestre (LST) do satélite para estimar a fração evaporativa em relação à evapotranspiração potencial. Dois sensores foram concatenados: MODIS (Terra) até 2017 e VIIRS (Suomi-NPP) de 2018 em diante. A sobreposição entre 2012 e 2023 foi harmonizada por média mensal. Já vem em mm/mês, sem necessidade de fator de escala adicional.

3.1.3 Armazenamento Hídrico Terrestre

3.1.3.1 GRACE Mascon v04 — Gravity Recovery and Climate Experiment (JPL RL06)

O GRACE e sua sucessora GRACE-FO medem variações no campo gravitacional terrestre causadas por redistribuição de massa de água. A Anomalia de Armazenamento Terrestre de Água (TWSA) representa a diferença entre o armazenamento atual e uma linha de base. TWSA positiva = mais água que o usual; negativa = menos. O produto integra toda a água: solo, aquíferos, rios, lagos. Os dados vêm em centímetros de equivalente de água e são multiplicados por 10 para obter milímetros. Há lacuna de ~11 meses entre 2017–2018 (transição GRACE → GRACE-FO), tratada excluindo meses sem dados.

3.2 Dados Fluviométricos (ANA — Hidroweb)

Estações fluviométricas utilizadas.
Estacao Codigo Rio Area km² Papel
Porto Murtinho 67100000 Paraguai 576.000 Variável resposta
Cáceres 66070004 Paraguai 32.400 Tributária
Cuiabá 66260002 Cuiabá 22.800 Tributária
Ac. Corumbá Grande 66460000 Taquari 23.000 Tributária
Coxim 66870000 Taquari 27.600 Tributária
Aquidauana 66945000 Aquidauana 15.700 Tributária
Miranda 66910000 Miranda 15.000 Tributária

Conversão de unidades: para realizar a comparação de vazão (m³/s) com P e ET (mm/mês) eu aplico a seguinte conversão nos dados de vazão: \[Q_{mm/mês} = \frac{Q_{m^3/s} \times 86{,}4 \times N_{dias}}{A_{bacia\,(km^2)}}\]

onde \(A_{bacia} = 576.000 \text{ km}^2\) (representa toda a área de drenagem em porto murtinho, os dados do balanço hidrico foram considerados apenas dentro da região do bioma do pantanal, porém utiliza a área da bacia completa para fazer a conversão da vazão observada)e \(N_{dias}\) é o número de dias do mês.

4 Séries Históricas das Variáveis de Entrada

# ==============================================================================
# Leitura dos Dados Brutos
# ==============================================================================

chirps_raw   <- read_csv("Dados/P/Serie_Historica_CHIRPS_Pantanal.csv", show_col_types = FALSE)
era5_raw     <- read_csv("Dados/P/Serie_Historica_ERA5_Land_Precipitacao_Mensal_Pantanal.csv", show_col_types = FALSE)
mod16_raw    <- read_csv("Dados/ET/Serie_Historica_MOD16_Pantanal.csv", show_col_types = FALSE)
ssebop_modis <- read_csv("Dados/ET/Serie_Historica_SSEBop_MODIS_Mensal_Pantanal.csv", show_col_types = FALSE)
ssebop_viirs <- read_csv("Dados/ET/Serie_Historica_SSEBop_VIIRS_Mensal_Pantanal.csv", show_col_types = FALSE)
grace_raw    <- read_csv("Dados/Serie_Historica_GRACE_Mascon_V04_Pantanal.csv", show_col_types = FALSE)
vazao_raw    <- read_excel("Dados/Estações de Sedes Municipais/Estacao_67100000.xlsx")

# ==============================================================================
# Processamento e Agregação Mensal
# ==============================================================================

# Precipitação: CHIRPS
# Agrega os dados diários para o total da precipitação mensal (mm/mês)
chirps_m <- chirps_raw %>% 
  mutate(ano_mes = floor_date(as.Date(data), "month")) %>%
  group_by(ano_mes) %>% 
  summarise(P_CHIRPS = sum(valor, na.rm = TRUE))

# Precipitação: ERA5-Land
# O dado nativo vem em metros. Multiplicamos por 1000 para converter para mm e somamos para obter o total mensal.
era5_m <- era5_raw %>% 
  mutate(ano_mes = floor_date(as.Date(data), "month")) %>%
  group_by(ano_mes) %>% 
  summarise(P_ERA5 = sum(valor, na.rm = TRUE) * 1000)

# Evapotranspiração: MOD16
# Possui fator de escala nativo de 0.1. Multiplicamos por 0.1 para converter para mm e somamos para o total mensal.
mod16_m <- mod16_raw %>% 
  mutate(ano_mes = floor_date(as.Date(data), "month")) %>%
  group_by(ano_mes) %>% 
  summarise(ET_MOD16 = sum(valor, na.rm = TRUE) * 0.1)

# Evapotranspiração: SSEBop
# Empilha os sensores MODIS e VIIRS. A sobreposição entre eles é harmonizada extraindo a média (mean) do valor para aquele mês. O dado já está em mm.
ssebop_m <- bind_rows(ssebop_modis, ssebop_viirs) %>%
  mutate(ano_mes = floor_date(as.Date(data), "month")) %>%
  group_by(ano_mes) %>% 
  summarise(ET_SSEBop = mean(valor, na.rm = TRUE))

# Armazenamento Hídrico: GRACE (TWSA e dS)
# O dado vem em cm equivalentes de água. Multiplicamos por 10 para mm.
# A variação mensal (dS) é a diferença entre a TWSA do mês atual e do mês anterior.
grace_m <- grace_raw %>% 
  mutate(ano_mes = floor_date(as.Date(data), "month")) %>%
  group_by(ano_mes) %>% 
  summarise(TWSA = mean(valor, na.rm = TRUE) * 10) %>%
  arrange(ano_mes) %>% 
  mutate(dS = TWSA - lag(TWSA))

# Vazão Observada (Porto Murtinho)
# Converte a vazão média (m³/s) para lâmina d'água (mm/mês) distribuída na bacia.
# Fórmula: Q_mm = (Q_m3s * 86.4 * n_dias) / area_bacia_km2
vazao_m <- vazao_raw %>%
  mutate(
    DataHora = as.Date(DataHora),
    ano_mes  = floor_date(DataHora, "month")
  ) %>%
  group_by(ano_mes) %>%
  summarise(
    Q_obs_m3s = mean(Vazao, na.rm = TRUE),
    dias      = days_in_month(first(ano_mes))
  ) %>%
  mutate(Q_obs = (Q_obs_m3s * 86.4 * dias) / area_bacia) %>% 
  select(ano_mes, Q_obs)

4.1 Comparação entre produtos de precipitação

# ==============================================================================
# Gráfico de Comparação de Precipitação (CHIRPS vs ERA5-Land)
# ==============================================================================

p_p <- left_join(chirps_m, era5_m, by = "ano_mes") %>%
  pivot_longer(
    cols = c(P_CHIRPS, P_ERA5), 
    names_to = "Produto", 
    values_to = "P"
  ) %>%
  mutate(Produto = recode(Produto, P_CHIRPS = "CHIRPS", P_ERA5 = "ERA5-Land")) %>%
  
  ggplot(aes(x = ano_mes, y = P, color = Produto)) +
  
  geom_line(linewidth = 0.8, alpha = 0.85) +
  scale_color_manual(values = c("CHIRPS" = COR_P, "ERA5-Land" = COR_ET)) +
  theme_minimal(base_size = 11) + 
  theme(legend.position = "bottom") +
  labs(
    x = NULL, 
    y = "Precipitação (mm/mês)", 
    color = NULL
  )

ggplotly(p_p, tooltip = c("x", "y", "colour")) %>% 
  layout(legend = list(orientation = "h", y = -0.2))

Série histórica de precipitação mensal.

Como avaliar qual produto é melhor?:

CHIRPS e ERA5-Land mostram boa concordância sazonal, com pequenas diferenças nos valores absolutos. O qye pode gerar duvida na hora de selecionar qual produto vale mais a pena ser utilizado. será se compara com pluviômetros na região?

4.2 Comparação entre produtos de evapotranspiração

# ==============================================================================
# Gráfico de Comparação de Evapotranspiração (MOD16 vs SSEBop)
# ==============================================================================

p_et <- left_join(mod16_m, ssebop_m, by = "ano_mes") %>%
  pivot_longer(
    cols = c(ET_MOD16, ET_SSEBop), 
    names_to = "Produto", 
    values_to = "ET"
  ) %>%
  mutate(
    Produto = recode(Produto, ET_MOD16 = "MOD16", ET_SSEBop = "SSEBop (MODIS+VIIRS)")
  ) %>%
  
  ggplot(aes(x = ano_mes, y = ET, color = Produto)) +
  geom_line(linewidth = 0.8, alpha = 0.85) +
  scale_color_manual(values = c("MOD16" = COR_ET, "SSEBop (MODIS+VIIRS)" = COR_Q_M1)) +
  theme_minimal(base_size = 11) + 
  theme(legend.position = "bottom") +
  labs(
    x = NULL, 
    y = "ET (mm/mês)", 
    color = NULL
  )

ggplotly(p_et, tooltip = c("x", "y", "colour")) %>% 
  layout(legend = list(orientation = "h", y = -0.2))

Série histórica de ET mensal.

MOD16 e SSEBop (MODIS + VIIRS) apresentam disparidades até 2009 (MOD16 subestima a evapotranspiração). É perceptivel que o produto SSEBop possui valores medios dentro de uma faixa entre 50 e 100 mm/mês, enquanto MOD16 possui minimos de evapotranspiração mais baixos. Neste caso optei por selecionar o SSEBop como o melhor por se tratar da média de 2 produtos, mas posso estar completamente equivocado, já que a média pode justamente eliminar picos de maxima ou minima.

4.3 GRACE — Anomalia de armazenamento

# ==============================================================================
# Gráfico da Anomalia de Armazenamento Terrestre de Água (TWSA - GRACE)
# ==============================================================================

grace_m_completo <- grace_m %>%
  complete(ano_mes = seq(min(ano_mes), max(ano_mes), by = "month"))

# Inicializa o gráfico com a base contendo os NAs explícitos para visualização da quebra entre os períodos GRACE e GRACE-FO
p_grace <- ggplot(grace_m_completo, aes(x = ano_mes)) +
  
  geom_line(aes(y = TWSA), color = COR_TWSA, linewidth = 0.9) +
  geom_hline(yintercept = 0, linetype = "dashed", color = "gray50") +
  theme_minimal(base_size = 11) + 
  labs(
    x = NULL, 
    y = "TWSA (mm)"
  )

ggplotly(p_grace) %>% 
  layout(hovermode = "x unified")

Série histórica de TWSA mensal.

Lacuna 2017–2018 = transição GRACE → GRACE-FO. A série mostra algumas falhas também antes de 2017, talvez seja interessante utilizar series sintéticas do GRACE, procurar produtos corrigidos.

4.4 Vazão observada em Porto Murtinho

# ==============================================================================
# Gráfico da Série Histórica de Vazão Observada (Porto Murtinho)
# ==============================================================================

p_vaz <- ggplot(vazao_m, aes(x = ano_mes, y = Q_obs)) +
  
  geom_line(color = COR_OBS, linewidth = 0.9) +
  
  theme_minimal(base_size = 11) + 
  labs(
    x = NULL, 
    y = "Q observada (mm/mês)"
  )

ggplotly(p_vaz) %>% 
  layout(hovermode = "x unified")

Série histórica de vazão mensal em Porto Murtinho (mm/mês).

5 Etapa 1: O Balanço Hídrico Tradicional

5.1 A equação e sua lógica

\[\Delta S = P - ET - Q \quad \Longrightarrow \quad Q = P - ET - \Delta S\]

A lógica é simples: tudo que entra como chuva e não evapora nem fica armazenado deve sair como escoamento. O GRACE mede \(\Delta S\); se tivéssemos \(P\) e \(ET\) sem erros teríamos \(Q\) estimado igual ao observado.

5.2 O problema: propagação de erros

Moreira et al. (2019) avaliaram essa abordagem para o Alto Paraguai e encontraram correlação de apenas r = 0,26 e RMSE de 279%. \(P\) e \(ET\) de satélite carregam erros de 30–50%, e ao subtraí-los o erro resultante em \(Q\) domina completamente o sinal.

# ==============================================================================
# União dos Dados e Cálculo do Balanço Hídrico
# ==============================================================================

# Juntando as séries temporais dos produtos e da vazão observada,
# removendo valores ausentes (drop_na) e ordenando cronologicamente.
#
# Em seguida, calculamos o balanço hídrico em dois sentidos:
# 1. Estimando Q como resíduo (Q = P - ET - dS)
# 2. Estimando a variação de armazenamento (dS) como resíduo.

dados_bh <- era5_m %>%
  inner_join(ssebop_m, by = "ano_mes") %>% 
  inner_join(grace_m, by = "ano_mes") %>%
  inner_join(vazao_m, by = "ano_mes") %>% 
  drop_na() %>% 
  arrange(ano_mes) %>%
  mutate(
    Q_balanco  = dplyr::lag(P_ERA5, 0) - dplyr::lag(ET_SSEBop, 0) - dS,
    dS_balanco = dplyr::lag(P_ERA5, 0) - dplyr::lag(ET_SSEBop, 0) - Q_obs
  )

# ==============================================================================
# Avaliação de Desempenho (R²)
# ==============================================================================

# Função para calcular o R² entre a série observada (o) e simulada (s).
r2_fn <- function(o, s) {
  df <- data.frame(o = o, s = s) %>% 
    drop_na()
  
  SST <- sum((df$o - mean(df$o))^2) # Soma dos Quadrados Totais
  SSR <- sum((df$s - df$o)^2)       # Soma dos Quadrados dos Resíduos
  
  return(1 - (SSR / SST))
}

r2_q  <- round(r2_fn(dados_bh$Q_obs, dados_bh$Q_balanco), 3)
r2_dS <- round(r2_fn(dados_bh$dS, dados_bh$dS_balanco), 3)

rmse_q  <- round(sqrt(mean((dados_bh$Q_obs - dados_bh$Q_balanco)^2, na.rm = TRUE)), 3)
rmse_dS <- round(sqrt(mean((dados_bh$dS - dados_bh$dS_balanco)^2, na.rm = TRUE)), 3)

# print(paste0("R² para Q estimado: ", r2_q, "e R² para dS estimado: ", r2_dS))

5.3 Sentido 1 — Estimando Q pelo balanço

# ==============================================================================
# Gráfico do Balanço Hídrico: Sentido Q (Estimativa de Vazão)
# ==============================================================================

p1 <- ggplot(dados_bh, aes(x = ano_mes)) +
  
  geom_line(aes(y = Q_obs, color = "Q Observado"), linewidth = 1) +
  geom_line(aes(y = Q_balanco, color = "Q = P − ET − ΔS"), linewidth = 0.9) +
  scale_color_manual(values = c("Q Observado" = COR_OBS, "Q = P − ET − ΔS" = "#8e44ad")) +
  theme_minimal(base_size = 11) + 
  theme(legend.position = "bottom") +
  labs(
    x = NULL, 
    y = "Vazão (mm/mês)", 
    color = NULL,
    title = paste0("Q Estimado | R² = ", r2_q, " | RMSE = ", rmse_q)
  )

ggplotly(p1, tooltip = c("x", "y", "colour")) %>% 
  layout(
    legend = list(orientation = "h", y = -0.2),
    hovermode = "x unified"
  )

Balanço hídrico no sentido Q = P − ET − ΔS.

A estimativa não acompanha nem de perto a escala da vazão observada, o que pode indicar tanto possiveis erros na hora de aplicar o balanço quanto a propagação de erros presentes nas séries utilizadas.

5.4 Sentido 2 — Estimando \(\Delta S\) pelo balanço

# ==============================================================================
# Gráfico do Balanço Hídrico: Sentido dS (Variação de Armazenamento)
# ==============================================================================

p2 <- ggplot(dados_bh, aes(x = ano_mes)) +
  
  geom_line(aes(y = dS, color = "dS GRACE"), linewidth = 1) +
  geom_line(aes(y = dS_balanco, color = "dS = P − ET − Q_obs"), linewidth = 0.9) +
  scale_color_manual(values = c("dS GRACE" = COR_TWSA, "dS = P − ET − Q_obs" = COR_Q_M3)) +
  theme_minimal(base_size = 11) + 
  theme(legend.position = "bottom") +
  labs(
    x = NULL, 
    y = "dS (mm)", 
    color = NULL,
    title = paste0("Sentido dS | R² = ", r2_dS, " | RMSE = ", rmse_dS)
  )

ggplotly(p2, tooltip = c("x", "y", "colour")) %>% 
  layout(
    legend = list(orientation = "h", y = -0.2),
    hovermode = "x unified"
  )

Balanço hídrico no sentido dS = P − ET − Q_obs.

Quando a vazão observada é âncora, o dS estimado acompanha razoavelmente o GRACE, indicando que P e ET são coerentes — o problema está em estimar Q como resíduo.

Tentei aplicar CCF para encontrar o lag ótimo para P e ET, mas o resultado não melhorou o R² então retirei o trecho que escolhia o lag por Correlação cruzada no balanço hidrico tradicional, estranhamente os dados de P e ET sem lag possuem maior aderencia a serie (não são tão piores quanto aplicando o atraso de 2 a 4 meses que discutimos anteriormente).

Conclusão da Etapa 1: O balanço parece seguir a ordem de grandeza quando Q é âncora (R² ≈ 0.265), mas falha ao estimar a vazão como resíduo (R² ≈ -219.549). Isso motivou a busca por outro método.

6 Etapa 2: Relação Exponencial TWSA–Vazão

6.1 Fundamentação teórica

Kirchner (2009) demonstrou que quando o expoente \(b\) da lei de potência de recessão é 2, a solução da equação de sensibilidade resulta em uma relação exponencial entre armazenamento e descarga. Duvvuri e Beighley (2023) aplicaram isso com TWSA do GRACE em 2.870 estações nos EUA, com desempenho superior em 70% dos locais:

\[Q = \alpha \cdot e^{\beta \cdot TWSA}\]

Parâmetros:

  • α: descarga quando TWSA = 0 corresponde à vazão média anual da bacia.

  • β: sensibilidade da descarga ao armazenamento. β alto significa que o escoamento responde fortemente ao estado de armazenamento.

6.2 Calibração: linearização logarítmica com dados reais

\[\ln(Q) = \underbrace{\ln(\alpha)}_{\text{intercepto}} + \underbrace{\beta}_{\text{coeficiente}} \cdot TWSA\]

# ==============================================================================
# Preparação dos Dados
# ==============================================================================

dados_mb <- era5_m %>% 
  inner_join(grace_m, by = "ano_mes") %>%
  inner_join(ssebop_m, by = "ano_mes") %>% 
  inner_join(vazao_m, by = "ano_mes") %>%
  drop_na() %>% 
  arrange(ano_mes)

# Filtrando apenas meses com vazão estritamente maior que zero,
# passo necessário para viabilizar a transformação logarítmica (ln).
dados_fit_b <- dados_mb %>% 
  filter(Q_obs > 0)

# ==============================================================================
# Ajuste do Modelo Exponencial (M1: Base GRACE)
# ==============================================================================

# Ajustando o modelo linearizado: ln(Q) = ln(α) + β * TWSA
mod_b <- lm(log(Q_obs) ~ TWSA, data = dados_fit_b)

# Extraindo os parâmetros α e β.
# O parâmetro α (vazão de base) é obtido exponenciando o intercepto.
# O parâmetro β (sensibilidade) é o coeficiente associado ao TWSA.
a_r <- exp(coef(mod_b)[1])
b_r <- coef(mod_b)[2]

# R² do ajuste log-linear (arredondado para exibição nos gráficos)
r2_log <- round(summary(mod_b)$r.squared, 3)

# Gerando a série de vazão estimada ("bruta", sem lag).
# O predict retorna os valores logarítmicos, então usamos exp() para voltar à escala original (mm/mês).
dados_mb$Q_m1_raw <- exp(predict(mod_b, newdata = dados_mb))

# ==============================================================================
# Função de Otimização de Lag 
# ==============================================================================

# Função para encontrar o lag (atraso) ótimo entre a série observada 
# e a série simulada. O loop testa de 0 até 'maxl' meses, retornando o lag que maximiza a métrica R².
get_best_lag <- function(obs, sim_raw, maxl = 18) {
  
  best <- list(lag = 0, r2 = -Inf)
  
  for (l in 0:maxl) {
    # Aplica o atraso atual (l) na série simulada
    s <- dplyr::lag(as.numeric(sim_raw), l)
    
    # Cria um vetor lógico (v) para identificar os meses onde tem tanto o dado observado quanto o simulado
    v <- !is.na(obs) & !is.na(s)
    
    # Se houver menos de 5 pontos coincidentes, pula para o próximo lag
    if (sum(v) < 5) next
    
    # Cálculo do R² para o lag atual
    SST <- sum((obs[v] - mean(obs[v]))^2) # Soma dos Quadrados Totais
    SSR <- sum((s[v] - obs[v])^2)         # Soma dos Quadrados dos Resíduos
    r2  <- 1 - (SSR / SST)
    
    # 5. Se este R² for o maior até agora, salva na lista 'best'
    if (r2 > best$r2) {
      best <- list(lag = l, r2 = round(r2, 3))
    }
  }
  
  return(best)
}

# ==============================================================================
# Aplicação do Lag Ótimo
# ==============================================================================

# Calcula o lag ideal testando a vazão observada contra a bruta do M1
res1 <- get_best_lag(dados_mb$Q_obs, dados_mb$Q_m1_raw)

# Desloca a série bruta permanentemente usando o lag encontrado
dados_mb$Q_m1 <- dplyr::lag(dados_mb$Q_m1_raw, res1$lag)
# ==============================================================================
# Gráficos de Ajuste do Modelo: Original vs Log-Linear
# ==============================================================================

po <- ggplot(dados_fit_b, aes(x = TWSA, y = Q_obs)) +
  geom_point(color = COR_TWSA, alpha = 0.4, size = 1.5) +
  stat_function(
    fun = function(x) a_r * exp(b_r * x), 
    color = COR_Q_M3, 
    linewidth = 1.1
  ) +
  annotate(
    "text", 
    x = min(dados_fit_b$TWSA) + 5, 
    y = max(dados_fit_b$Q_obs) * 0.9,
    label = paste0("Q = ", round(a_r, 2), " · e^(", round(b_r, 4), " · TWSA)"),
    hjust = 0, 
    size = 3.2
  ) +
  
  theme_minimal(base_size = 10) +
  labs(
    title = "Escala original", 
    x = "TWSA (mm)", 
    y = "Q obs (mm/mês)"
  )

pl <- ggplot(dados_fit_b, aes(x = TWSA, y = log(Q_obs))) +
  geom_point(color = COR_TWSA, alpha = 0.4, size = 1.5) +
  
  geom_smooth(
    method = "lm", 
    color = COR_Q_M3, 
    linewidth = 1.1, 
    se = TRUE, 
    fill = "#f7c6c6"
  ) +
  
  annotate(
    "text", 
    x = min(dados_fit_b$TWSA) + 5, 
    y = max(log(dados_fit_b$Q_obs)) - 0.35,
    label = paste0("R² ajuste log-linear = ", r2_log), 
    hjust = 0, 
    size = 3.2
  ) +
  
  theme_minimal(base_size = 10) +
  labs(
    title = "Escala log-linear (Regressão OLS)", 
    x = "TWSA (mm)", 
    y = "ln(Q obs)"
  )

po + pl
Dados reais da pesquisa: escala original (esquerda) e log-linear (direita).

Dados reais da pesquisa: escala original (esquerda) e log-linear (direita).

7 Etapa 3: Modelos Estendidos com P e ET

Foi testada a inclusão de P e ET como preditores adicionais na regressão log-linear, resultando em quatro modelos para comparação:

7.1 Os quatro modelos

\[M1: \ln(Q) = \ln\alpha + \beta \cdot TWSA\] \[M2: \ln(Q) = \ln\alpha + \beta_1 TWSA + \beta_2 P\] \[M3: \ln(Q) = \ln\alpha + \beta_1 TWSA + \beta_2 ET\] \[M4: \ln(Q) = \ln\alpha + \beta_1 TWSA + \beta_2 P + \beta_3 ET\]

# ==============================================================================
# Ajuste dos Modelos Estendidos (M2, M3 e M4)
# ==============================================================================

# Usamos o mesmo dataframe filtrado (sem NAs e Q_obs > 0) da etapa anterior
df_f <- dados_fit_b

# M2: Adiciona Precipitação (ERA5) como preditor
mod2 <- lm(log(Q_obs) ~ TWSA + P_ERA5, data = df_f)

# M3: Adiciona Evapotranspiração (SSEBop) como preditor
mod3 <- lm(log(Q_obs) ~ TWSA + ET_SSEBop, data = df_f)

# M4: Modelo completo com TWSA, Precipitação e Evapotranspiração
mod4 <- lm(log(Q_obs) ~ TWSA + P_ERA5 + ET_SSEBop, data = df_f)

# ==============================================================================
# Geração das Séries Estimadas (Sem Lag)
# ==============================================================================

# O modelo prevê na escala logarítmica (ln). 
# Usamos exp() para retornar os valores à escala original de vazão (mm/mês).
dados_mb$Q_m2_raw <- exp(predict(mod2, newdata = dados_mb))
dados_mb$Q_m3_raw <- exp(predict(mod3, newdata = dados_mb))
dados_mb$Q_m4_raw <- exp(predict(mod4, newdata = dados_mb))

# ==============================================================================
# Otimização e Aplicação do Lag Ótimo
# ==============================================================================

# Encontra o melhor atraso (lag que maximiza o R²) para cada modelo específico
res2 <- get_best_lag(dados_mb$Q_obs, dados_mb$Q_m2_raw)
res3 <- get_best_lag(dados_mb$Q_obs, dados_mb$Q_m3_raw)
res4 <- get_best_lag(dados_mb$Q_obs, dados_mb$Q_m4_raw)

# Desloca as séries temporais permanentemente aplicando o lag encontrado
dados_mb$Q_m2 <- dplyr::lag(dados_mb$Q_m2_raw, res2$lag)
dados_mb$Q_m3 <- dplyr::lag(dados_mb$Q_m3_raw, res3$lag)
dados_mb$Q_m4 <- dplyr::lag(dados_mb$Q_m4_raw, res4$lag)

# ==============================================================================
# Cálculo do RMSE (Root Mean Square Error)
# ==============================================================================

# Avalia o erro absoluto médio em mm/mês para os quatro modelos.
# O 'na.rm = TRUE' garante que o cálculo utilize apenas o período de sobreposição válido.
rmse1 <- sqrt(mean((dados_mb$Q_obs - dados_mb$Q_m1)^2, na.rm = TRUE))
rmse2 <- sqrt(mean((dados_mb$Q_obs - dados_mb$Q_m2)^2, na.rm = TRUE))
rmse3 <- sqrt(mean((dados_mb$Q_obs - dados_mb$Q_m3)^2, na.rm = TRUE))
rmse4 <- sqrt(mean((dados_mb$Q_obs - dados_mb$Q_m4)^2, na.rm = TRUE))
Desempenho dos quatro modelos testados
Modelo R2 RMSE Lag
M1-TWSA 0.474 2.677 3
M2-TWSA+P 0.534 2.511 1
M3-TWSA+ET 0.529 2.525 1
M4-TWSA+P+ET 0.537 2.502 1

A inclusão de P e ET melhora o R², mas o ganho é modesto. O modelo M4 (com ambos os preditores) tem o melhor desempenho, mas a diferença em relação ao M2 e M3 é pequena. Isso sugere que a informação de P e ET pode estar sendo capturada indiretamente.

7.2 Comparação das séries dos quatro modelos

# ==============================================================================
# Gráfico de Comparação: Modelos Iniciais vs. Vazão Observada
# ==============================================================================

nms <- c(
  "Observado",
  paste0("M1-TWSA (R²=", res1$r2, ")"),
  paste0("M2-TWSA+P (R²=", res2$r2, ")"),
  paste0("M3-TWSA+ET (R²=", res3$r2, ")"),
  paste0("M4-TWSA+P+ET (R²=", res4$r2, ")")
)

pm_i <- ggplot(dados_mb, aes(x = ano_mes)) +
  
  geom_line(aes(y = Q_obs, color = nms[1]), linewidth = 1.1) +
  
  geom_line(aes(y = Q_m1, color = nms[2]), linewidth = 0.7, alpha = 0.9) +
  geom_line(aes(y = Q_m2, color = nms[3]), linewidth = 0.7, alpha = 0.9) +
  geom_line(aes(y = Q_m3, color = nms[4]), linewidth = 0.7, alpha = 0.9) +
  geom_line(aes(y = Q_m4, color = nms[5]), linewidth = 0.7, alpha = 0.9) +
  
  scale_color_manual(
    values = setNames(
      c(COR_OBS, COR_Q_M1, COR_Q_M2, COR_Q_M3, COR_Q_M4), 
      nms
    )
  ) +
  
  theme_minimal(base_size = 11) + 
  theme(legend.position = "bottom") +
  guides(color = guide_legend(nrow = 3)) + 
  labs(
    x = NULL, 
    y = "Vazão (mm/mês)", 
    color = NULL
  )

ggplotly(pm_i, tooltip = c("x", "y", "colour")) %>%
  layout(
    legend = list(orientation = "h", y = -0.25), 
    hovermode = "x unified"
  )

Comparação entre a série observada e os quatro modelos com lag ótimo.

Como eu havia falado, aqui os modelos 2 e 4 possuem uma sobreposição quase perfeita, indicando que a inclusão de ET não traz ganho adicional quando P já está presente ou que a inclusão de P não traz ganho adicional quando ET já está presente. Isso sugere que a informação de um dos dois preditores pode estar sendo capturada indiretamente. O que leva a necessidade de entender melhor a interação entre esses fatores.

Algumas possibilidades de análise podem incluir futuros modelos que explorem interações entre os preditores, como TWSA:P, TWSA:ET e P:ET. Essas interações podem revelar mecanismos hidrológicos mais complexos, como o excesso de saturação, que é particularmente relevante para o Pantanal.

Interpretação física de cada termo:

Termo Significado hidrológico
TWSA Estado de armazenamento da bacia (solo, aquíferos, planície). TWSA positiva → bacia saturada → mais escoamento.
P Forçamento pluviométrico recente. Chuva alta aumenta a descarga diretamente.
ET Perda evaporativa. ET alta reduz o escoamento ao drenar água do solo antes que chegue ao rio.
TWSA:P Interação armazenamento–chuva: a resposta da bacia à precipitação depende do estado de saturação. Com TWSA positiva (bacia cheia), uma chuva gera muito mais escoamento que com TWSA negativa (bacia seca). Este é o mecanismo de excesso de saturação
TWSA:ET A eficiência da ET depende da disponibilidade hídrica. Com TWSA negativa, as plantas sofrem estresse e a ET reduz; com TWSA positiva, a ET opera em plena capacidade.
P:ET Dias chuvosos reduzem a demanda atmosférica (menor radiação, maior umidade do ar), diminuindo a ET e aumentando o aproveitamento da chuva para escoamento.

8 Etapa 4: Análise das Vazões a Montante

Por que ainda limitado? A variável-resposta inclui a contribuição das sub-bacias do Planalto (Cáceres, Cuiabá, etc.), que têm dinâmica diferente do Pantanal.

Para avaliar as vazões a montante e seu impacto na série de Porto Murtinho, foi realizada uma análise de correlação cruzada entre as estações tributárias e a estação principal. O objetivo é identificar o lag de maior correlação, que indica o tempo de viagem da onda de cheia do tributário até Porto Murtinho. Isso se faz necessário para que possa ser separado a vazão vinda das estações a montante da vazão incremental gerada no proprio bioma do pantanal, que é a variavel que queremos entender pois os nossos dados correspondem a essa região.

8.1 Correlação cruzada — todas as estações

A CCF (Cross-Correlation Function) mede a correlação entre dois sinais em diferentes deslocamentos temporais. O lag de maior correlação nos valores negativos indica o tempo de viagem da onda de cheia do tributário até Porto Murtinho.

# ==============================================================================
# Leitura dos Dados das Estações Fluviométricas
# ==============================================================================

ps <- "Dados/Estações de Sedes Municipais"
pa <- "Dados/Estações Auxiliares"

pm_ <- read_excel(file.path(ps, "Estacao_67100000.xlsx")) %>% rename(Porto_Murtinho  = Vazao)
cac <- read_excel(file.path(ps, "Estacao_66070004.xlsx")) %>% rename(Caceres         = Vazao)
cui <- read_excel(file.path(ps, "Estacao_66260002.xlsx")) %>% rename(Cuiaba          = Vazao)
acg <- read_excel(file.path(pa, "Estacao_66460000.xlsx")) %>% rename(AC_Grande       = Vazao)
cox <- read_excel(file.path(ps, "Estacao_66870000.xlsx")) %>% rename(Coxim           = Vazao)
aqu <- read_excel(file.path(ps, "Estacao_66945000.xlsx")) %>% rename(Aquidauana      = Vazao)
mir <- read_excel(file.path(ps, "Estacao_66910000.xlsx")) %>% rename(Miranda         = Vazao)

# ==============================================================================
# União das Séries Temporais Diárias
# ==============================================================================

# Juntando todas as estações pela data e hora (DataHora)
df_vz <- pm_ %>% 
  full_join(cac, by = "DataHora") %>% 
  full_join(cui, by = "DataHora") %>%
  full_join(acg, by = "DataHora") %>% 
  full_join(cox, by = "DataHora") %>%
  full_join(aqu, by = "DataHora") %>% 
  full_join(mir, by = "DataHora") %>% 
  arrange(DataHora)

# ==============================================================================
# Agregação Mensal e Formatação da Tabela
# ==============================================================================

# Transforma os dados para o formato longo (pivot_longer), calcula a média mensal por estação e retorna para o formato largo (pivot_wider).
df_larg <- df_vz %>%
  pivot_longer(
    cols = -DataHora, 
    names_to = "est", 
    values_to = "Q"
  ) %>%
  mutate(ano_mes = floor_date(DataHora, "month")) %>%
  group_by(est, ano_mes) %>% 
  summarise(
    Q = mean(Q, na.rm = TRUE), 
    .groups = "drop"
  ) %>%
  pivot_wider(
    names_from = est, 
    values_from = Q
  ) %>% 
  arrange(ano_mes)
# ==============================================================================
# Definição dos Tributários e Preparação da CCF
# ==============================================================================

# Vetor com os nomes das colunas das estações a montante
ests_mont <- c(
  "Caceres", "Cuiaba", "AC_Grande", 
  "Coxim", "Aquidauana", "Miranda"
)

# Nomes formatados para exibição nos gráficos
nomes_disp <- c(
  Caceres    = "Cáceres", 
  Cuiaba     = "Cuiabá", 
  AC_Grande  = "Ac. C. Grande",
  Coxim      = "Coxim", 
  Aquidauana = "Aquidauana", 
  Miranda    = "Miranda"
)

# Inicialização das variáveis para armazenar resultados
lags_rios <- c()
ccf_list  <- list()

# ==============================================================================
# Cálculo da Correlação Cruzada (Iterando por Tributário)
# ==============================================================================

for (e in ests_mont) {
  
  # Verifica se a estação existe no dataframe
  if (!e %in% names(df_larg)) next
  
  # Identifica os meses (linhas) que possuem dados válidos tanto para o tributário quanto para Porto Murtinho
  v <- complete.cases(df_larg[[e]], df_larg$Porto_Murtinho)
  
  # Calcula a CCF (Cross-Correlation Function)
  rc <- ccf(
    x = df_larg[[e]][v], 
    y = df_larg$Porto_Murtinho[v], 
    lag.max = 12, 
    plot = FALSE, 
    na.action = na.pass
  )
  
  # Filtra apenas os lags <= 0 
  # (O tributário a montante influencia Porto Murtinho depois, não antes)
  lv <- rc$lag[rc$lag <= 0]
  av <- rc$acf[rc$lag <= 0]
  
  # Encontra o lag (em módulo) que possui a maior correlação positiva
  bl <- abs(lv[which.max(av)])
  
  # Armazena o lag ótimo encontrado no vetor nomeado
  lags_rios <- c(lags_rios, setNames(bl, e))
  
  ccf_list[[e]] <- data.frame(
    est = nomes_disp[e], 
    Lag = rc$lag, 
    Cor = rc$acf, 
    lo  = -bl
  )
}

# ==============================================================================
# Construção do Gráfico (ggplot2)
# ==============================================================================

df_ccf_all <- bind_rows(ccf_list)

lag_lbs <- data.frame(
  est = nomes_disp[names(lags_rios)],
  lo  = -as.integer(lags_rios),
  lb  = paste0("lag=", lags_rios, "m")
)

ggplot(df_ccf_all, aes(x = Lag, y = Cor)) +
  
  geom_col(fill = COR_TWSA, alpha = 0.65, width = 0.3) +
  geom_hline(yintercept = 0, color = "black", linewidth = 0.35) +
  geom_vline(
    data = lag_lbs, aes(xintercept = lo), 
    color = "red", linetype = "dashed", linewidth = 0.9
  ) +
  geom_text(
    data = lag_lbs, aes(x = lo - 0.3, y = 0.88, label = lb), 
    hjust = 1, size = 3, color = "red"
  ) +
  facet_wrap(~est, ncol = 3) + 
  theme_minimal(base_size = 10) +
  labs(
    x = "Defasagem (meses)", 
    y = "Correlação", 
    title = "CCF: tributários vs Porto Murtinho (dados reais)"
  )
Correlação cruzada entre cada tributário e Porto Murtinho usando os dados reais da pesquisa. A linha vermelha tracejada marca o lag ótimo selecionado automaticamente.

Correlação cruzada entre cada tributário e Porto Murtinho usando os dados reais da pesquisa. A linha vermelha tracejada marca o lag ótimo selecionado automaticamente.

Não fez muito sentido levando em consideração o diagrama unifilar:

Exemplos:

  • carceres tem um tempo de resposta para Porto Murtinho de 65 dias, mas o CCF estimou 4 meses;

  • Miranda esta relativamente proximo de Porto Murtinho mas também apresentou 4 meses.

# ==============================================================================
# Preparação
# ==============================================================================

# Copia do dataframe original para não alterá-lo diretamente
df_calc <- df_larg

# Vetor para armazenar os nomes das novas colunas defasadas
cols_lag_v <- c()

# Cria uma tabela auxiliar com o número exato de dias de cada mês
# (necessário para a conversão de m³/s para mm/mês)
dias_df <- df_vz %>% 
  pivot_longer(
    cols = -DataHora, 
    names_to = "e", 
    values_to = "v"
  ) %>%
  mutate(ano_mes = floor_date(DataHora, "month")) %>%
  group_by(ano_mes) %>% 
  summarise(
    dias = days_in_month(first(ano_mes)), 
    .groups = "drop"
  )

# ==============================================================================
# Aplicação dos Lags por Tributário
# ==============================================================================

# Iteramos sobre as estações a montante para defasar cada uma com seu respectivo lag
for (e in ests_mont) {
  
  if (!e %in% names(df_larg)) next
  
  # Define o nome da nova coluna defasada (ex: Caceres_lag)
  nm <- paste0(e, "_lag")
  
  # Aplica o lag específico daquele rio (lags_rios foi gerado na CCF)
  df_calc[[nm]] <- dplyr::lag(df_calc[[e]], lags_rios[e])
  
  # Guarda o nome da coluna para podermos somar todas elas facilmente no próximo passo
  cols_lag_v <- c(cols_lag_v, nm)
}

# ==============================================================================
# Balanço a Montante e Cálculo da Vazão Incremental
# ==============================================================================

df_bal <- df_calc %>% 
  left_join(dias_df, by = "ano_mes") %>%
  
  # Agrupamento e soma sum() das vazões a montante (defasadas) para cada mês. 
  rowwise() %>% 
  mutate(Soma_Mont = sum(c_across(all_of(cols_lag_v)), na.rm = TRUE)) %>% 
  ungroup() %>%
  
  # Conversão de unidades e cálculo final da incremental
  mutate(
    # Fator de conversão de m³/s para mm/mês
    fator   = (86.4 * dias) / area_bacia,
    
    # Converte as vazões para mm/mês
    PM_mm   = Porto_Murtinho * fator,
    Mont_mm = Soma_Mont * fator,
    
    # Calcula a vazão incremental gerada dentro do Pantanal.
    # O pmax(..., 0) garante que meses onde o Pantanal apenas reteve água (incremental "negativa") sejam tratados fisicamente como zero escoamento direto.
    Q_inc   = pmax(PM_mm - Mont_mm, 0)
  )
# ==============================================================================
# Gráfico da Decomposição da Vazão (Total, Montante e Incremental)
# ==============================================================================

df_plot_inc <- df_bal %>% 
  select(ano_mes, PM_mm, Mont_mm, Q_inc) %>%
  pivot_longer(
    cols = -ano_mes, 
    names_to = "var", 
    values_to = "Q"
  ) %>%
  mutate(
    var = recode(
      var, 
      PM_mm   = "Vazão Total (PM)",
      Mont_mm = "Contribuição Montante",
      Q_inc   = "Vazão Incremental (Pantanal)"
    )
  )

p_inc <- ggplot(df_plot_inc, aes(x = ano_mes, y = Q, color = var)) +
  
  geom_line(linewidth = 0.85, alpha = 0.9) +
  geom_hline(yintercept = 0, linetype = "dashed", color = "gray50") +
  scale_color_manual(
    values = c(
      "Vazão Total (PM)"             = COR_OBS,
      "Contribuição Montante"        = COR_TWSA,
      "Vazão Incremental (Pantanal)" = COR_Q_M3
    )
  ) +
  theme_minimal(base_size = 11) +
  theme(legend.position = "bottom") +
  labs(
    x = NULL, 
    y = "Vazão (mm/mês)", 
    color = NULL
  )
ggplotly(p_inc, tooltip = c("x", "y", "colour")) %>% 
  layout(
    legend = list(orientation = "h", y = -0.2), 
    hovermode = "x unified"
  )

Decomposição da vazão total: contribuição das sub-bacias do Planalto (azul) e vazão incremental gerada dentro do Pantanal (vermelho).

As vazões a montante de contribuição dos tributários do Planalto representa boa parte do que é observado em Porto Murtinho e possuem um tempo de viagem de 4-5 meses até a exutoria da bacia.

Em alguns trechos da série é possivel ver que a vazão gerada no pantanal tende a ser negativa, isso implica que a planicie reteve à agua que gerou, podendo inclusive reter parte da vazão de contribuição a montante.

Isso dificulta a ánalise pois os modelos são log lineares e não conseguem lidar com valores negativos, o que levou a necessidade de filtrar apenas os meses com vazão incremental estritamente maior que zero para a etapa seguinte de modelagem.

8.2 Aplicação dos modelos log-lineares à série incremental

dados_inc <- df_bal %>% 
  select(ano_mes, PM_mm, Mont_mm, Q_inc) %>%
  inner_join(grace_m, by = "ano_mes") %>%
  inner_join(era5_m, by = "ano_mes") %>%
  inner_join(ssebop_m, by = "ano_mes") %>%
  drop_na() %>% 
  arrange(ano_mes)

df_f_inc <- dados_inc %>% filter(Q_inc > 0)

# ==============================================================================
# Ajuste dos Modelos (Treinados na Incremental)
# ==============================================================================

mod1i <- lm(log(Q_inc) ~ TWSA, data = df_f_inc)
mod2i <- lm(log(Q_inc) ~ TWSA + P_ERA5, data = df_f_inc)
mod3i <- lm(log(Q_inc) ~ TWSA + ET_SSEBop, data = df_f_inc)
mod4i <- lm(log(Q_inc) ~ TWSA + P_ERA5 + ET_SSEBop, data = df_f_inc)

# Previsão bruta da incremental (escala original em mm/mês)
dados_inc$Q_m1i_raw <- exp(predict(mod1i, newdata = dados_inc))
dados_inc$Q_m2i_raw <- exp(predict(mod2i, newdata = dados_inc))
dados_inc$Q_m3i_raw <- exp(predict(mod3i, newdata = dados_inc))
dados_inc$Q_m4i_raw <- exp(predict(mod4i, newdata = dados_inc))

# ==============================================================================
# Otimização do Lag e Reconstrução da Vazão Total
# ==============================================================================

# Encontra o lag ótimo comparando a incremental simulada com a incremental observada
res1i <- get_best_lag(dados_inc$Q_inc, dados_inc$Q_m1i_raw)
res2i <- get_best_lag(dados_inc$Q_inc, dados_inc$Q_m2i_raw)
res3i <- get_best_lag(dados_inc$Q_inc, dados_inc$Q_m3i_raw)
res4i <- get_best_lag(dados_inc$Q_inc, dados_inc$Q_m4i_raw)

# Aplica o lag
dados_inc$Q_m1i <- dplyr::lag(dados_inc$Q_m1i_raw, res1i$lag)
dados_inc$Q_m2i <- dplyr::lag(dados_inc$Q_m2i_raw, res2i$lag)
dados_inc$Q_m3i <- dplyr::lag(dados_inc$Q_m3i_raw, res3i$lag)
dados_inc$Q_m4i <- dplyr::lag(dados_inc$Q_m4i_raw, res4i$lag)

# Reconstrução: Vazão Total = Vazão Montante (Real) + Vazão Incremental (Simulada)
dados_inc$PM_m1i <- dados_inc$Mont_mm + dados_inc$Q_m1i
dados_inc$PM_m2i <- dados_inc$Mont_mm + dados_inc$Q_m2i
dados_inc$PM_m3i <- dados_inc$Mont_mm + dados_inc$Q_m3i
dados_inc$PM_m4i <- dados_inc$Mont_mm + dados_inc$Q_m4i

# ==============================================================================
# Cálculo de Métricas (Em Relação à Vazão Total PM_mm)
# ==============================================================================

# Função auxiliar para recalcular o R² baseado na vazão TOTAL reconstruída
calc_r2 <- function(obs, sim) {
  df <- data.frame(o = obs, s = sim) %>% drop_na()
  round(1 - sum((df$s - df$o)^2) / sum((df$o - mean(df$o))^2), 3)
}

# R² da Reconstrução
r2_1i <- calc_r2(dados_inc$PM_mm, dados_inc$PM_m1i)
r2_2i <- calc_r2(dados_inc$PM_mm, dados_inc$PM_m2i)
r2_3i <- calc_r2(dados_inc$PM_mm, dados_inc$PM_m3i)
r2_4i <- calc_r2(dados_inc$PM_mm, dados_inc$PM_m4i)

# RMSE da Reconstrução
rmse1i <- sqrt(mean((dados_inc$PM_mm - dados_inc$PM_m1i)^2, na.rm = TRUE))
rmse2i <- sqrt(mean((dados_inc$PM_mm - dados_inc$PM_m2i)^2, na.rm = TRUE))
rmse3i <- sqrt(mean((dados_inc$PM_mm - dados_inc$PM_m3i)^2, na.rm = TRUE))
rmse4i <- sqrt(mean((dados_inc$PM_mm - dados_inc$PM_m4i)^2, na.rm = TRUE))
Desempenho dos quatro modelos estimando a vazão incremental, reconstruídos para a Vazão Total de Porto Murtinho.
Modelo R2 RMSE Lag
M1i-TWSA 0.597 2.335 0
M2i-TWSA+P 0.601 2.323 0
M3i-TWSA+ET 0.633 2.226 0
M4i-TWSA+P+ET 0.646 2.189 0

Isso é um pouco confuso, mas eu entendi que o lag ótimo calculado foi a 0 meses pois foi feito com os dados da vazão incremental, talvez tenha perdido a sensibilidade devido acumulo de erros na operação ou talvez os dados do balanço hidrico (P, TWSA, ET) realmente respondem mais rapido quando se trata do que é gerado somente na planicie.

Eu chuto que tem acumulo de erro devido a elasticidade do atraso dos fatores, ao invês de um lag ótimo para toda a série ele seria uma variavel ao longo do tempo que depende do estado de armazenamento do pantanal, porém não sei como aplicar essa lógica.

Um problema que surge disso também é a falta de estoque de dados para predições ( necessáriamente os dados do mês atual melhor se relacionam com a previsão da vazão do mesmo mês)

# ==============================================================================
# Gráfico de Comparação: Vazão Total Reconstruída
# ==============================================================================

nms_i <- c(
  "Vazão Total Observada",
  paste0("M1i-TWSA (R²=", r2_1i, ")"),
  paste0("M2i-TWSA+P (R²=", r2_2i, ")"),
  paste0("M3i-TWSA+ET (R²=", r2_3i, ")"),
  paste0("M4i-TWSA+P+ET (R²=", r2_4i, ")")
)

pm_i <- ggplot(dados_inc, aes(x = ano_mes)) +
  
  geom_line(aes(y = PM_mm, color = nms_i[1]),  linewidth = 1.1) +
  geom_line(aes(y = PM_m1i, color = nms_i[2]), linewidth = 0.7, alpha = 0.9) +
  geom_line(aes(y = PM_m2i, color = nms_i[3]), linewidth = 0.7, alpha = 0.9) +
  geom_line(aes(y = PM_m3i, color = nms_i[4]), linewidth = 0.7, alpha = 0.9) +
  geom_line(aes(y = PM_m4i, color = nms_i[5]), linewidth = 0.7, alpha = 0.9) +
  
  scale_color_manual(
    values = setNames(
      c(COR_OBS, COR_Q_M1, COR_Q_M2, COR_Q_M3, COR_Q_M4), 
      nms_i
    )
  ) +
  theme_minimal(base_size = 11) + 
  theme(legend.position = "bottom") +
  guides(color = guide_legend(nrow = 3)) + 
  labs(
    x = NULL, 
    y = "Vazão Total (mm/mês)", 
    color = NULL
  )

ggplotly(pm_i, tooltip = c("x", "y", "colour")) %>%
  layout(
    legend = list(orientation = "h", y = -0.25), 
    hovermode = "x unified"
  )

Comparação entre a Vazão Total Observada e a reconstrução feita pelos quatro modelos modelados via vazão incremental.

Os modelos com a vazão incremental apresentaram um desempenho melhor, mas ainda estão longe de capturar a dinâmica total da vazão em Porto Murtinho. Porem foi o melhor resultado obtido até o presente momento.

8.3 Teste do modelo em prever a vazão observada a partir da vazão incremental simulada

Para realizar a predição o modelo será calibrado somente com dados antes do periodo de testes (jan/2022) e depois será avaliado o desempenho do modelo em prever a vazão observada a partir da vazão incremental simulada para os meses de teste (jan/2022 em diante). O modelo de teste será o modelo 4i, que é o modelo mais completo e apresentou o melhor desempenho na etapa anterior.

O corte foi selecionado em 2022 devido a seca extrema que ocorreu no pantanal, se o modelo não tivesse pelo menos alguns dados desse periodo de minima hda serie historica disponivel para modelagem não teria como aprender a dinâmica de seca extrema. O que me leva a pensar no uso de séries sinteticas do GRACE, implementação de series historicas de pluviometros no lugar de satélites e o calculo de ET a partir da cobertura vegetal obtida por imagens de satélite. (sim, muda completamente a proposta da pesquisa, que é trabalhar com sensoriamento remoto rsrsrsrs.)

Os lags dos satélites são calculados apenas com os dados de treino, evitando contaminação do período de teste.

9 Execução

# ==============================================================================
# Preparação da Base Única
# ==============================================================================

dados_br <- df_bal %>% 
  select(ano_mes, PM_mm, Mont_mm, Q_inc) %>%
  inner_join(grace_m  %>% select(ano_mes, TWSA),      by = "ano_mes") %>%
  inner_join(ssebop_m %>% select(ano_mes, ET_SSEBop), by = "ano_mes") %>%
  inner_join(era5_m   %>% select(ano_mes, P_ERA5),    by = "ano_mes") %>%
  arrange(ano_mes)

# Filtro de treino para cálculo de Lag
tr_lag <- dados_br %>% 
  filter(ano_mes < data_corte) %>% 
  drop_na()

# ==============================================================================
# Otimização dos Lags (CCF)
# ==============================================================================

encontrar_lag_ccf <- function(preditor, resposta) {
  res <- ccf(preditor, resposta, lag.max = 18, plot = FALSE, na.action = na.pass)
  lv <- res$lag[res$lag <= 0]
  av <- res$acf[res$lag <= 0]
  return(abs(lv[which.max(abs(av))]))
}

lag_P    <- encontrar_lag_ccf(tr_lag$P_ERA5,    tr_lag$Q_inc)
lag_ET   <- encontrar_lag_ccf(tr_lag$ET_SSEBop, tr_lag$Q_inc)
lag_TWSA <- encontrar_lag_ccf(tr_lag$TWSA,      tr_lag$Q_inc)

# ==============================================================================
# Aplicação dos Lags e Divisão de Dados
# ==============================================================================

dados <- dados_br %>%
  mutate(
    P_lag    = dplyr::lag(P_ERA5, lag_P),
    ET_lag   = dplyr::lag(ET_SSEBop, lag_ET),
    TWSA_lag = dplyr::lag(TWSA, lag_TWSA)
  ) %>% 
  drop_na()

# Divisão Treino/Teste
dados_treino <- dados %>% filter(ano_mes < data_corte)
dados_teste  <- dados %>% filter(ano_mes >= data_corte)

# Para o Modelo 4 (log-linear), precisamos filtrar Q_inc > 0 no treino
dados_treino_fit <- dados_treino %>% filter(Q_inc > 0)

# ==============================================================================
# Modelagem: Modelo 4 (Log-Linear na Incremental)
# ==============================================================================

# Ajuste do Modelo 4 usando as variáveis defasadas
mod4i_final <- lm(log(Q_inc) ~ TWSA_lag + P_lag + ET_lag, data = dados_treino_fit)

# ==============================================================================
# Previsão e Reconstrução da Vazão Total (PM)
# ==============================================================================

# Prever a vazão incremental (exp para voltar da escala log)
dados$Q_inc_prev <- exp(predict(mod4i_final, newdata = dados))

# Reconstruir a Vazão Total: Montante Real + Incremental Prevista
dados$PM_prev <- dados$Mont_mm + dados$Q_inc_prev

# ==============================================================================
# Avaliação de Métricas (Comparando com a Vazão Total Observada)
# ==============================================================================

calcular_metricas <- function(obs, sim, label_periodo) {
  df <- data.frame(o = obs, s = sim) %>% drop_na()
  SST <- sum((df$o - mean(df$o))^2)
  SSR <- sum((df$s - df$o)^2)
  
  data.frame(
    Periodo = label_periodo,
    R2      = round(1 - (SSR / SST), 3),
    RMSE    = round(sqrt(mean((df$s - df$o)^2)), 2)
  )
}

met_tr_m4 <- calcular_metricas(
  obs = dados$PM_mm[dados$ano_mes < data_corte], 
  sim = dados$PM_prev[dados$ano_mes < data_corte], 
  "Treino (M4i)"
)

met_te_m4 <- calcular_metricas(
  obs = dados$PM_mm[dados$ano_mes >= data_corte], 
  sim = dados$PM_prev[dados$ano_mes >= data_corte], 
  "Teste (M4i)"
)

met_gl_m4 <- calcular_metricas(
  obs = dados$PM_mm, 
  sim = dados$PM_prev, 
  "Global (M4i)"
)

# ==============================================================================
# Tabela Consolidada de Resultados (Treino, Teste e Global)
# ==============================================================================

# Unimos as métricas geradas no bloco de modelagem e formatamos a tabela
bind_rows(met_tr_m4, met_te_m4, met_gl_m4) %>%
  kable(
    caption = "Métricas finais de desempenho do Modelo 4 (Vazão Incremental Reconstruída).",
    col.names = c("Período", "R²", "RMSE (mm/mês)"),
    digits = 3
  ) %>%
  kable_styling(
    bootstrap_options = c("striped", "hover", "condensed")
  ) %>%
  # Destaca a linha Global (linha 3)
  row_spec(3, bold = TRUE, background = "#e8f4fd")
Métricas finais de desempenho do Modelo 4 (Vazão Incremental Reconstruída).
Período R² RMSE (mm/mês)
Treino (M4i) 0.604 2.29
Teste (M4i) 0.465 1.98
Global (M4i) 0.629 2.24 
# ==============================================================================
# Gráfico de Série Temporal: Observado vs. Treino vs. Teste
# ==============================================================================

# Preparação: Separando a previsão em duas colunas para plotar cores diferentes
dados_plot <- dados %>%
  mutate(
    # A linha de treino inclui o mês de corte apenas para conectar visualmente
    PM_prev_tr = ifelse(ano_mes <= data_corte, PM_prev, NA),
    # A linha de teste começa a partir da data de corte
    PM_prev_te = ifelse(ano_mes >= data_corte, PM_prev, NA)
  )

# Construção do Gráfico (ggplot2)
p_sf <- ggplot(dados_plot, aes(x = ano_mes)) +
  
  # Fundo sombreado para a área de teste
  annotate("rect", xmin = data_corte, xmax = max(dados_plot$ano_mes),
           ymin = -Inf, ymax = Inf, alpha = 0.07, fill = "gray40") +
  
  geom_vline(xintercept = as.Date(data_corte), linetype = "dashed", color = "gray50") +
  geom_line(aes(y = PM_mm, color = "Observado"), linewidth = 1) +
  geom_line(aes(y = PM_prev_tr, color = "Previsto (Treino)"), linewidth = 0.85, alpha = 0.9) +
  geom_line(aes(y = PM_prev_te, color = "Previsto (Teste)"), linewidth = 0.9) +
  annotate("text", x = as.Date("2010-01-01"), y = max(dados_plot$PM_mm, na.rm = TRUE) * 0.93,
           label = "TREINO", fontface = "bold", color = "gray40", size = 4) +
  annotate("text", x = data_corte + 365, y = max(dados_plot$PM_mm, na.rm = TRUE) * 0.93,
           label = "TESTE", fontface = "bold", color = "black", size = 4) +
  
  # Configuração de Cores
  scale_color_manual(
    values = c(
      "Observado"         = COR_OBS,
      "Previsto (Treino)" = "#377eb8", # Azul
      "Previsto (Teste)"  = COR_Q_M3   # Vermelho padrão dos seus gráficos
    )
  ) +
  
  theme_minimal(base_size = 11) +
  theme(legend.position = "bottom") +
  labs(x = NULL, y = "Vazão Total (mm/mês)", color = NULL)

# 3. Conversão para Plotly
ggplotly(p_sf, tooltip = c("x", "y", "colour")) %>%
  layout(
    hovermode = "x unified", 
    legend = list(orientation = "h", y = -0.15)
  )

Gráfico de previsão de vazão

10 Resultados do Modelo Final

# ==============================================================================
# Gráfico de Dispersão e Aderência do Modelo 4
# ==============================================================================

# Preparação dos dados para a legenda dinâmica
dd <- dados %>% 
  mutate(
    Fase = ifelse(ano_mes < data_corte, "Treino", "Teste"),
    ds   = format(ano_mes, "%b/%Y") # Formatação de data (ex: Jan/2015) para o tooltip
  )

p <- ggplot(dd, aes(
      x = PM_mm, 
      y = PM_prev, 
      color = Fase,
      text = paste0(ds, "<br>Obs: ", round(PM_mm, 2), " mm<br>Prev: ", round(PM_prev, 2))
    )) +
  
  geom_abline(slope = 1, intercept = 0, linetype = "dashed", color = "black", linewidth = 0.9) +
  
  geom_point(size = 2.5, alpha = 0.75) +
  geom_smooth(method = "lm", se = FALSE, linewidth = 0.8) +
  scale_color_manual(
    values = c("Treino" = "#377eb8", "Teste" = COR_Q_M3)
  ) +
  
  theme_minimal(base_size = 11) +
  theme(legend.position = "bottom") +
  labs(
    x = "Observado (mm/mês)",
    y = "Previsto (mm/mês)",
    color = NULL,
    title = paste0("R² Global: ", met_gl_m4$R2, " | R² Teste: ", met_te_m4$R2)
  )

ggplotly(p, tooltip = "text") %>%
  layout(legend = list(orientation = "h", y = -0.15))

Gráfico de Dispersão Observado vs. Previsto para o Modelo 4i. A linha tracejada representa a linha 1:1 (previsão perfeita). Os pontos azuis correspondem ao período de treino e os vermelhos ao teste. Passe o mouse para ver os valores individuais.

11 Discussão e Próximos Passos

O modelo atual melhorou consideravelmente ao utilizar a vazão incremental, porém perdeu o tempo de transporte dos dados observados no pantanal até porto murtinho (lag ótimo = 0). A pergunta da pesquisa até agora segue com uma parcial de que: É possivel modelar a vazão observada em Porto Murtinho utilizando dados de sensoriamento remoto, entretando ao incrementar dados de vazão observado em estações fluviometricas a montante da planicie os dados melhoram consideravelmente, o que pode dificultar o uso restrito de dados de sensoriamento remoto para tal.

  • Porque no Balanço Tradicional o lag otimo é igual a 0 (dS = P - ET - Q)
  • Porque não há ganho significativo no modelo 4 em relação ao modelo 2 log-linear?
  • Como lidar com o lag igual a 0 para construição de modelos de predição?

Referências

Duvvuri, Bhavya, e Edward Beighley. 2023. “Estimating Monthly River Discharges from GRACE/GRACE-FO Terrestrial Water Storage Anomalies”. Remote Sensing 15: 4516. https://doi.org/10.3390/rs15184516.
Hamilton, Stuart K., Steven J. Sippel, e John M. Melack. 1996. “Inundation patterns in the Pantanal wetland of South America determined from passive microwave remote sensing”. Archiv für Hydrobiologie 137: 1–23.
Kirchner, James W. 2009. “Catchments as simple dynamical systems”. Water Resources Research 45. https://doi.org/10.1029/2008WR006912.
Moreira, Adriana Aparecida, Alice Fassoni-Andrade, Anderson Ruhoff, e Rodrigo Dias de Paiva. 2019. “Water Balance Based on Remote Sensing Data in Pantanal”. Revista Ra’e Ga 46: 20–32. https://doi.org/10.5380/raega.