Modelagem numérica para o processo da evaporação da água do solo

Numeric modeling to simulate the soil water evaporation process

SEÇÃO I - FÍSICA DO SOLO

Modelagem numérica para o processo da evaporação da água do solo

Numeric modeling to simulate the soil water evaporation process

Celso Luiz Prevedello^I; Jocely Maria Thomazoni Loyola^II

^IProfessor do Departamento de Solos, Universidade Federal do Paraná - UFPR. R. dos Funcionários 1540, Cabral, CEP 80.035-050 Curitiba (PR). E-mail: [e-mail]

^IIProfessora do Departamento de Construção Civil, UFPR. E-mail: [e-mail]

RESUMO

Foram realizados vários ensaios laboratoriais para avaliar o desempenho de um modelo numérico em simular o processo unidimensional da evaporação da água do solo. Este modelo foi desenvolvido a partir da linearização da equação de Richards no espaço e do uso da técnica iterativa de Newton-Raphson, para resolver essa equação não-linear no tempo, e leva em conta o fluxo de vapor dentro do perfil de solo e a taxa de evaporação na sua superfície. Os resultados mostraram que o modelo foi capaz de simular satisfatoriamente o processo, tanto em meios porosos homogêneos quanto estratificados.

Termos de indexação: meios porosos homogêneos, meios porosos estratificados, equação de Richards, método de Newton-Raphson.

SUMMARY

Several laboratory tests were carried out with the objective of evaluating the performance of a numerical model to simulate the unidimensional soil water evaporation process. This model was created based on Richard's equation linearized in space and by an iterative Newton-Raphson method to solve this unlinearized equation in time. The model considers the vapor flow within the soil profile and its evaporation rate on the surface. Results showed that the model simulated the soil water evaporation process satisfactorily for both homogeneous and stratified porous media.

Index terms: homogeneous porous media, layered porous media, Richards equation, Newton-Raphson's.

INTRODUÇÃO

Todos os processos dinâmicos da água no solo, como infiltração, redistribuição e evaporação, são tipicamente transientes, e a solução analítica da equação diferencial parcial que descreve esses processos (equação de Richards), quando possível, resulta em grande grau de dificuldade. Isso, do ponto de vista matemático, decorre do fato de a equação de Richards não ser linear, uma vez que a condutividade hidráulica K na equação é função da variável dependente umidade q ou da pressão da água no solo p/rg que, por sua vez, é função da posição no espaço e do tempo. E a relação entre a condutividade hidráulica e a umidade (ou p/rg) é tal que a equação de Richards resulta altamente não-linear. Uma equação assim, via de regra, só admite solução analítica quando linearizada. Nem sempre, porém, uma equação diferencial pode ser linearizada, a menos que erros até significativos sejam introduzidos. Atualmente, com o advento dos microcomputadores, as soluções desses problemas podem ser encontradas por meio de métodos numéricos, com grandes vantagens.

Este trabalho teve por objetivo avaliar experimentalmente o desempenho de um modelo numérico em simular o processo unidimensional da evaporação da água do solo, desenvolvido a partir da linearização da equação de Richards no espaço e do uso da técnica iterativa de Newton-Raphson para resolver essa equação não-linear no tempo.

CONSIDERAÇÕES TEÓRICAS

Com base na teoria desenvolvida por Mualem (1976), nas equações combinadas de Van Genuchten (1980) para as funções q(p/rg), K(p/rg) (sendo K a condutividade hidráulica, q a umidade, e p/rg a pressão da água no solo), e no argumento de que a forma integrada da lei de Buckingham-Darcy pode produzir uma expressão analítica para a densidade de fluxo mátrico em regime de fluxo não-transiente, então é possível linearizar a equação de Richards no espaço. Com isso, o que restará será uma equação não-linear no tempo, cuja solução poderá ser obtida mediante uma técnica iterativa, como a de Newton-Raphson, por exemplo. As considerações que se seguem visam ao desenvolvimento analítico deste postulado.

Propriedades hidráulicas do solo e suas funções combinadas

Desde que muitas das propriedades hidráulicas do solo são difíceis de ser obtidas em condições de campo, seja pela variabilidade a que estão sujeitas, seja pelo custo e, ou, demanda de tempo dessas medidas, muitos autores têm procurado estabelecer, modificar ou substabelecer modelos para predizer a condutividade hidráulica como uma função da umidade a partir de outras propriedades hidráulicas mais facilmente mensuráveis, como, por exemplo, da curva de retenção de água do solo e da condutividade hidráulica na saturação (Mualem, 1976; Van Genuchten, 1980).

A teoria desenvolvida por Mualem (1976) mostra-se particularmente interessante porque permite a derivação analítica de uma série de importantes funções hidráulicas do solo, desde que a dependência p/pg versus q (curva de retenção de água no solo) também seja conhecida analiticamente. Um modelo empírico bastante versátil e que tem sido bastante utilizado (Marion et al., 1994; Prevedello et al., 1995; Simunek et al., 1998) para representar analiticamente a curva de retenção de água no solo foi proposto por Van Genuchten (1980), ou seja,

sendo Q uma função que pode ser expressa tanto em termos da umidade quanto da pressão da água no solo, ou seja, em que a, n, q_s (umidade na saturação), q_r (umidade residual) e m (m = 1-1/n) são parâmetros independentes que devem ser estimados pelo ajuste da equação (1) aos dados observados da curva de retenção de água no solo. Esses parâmetros podem ser simultaneamente estimados pela técnica dos mínimos quadrados para ajuste de funções não-lineares (Boratto, 1984), dentre outras.

A equação (1), quando introduzida na teoria de Mualem (1976), resulta que a condutividade hidráulica pode ser expressa como uma função de q ou p/rg, por meio da seguinte expressão:

em que K_s é a condutividade hidráulica do solo saturado e Q e m os parâmetros já definidos na equação (1).

Solução numérica da equação de Richards

Para expressar os processos unidimensionais da infiltração, redistribuição e evaporação da água do solo, Prevedello (1996) obteve a seguinte expressão:

em que Dt é o intervalo de tempo em que ocorre a variação q_i^j+1-q_i^j , sendo q_i^j+1 e q_i^j, as umidades na profundidade i nos tempos j+1 e j, respectivamente. O primeiro membro da equação 3 assume que a umidade na camada z_i-z_i-1 é representada pela umidade da profundidade i. Esse membro representa a variação de armazenamento de água ocorrida na camada z_i-z_i-1 durante o intervalo de tempo Dt. Num perfil de solo com i profundidades, é necessário encontrar os valores de p/rg que satisfaçam a igualdade da (3) para assegurar o princípio da conservação da matéria em todas essas profundidades, desde i = 1 até N, em que N é o número total de profundidades consideradas. E esses valores de p/rg que asseguram a igualdade da (3) podem ser conhecidos, já que as funções K(p/rg) e q(p/rg) são conhecidas (equações 2 e 1, respectivamente). Dessa forma, a equação de Richards, na qual a condutividade hidráulica é função da variável dependente (q ou p/rg), passa a ser linearizada no espaço. E o problema agora consiste em resolver uma equação não-linear no tempo. E é nesse aspecto que o método iterativo de Newton-Raphson apresenta-se particularmente interessante. Para aplicar esse método na (3), é necessário derivá-la com respeito a p_i/rg, o que resulta:

As equações (3) e (4) devem ser resolvidas para se encontrar p_i/rg no final de cada tempo, juntamente com as condições de fronteira. Essas condições, via de regra, são a definição de um fluxo ou de um valor de p/rg no topo ou no fundo do perfil (ou de uma coluna) de solo. Se existe lençol freático em z = i, por exemplo, há que se impor p_i/rg = 0. Em caso contrário, p_i/rg pode assumir um valor constante, previamente estimado ou admitido.

As equações (3) e (4) devem ser utilizadas para avaliar p^j/rg e encontrar p^j+1/rg para cada profundidade de interesse no solo e elas deverão ser novamente utilizadas para reavaliar esses termos quantas vezes forem necessárias (iterações) até que ocorra a convergência, ou seja, quando os dois membros da (3) se igualem, a fim de que o princípio da conservação da matéria seja assegurado.

Com base nas equações (3) e (4), Prevedello (1996) elaborou o programa abaixo, escrito em linguagem BASIC (Beginner's All-purpose Symbolic Instructional Code) e que leva em conta o fluxo de vapor dentro do perfil de solo e a taxa de evaporação na sua superfície. Nesse programa, a condutividade do solo ao vapor de água e a condutividade hidráulica foram somadas para se obter a densidade de fluxo total de água (na forma líquida e de vapor). A densidade de fluxo de vapor dentro do solo é calculada na linha 270 e adicionada na equação do balanço, linha 320. A derivada da densidade de fluxo de vapor também é calculada na linha 270 e está presente na linha 310, para assegurar a convergência nos cálculos dos novos valores de p/rg. A densidade de fluxo de vapor e sua derivada na superfície do solo é calculada na linha 240. Para meios isotérmicos, a densidade de fluxo de vapor (q_v) pode ser calculada (Prevedello, 1996) pela equação:

em que D_v é a difusividade do solo ao vapor (m²/s), c'_v é a concentração de saturação por vapor, UR_es é a umidade relativa do ar na superfície do solo, e z é a coordenada de posição. Por isso, a umidade do solo em cada incremento de profundidade deve ser conhecida. Esse cálculo é inicialmente realizado na linha 150 e, posteriormente, na linha 460. Os resultados dos perfis de umidade e de pressão da água no solo, para cada intervalo de hora, são emitidos pelas linhas 490 a 510. Também para cada intervalo de tempo é calculada, na linha 520, a taxa de evaporação na superfície do solo e impressa, pelo comando da linha 530, em mm/d. A relação evaporação atual/evaporação potencial também é impressa nessa linha.

Os principais passos utilizados na elaboração desse programa são comentados a seguir. O programa simula o processo da evaporação da água do solo com determinada umidade inicial e fornece os perfis de umidade e de pressão da água no solo em intervalos de tempo previamente solicitados, bem como a relação evaporação atual/potencial. Os dados do solo são solicitados da linha 40 até 87. Na linha 100, é solicitado o intervalo de tempo desejado entre os perfis de umidade e de pressão. As linhas 103 e 104 definem as profundidades do solo para as quais os valores de umidade e de pressão serão calculados. As linhas 130 e 150 definem os valores de umidade e de pressão no tempo inicial (t = 0) e a linha 235 calcula os respectivos valores de K. A equação (3) é calculada na linha 320, para cada profundidade, e os valores da função derivada (equação 4), nas respectivas profundidades, são calculados na linha 310. Na linha 330, são acumulados os valores das diferenças entre o primeiro e o segundo membro da (3), para todas as profundidades consideradas. As linhas 350 e 360 calculam os incrementos de pressão que devem ser somados aos valores de pressão do último intervalo de tempo e a linha 465 calcula os valores de umidade do próximo perfil, por meio da equação (1), se a soma acumulada dos valores das diferenças entre o primeiro e o segundo membro da (3) de todas as profundidades não se anula dentro da precisão imposta na linha 200 e verificada pelo comando da linha 470. Se, por outro lado, a condição de convergência for satisfeita, os resultados serão impressos pelo comando das linhas 460 a 485, bem como a taxa de evaporação atual/potencial, a qual é calculada nas linhas 490 e 500. Finalmente, tudo se repete para o novo intervalo de tempo pelo comando da linha 540.

MATERIAL E MÉTODOS

Avaliação experimental do modelo numérico

Para uma avaliação do modelo numérico apresentado no item anterior, foram realizados vários ensaios laboratoriais por meio do processo da evaporação da água do solo, tanto em colunas constituídas de um único material poroso (arenoso ou argiloso), como em colunas formadas por camadas de espessuras e propriedades hidráulicas distintas, alternando-se esses materiais, a fim de que se pudesse avaliar a versatilidade do modelo numérico frente a casos de estratificação e sua exeqüibilidade frente às situações reais de campo.

Todos os ensaios foram efetuados em condições controladas no Laboratório de Física de Solos da UFPR. Foram utilizados dois materiais porosos fisicamente distintos: uma Areia Marinha, do balneário Caiobá, município de Matinhos (PR), e um Latossolo Vermelho-Escuro textura argilosa, do município de Santa Terezinha de Itaipu (PR). A profundidade de coleta esteve entre 0 e 0,25 m. Após secagem ao ar, ambos os materiais foram peneirados em malha de 2 mm e acondicionados em colunas de acrílico com 0,45 m de altura e 50 mm de diâmetro interno. Para minimizar o efeito da segregação de partículas por ocasião do preenchimento das colunas, utilizou-se um tubo extensor munido de um funil pelo qual o material poroso era continuamente depositado, de tal forma que o tubo ficasse permanentemente completo, evitando-se, dessa forma, que as partículas caíssem em queda livre. Foram montadas quatro a cinco colunas para cada tipo de material poroso.

As propriedades hidráulicas do solo foram obtidas a partir da curva de retenção de água por meio de mesa de tensão, empregando-se papel "mata-borrão" de 1,0 m de coluna de H₂O de pressão de borbulhamento e câmara de Richards. O solo foi acondicionado em anéis de 36 mm de diâmetro e 25 mm de altura e submetido a tensões entre 0,01 e 0,81 m de coluna de H₂O, na mesa de tensão, e entre 1,0 e 150 m de coluna de H₂O, na câmara de Richards. Para cada tensão foram utilizadas seis repetições. Os valores de umidade com base na massa seca (U) foram convertidos à base de volume (q) (q = U.r_s/r_água, em que r_s e r_água são, respectivamente, a massa específica do solo seco e da água) e, posteriormente, ajustados ao modelo de Van Genuchten (1980), pela equação (1). Foi utilizado um programa de regressão não-linear para ajustar os resultados experimentais à equação (1), conforme Boratto (1984). Desse ajuste, isto é, com os parâmetros independentes a, q_r, n e m conhecidos, além de K_s, a função K(q) foi obtida por meio da equação (2). A condutividade hidráulica saturada (K_s) foi obtida mediante aplicação da equação de Darcy (1856) nas amostras acondicionadas nos anéis anteriormente mencionados.

Nos experimentos de evaporação, foram obtidos quatorze perfis de umidade, correspondentes às seguintes situações: (a) cinco colunas de material homogêneo de Areia Marinha nos tempos de evaporação de 0; 12; 36; 72 e 240 min; (b) quatro colunas de material homogêneo de Latossolo nos tempos de 0; 12; 36; 72 e 240 min; e (c) quatro colunas estratificadas (com camadas alternadas de ambos os materiais) nos tempos de 0; 12; 36 e 240 min. Cada coluna foi saturada com água e o excesso deixado drenar livremente, quando então a extremidade inferior da coluna era vedada com um filme plástico e iniciada a contagem do tempo, sujeitando-a a uma intensidade de evaporação aproximadamente constante e regulada pela velocidade máxima de um ventilador. A taxa de evaporação alcançada por esse artifício foi medida com auxílio de um evaporímetro Piché, com a superfície evaporante na altura da superfície da coluna. Transcorrido cada tempo destinado à evaporação, a coluna era posicionada horizontalmente e amostras eram coletadas nas profundidades de 0, 0,03, 0,08; 0,10; 0,15; 0,20; 0,25; 0,30 e 0,35 m para determinação da umidade volumétrica. Esses perfis experimentais de umidade foram então comparados com os respectivos perfis obtidos pela modelagem numérica e confrontados estatisticamente pelo teste de "Kolmogorov-Smirnov" (K-S) do programa estatístico Statgraphics. Este procedimento calcula o máximo desvio entre as funções de distribuição acumulada de duas amostras. Se este desvio for suficientemente grande, a hipótese de que as distribuições são as mesmas é rejeitada. Neste procedimento, se o nível de significância cai abaixo de 0,05, então as duas distribuições são significativamente diferentes a 5 %. Segundo Mendenhall (1985), as pesquisas têm mostrado que os testes não-paramétricos são praticamente tão úteis para detectar diferenças entre populações quanto os métodos paramétricos, desde que para esses últimos seja satisfeita a suposição da distribuição normal. Segundo esse autor, os testes não-paramétricos são, freqüentemente, muito úteis e mais poderosos do que os paramétricos, quanto à detecção de diferenças entre populações, razão por que muitos estatísticos advogam o emprego dos testes não-paramétricos preferencialmente aos procedimentos paramétricos.

RESULTADOS E DISCUSSÃO

Propriedades hidráulicas dos meios porosos utilizados

O quadro 1 reúne os valores de condutividade hidráulica saturada e dos parâmetros independentes obtidos por ajuste do modelo de Van Genuchten (1980) aos valores experimentais da curva de retenção. A figura 1 mostra os pontos experimentais e a curva de retenção ajustada por esse modelo para o Latossolo, e a figura 2 para a Areia Marinha.

Como se observa, o modelo de Van Genuchten (1980) ajustou-se bem aos valores experimentais dos dois meios porosos, com coeficiente de determinação acima de 0,988.

Os valores do quadro 1 foram então introduzidos no programa para simulação do processo da evaporação, para cada um dos casos formulados. Os resultados dos perfis de umidade simulados e observados da evaporação na Areia Marinha e no Latossolo (colunas homogêneas) são mostrados nas figura 3 e 4, respectivamente. No caso da Areia Marinha (Figura 3), percebe-se que os perfis se transladaram bem no tempo, aproximando-se muito bem dos valores experimentais ao longo de toda a coluna, não diferindo estatisticamente, pelo teste K-S, a 5 %. No Latossolo (Figura 4), apesar de os perfis simulados se aproximarem relativamente bem dos valores medidos abaixo de 0,10 m na camada superficial eles ficaram prejudicados. A causa desse pobre sajuste na camada de 0-0,10 m pode ser devida a uma limitação experimental inerente do método de obtenção da curva de retenção de água pelas câmaras de Richards, em que os máximos valores de tensão são da ordem de 1.500 kPa (» 150 m de coluna de H₂O), valor esse que nem sempre extrai toda a água das amostras, principalmente dos meios porosos argilosos.

Como conseqüência, a umidade residual q_r para esses meios assume valor relativamente alto. Acontece, entretanto, que, no processo da evaporação, a superfície do solo seca muito e a umidade atinge valores próximos de zero, enquanto, por simulação, o menor valor de umidade que o solo pode atingir nunca poderá ser inferior a q_r. De fato, os valores medidos de umidade na superfície das colunas chegaram a acusar 1 a 2 % (em volume) para tempos de 72 e 240 h, enquanto se observou (Quadro 1) que o valor da umidade residual q_r da Areia Marinha foi inferior a 2 % e para o Latossolo foi superior a 19 %. Mesmo assim, o teste K-S não acusou diferença significativa, a 5 %, em nenhum dos perfis.

A figura 5 mostra os resultados simulados e observados nos diversos períodos de tempo de evaporação em colunas estratificadas, constituídas de 2 camadas: a primeira, de areia, da superfície até os primeiros 120 mm, e a segunda, de material argiloso, de 120 a 450 mm. Verificou-se que, de modo geral, os resultados foram satisfatórios em todos os casos simulados, como também ficou confirmado pelo teste estatístico K-S.

Kolasew (1941) e Lemon (1956) foram os primeiros a reconhecerem que o processo da evaporação de um solo com alta umidade e sem vegetação passa por três estádios distintos. No primeiro estádio, a velocidade de evaporação é constante e independente da umidade do perfil (ou da coluna) de solo. Nesse estádio, a evaporação depende das condições reinantes externamente, tais como: velocidade do vento, temperatura, umidade do ar e energia radiante. Esse estádio termina quando se estabelece uma resistência ao fluxo da água na superfície do solo e a velocidade de evaporação decresce. Inicia-se, então, o segundo estádio, no qual a velocidade de evaporação é função linear da umidade do solo e as condições reinantes externamente não são mais de importância fundamental, enquanto as condições intrínsecas do solo governam o transporte da água dentro do solo e, conseqüentemente, a velocidade de evaporação. Quando a função que correlaciona a velocidade de evaporação com a umidade começa a perder linearidade, inicia-se o terceiro estádio, que se caracteriza por um movimento bastante lento da água dentro do solo, decorrente de forças de adsorção que aparecem entre a água e as partículas sólidas do solo e de uma condutividade hidráulica muito baixa. Durante o primeiro estádio, a evaporação é máxima (potencial).

A figura 6 mostra a velocidade de evaporação simulada para as duas situações consideradas anteriormente: coluna homogênea e estratificada, ambas partindo da condição de umidade próxima da saturação (0,05 m de coluna de H₂O de tensão). Nota-se que a duração do primeiro estádio foi maior para a condição de coluna homogênea. A areia sobre material argiloso, nas colunas estratificadas, fez com que essas colunas não percebessem o primeiro estádio da evaporação. Semelhante situação foi observada nos experimentos de evaporação em colunas realizados por Reichardt (1972).

A rápida passagem do primeiro para o segundo estádio nas colunas estratificadas é decorrente da rápida secagem da capa de areia (120 mm), que, já nas primeiras 12 h de evaporação, teve a umidade decrescida para valores inferiores a 0,10 m³ m^-³, comportando-se como uma capa isolante ("mulching"), iniciando-se, então, uma forte resistência ao fluxo da água na superfície dessas colunas. O teste estatístico "Kolmogorov-Smirnov" revelou não haver diferença significativa entre os valores estimados e observados, a 5 %, em todos os tempos de evaporação considerados.

CONCLUSÃO

O modelo proposto, resumido na forma das equações (3) e (4) ou do programa computacional anexo, foi capaz de simular o processo unidimensional da evaporação da água do solo, tanto em meios porosos homogêneos quanto estratificados.

AGRADECIMENTO

Ao CNPq, pela concessão de bolsa de auxílio ao primeiro autor.

LITERATURA CITADA

Recebido para publicação em março de 2004 e aprovado em março de 2005.

• BORATTO, F. BASIC para engenheiros e cientistas. 2.ed. Rio de Janeiro, Livros Técnicos e Científicos, 1984. 120p.
• KOLASEW, F.E. Ways of suppressing evaporation of soil moisture. Sborn. Rab. Agron. Fiz., 3:67, 1941.
• LEMON, E.R. The potencialities for decreasing soil moisture evaporation loss. Soil Sci. Soc. Am. Proc., 20:120-125, 1956.
• MARION, J.M.; OR, D. & ROLSTON, D.E. Evaluation of methods for determining soil-water retentivity and unsaturated hydraulic conductivity. Soil. Sci., 18:1-13, 1994.
• MENDENHALL, W. Probabilidade e estatística. Vol.2. Tradução de José Fabiano da Rocha. Rio de Janeiro, Campus, 1985. 290p.
• MUALEM, Y. A new model for predicting the hydraulic conductivity of unsaturated porous media. Water Res. Res., 12:513-522, 1976.
• PREVEDELLO, C.L.; KOBIYAMA, M.; JACOBS, G.A. & DIVARDIN, C.R. Comparação dos métodos do perfil instantâneo e de van Genuchten na obtenção da condutividade hidráulica de uma areia marinha. R. Bras. Ci. Solo, 19:1-5, 1995.
• PREVEDELLO, C.L. Física do solo com problemas resolvidos. Curitiba, Salesward-Discovery, 1996. 446p.
• REICHARDT, K. Evaporação da água do solo. Piracicaba, CENA, 1972. 30p. (Boletim Científico do CENA)
• SIMUNEK, J.; VAN GENUCHTEN, M.T.; GRIBB, M.M. & HOPMANS, J. Parameter estimation of unsaturated soil hydraulic properties from transient flow processes. Soil Tillage Res., 47:27-36, 1998.
• VAN GENUCHTEN, M. Th. A closed-form equation for predicting the hidraulic conductivity of unsaturated soils. Soil Sci. Soc. Am. J., 44:892-898, 1980.

Datas de Publicação

• Publicação nesta coleção
  03 Mar 2006
• Data do Fascículo
  Out 2005