O TOPMODEL e o TWI

Modelos hidrológicos objetivam melhorar a tomada de decisão em relação a problemas de recursos hídricos. Como qualquer modelo, são representações imperfeitas da realidade, ainda que úteis.

Um modelo de especial relevância na hidrologia é o TOPMODEL, apresentado por Beven e Kirkby em um artigo de 1979. Nessa época, os Sistemas de Informações Geográficas já exisitiam em computadores centralizados, mas ainda estavam para ser distribuídos em PCs. Inclusive, segundo Beven em um entrevista, o manuscrito foi rejeitado inicialmente por uma revista em razão de sua abordagem “geográfica” demais.

Entre as particularidades do TOPMODEL é que ele permite o mapeamento dos processos hidrológicos na bacia hidrográfica mesmo sem ser um modelo verdadeiramente distribuído. Na terminologia vigente, é um modelo semi-distribuído.

O mapeamento é possível por que existe um índice topográfico que reflete a tendência da água subterrânea aflorar em um ponto qualquer do terreno. Esse índice chama-se TWI, de “topographical wetness index” e pode ser derivado de um modelo digital de elevação no QGIS.

A área de contribuição variável (VSA)

Imagine que uma bacia hidrográfica enfrenta um longo período de estiagem. Simplesmente não chove. Por meses. Durante esse período, toda a água presente no solo ou é drenada para a zona saturada da água subterrânea, ou é sequestrada pelas raízes das plantas, que já estão sofrendo com a falta de água. O solo está seco. A água subterrânea, por sua vez, está sendo completamente drenada para o fundos dos vales, contribuindo assim para um ínfimo escoamento de base dos riachos que afluem para o rio principal. Nesse ponto, a zona saturada está a ponto de secar também. O rio vai secar se não chover.

Como que por um milagre, as chuvas voltam. As águas passam a se infiltrar no solo e a contribuir para a recarga da zona saturada. Com o aumento de carga, a água começa a aflorar. No início, esse afloramento ocorre em pontos mais baixos do terreno, em fundos de vale. Áreas de convergência topográfica, antes secas, se transformam em banhados e pântanos.

Mas as chuvas não param. A zona saturada no solo começa então a crescer em volume, aflorando a água paulatinamente em diversas manchas cada vez maiores e mais distantes dos fundos dos vales até que sobem o terreno em direção aos talvegues das encostas íngremes.

Nesse ponto, se as chuvas não pararem, as coisas podem ficar perigosas. A chuva que antige as áreas de solo saturado simplesmente não se infiltra, pois não há mais poros vazios para absorvê-la. Assim, pulsos de escoamento superficial passam a ocorrer com mais intensidade e frequência. Isso amplia os picos de vazão e as taxas de erosão. Com mais chuvas ainda, as encostas mais íngremes não suportam mais seu próprio peso de tanta água e colapsam em violentos movimentos de massa.

Esse cenário imaginário vai de um extremo até outro em termos de saturação da bacia hidrográfica. Para lugares com climas úmidos e solos relativamente rasos esse retrato pode ser bastante fidedigno, como atesta o histórico dos movimentos de massa na Serra do Mar.

O conceito chave por trás desse modelo perceptual é a área de contribuição variável (VSA – variable source area): a mancha de solo saturada pelo afloramento da água subterrânea varia de forma previsível em função da saturação da bacia, iniciando nos fundos dos vales e indo em direção aos talvegues das encostas. A chuva incidente sobre essa área não tem como se infiltrar, sendo uma área ativa de produção de escoamento superficial pelo excesso de saturação (também chamado processo Dunneano).

O TWI e a VSA

O papel do TWI no TOPMODEL é determinar a dinâmica da zona saturada, representada pelo déficit local de água. O mapa da VSA, por sua vez, é onde o déficit local for nulo ou negativo.

Com os conceitos do TOPMODEL, a VSA para derterminada condição de saturação da bacia pode ser mapeada em uma camada raster da seguinte forma:

Primeiramente, é preciso determinar o déficit local de água (isto é, nos píxeis da camada raster):

Di = D + m * (TWIm – TWIi)

Onde D é o deficit global de água na bacia, TWIm é o TWI médio da bacia, TWIi é o TWI local e m é um parâmetro que escalona a relação (o m também tem uma interpretação hidrológica, refletindo a transmissividade do solo da bacia).

O mapa da VSA será onde o solo está saturado, ou seja, o déficit local é nulo ou negativo:

VSA = Di >= 0

Se você implementar um modelo hidrológico para computar o déficit global da bacia, seja esse modelo o TOPMODEL ou não, você poderá visualizar a distribuição espacial do déficit local na bacia hidrográfica. Com o TOPMODEL, foi exatamente isso que eu fiz na animação a seguir. Se você prestar atenção, vais perceber a dinâmica da zona saturada à medida que a bacia é submetida a perídos de estiagem e vice-versa.

Mas então como então calculamos o TWI?

Derivando o TWI de um MDE

Os pacotes do SAGA apresentam uma ferramenta específica para determinar o TWI de um modelo digital de elevação (MDE). Se você está com pressa, corra lá e use essa ferramenta. Boa sorte.

Agora se você não confia muito na ferramenta do SAGA (principalmente em razão de existir pouca documentação sobre ela), vou apresentar abaixo os passos para produzir o TWI de forma artesanal.

Para qualquer ponto i do espaço, a fórmula do TWI é a seguinte:

TWI = LN( a / Tan(B))

Onde a é a área de drenagem acumulada por unidade de contorno no ponto i e B é a declividade do terreno no ponto i (em radianos).

Como sempre, o fluxograma:

Passo -1: obter um MDE

Se você caiu de pára-quedas aqui (seja bem-vind@!), tudo que vamos fazer é extraído de um Modelo Digital de Elevação (MDE). Nesse post aqui eu dou dicas de como obter um MDE do zero, diretamente dos dados do USGS.

Passo 0: obter um MDE sem depressões

De novo, se você apareceu aqui do nada (seja bem-vind@!), saiba que precisamos derivar um MDE sem buracos no caminho da água. Nesse post aqui eu dou dicas de como fazer isso usando o QGIS.

Passo 1: obter a Área de Drenagem

Essa etapa é muito conhecida como produzir o acúmulo de fluxo (flow accumulation). Nos pacotes do SAGA essa ferramenta é chamada de Catchment Area (disponível em Terrain Analysis – Hydrology), que faz muito mais sentido na minha opinião.

A ferramenta requer apenas que você forneça o MDE (sem depressões) e espere o processamento acabar. Dependendo do tamanho do MDE, esperar bastante.

Passo 2: obter a Declividade

Esse passo é para ser bem mas rápido. Uma opção é ir em GDAL/OGR – [GDAL] Analysis e selecionar a ferramenta Slope. Só precisa indicar o MDE.  

Passo 3: obter o TWI pela calculadora raster, não esquecendo dois detalhes

Vou repetir a fórmula do TWI:

TWI = LN( a / Tan(B))

Onde a é a área de drenagem acumulada por unidade de contorno no ponto i e B é a declividade do terreno no ponto i (em radianos).

Ou seja: o valor de a não é a Área de Drenagem que derivamos! Para o obter o valor de a é preciso dividir a Área de Drenagem pelo comprimento do píxel, isto é, pela resolução da camada raster. Por exemplo, se estiver derivando os dados com o MDE SRTM de 30 metros, dividir por 30.

Outra coisa: o valor de B é em radianos! A declividade produzida pelo Slope foi computada em graus (confira você mesmo).

Portanto, a fómula na calculadora raster para derivar o TWI deve incluir esses dois detalhes:

ln ( "catchment_area@1"/ (res * tan ( 2 * 3.141592* "slope@1" / 360 ))

Onde “catchment_area@1” é a camada raster da área de drenagem, res é o comprimento da resolução da camada raster e  “slope@1” é a camada raster de declividade em graus. Os argumentos dentro da função tangente são os conversores de graus para radianos.

Referência original do TOPMODEL

Beven and Kirkby. (1979). A physically based, variable contributing area model of basin hydrology. In Hydrological Sciences.