Uma Aplicação dos Algoritmos Genéticos e do Método dos Volumes Finitos para Simular o Fluxo da Água na Zona Vadosa

RESUMO

O fluxo da água em meios porosos é governado por uma equação diferencial parcial (a equação de Richards), cuja forma mista envolve as variáveis umidade do solo e potencial mátrico. A curva de retenção é uma relação não linear entre estas variáveis, fundamental no estudo da dinâmica da água no solo na zona vadosa. Um dos objetivos deste trabalho é avaliar o modelo de curva de retenção da água no solo que melhor ajusta os dados mensurados de umidade e potencial mátrico, obtidos a partir de um experimento realizado na Universidade Federal Rural do Rio de Janeiro. No presente artigo, um Algoritmo Genético (AG) é proposto de forma a buscar os parâmetros de ajuste que maximizam o coeficiente de determinação, considerando três horizontes de solo e os seguintes modelos de curvas de retenção: van Genuchten, Brooks-Corey e Haverkamp. A performance do AG implementado é avaliada comparando os resultados obtidos com o programa SWRC Fit, que usa o método determin´ıstico de Levenberg-Marquardt. A partir dos resultados obtidos pode-se observar que o AG ajustou com maior precisão os dados mensurados e o modelo Haverkamp apresentou o maior coeficiente de determinação. Além disso, o esquema de discretização da equação de Richards proposto se apresentou mais estável usando o modelo Haverkamp como curva de retenção da água no solo.

Palavras-chave:
ensaios laboratoriais; equação de Richards; regressão não linear

ABSTRACT

The water flux in porous media obeys a nonlinear partial differential equation (the Richards equation), whose the mixed form contains the variables volumetric water content and matric potential. The soil water retention curve is the nonlinear relation between the volumetric water content and matric potential, essential for analysing the soil water dynamics in the vadose zone. The main purpose of this paper is to evaluate the best fit model for the retention curve considering experimental analysis realized at Universidade Federal Rural do Rio de Janeiro. In this work, a genetic algorithm (GA) is proposed in order to search the design variables that optimize the coefficient of determination (R Squared), considering three soil horizons and the van Genuchten, Brooks-Corey and Haverkamp models to describe the retention curve. The GA proposed was evaluate comparing the results with the software SWRC Fit , which performs nonlinear fitting of soil water retention curves by Levenberg-Marquardt method. GA adjusted the measured data more precisely and the Haverkamp model presented the highest coefficient of determination. In addition, the discretization for the equation Richards was more stable using the Haverkamp model as a soil water retention curve.

Keywords:
laboratory tests; Richards equation; non-linear regression

Na agricultura, a importância do solo se deve à grande diversidade de culturas que se desenvolve sobre ele. Vale salientar que boa parte dessas culturas agrícolas consolidam seu desenvolvimento na zona não saturada do solo ²¹ . Por este motivo, o desenvolvimento de modelos que sejam capazes de descrever e quantificar o fluxo de água através do solo corrobora na otimização de práticas agrícolas.

A equação de Richards governa o movimento da água em um meio poroso saturado/insaturado e, em geral, não possui solução analítica. O método das diferenças finitas ⁵ , o método dos elementos finitos ⁸ e, mais recentemente, o método dos volumes finitos ^{7), (29), (35} são os mais utilizados para resolver de forma numérica a equação de Richards. No entanto, uma solução numérica robusta, capaz de propor acurácia e eficiência para problemas multidimensionais, continua sendo um desafio, devido ao comportamento não-linear da equação ¹⁵ . Em geral, os aspectos numéricos analisados incluem: conservação de massa, estabilidade, acurácia, eficiência, continuidade a partir do regime insaturado para o saturado (saturação variável) e a não linearidade da condutividade hidráulica ⁴ .

Uma das ferramentas mais utilizadas para descrever e predizer o processo de infiltração da água no solo é o software Hydrus , desenvolvido por Simunek, van Genuchten & Sejna ²⁸ . Tal programa é um modelo de elementos finitos, que numericamente resolve a equação de Richards para o fluxo saturado/insaturado da água no solo e a equação de convecção-dispersão para o transporte de solutos. Uma abordagem compartimental alternativa em que para a redistribuição da água no solo é feita uma analogia com os modelos de competição em dinâmica de populações foi apresentada em ²⁴ , cujos resultados, apesar de preliminares, se mostraram promissores.

A determinação e a análise das relações entre as variáveis umidade e potencial mátrico tem significativa importância, visto que decorre delas a não linearidade da equação de Richards. As funções empíricas de van Genuchten ³¹ e Brooks-Corey ³ são comumente usadas para representar analiticamente a curva de retenção da água no solo e a condutividade hidráulica. As expressões analíticas apresentadas em ¹² para as curvas de retenção e de condutividade podem ser vistas como uma versão mais simplificada do modelo de van Genuchten, mas que ainda preservam as características de uma curva logística ( sigmoid function ).

A estimativa dos parâmetros que caracterizam os diferentes modelos de curvas de retenção da água no solo é realizada, tradicionalmente, a partir da câmara de pressão de Richards ¹⁶ , que vai mensurar umidade em função do potencial mátrico, considerando amostras de solo inicialmente saturadas. Em seguida, um problema de quadrados mínimos não linear pode ser modelado e resolvido a fim de buscar os parâmetros que melhor ajustam as curvas de retenção da água no solo. No presente artigo, a qualidade da solução obtida por um algoritmo genético (AG) ¹⁰ foi avaliada, comparando os resultados com o programa já altamente testado SWRC Fit²⁶ , que usa o método determinístico de Levenberg-Marquardt ²⁰ . AGs pertencem à classe das técnicas evolucionárias, que não requer nenhum cálculo de derivadas, cuja otimização é obtida usando apenas a função objetivo. Os procedimentos envolvidos nos AGs são probabilísticos, baseados em métodos de busca global, o que pode tornar a convergência mais lenta quando comparados com os métodos determinísticos. Além disso, nos AGs, as variáveis de projeto estão sujeitas a restrições de domínio. Seja de representação binária, real ou de forma híbrida, os AGs têm sido aplicados com sucesso em problemas de otimização e análise inversa nos diferentes campos de pesquisa, tais como danos estruturais ²² , proteção catódica ²³ e controle populacional ²⁷ .

O objetivo deste trabalho é verificar se um modelo mais simples e pouco utilizado de curvas de retenção pode ser usado para ajustar dados experimentais de umidade e potencial mátrico a partir da implementação de um AG eficiente. Apesar de ser o menos citado na literatura, o modelo Haverkamp apresentou o melhor ajuste quando comparado com os modelos mais populares, considerando os dados mensurados em laboratório e o AG implementado. Simulações numéricas para a equação de Richards foram realizadas, usando um esquema de discretização explícito no tempo e o método dos volumes finitos para a discretização espacial. Também foi possível utilizar um passo de tempo maior com o uso do modelo Haverkamp, de forma que ainda garantisse a taxa de convergência esperada para o método de Euler explícito, baseada no erro de truncamento do esquema de discretização.

1 MODELO MATEMÁTICO

A equação de Richards descreve o balanço hídrico no solo e define a quantidade relativa de fluido (umidade) como função de alguns parâmetros, tais como condutividade e potencial mátrico. Basicamente, ao despejar uma certa quantidade de água no solo, por diferença de potencial, este fluido se dispersará com maior (ou menor) dificuldade dependendo de propriedades específicas do solo estudado.

Sejam $\theta =\theta (t,z)[L^3/L^3]$ a umidade volumétrica, definida como o volume de água que atravessa um volume total de solo, função do tempo t [ s ] e do espaço z [ L ], $K=K\left(h\right)[L/T]$ a condutividade hidráulica do solo e $h=h\left(\theta \right)\left[L\right]$ o potencial mátrico em dimensões de altura. Como é de praxe, combinando a equação do fluxo, no caso a equação de Darcy-Buckingham para o fluxo da água, com a equação da continuidade, pode-se chegar na seguinte equação de Richards, na direção z ,

$\frac{\partial \theta }{\partial t}=\frac{\partial }{\partial z}\left[K\left(h\right)\left(1+\frac{\partial h}{\partial z}\right)\right],$ (1.1)

também conhecida como forma mista pois envolve tanto θ quanto h em sua formulação. Esta costuma ser a formulação mais utilizada por aliar, na teoria, a conservação de massa da forma θ e a continuidade a partir do regime insaturado para o saturado (saturação variável) da forma h .

A forma θ da equação de Richards é dada como segue

$\frac{\partial \theta }{\partial t}=\frac{\partial }{\partial z}\left[K\left(\theta \right)+D\left(\theta \right)\frac{\partial \theta }{\partial z}\right],$ (1.2)

em que $D\left(\theta \right)=K\left(\theta \right)\frac{\partial h}{\partial \theta }[L^2/T]$ representa a difusividade da água no solo, que pode ser obtida aplicando-se a regra da cadeia ao termo $\frac{\partial h}{\partial z}$ em (1.1) .

A equação (1.1) em sua forma h é escrita como

$C\left(h\right)\frac{\partial h}{\partial t}=\frac{\partial }{\partial z}\left[K\left(h\right)\left(1+\frac{\partial h}{\partial z}\right)\right],$ (1.3)

onde $C\left(h\right)=\frac{\partial \theta }{\partial h}[1/L]$ é chamada capacidade hídrica e pode ser obtida desenvolvendo o termo $\frac{\partial \theta }{\partial t}$ através da regra da cadeia, pois também $\theta =\theta \left(h\right)$ . Neste trabalho, apenas a forma mista da equação de Richards é avaliada. As formas θ e h serão abordadas do ponto de vista numérico, para fins comparativos, em trabalhos futuros.

2 DISCRETIZAÇÃO DA EQUAÇÃO DE RICHARDS

Neste artigo, a discretização temporal é realizada através do método de Euler explícito e o método dos volumes finitos (MVF) é usado para a discretização espacial. O MVF inicia com a equação diferencial na forma conservativa, e tem por objetivo integrá-la sobre o volume elementar, no espaço e no tempo. É importante observar que a forma conservativa de uma equação é aquela onde os fluxos aparecem dentro do sinal da derivada ¹⁸ .

No problema abordado, a equação diferencial (1.1) será utilizada. Devido ao balanço de massa usado na concepção da equação de Richards, esta já se encontra na forma conservativa. O volume elementar considerado, não será exatamente um volume, mas dois pontos, pois o objetivo é fazer um cálculo unidimensional, ou seja, ao longo do segmento de reta que inicia no ponto z = 0 e se prolonga até o limite do perfil de solo em questão, onde z = L .

Demais esquemas de discretização temporal, implícitos e semi-implícitos por exemplo, também serão abordados em trabalhos futuros.

2.1 Discretização Espacial

Dividindo o domínio em N células não sobrepostas $[z_i,z_{i+1}]$ e integrando a forma mista da Equação de Richards (1.1) sobre a célula i tem-se

${\int }_{\Delta z}\frac{\partial \theta }{\partial t}dz={\int }_{\Delta z}\frac{\partial }{\partial z}\left[K\left(\theta \right)\left(1+\frac{\partial h}{\partial z}\right)\right]dz,$ (2.1)

onde $\Delta z=z_{i+1}-z_i$ é o tamanho do passo da malha espacial. O desenho da célula i é mostrado na Figura 1 .

Figura 1
Esquema de volumes finitos.

Ao lado esquerdo da equação (2.1) aplica-se o Teorema do Valor Médio para Integrais, obtendo a derivada $\frac{\partial \theta }{\partial t}$ que será aproximada pelo método explícito de Euler e, além disso, desconsidera-se a dependência de K e θ em h para que, aplicando o Teorema Fundamental do Cálculo para o termo à direita da equação, tenha-se:

(2.2) $\begin{array}{ccc}\frac{\partial {\theta }_i}{\partial t}{\int }_{\Delta z}dz& =& K\left(\theta \right)\left(1+\frac{\partial h}{\partial z}\right)\overline{)\begin{array}{c}{z}_{i+1/2}\\ z_{i-1/2}\end{array}}\\ & =& K_{i+1/2}\left(1+\frac{\partial h}{\partial z}\left(t,z_{i+1/2}\right)\right)-K_{i-1/2}\left(1+\frac{\partial h}{\partial z}\left(t,z_{i-1/2}\right)\right),\end{array}$

onde $K_{i-1/2}$ representa a condutividade hidráulica no ponto médio da interface que a célula i − 1 compartilha com a célula i . De igual modo, $K_{i+1/2}$ é a condutividade no ponto médio da interface entre as células i e i +1.

Para simplificar a notação, será usado:

$\frac{\partial h}{\partial z}(t,z_{i\pm 1/2})=\frac{\partial h_{i\pm 1/2}}{\partial z}.$ (2.3)

Assim, obtém-se

(2.4) $$\begin{gathered} \frac{\partial {\theta }_i}{\partial t}\Delta z=K_{i+1/2}+K_{i+1/2}\frac{h_{i+1}-h_i}{z_{i+1}-z_i}-K_{i-1/2}-K_{i-1/2}\frac{h_i-h_{i-1}}{z_i-z_{i-1}} \\ =K_{i+1/2}+K_{i+1/2}\frac{h_{i+1}-h_i}{\Delta z}-K_{i-1/2}-K_{i-1/2}\frac{h_i-h_{i-1}}{\Delta z} \\ =\frac{K_{i+1/2}{h}_{i+1}-(K_{i+1/2}+K_{i-1/2})h_i+K_{i-1/2}{h}_{i-1}}{\Delta z}+K_{i+1/2}-K_{i-1/2}. \end{gathered}$$

Finalmente, a discretização espacial para a forma mista da Equação de Richards pode ser escrita como

(2.5) $\frac{\partial {\theta }_i}{\partial t}=\frac{1}{\Delta z^2}\left[K_{i+1/2}{h}_{i+1}-(K_{i+1/2}+K_{i-1/2})h_i+K_{i-1/2}{h}_{i-1}\right]+\frac{K_{i+1/2}-K_{i-1/2}}{\Delta z}.$

A literatura relata várias formas de calcular os valores das variáveis $K_{i-1/2}e{K}_{i+1/2}$ , dentre elas destacam-se as médias aritmética, geométrica, harmônica e a ascendente, como pode ser visto nas referências ^{1), (2), (9), (11), (13), (30), (32} . Neste trabalho, optou-se por realizar o cálculo das condutividades através da média aritmética, devido à sua ampla utilização. Dessa forma,

$K_{i-1/2}=\frac{K_i-K_{i-1}}{2},$ (2.6)

$K_{i+1/2}=\frac{K_{i+1}-K_i}{2}.$ (2.7)

2.2 Discretização Temporal

O método de Euler avançado faz uso da seguinte aproximação para a derivada da umidade em relação ao tempo

$\frac{\partial {\theta }_i}{\partial t}\approx \frac{{\theta }_i^{n+1}-{\theta }_i^n}{\Delta t},$ (2.8)

onde Δ t é o tamanho do passo de tempo.

Substituindo (2.8) na equação (2.5) , obtém-se

(2.9) $\frac{{\theta }_i^{n+1}-{\theta }_i^n}{\Delta t}=\frac{1}{\Delta z^2}\left[K_{i+1/2}^n{h}_{i+1}^n-(K_{i+1/2}^n+K_{i-1/2}^n)h_i^n+K_{i-1/2}^n{h}_{i-1}^n\right]+\frac{K_{i+1/2}^n-K_{i-1/2}^n}{\Delta z},$

e isolando o termo chega-se ao esquema explícito de discretização temporal para a forma mista da equação de Richards.

(2.10) ${\theta }_i^{n+1}=\frac{\Delta t}{\Delta z^2}\left[K_{i+1/2}^n{h}_{i+1}^n-(K_{i+1/2}^n+K_{i-1/2}^n)h_i^n+K_{i-1/2}^n{h}_{i-1}^n\right]+\frac{\Delta t}{\Delta z}[K_{i+1/2}^n-K_{i-1/2}^n]+{\theta }_i^n.$

3 MODELOS PARA A CONDUTIVIDADE HIDRÁULICA

Não é possível resolver a equação de Richards sem conhecer previamente a curva de retenção da água no solo $\theta =\theta \left(h\right)$ e a condutividade hidráulica $K=K\left(\theta \right(h\left)\right)$ . Uma das alternativas é o cálculo teórico da condutividade hidráulica a partir de dados da curva de retenção da água no solo, mais facilmente medidos no campo ou no laboratório ¹⁸ .

Para confeccionar as curvas de retenção em laboratório usou-se o método da câmara de pressão de Richards. Este consiste em colocar amostras de solo em um recipiente específico que permite aplicar diversas pressões a essas amostras até que não saia mais água do equipamento. Neste ponto, obtém-se um equilíbrio entre o potencial matricial das amostras e a pressão aplicada. Em seguida, mede-se o teor de umidade do solo. Esse procedimento foi repetido para pressões entre 0 e 15 atm, visto que o primeiro equivale ao ponto em que o solo está saturado e o segundo é considerado o ponto em que as amostras atingem a umidade residual. Os valores de umidade obtidos em função das pressões entre 0 e 15 atm formam a curva (discreta) de retenção da água no solo ¹⁴ .

O experimento da panela de pressão de Richards foi conduzido nos laboratórios de física do solo da Universidade Federal Rural do Rio de Janeiro. O solo utilizado foi um Planossolo Háplico provenientes do município de Seropédica - RJ. Foram coletadas três amostras deformadas referentes aos horizontes A (0-22 cm), E (22-69 cm) e B (69-98 cm) do perfil de solo considerado, as quais foram submetidas ao procedimento para a obtenção de ( h ,θ).

Na prática da física do solo, diversos modelos para curva de retenção são encontrados na literatura ^{25), (33), (34} . Neste trabalho, foram utilizados os modelos de Haverkamp, van Genuchten e Brooks-Corey dados, respectivamente pelas equações (3.1) , (3.2) , (3.3) , como segue.

$\theta \left(h\right)=\frac{A({\theta }_s-{\theta }_r)}{A+|h{|}^{\beta }}+{\theta }_r,$ (3.1)

$\theta \left(h\right)=\frac{{\theta }_s-{\theta }_r}{[1+(\alpha \left|h\right|{)}^n{]}^m}+{\theta }_r,$ (3.2)

em que $m=1-\frac{1}{n}$ ,

$\theta \left(h\right)=\left\{\begin{array}{l}\frac{{\theta }_s-{\theta }_r}{{\left({\displaystyle \frac{h}{h_b}}\right)}^{\lambda }},h>h_b,\\ {\theta }_s,h\le h_b,\end{array}\right.$ (3.3)

Nas expressões das curvas de retenção existem parâmetros comuns tais como θ _s e θ _r , chamados de umidade de saturação e umidade residual, respectivamente. Estes são relativos a cada solo e podem ser obtidos experimentalmente. Os demais são considerados parâmetros de ajuste da curva, que também armazenam informações acerca das características físicas do solo. Em (3.2) , $\alpha >0\left[L^{-1}\right]$ é considerado o inverso da densidade de poros e n ≥ 1 é adimensional. Em (3.3) , h _b [ L ] é a pressão de borbulhamento, isto é, a pressão aplicada para retirar as moléculas de água dos maiores poros do solo. Neste trabalho, todos os parâmetros são obtidos numericamente a partir da solução de um problema de mínimos quadrados não linear, descrito na próxima seção.

Cada uma das curvas apresentadas acima possui expressão para representar K ( θ ). Para (3.2) , a seguinte equação pode ser obtida de forma teórica ¹⁹

$K\left(\theta \right)=K_{sat}{S}^{1/2}{\left[1-{\left(1-S^{1/m}\right)}^m\right]}^2,$ (3.4)

onde $K_{sat}=K\left({\theta }_s\right)$ é a condutividade de saturação e <mml:math><mml:mi>S</mml:mi><mml:mo>=</mml:mo><mml:mi>S</mml:mi><mml:mo>(</mml:mo><mml:mi>θ</mml:mi><mml:mo>)</mml:mo></mml:math> é dado por

$S\left(\theta \right)=\frac{\theta -{\theta }_r}{{\theta }_s-{\theta }_r}.$ (3.5)

O modelo apresentado na equação (3.1) estabelece a seguinte relação entre K e h

$K\left(h\right)=K_{sat}\frac{A}{A+|h{|}^{\beta }}.$ (3.6)

Já para o modelo apresentado em (3.3) , a equação para a condutividade hidráulica fica

$K\left(h\right)=K_{sat}{\left(\frac{h_b}{h}\right)}^{\rho },$ (3.7)

sendo $\rho =2+3\lambda$ .

4 AGS E O PROBLEMA DE QUADRADOS MÍNIMOS NÃO LINEAR

Algoritmos genéticos (AGs) fazem parte do ramo dos algoritmos evolucionários e, como tais, podem ser definidos como uma técnica de busca baseada em uma metáfora do processo ecológico de evolução natural ⁶ .

Nos AGs, populações de indivíduos (cromossomos) são criadas e submetidas aos operadores genéticos: seleção, recombinação ( crossover ) e mutação. Estes operadores utilizam uma caracterização da qualidade de cada indivíduo como solução do problema em questão através de uma função de avaliação ou da função objetivo. Com isto, vão gerar um processo de evolução natural destes indivíduos que, eventualmente, deverá gerar um indivíduo que caracterizará uma boa solução ao problema.

Em problemas de quadrados mínimos, a função objetivo f possui a seguinte forma especial

$f\left(\mathbf{x}\right)=\frac{1}{2}\sum _{j=1}^m{r}_j^2\left(x\right),$ (4.1)

onde cada r _j ( x ) é uma função resíduo que leva ${\mathfrak{R}}^{n}a\mathfrak{R}$ . Considera-se m ≥ n .

Aqui, o objetivo é estimar os n parâmetros que caracterizam as curvas de retenção $\theta =\theta \left(h\right)$ de forma que ajustem os dados observados em laboratório, ou seja, um conjunto de pares da forma $(h^j,{\theta }^j),j=1,\cdots ,m$ , onde m é o número de dados observados. Tais parâmetros são as denominadas variáveis de projeto do problema de otimização, no caso do modelo Haverkamp definidas pelo vetor $x=(\alpha ,\beta ,{\theta }_r,{\theta }_s)$ . Assim, o resíduo fica definido como

$r_j\left(\mathbf{x}\right)=\left({\theta }_j^{OBS}-{\theta }_j^{AG}\right),$ (4.2)

onde θ ^OBS são os valores de umidade observados e θ ^AG é a umidade estimada em cada cromossomo na população do AG.

Fundamentalmente, o que liga o problema de otimização ao AG é a função objetivo. Desta maneira, é proposta uma nova metodologia de tal forma que o problema de minimização acima se torna em um problema de maximização, cuja função objetivo a ser tratada é o coeficiente R ², de valor máximo igual a 1, podendo ser escrito pela seguinte equação

$R^2\left(\mathbf{x}\right)=1-\frac{{\displaystyle \sum _{i=1}^N({\theta }_i^{OBS}-{\theta }^{AG}{)}^2}}{{\displaystyle \sum _{i=1}^N({\theta }_i^{OBS}-\overline{\theta }{)}^2}},$ (4.3)

onde $\overline{\theta }$ é a média dos dados observados.

Neste artigo, a otimização da função objetivo (4.3) é alcançada pela implementação em JAVA de um AG com representação binária que busca maximizadores, como apresentado em ¹⁷ . Contudo, algumas características foram incluídas como o elitismo, o cruzamento uniforme e a probabilidade de mutação variar linearmente ao longo das gerações. Além disso, pode-se mostrar que os maximizadores de uma função g são equivalentes aos minimizadores da função h quando h =− g . Muitas vezes, ainda, pode ser necessário trabalhar com uma função objetivo toda positiva. Em casos deste tipo, pode-se adicionar uma constante positiva C a função objetivo, que seja suficiente para torná-la toda positiva, visto que $maxg\left(x\right)=max\left\{g\left(x\right)+C\right\}$ .

5 RESULTADOS E DISCUSSÃO

No primeiro exemplo, os dados mensurados ( h , θ ) são ajustados pelo AG proposto e por um método determinístico, considerando os três modelos de curvas já apresentados durante o texto. No segundo exemplo, propõe-se uma simulação numérica da equação de Richards considerando os dados mensurados e os esquemas de discretização também já descritos.

5.1 Exemplo 1

As Figuras 2 , 3 e 4 mostram os valores obtidos em laboratório de ( h , θ ), assim como as curvas de retenção, para os três horizontes de solo analisados, estimadas pelo AG após 100 gerações e considerando os modelos van Genuchten (VG), Brooks-Corey (BC) e Haverkamp (HK). Por um problema técnico na análise laboratorial, não foi possível medir a umidade referente ao potencial mátrico $h=1000cm{H}_{2}O$ no Horizonte B. Dessa forma, optou-se em retirar este ponto nos outros dois horizontes para a estimativa das variáveis de projeto. Após o ajuste das curvas, estes pontos foram inseridos nos gráficos dos Horizontes A e E, verificando-se que o problema técnico relatado não causou prejuízos no ajuste realizado. A partir da Tabela 1 e o coeficiente de determinação R ², pode-se verificar que o uso do modelo Haverkamp (HK) com o AG apresentou o melhor ajuste em relação aos valores observados, visto que se demonstrou mais eficiente que o SWRC Fit e os modelos van Genuchten (VG) e Brooks-Corey (BC), que estão disponíveis no programa.

Figure 2
Curvas de retenção da água no horizonte A usando o AG.

Figure 3
Curva de retenção da água no horizonte E usando o AG.

Figure 4
Curva de retenção da água no horizonte B usando o AG.

Table 1
Análise dos coeficientes de determinação do AG com o SWRC Fit .

Horizonte A	VG	BC	HK
AG	R2 = 0;99041	R2 = 0;99044	R2 = 0;99208
SWRC Fit	R2 = 0;98667	R2 = 0;98917
Horizonte E	VG	BC	HK
AG	R2 = 0;99530	R2 = 0;99531	R2 = 0;99617
SWRC Fit	R2 = 0;99485	R2 = 0;99487
Horizonte B	VG	BC	HK
AG	R2 = 0;93624	R2 = 0;90318	R2 = 0;97798
SWRC Fit	R2 = 0;88933	R2 = 0;83278

5.2 Exemplo 2

Neste exemplo, simulações numéricas da equação de Richards são propostas considerando as características encontradas no Horizonte E (maior R ²), composto principalmente de areia, com profundidade de 22 a 69 cm. Foram prescritas condições de contorno de Dirichlet para $h=h(z,t)$ , com potenciais prescritos $h(22,t)=-0,001cmeh(69,t)=-350cm$ . A condição inicial para o Horizonte E foi de $h(z,0)=-350cm$ . A condutividade de saturação para o Horizonte E é de $K_{sat}=0,0067cm·s^{-1}$ .

Primeiramente, a simulação numérica é proposta considerando o modelo van Genuchten e os parâmetros obtidos pelo AG no exemplo anterior, para as curvas de retenção da água no solo e da condutividade hidráulica. As variáveis de projeto estimadas pelo AG após 100 gerações para o modelo van Genuchten foram ${\theta }_s=0,399c{m}^3·c{m}^{-3},{\theta }_r=3,246\times {10}^{-4}c{m}^3·c{m}^{-3},\alpha =1,347c{m}^{-1}en=1,228$ . O teor de umidade para o Horizonte E, obtido usando um passo de tempo ${\Delta }_t=0,0006s$ após 300 s é apresentado no gráfico na Figura 5 . O tamanho da malha espacial é de $\Delta z=1cm$ .

Figure 5
Teor de umidade para o horizonte E.

Visto que a consistência decorre dos desenvolvimentos de Taylor usados para as discretizações espacial e temporal, e que a estabilidade é localmente verificada a cada passo de tempo considerado para simulação, a convergência é verificada usando métodos numéricos, isto é, repetindo, sucessivamente, o cálculo para uma série de malhas temporais refinadas. Dessa forma, fixando-se a malha espacial com $\Delta z=1cm$ , variou-se os passos de tempo calculando erro relativo cometido em cada iteração do método da seguinte forma

$\left|\right|\overline{e}|{|}_{L_2}=\sqrt{\frac{{\displaystyle \sum _{i=1}^n({\theta }_i^{j+1}-{\theta }_i^j{)}^2}}{{\displaystyle \sum _{i=1}^n({\theta }_i^{j+1}{)}^2}}}.$ (5.1)

Os pontos $(\Delta t,|\left|\overline{e}\right|{|}_{L_2})$ foram calculados, onde as informações desta análise na escala logarítmica estão ilustradas na Figura 6 . A linha de tendência com coeficiente angular igual a 1 mostra a ordem do erro de truncamento esperada para o esquema de discretização temporal proposto, $O\left(\Delta t\right)=1$ , onde maior valor de Δ t foi de $0,007s\left(\right|\left|\overline{e}\right|{|}_{L_2}=8,40\times {10}^{-7})$ .

Figure 6:
Erro relativo considerando o modelo van Genuchten.

Finalmente, a análise numérica anterior é novamente realizada, considerando, agora, o modelo Haverkamp e os parâmetros obtidos pelo AG, a saber: ${\theta }_s=0,399c{m}^3·c{m}^{-3},{\theta }_r=8,224\times {10}^{-4}c{m}^3·c{m}^{-3},A=1.693c{m}^{-1}\mathrm{e}\beta =0.359$ . O teor de umidade para o Horizonte E obtido nesta simulação pode ser visto na Figura 7 . A Figura 8 mostra a linha de tendência para o erro relativo, onde foi possível manter a ordem $O\left(\Delta t\right)=1$ com uma malha temporal menos refinada em relação à simulação anterior, a partir de $\Delta t=0,02s\left(\right|\left|\overline{e}\right|{|}_{L_2}=1,45\times {10}^{-6})$ .

Figure 7
Teor de umidade para o horizonte E.

Figure 8
Erro relativo considerando o modelo Haverkamp.

6 CONCLUSÕES

O presente artigo teve como objetivo propor uma abordagem não convencional para o ajuste das curvas de retenção da água no solo, considerando amostras de três horizontes de solo, coletadas no entorno da Universidade Federal Rural do Rio Janeiro. A abordagem se mostrou eficaz, visto que os parâmetros de ajustes alcançados pelo AG obtiveram os maiores coeficientes de determinação R ², quando comparados aos resultados obtidos pelo programa SWRC-Fit, que usa o método determinístico Levenberg-Marquardt. Tal fato ocorreu com todos os modelos de curvas de retenção propostos que também estão disponíveis no software . O sucesso dos resultados pode ser explicado pela estratégia de considerar como problema de otimização para o AG a busca dos maximizadores do coeficiente de determinação, que se comportou como a função objetivo do problema. Em relação aos dados analisados neste trabalho, o modelo Haverkamp para a curva de retenção apresentou o maior coeficiente R ² quando comparado aos modelos van Genuchten e Brooks-Corey. Por fim, foram realizadas simulações numéricas da equação de Richards considerando um simples esquema discretização temporal explícito. Além do melhor ajuste, foi usando o modelo Haverkamp para a curvas de retenção da água no solo e de condutividade hidráulica na equação de Richards que uma maior relaxação no tamanho de passo de tempo foi possível, como pode ser observado nas Figuras 6 e 8 . A taxa de convergência calculada mantendo uma malha espacial com $\Delta z=1cm$ e refinando a malha temporal se mostrou consistente com a esperada, considerando o erro de truncamento do método de Euler explícito, a partir de $\Delta t=0,02s\left(\right|\left|\overline{e}\right|{|}_{L_2}=1,45\times {10}^{-6})$ para o modelo Haverkamp e apenas a partir de $\Delta t=0,007s\left(\right|\left|\overline{e}\right|{|}_{L_2}=8,40\times {10}^{-7})$ para o modelo van Genuchten.

Novas simulações numéricas estão sendo realizadas considerando outros exemplos e esquemas de discretização a fim de avaliar a importância da escolha dos modelos de curvas de retenção e de condutividade hidráulica, aliando ajuste de curvas e estabilidade numérica. Novos experimentos também estão sendo finalizados a fim de comparar os valores de umidade mensurados em campo controlado com os estimados pela simulação numérica da equação de Richards, usando as curvas de retenção já obtidas em laboratório. Estes serão temas de futuras publicações.

ACKNOWLEDGEMENTS

Uma versão preliminar deste trabalho foi apresentada no XXXIX CNMAC 2019.

REFERÊNCIAS

• ¹
  D.L. Baker. A Darcian integral approximation to interblock hydraulic conductivity means in vertical infiltration. Computers & Geosciences, 26(5) (2000), 581-590.
• ²
  B. Belfort & F. Lehmann. Comparison of equivalent conductivities for numerical simulation of onedimensional unsaturated flow. Vadose Zone Journal, 4(4) (2005), 1191-1200.
• ³
  R.H. Brooks & A.T. Corey. Properties of porous media affecting fluid flow. Journal of the irrigation and drainage division, 92(2) (1966), 61-90.
• ⁴
  D. Caviedes-Voullième, P. Garcı, J. Murillo et al. Verification, conservation, stability and efficiency of a finite volume method for the 1D Richards equation. Journal of hydrology, 480 (2013), 69-84.
• ⁵
  M.A. Celia, E.T. Bouloutas & R.L. Zarba. A general mass-conservative numerical solution for the unsaturated flow equation. Water resources research, 26(7) (1990), 1483-1496.
• ⁶
  A.E. Eiben & J.H. Smith. “Introduction to Evolutionary Computing”. Springer (2003).
• ⁷
  R. Eymard, M. Gutnic & D. Hilhorst. The finite volume method for Richards equation. Computational Geosciences, 3(3-4) (1999), 259-294.
• ⁸
  P.A. Forsyth, Y. Wu & K. Pruess. Robust numerical methods for saturated-unsaturated flow with dry initial conditions in heterogeneous media. Advances in Water Resources, 18(1) (1995), 25-38.
• ⁹
  J.M. Gastó, J. Grifoll & Y. Cohen. Estimation of internodal permeabilities for numerical simulation of unsaturated flows. Water Resources Research, 38(12) (2002), 62-1.
• ¹⁰
  D.E. Godberg. Genetic algorithms in search. Optimization, and Machine Learning, (1989).
• ¹¹
  R. Haverkamp & M. Vauclin. A note on estimating finite difference interblock hydraulic conductivity values for transient unsaturated flow problems. Water Resources Research, 15(1) (1979), 181-187.
• ¹²
  R. Haverkamp, M. Vauclin, J. Touma, P. Wierenga & G. Vachaud. A comparison of numerical simulation models for one-dimensional infiltration 1. Soil Science Society of America Journal, 41(2) (1977), 285-294.
• ¹³
  U. Hornung & W. Messing. Truncation errors in the numerical solution of horizontal diffusion in saturated/unsaturated media. Advances in Water Resources, 6(3) (1983), 165-168.
• ¹⁴
  A.E. Klar. “A água no sistema solo-planta-atmosfera”, volume 385. Nobel São Paulo (1984).
• ¹⁵
  W. Lai & F.L. Ogden. A mass-conservative finite volume predictor-corrector solution of the 1D Richards’ equation. Journal of Hydrology, 523 (2015), 119-127.
• ¹⁶
  P.L. Libardi. “Dinâmica de Água no Solo”. Edusp (2005).
• ¹⁷
  R. Linden. “Algoritmos genéticos - Uma importante ferramenta da inteligência computacional”. Brasport, Rio de Janeiro (2008).
• ¹⁸
  C.R. Maliska. “Transferência de calor e mecânica dos fluidos computacional.”. Grupo Gen-LTC (2017).
• ¹⁹
  Y. Mualem. A new model for predicting the hydraulic conductivity of unsaturated porous media. Water resources research, 12(3) (1976), 513-522.
• ²⁰
  J. Nocedal & S. Wright. “Numerical optimization”. Springer Science & Business Media (2006).
• ²¹
  K. Reichardt & L.C. Timm. “Solo, planta e atmosfera: conceitos, processos e aplicações”. Manole Barueri (2004).
• ²²
  L. Santos, H. Campos-Velho & L. Chiwiacowsky. Uma l lise de robustez do método híbrido de estimativa c c a o dano estrutural. TEMA (São Carlos), 12(3) (2011), 245-252.
• ²³
  W. Santos, J. Santiago & J. Telles. Optimal positioning of anodes and virtual sources in the design of cathodic protection systems using the method of fundamental solutions. Engineering Analysis with Boundary Elements, 46 (2014), 67-74.
• ²⁴
  W.J. Santos, R.F. Oliveira, M.B. Ceddia & J.L. Almeida. Conteúdo volumétrico da água no solo via modelos de competição interespecífica. Pesquisa e Ensino em Ciências Exatas e da Natureza, 2 (2018), 30-39.
• ²⁵
  M.G. Schaap & M.T. Van Genuchten. A modified Muallen-van Genuchten formulation for improved description of the hydraulic conductivity near saturation. Vadose Sone Journal, 5(1) (2006), 27-34.
• ²⁶
  K. Seki. SWRC fit-a nonlinear fitting program with a water retention curve for soils having unimodal and bimodal pore structure. Hydrology and Earth System Sciences Discussions, 4(1) (2007), 407-437.
• ²⁷
  L.S.B. Silva, A. Vasconcelos, A.L. Sanches, R.T.N. Cardoso, J.L.A. Fernandes & A. Eiras. Otimização Mono-objetivo no Controle do Mosquito Aedes aegypti por meio de um Modelo de Duas Populações com Influência da Precipitação. TEMA (São Carlos), (1) (2019).
• ²⁸
  J. Simunek, M. Van Genuchten & M. Sejna. Development and applications of the HYDRUS and STANMOD software packages and related codes. Vadose Zone Journal, 7(2) (2008), 587-600.
• ²⁹
  J. Takeuchi, T. Kawachi, C. Imagawa, N. Buma, K. Unami & S. Maeda. A physically based FVM watershed model fully coupling surface and subsurface water flows. Paddy and Water Environment, 8(2) (2010), 145-156.
• ³⁰
  J. Van Dam & R. Feddes. Numerical simulation of infiltration, evaporation and shallow groundwater levels with the Richards equation. Journal of Hydrology, 233(1-4) (2000), 72-85.
• ³¹
  M.T. Van Genuchten. A closed-form equation for predicting the hydraulic conductivity of unsaturated soils 1. Soil Science Society of America Journal, 44(5) (1980), 892-898.
• ³²
  J. Vanderborght, R. Kasteel, M. Herbst, M. Javaux, D. Thiéry, M. Vanclooster, C. Mouvet & H. Vereecken. A set of analytical benchmarks to test numerical models of flow and transport in soils. Vadose Zone Journal, 4(1) (2005), 206-221.
• ³³
  T. Vogel & M. Cisrelova. on the reliability of unsaturated hydraulic conductivity calculated from the moisture retention curve. Transport in porous media, 3(1) (1988), 1-15.
• ³⁴
  T. Vogel, M.T. Van Genuchten & M. Cislerova. Effect of the shape of the soil hydraulic functions near saturation on variably-saturated flow predictions. Advances in Water Resources, 24(2) (2000), 133-144.
• ³⁵
  C. Zambra, M. Dumbser, E. Toro & N. Moraga. A novel numerical method of high-order accuracy for flow in unsaturated porous media. International Journal for Numerical Methods in Engineering, 89(2) (2012), 227-240.

Datas de Publicação

• Publicação nesta coleção
  05 Abr 2021
• Data do Fascículo
  Jan-Mar 2021

Histórico

• Recebido
  27 Jan 2020
• Aceito
  01 Set 2020