Geoprocessamento para a Solução Fraca do Transporte de Contaminantes Acoplado ao Fluxo de Água Subterrânea ^†

RESUMO

A solução automatizada das equações do transporte de contaminantes no interior de um aquífero depende de parâmetros físicos e químico e, principalmente, do campo de velocidades de suas águas subterrâneas. Este campo de velocidades pode ser obtido da solução numérica do modelo matemático que descreve o potencial hidráulico da região de fluxo saturado. Na fase de pré-processamento da solução, o uso de técnicas de geoprocessamento foi aplicado na geração de malhas computacionais com a imposição geoespacial das condições iniciais e de fronteira para simular o comportamento do lençol freático na localidade do contaminante. Neste trabalho, após a transferência do potencial hidráulico, em termos de campo de velocidades para a equação do transporte, um código numérico desenvolvido em Python foi capaz de prever a migração de um contaminante no interior de um aquífero descrevendo a quantidade do poluente através das isolinhas de concentração. Outra contribuição importante das técnicas de geoprocessamento foi observado na padronização e especificação das informações georreferenciadas das várias fontes que forneceram as bases de dados geoespacias para a descrição hidrogeológica da região de estudos.

Palavras-chave:
geocomputação; nível freático; aquíferos; dolfin; modelo ADR

ABSTRACT

The automated solution of contaminant transport equations inside an aquifer depends on physical and chemical parameters and, most importantly, on the groundwater velocity field. Such a field can be obtained from the numerical solution of a mathematical model describing the hydraulic potential of the saturated flux region. In the pre-processing solution step, geoprocessing techniques were applied to produce computational meshes respecting the geospatial constraints of the initial and boundary conditions to simulate the behavior of the water table around the contaminant. In this work, after the hydraulic potential transference - in terms of velocities field for the transport equation -, a Python numerical code was developed to predict the contaminant migration within an aquifer, describing the amount of pollutant through concentration isolines. Another important contribution of the geoprocessing techniques was observed in the standardization and specification of the georeferenced information from several sources, which provided the geospatial database for the hydrogeologic description of the studied area.

Keywords:
geocomputation; water table; aquifers; dolfin code; ADR model

1 INTRODUÇÃO

O campo de velocidades das águas subterrâneas é um parâmetro importante na descrição apropriada do componente advectivo da equação do transporte de contaminantes no interior de um aquífero. A direção e a magnitude do campo de velocidades também são aplicadas na deter minação dos coeficientes da matriz de dispersibilidade para uma descrição da contribuição do transporte dispersivo.

Este comportamento hidráulico, que rege o movimento da água subterrânea de uma região de interesse, fornece o campo de velocidades pela Lei de Darcy e poderá ser obtido através da resolução automatizada das equações diferenciais do modelo matemático que descrevem as linhas de fluxo em solo saturado. A delimitação do domínio computacional para representar o contorno dos limites da extensão do aquífero, a imposição das condições de fronteira e da condição inicial do problema de fluxo e a representação georreferenciada de uma fonte pontual de contaminação de água subterrânea são informações úteis disponibilizadas pelas técnicas de geoprocessamento e de suas ferramentas avançadas.

É importante ressaltar neste momento que, de acordo com⁶, existe uma preocupação constante na comunidade científica que buscam contribuições inovadoras para tornarem os Sistemas de Informações Geográficas (SIG) uma ferramenta útil para auxiliar no planejamento, modelagem e sínteses dos recursos hídricos através da geocomputação. Outra justificativa a ser mencionada é que as específicas informações relacionadas aos recursos de águas subterrâneas são escassas e, frequentemente, tornam-se inacessíveis ou possuem suas aquisições de dados através de elevados custos de exploração⁴. Neste sentido, a aplicação das técnicas de geoprocessamento torna-se fundamental para a exibição e análise dos dados da representação espacial que, mesmo quando coletados de diferentes fontes, são capazes de compor uma única geodatabase da região de interesse.

Neste artigo, o uso das técnicas de geoprocessamento das aplicações hidrogeológicas com as soluções automatizadas dos modelos de simulação visam, além da integração de uma grande quantidade de informações georreferenciadas provenientes das mais variadas bases de dados, a compreensão e o manejo dos recursos hídricos para prever e interpretar o fluxo de água e o transporte de contaminantes no interior de um aquífero. Assim, o objetivo deste trabalho é a implementação conjunta de uma geodatabase de estrutura SIG acoplada a um código de programação Python para obter a solução automatizada, via formulação fraca e método dos elementos finitos, desenvolvidas com o uso de bibliotecas Dolfin do Projeto FEniCS¹¹. A contribuição relevante do trabalho pode ser entendida como a apresentação de uma metodologia geocomputacional que obtém o campo de velocidades de água subterrânea de uma região geor referenciada e, no mesmo código numérico, implementa os resultados numéricos nos termos advectivos e dispersivos visando a resolução automatizada do transporte de contaminantes em aquíferos.

Para isto, o artigo organiza-se da seguinte forma: a Seção 2 descreve as equações governantes da equação do fluxo em um aquífero confinado e as equações do transporte de contaminantes em meio saturado, visando a automação da resolução destas equações através de solução fraca pelo método dos elementos finitos; a Seção 3 fornece a implementação Pyhton da respectiva formulação fraca das equações; a Seção 4, foi subdividida na subseção 4.1 que apresenta as ferramentas de geoprocessamento usadas no artigo para aplicação da metodologia proposta, e na subseção 4.2 que identifica na cidade de Bauru-SP uma fonte de contaminação pontual do lençol freático. Outros exemplos de SIG que são aplicados na resolução de problemas envolvendo águas subterrâneas são encontrados no pacote ArcHydro Groundwater do software ArcGis¹⁵. No entanto, estes exemplos não apresentam vínculo entre as soluções numéricas das equações envolvidas.

2 EQUAÇÕES GOVERNANTES

O modelo matemático para a distribuição de cargas hidráulicas que regem o fluxo de água subterrânea, obtido pela lei de Darcy e pelo princípio da conservação de massa em um volume elementar representativo (REV) de um aquífero é dado por⁵ como sendo a equação linear elíptica:

$\frac{\partial }{\partial x}\left[K_{xx}\frac{\partial h}{\partial x}\right]+\frac{\partial }{\partial y}\left[K_{yy}\frac{\partial h}{\partial y}\right]+\frac{\partial }{\partial z}\left[K_{zz}\frac{\partial h}{\partial z}\right]+W\left(x,y,z,t\right)=S_e\frac{\partial h}{\partial t}$ (2.1)

sendo que h = h(x, y, z, t) é a carga hidráulica total [L], _^Kxx , _^Kyy e _^Kzz são os componentes principais do tensor de condutividade hidráulica $\left[\frac{L}{T}\right]$, _^Ss é o coeficiente de armazenamento específico $\left[\frac{1}{L}\right]$ e W o termo de fonte ou sorvedouro no interior do aquífero $\left[\frac{1}{L}\right]$. Para o caso tridimensional W poderia ser representado por:

$W={\displaystyle \sum }_i{Q}_i\delta \left(x-x_i\right)\delta \left(y-y_i\right)\delta \left(z-z_i\right)$ (2.2)

sendo $Q_i\left[\frac{M^3}{T}\right]$ a taxa de bombeamento ou injeção do poço i em (_{^{xi, yi, zi}} ) e δ é a função delta de Dirac com unidades $\frac{1}{L}$. Considerando desprezíveis as variações de carga ao longo da dimensão vertical (hipótese de Dupuit), e que a dimensão horizontal dos aquíferos em escalas regionais pode ser da ordem de dezenas de quilômetros, o fluxo pode ser modelado por uma equação bidimensional em x e y dada por⁵:

$\frac{\partial }{\partial x}\left[T_{xx}\frac{\partial h}{\partial x}\right]+\frac{\partial }{\partial y}\left[T_{yy}\frac{\partial h}{\partial y}\right]+W\left(x,y,t\right)=S\frac{\partial h}{\partial t}$ (2.3)

sendo

$T_{xx}=b·K_{xx}e{T}_{yy}=b·K_{yy}\left[\frac{L^2}{T}\right]$

as transmissividades nas direções x e y do aquífero confinado, b a espessura do aquífero [L] e S = _{^{Ss · b}} , o coeficiente de armazenamento [adimensional]. A função W representa os termos de drenança dos fluxos verticais remanescentes em z = 0 (camada confinante inferior) e em z = b (camada confinante superior) acrescidos da atividade de um poço i na posição (_{^{xi, yi}} ). Assim, a formulação contínua da equação transiente do fluxo em aquífero confinado, acrescidas das suas condições iniciais e de fronteira, é dada por:

$$\{\begin{array}{l}\frac{\partial }{\partial x}\left[T_{xx}\frac{\partial h}{\partial x}\right]+\frac{\partial }{\partial y}\left[T_{yy}\frac{\partial h}{\partial y}\right]+W\left(x,y,t\right)=S\frac{\partial h}{\partial t}em\Omega \times \left(\mathrm{0,}{t}_f\right]\hfill \\ \begin{array}{l}\begin{array}{cccccccc}& & & & & & & \end{array}\begin{array}{cc}& \end{array}h=h_{D}sobre{\Gamma }_D\times \left(\mathrm{0,}{t}_f\right]\\ \begin{array}{cccccccc}& & & & & & & \end{array}n·D\nabla h=gsobre{\Gamma }_N\times \left(\mathrm{0,}{t}_f\right]\\ \begin{array}{cccccccc}& & & & & & & \end{array}h\left(·,t_0\right)=h_{0}em\Omega parat=0\end{array}\hfill \end{array}\text{}$$ (2.4)

sendo Ω ⊂ ℝ² o domínio poligonal limitado e com fronteira Lipschitz¹ Γ que consiste de duas partes disjuntas, Γ_D a fronteira de Dirichlet e Γ_N a fronteira de Neumann, tais que Γ_D ∪ Γ_N = ∂Ω. O tempo final T é arbitrário, no entanto, precisa ser especificado.

O modelo de simulação do transporte de contaminantes dissolvidos em água subterrânea será descrito por equações de transporte de soluto que descrevem os fenômenos de advecção-difusão (AD) em solo saturado. A equação diferencial parcial AD deriva do princípio da conservação de massa - a equação da continuidade, afirmando que a taxa de variação da massa de soluto dentro de um volume de controle do meio poroso é igual a diferença entre o fluxo de soluto que entra e o fluxo que sai do volume, ajustados por perdas ou ganhos da massa de soluto devido a reações químicas ou degradação biológicas. Assim, a taxa de soluto entrando no volume de controle será caracterizado por dois processos principais: a advecção e a dispersão hidrodinâmica². E ainda, conforme modelado na equação (2.5), esse segundo processo físico será representado pelo efeito combinado de outros dois processos: a difusão molecular _^Dm e a dispersão mecânica _^Dd .

A equação governante do transporte de contaminantes em águas subterrâneas será modelada pela equação parabólica linear e de segunda ordem (2.5) com as seguintes condições iniciais e de fronteira, conforme fornecidas em¹⁴:

$\{\begin{array}{cc}{\partial }_{t}C-div\left(D\nabla C\right)+v·\nabla C+\lambda C=f& \text{in}\Omega \times \left(\mathrm{0,}{t}_f\right]\\ C=C_D& \text{em}{\Gamma }_D\times \left(\mathrm{0,}{t}_f\right]\\ n·D\nabla C=G& \text{sobre}{\Gamma }_N\times \left(\mathrm{0,}{t}_f\right]\\ C\left(·,t_0\right)=C_0& \text{em}\Omega \text{para}t=0\end{array}$ (2.5)

sendo C a concentração de soluto $\left[\frac{M}{L^3}\right]$, ?C o gradiente de concentração $\left[\frac{M}{L^4}\right]$, _{^{D = Dm + Dd}} o tensor da dispersão hidrodinâmica $\left[\frac{L^2}{T}\right]$, v vetor de velocidade real $\left[\frac{L}{T}\right]$, λ o coeficiente de decaimento de 1^a ordem $\left[\frac{1}{T}\right]$, f o termo de fonte ou sorvedouro $\left[\frac{M}{L^{2}T}\right]$ e o domínio poligonal limitado Ω ⊂ ℝ² com fronteira Γ = Γ_D ∪ Γ_N Lipschitz.

Para que as informações sobre o campo de velocidades provenientes da equação do fluxo (2.4) sejam aplicadas na equação do transporte de contaminantes (2.5), algumas hipóteses físicas e do meio precisam ser estipuladas. Ou seja, será considerado que o contaminante não influenciará na velocidade do fluxo da água subterrânea, caracterizando o aquífero em estudo como um campo conservativo; a temperatura média no meio poroso não sofrerá alterações significantes durante o tempo observado; a matriz do solo saturado em estudo é rígida, ou seja, não sofrerá deformações ou relaxamento e seus poros estarão todos conectados para a ocorrência do transporte de contaminante.

3 SOLUÇÃO FRACA E IMPLEMENTAÇÃO DAS EQUAÇÕES

A resolução numérica automatizada neste trabalho considera a formulação fraca da equação do fluxo (2.1) dada por¹⁴. As linhas de comando Python que implementam soluções variacionais, através de programação computacional utilizando a biblioteca Dolfin do projeto FEniCS¹¹, são:

u = TrialFunction(V) v = TestFunction(V) f2 = Constant(0.0)

a = S_e*b*u*v*dx + theta*dt*K*b*inner(nabla_grad(u), nabla_grad(v))*dx

L = (S_e*b*u_1*v + dt*f2*v -(1.0-theta)*dt*T*inner(nabla_grad(u_1),

nabla_grad(v)))*dx

A = assemble(a) b = None # variable used for memory savings in

assemble calls

u = Function(V), t = dt # unknown at the current time level

while t <= t_stop:

b = assemble(L, tensor=b)

for bc in bcs:

bc.apply(A,b)

solve(A, u.vector(), b), t += dt

u_1.assign(u)

O campo de velocidades v(x, y) para implementação do termo advectivo e dos coeficientes dispersivos da equação do transporte são obtidos através do cálculo dos gradientes de potenciais hidráulicos e da Lei de Darcy aplicando, de acordo com¹, as seguintes linhas de comando Python:

degree = 1

Vg = VectorFunctionSpace(mesh, "Lagrange", degree)

w = TrialFunction(Vg)

v1 = TestFunction(Vg)

a = inner(w, v1)*dx

L = inner(grad(h), v1)*dx

gradh = Function(Vg)

solve(a == L, gradh)

parameters['allow_extrapolation'] = True

plot(-K*gradh,interactive=True,title='Real Velocity of Groundwater')

gradx, grady = gradh.split(deepcopy=True)

# CLASSES TO DEFINE VELOCITY FIELD IN COORDINATE DIRECTIONS AND USE

IT ON DISPERSION TENSOR

class vx(Expression):

def eval(self, values, x):

values[0]=(-1)*(kxx/n_e)*gradx(x[0],x[1])

vx = vx(element=V.ufl_element());

vx = interpolate(vx,V)

class vy(Expression):

def eval(self, values, x):

values[0]=(-1)*(kyy/n_e)*grady(x[0],x[1])

vy = vy(element=V.ufl_element());

vy = interpolate(vy,V)

Do processo de discretização θA-esquema estável ¹⁶, são fornecidos:

$\begin{array}{l}a\left(C^n,w\right)=\underset{\Omega }{{\displaystyle \int }}\frac{1}{{\tau }_n}{C}^n{w}_{n}d\Omega +\underset{\Omega }{{\displaystyle \int }}\left(\theta \nabla C^n\right)·D_h^{n,\theta }\nabla w_{n}d\Omega \text{}\\ \begin{array}{cccc}& & & \end{array}+\underset{\Omega }{{\displaystyle \int }}{v}_h^{n,\theta }·\nabla \left(\theta C^n\right)w_{n}d\Omega +\underset{\Omega }{{\displaystyle \int }}{\lambda }_h^{n,\theta }\left(\theta C^n\right)w_{n}d\Omega \text{}\\ \begin{array}{cccc}& & & \end{array}-\underset{{\Gamma }_N}{{\displaystyle \int }}\theta n·D_h^n\nabla C^n{w}_{n}dS\end{array}$ (3.1)

$$\begin{gathered} L\left(w_n\right){\int }_{\Omega }\frac{1}{{\tau }_n}{C}^{n-1}{w}_{n}d\Omega +{\int }_{\Omega }\left(\left(\theta -1\right)\nabla C^{n-1}\right)·D_h^{n,\theta }\nabla w_{n}d\Omega \\ +{\int }_{\Omega }{v}_h^{n,\theta }·\nabla \left(\left(\theta -1\right)C^{n-1}\right)w_{n}d\Omega \\ +{\int }_{\Omega }{\lambda }_h^{n,\theta }\left(\left(\theta -1\right)C^{n-1}\right)w_{n}d\Omega \\ +{\int }_{\Omega }\left(\theta f_h^n+\left(1-\theta \right)f_h^{n-1}\right)w_{n}d\Omega \\ +{\int }_{{\Gamma }_N}\left(1-\theta \right)\mathbf{n}·D_h^{n-1}\nabla C^{n-1}{w}_{n}dS. \end{gathered}$$ (3.2)

Esta formulação fraca representa a parte central da solução fraca da equação do modelo ADR do transporte. Seus respectivos termos bilineares (3.1) e lineares (3.2), foram desenvolvidos em ¹⁴. Os sobrescritos n e n ? 1 representam o n-ésimo (ou (n ? 1)-ésimo) passo de tempo na discretização temporal e o subescrito _^Tn representa, no passo de tempo n, a triangulação da discretização espacial do θA-esquema estável¹⁶. A implementação Python dos termos (3.1) e (3.2) na solução numérica que utiliza o método de elementos finitos é descrita pelo seguinte código computacional:

# the bilinear form for finite element formulation

a = (1.0/dt)*c*w*dx + theta*inner(D*nabla_grad(c),nabla_grad(w))*dx+\

theta*(inner(grad(H1), nabla_grad(c)))*w*dx + theta*Lambda*c*w*dx-\

theta*inner(D*grad(c),n)*w*ds(0)-theta*inner(D*grad(c),n)*w*ds(1)

# the linear form for finite element formulation

L = ((1.0/dt)*c1*w+(theta-1.0)*inner(D*nabla_grad(c1),

nabla_grad(w)))*dx+\

(theta-1.0)*inner(grad(H1), grad(c1))*w*dx+(theta-1.0)*Lambda*c1*w*dx

+(theta*f1+(1.0-theta)*f)*w*dx+(1.0-theta)*inner(D*grad(c1),n)*w*ds(0)

+(1.0-theta)*inner(D*grad(c1),n)*w*ds(1)

A escolha da linguagem de programação open source Python para a resolução numérica da equação do transporte (2.5), através da formulação variacional definida pela forma bilinear (3.1) e forma linear (3.2), foi em razão da disponibilidade de uma extensa biblioteca de computação científica que empregam o método dos elementos finitos^{9), (10}. Desta forma, pode-se citar outra contribuição do artigo como sendo a implementação Python da formulação fraca das equações governantes (2.4) e (2.5) com condições inicias e de contorno automatizadas pelo emprego de ferramentas de geoprocessamento.

4 GEOPROCESSAMENTO DO DOMÍNIO COMPUTACIONAL

Ferramentas de geoprocessamento foram aplicadas inicialmente na adequação de metadadosreferentes à padronização em seus georreferenciamentos para um único sistema de coordenadas projetadas-UTM (EPSG 31988-SIRGAS 2000/UTM zone 22S), pois, muitas vezes são obtidos de diferentes websites descritos através de coordenadas geográficas. Após a identificação da região de interesse, realiza-se a união de polylines que descrevem a fronteira do domínio computacional. Desta fronteira, um shapefile de pontos é criado contendo a orientação dos nós representados pelos vértices dos polylines para importação e geração das malhas Dolfin. Após identificação da região de estudos e importação do correspondente Modelo de Elevação Digital (DEM) uma interpolação IDW³, definida sobre uma máscara da região interna, conduzem o preenchimento de tabelas dinâmicas. Através de informações de localização espacial e de valores dos níveis freáticos dos vários poços tubulares, o passo seguinte passa a ser a definição de condições iniciais e de fronteira para a equação do fluxo.

No geoprocessamento da região de interesse também são implementados formatos digitais de corpos d’águas superficiais, tais como pequenas lagoas ou rios que cruzam o domínio, para importar os valores de suas cotas hidráulicas e representar a interação com aquífero em estudo¹³.

4.1 Região de Estudos e o Modelo de Fluxo

Para simulação numérica do transporte de contaminantes numa situação realística, o domínio computacional escolhido representará o aquífero abaixo da região urbana da cidade de Bauru-SP (Fig. 1), sujeito à contaminação pontual a partir das coordenadas geográficas 49,008W e 22,327S, conforme informações extraídas do Relatório Anual de Áreas Contaminadas e Reabilitadas no Estado de São Paulo.

Figura 1
Região urbana de Bauru-SP sobre divisas da microbacia e geopolítica, e ampliação da vizinhança populacional da fonte poluente de coordenadas 49,008W e 22,327S.

Esta fonte de contaminação por infiltração, nos dias atuais, encontra-se em situação de remediação. No entanto, antes da interdição, houve um impacto no meio ambiente com agentes poluentes atingindo o solo superficial, o subsolo, as águas subterrâneas, os sedimentos e a biota³.

Nesta fase, a metodologia foi subdividida em duas partes principais estabelecendo-se, inicialmente, uma geodatabase da região de estudo com shapefiles que representam o município de Bauru-SP (polígono), o domínio limitado pela microbacia a oeste e pela divisa geopolítica atravessada por um canal urbano a leste (linhas), a disposição de poços tubulares e dos vértices da fronteira do domínio (pontos), e um modelo de elevação digital (DEM). Posteriormente, umamalha de elementos finitos com condições iniciais e de fronteira para o problema do fluxo importa valores de coordenadas projetadas e de níveis estáticos de tabelas dinâmicas da geodatabase, conforme representado na Figura 2.

Figura 2
Geodatabase da região de estudo, extração de um shapefile representando o domínio computacional e a correspondente malha de elementos finitos no formato Dolfin para simulação no código Python.

A resolução automatizada do modelo matemático do fluxo e do transporte foi obtida com a programação Python usando bibliotecas Dolfin do projeto FEniCS. As fases do geoprocessamento aplicadas na identificação da fronteira da região de interesse e a obtenção da correspondente malha Dolfin (.xml) seguem apresentadas na Figura 2. Ainda na execução desta segunda etapa, um shapefile com características Polyline precisou ser criado através da aplicação da ferramenta clip para extrair o segmento fronteiriço da microbacia sobre a porção homogênea que representa o município de Bauru-SP na Figura 2 e de outra ferramenta do ArcGis para a reunião deste com o segmento da divisa geopolítica do município fora da mesma microbacia.

O polyline resultante, cujos vértices representam nós orientados e de coordenadas projetadas EPSG:31982 (SIRGAS 2000/UTM), alimentou o script Python para construir a malha de elementos finitos no formato Dolfin com a restrição de 10 m² para o valor máximo da área de cada elemento triangular.

O polígono interior do domínio que circunda a fonte poluente é apenas figurativo e não faz parte da fronteira. Na Figura 2 observa-se também a importância obtida com a padronização do sistema de referência em coordenadas projetadas UTM das diversas fontes de dados digitais para polígonos, polylines, nós e modelo de elevação digital.

Se considerar que a cidade de Bauru-SP pertence ao Sistema de Aquífero Guarani (SAG), então o aquífero sob a micro bacia hidrográfica da região de estudos possui características de aquífero livre¹². Assim, o modelo matemático do fluxo de águas subterrâneas seria regido por uma equação não linear. No entanto, em virtude da escala regional do aquífero e da pequena variação na altura piezométrica com relação à sua magnitude, considerou-se um processo de linearização nas equações do aquífero livre para aplicar o código numérico Python que simula o fluxo em aquíferos confinados⁷, ou seja, o modelo matemático regido pelas equações lineares (2.4).

Em relação as tabelas dinâmicas da geodatabase, as informações georreferenciadas para locali zação dos poços tubulares e as respectivas medições dos níveis estáticos (_^Ne ), disponibilizados pela SIAGAS⁴ para o interior da região de estudo, adotou-se como estratégia a implantação no código Python um arquivo .csv contendo as condições inicias para resolução numérica do problema do fluxo subterrâneo. Informações de características hidrogeológicas da região de estudo foram baseadas em².

Para implementar as condições de fronteira das equações do fluxo (2.4), uma ferramenta de interpolação IDW do ArcGis, que aplica em nós da malha o método do inverso do quadrado da distância, foi responsável para estender as estimativas dos valores que delimitam o contorno oeste do domínio computacional. Para os nós da linha que descrevem o trecho do canal na região urbana de Bauru-SP, aplicou-se os mesmos valores correspondentes às suas cotas topográficas e no contorno geopolítico leste não foi atribuído nenhuma condição de fronteira. Após a imposição destas condições iniciais e de contorno sobre os nós da malha de elementos finitos no formato Dolfin, uma simulação do comportamento do lençol freático e o correspondente fluxo de velocidade das águas subterrâneas na proximidade da fonte poluidora pode ser visualizada na Figura 3 para um determinado passo de tempo t.

Figura 3
Geoprocessamento aplicado na delimitação do domínio computacional, na imposição das condições iniciais e de fronteira, na simulação do comportamento do lençol freático e na determinação de linha de fluxo das águas subterrâneas.

De posse do campo de velocidades v = v(x, y, t), previamente calculado conforme a descrição da Seção 3, avaliou-se tanto a matriz de dispersão como dependente da posição (x, y), e de suas componentes _^vx e _^vy , quanto o coeficiente variável para o termo de advecção da equação do transporte (2.5).

4.2 Fluxo acoplado ao Transporte de Contaminantes

Para o problema do transporte, considerou-se que a fonte em estudo descartou uma certa quantidade $C_0\left[\frac{mg}{litro}\right]$ de um poluente caracterizado como sendo um fluido de baixa densidade e miscível. Dessa forma, após percolação no subsolo, o contaminante após atingir as águas subterrâneas, se propagou pelo lençol freático.

A condição inicial para resolução numérica da equação do transporte foi imposta com o emprego de uma função gaussiana de suporte compacto para descrever a distribuição normalizada $\frac{C}{C_0}$ do poluente em uma distância radial r = 0,25 km da fonte poluidora. Para condições de contorno, considerou-se na fronteira oeste as condições de Dirichlet com C = 0 e na fronteira leste as condições de Neumann¹⁴.

Uma vez que a fonte poluidora se encontra, até o momento, em processo de descontaminação, a condição inicial dada para a simulação numérica foi do tipo pulso de contaminação a partir de um tempo t ₀, ou seja, considerou-se que a fonte poluidora não mais retroalimentava a região de estudos com outras quantidades de poluente para t > t ₀. As isolíneas de concentração de duas soluções gráficas apresentadas também são visualizadas na Figura 4 indicando que a frente de contaminação, com sua forma singular, seguiu em direção norte no rumo do canal da cidade e depois, em função de mudanças na direção do campo de velocidade na região, a frente adotou outro caminho preferencial dirigindo-se no sentido leste da cidade de Bauru-SP.

Figura 4
Solução numérica Python para o problema do transporte de contaminantes em aquíferos acoplada ao campo de velocidades das águas subterrâneas obtidas com o uso das técnicas de geoprocessamento do ArcGis.

A Figura 4 ilustra o comportamento do lençol freático, a região inicial de contaminação difusa na proximidade da fonte poluidora e duas soluções gráficas da frente de contaminação, uma em um tempo intermediário e a outra no tempo final da simulação. Verifica-se na simulação que esta frente de contaminação acompanhou o campo de velocidades obtido na segunda parte da implementação do código Python para a equação do fluxo subterrâneo. Para a simulação computacional do transporte considerou-se um período de 120 anos. Assim, o contaminante miscível percorreu sobre o lençol freático uma extensão total de 4,2 km, apresentando no final da trajetória uma persistência de baixos porcentuais do poluente em relação ao valor inicial normalizado, conforme valores indicados em suas isolíneas de concentração.

A validação desses resultados numéricos do fluxo de água subterrânea e do transporte de contaminantes é assegurada por estimadores de erro desenvolvidos em⁸ para resolução numérica das respectivas equações (2.4) e (2.5) e pelos métodos numéricos de resolução das equações algébricas de grande porte propostos por¹⁴.

5 CONSIDERAÇÕES FINAIS

Ferramentas de georreferenciamento para a área da hidrogeologia são utilizadas na descrição de seus sistemas hídricos naturais e suas conectividades entre os principais corpos d’água superficiais de uma bacia hidrográfica (nascentes, rios e lagoas), os poços tubulares no interior de um domínio e o aquífero simulado. Desta interação, a representação numérica que descreve a direção das linhas de fluxo de água subterrânea estabelece uma região de influência na simulação geocomputacional do transporte de contaminantes em aquíferos.

A simulação numérica apresentada neste trabalho considerou simplificações com relação à estra tificação geológica (morfologia, consistência e porosidade) da região de Bauru-SP. No entanto, estas descrições e propriedade hidráulicas são facilmente implementadas no código numérico. Os valores dos níveis estáticos obtidas dos poços tubulares apenas serviram de referências para estimativas no comportamento do lençol freático na região de estudos. Neste sentido, não foram considerados poços com bombeamento, ou seja, não foi retirado do solo saturado nenhuma quantidade de contaminante que estivesse na região de influência das conhecidas curvas de rebaixamento. Tal situação, uma vez que as informações estejam disponíveis, também possui fácil implementação no termo de fonte/sumidouro da equação do transporte. A versatilidade das ferramentas de geoprocessamento de uma região de interesse, acompanhadas de suas condições iniciais e de fronteiras, contribuíram para a modelagem computacional no sentido de padronizar seus dados geoespaciais e permitir o acoplamento da solução numérica, pelo método de elementos finitos, do fluxo de água subterrânea na solução do problema de transporte de contaminantes em meio saturado. Essas soluções gráficas sobre domínios georreferenciados possibilitam a criação de mapas digitais para representação visual do campo de velocidades de águas subterrâneas e das isolinhas de concentração com sobreposição apropriada do domínio geográfico previamente estabelecido.

AGRADECIMENTOS

Os autores agradecem o grupo de pesquisa WARSA (Water Resource System Analysis) pelo acolhimento científico e acadêmico que resultaram nos avanços da implementação.

REFERÊNCIAS

• ¹
  J. Bear & A.H. Cheng. Modeling Groundwater Flow and Contaminant Transport, Haifa, Springer, (2010).
• ²
  N. Cavaguti & F.P. Silva. Gestão dos recursos hídricos subterrâneos na cidade de Bauru-SP, face as características hidrogeológicas especiais da região, em “Congresso Brasileiro de Águas Subterrâneas”, Belo Horizonte-MG. Anais. ABAS, p. 74-79. (1992).
• ³
  CETESB Relatório Anual de Áreas Contaminadas e Reabilitadas no Estado de São Paulo, 3 CETESB Relatório Anual de Áreas Contaminadas e Reabilitadas no Estado de São Paulo, http://areascontaminadas.cetesb.sp.gov.br/ , acesso em 17/03/2015.
  » http://areascontaminadas.cetesb.sp.gov.br/
• ⁴
  R. Chesnaux et al. Building a geodatabase for mapping hydrogeological features and 3D modeling of groundwater systems: Application to the Saguenay-Lac-St-Jean Region. Canada, Computers & Geosciences, (2011).
• ⁵
  R. Cleary. Águas Subterrâneas, Clean Environment Brasil e Princenton Groundwater, (2008).
• ⁶
  B. Dixon, V. Uddameri & C. Ray. GIS and Geocomputation for Water Resource Science and Engineering, Wiley Publishers, (2015).
• ⁷
  J. Donea & A. Huerta. Finite Element Methods for Flow Problems, John Wiley & Sons, (2004).
• ⁸
  A. Firmiano & E. Wendland. Solução numérica adaptativa da equação do transporte de contaminantes. Biblioteca24horas, São Paulo-SP, (2012).
• ⁹
  M. Gockenbach. Understanding and Implementing the Finite Element Method. SIAM, (2006).
• ¹⁰
  H.P. Langtangen. A Primer on Scientific Programming with Python. Texts in Computational Science and Engineering, vol. 6. Springer, (2011).
• ¹¹
  A. Logg, K.A. Mardal & G. Wells. FEniCS Project: Lecture Notes in Computational Science and Engineering, Springer. New York, (2010). DOI 10.1007/978.3.642.23099.8.
  » https://doi.org/10.1007/978.3.642.23099.8
• ¹²
  A. Machado, E. Wendland & P. Krause. Hydrologic Simulation for Water Balance Improvement in an Outcrop Area of the Guarani Aquifer System. Environmental Processes, 3 (2016), 1-20.
• ¹³
  M. Morio, M. Finkel & E. Martac. Flow guided interpolation - A GIS-based method to represent contaminant concentration distributions in groundwater. Environmental Modelling & Software, 25(12) (2010), 1769-1780.
• ¹⁴
  J.P.M. Santos. Método Multigrid Algébrico: A Reutilização das Estruturas Multigrid no Transporte de Contaminantes. Tese de Doutorado, EESC, USP, São Carlos, SP, (2015).
• ¹⁵
  G. Strassberg et al. Arc Hydro Groundwater: GIS for Hydrogeology, ESRI Incorporated, (2011).
• ¹⁶
  R. Verfürth. Adaptive Finite Element Methods Lecture Notes Winter Term, Fakultät für Mathematik, Ruhr-Universität Bochum, Deutschland, (2008).

• †
  Trabalho apresentado no XXXVI CNMAC (Congresso Nacional de Matem´atica Aplicada e Computacional).

• 1
  O termo é devido ao matemático alemão Rudolf Lipschitz (1832-1903). Um domínio Lipschitz, ou domínio com fronteira Lipschitz, é um domínio no espaço Euclidiano cuja fronteira é suficientemente regular no sentido que ela pode ser localmente representada como o gráfico de uma função f : X → Y Lipschitz contínua (ie, existe K ≥ 0 tal que _^dy (f(x ₁), f(x ₂)) ≤ _^Kdx (x ₁, x ₂) para todo x ₁ e x ₂ em X).

• 2
  É a dispersão de uma solução que flui através de um meio poroso devido à difusão molecular e à não homogeneidade de velocidades microscópicas. Embora os respectivos parâmetros da dispersão mecânica e da difusão molecular, Dm e Dd, sejam provenientes de diferentes processos físicos, eles podem ser adicionados uma vez que possuem o mesmo efeito físico.

• 3
  Acrônimo do termo Inverse Distance Weighting.

• 4
  Sistema de Informações de Águas Subterrâneas desenvolvido pelo Serviço Geológico do Brasil – SGB, serviço disponível em http://siagasweb.cprm.gov.br/layout/visualizar_mapa.php

Datas de Publicação

• Publicação nesta coleção
  May-Aug 2017

Histórico

• Recebido
  18 Nov 2016
• Aceito
  30 Mar 2017