A análise do escoamento das águas subterrâneas exige uma compreensão de elementos da Mecânica dos fluidos. Albertson & Simons (1964) fornecem uma breve e útil revisão desses elementos; Streeter (1962) e Yennard (1961) são textos clássicos. Nosso objetivo aqui é simplesmente introduzir as propriedades básicas dos fluidos: densidade, peso específico, compressibilidade e viscosidade; e examinar os conceitos de pressão do fluido e carga de pressão.

Um exame dos princípios da mecânica dos fluidos deve começar, em geral, com uma revisão da mecânica dos materiais. A Tabela A1.1 fornece uma lista das propriedades mecânicas da matéria, juntamente com suas unidades dimensionais no sistema métrico SI, que possui dimensões básicas como: massa, comprimento e tempo; e unidades básicas no SI: quilograma (kg), metro (m) e segundo (s). Todas as outras propriedades são medidas em unidades derivadas deste conjunto básico. Algumas dessas propriedades são tão amplamente encontradas, e suas dimensões básicas tão complexas que o SI cunhou nomes especiais para suas unidades derivadas. Conforme apontado na Tabela A1.1, força e peso são medidos em Newton (N), pressão e tensão em N/m² ou Pascal (Pa), e trabalho e energia em Joule (J).

A Tabela A1.2 fornece uma análise do sistema internacional (SI) de algumas propriedades de fluidos e termos de águas subterrâneas que ocorrem neste texto. Cada uma delas é descrita mais detalhadamente no Capítulo 2.

Grande parte da tecnologia associada ao desenvolvimento de recursos hídricos subterrâneos na América do Norte ainda utiliza o sistema FPS (pé-libra-segundo).

A Tabela A1.3 fornece um conjunto de fatores para converter unidades FPS em unidades SI.

A massa específica (ou simplesmente, densidade absoluta) ρ de um fluido é definida como sua massa por unidade de volume (Tabela A1.1). A densidade de peso (ou peso específico, ou peso unitário) γ de um fluido é definido como o peso por unidade de volume. Os dois parâmetros estão relacionados por:

(A1.1)

Para a água, ρ = 1,0 g/cm³ = 1.000 kg/m³; γ = 9,8 × 10³ N/m³. No sistema FPS, γ = 62,4 lb_f/ft³.

A densidade relativa G de qualquer material é a relação entre sua densidade (ou peso específico) e da água. Para a água, G = 1,0; para a maioria dos solos e rochas, G ≈ 2,65.

A viscosidade de um fluido é a propriedade que permite que os fluidos resistam ao movimento relativo e à deformação por cisalhamento durante o fluxo. Quanto mais viscoso o fluido, maior a tensão de cisalhamento em qualquer gradiente de velocidade. De acordo com a lei da viscosidade de Newton,

(A1.2)

onde τ é a tensão de cisalhamento, dv/dy é o gradiente de velocidade e μ a viscosidade ou viscosidade dinâmica. A viscosidade cinemática v é dada por

(A1.3)

onde ρ é a densidade do fluido.

A compressibilidade de um fluido reflete suas propriedades tensão-deformação. A tensão é a resposta interna de um material a uma pressão externa. Nos fluidos, a tensão é transmitida pela pressão ao qual o fluido está submetido. A deformação é a medida da deformação linear ou volumétrica de um material tensionado. Nos fluidos, a deformação toma a forma do volume reduzido (e com densidade aumentada) sob pressões de fluido crescentes. A compressibilidade da água β foi amplamente discutida na Seção 2.9 e é definida pela Eq. (2.44).

A densidade, viscosidade e compressibilidade da água são funções da temperatura e pressão (Dorsey, 1940; Weast, 1972). Em geral, sua variação não é grande, e para a faixa de pressões e temperaturas que ocorrem na maioria das aplicações de águas subterrâneas, é comum considerá-las constantes. A 15,5 °C, ρ = 1,0 g/cm³, μ = 1,124 cP, e β = 4,4 × 10^–2 m²/N.

A pressão do fluido p em qualquer corpo imerso em água é a força por unidade de área que atua nesse ponto. Sob condições hidrostáticas, a pressão do fluido em um ponto reflete o peso da coluna d’água sobre a área da seção transversal em torno desse ponto. Assim, é possível expressar a pressão em relação à pressão absoluta mas, geralmente, é expressa em relação à pressão atmosférica. Neste último caso, a pressão é chamada de pressão relativa, pois esta é a pressão de leitura que é obtida em medidores zerados para a atmosfera.

A carga de pressão ψ em um determinado ponto do fluido é a altura que uma coluna d’água alcançaria em um manômetro instalado nesse ponto. Em um corpo d’água em repouso, ψ é igual à profundidade do ponto de medição abaixo da superfície. Se p é definido como pressão relativa, ψ é definida pela relação

(A1.4)

Em outras palavras, a carga de pressão ψ é uma medida da pressão p do fluido.

As pressões do fluido também são desenvolvidas nas águas subterrâneas que fluem nos solos através dos poros nas formações geológicas. Na Seção 2.2, os elementos da mecânica dos fluidos apresentados neste apêndice são aplicados no desenvolvimento da teoria do escoamento das águas subterrâneas.

Um desenvolvimento rigoroso da equação para o escoamento transiente em meios porosos saturados deve reconhecer o fato que mudanças transientes de pressão no fluido resultam em deformações da estrutura granular de um meio poroso e essas transformações implicam que o meio, assim como a água, está em movimento. Esta percepção cria a necessidade de dois refinamentos na clássica derivação apresentada por Jacob (1940) e apresentada na Seção 2.11 deste texto. Primeiramente, como reconhecido por Biot (1955), é necessário reescrever a lei de Darcy em termos da velocidade relativa do fluido em relação aos grãos. Em segundo lugar, como reconhecido por Cooper (1966), é necessário considerar a conservação da massa para o meio, bem como para o fluido no volume de controle elementar. Podem ser desenvolvidas relações de continuidade de uma das três seguintes maneiras: (1) pela consideração de um volume de controle deformável em coordenadas deformáveis; (2) pela consideração de um volume elementar deformável em coordenadas fixas; ou (3) pela consideração de um volume elementar fixo em coordenadas fixas. Conforme Gambolati & Freeze (1973), usaremos um volume elementar fixo em coordenadas fixas. A abordagem requer o uso de notação vetorial e derivada material (derivada total, derivada substantiva). Se estes conceitos não são familiares, Aris (1962) e Wills (1958) apresentam tratamentos introdutórios. O desenvolvimento será apresentado aqui para um meio isotrópico, homogêneo, com condutividade hidráulica K, porosidade n, e compressibilidade α. A mesma abordagem é facilmente adaptada para meios heterogêneos e anisotrópicos.

Em notação vetorial, a forma tridimensional da equação lei de Darcy (Eq. (2.34) é:

(A2.1)

onde  é a velocidade relativa do fluido em relação aos grãos e ∇h = (∂h/∂x, ∂h/∂y, ∂h/∂z) é o gradiente hidráulico. Podemos expandir o vetor  como:

(A2.2)

onde é a velocidade do fluido e  é a velocidade do meio deformável.

A equação de estado para a água [Eq. (2.47)] é

(A2.3)

e para os grãos de solo, que são incompressíveis, é:

ρ_s = constante (A2.4)

A equação da continuidade para a água é

(A2.5)

e para o solo é

(A2.6)

Nestas equações, é o operador divergente:

Expandindo a Eq. (A2.5), chegamos a

(A2.7)

Cancelando ρ_sda Eq. (A2.6) e reescrevendo aquela equação, obtemos uma expressão para ∂n/∂t. Essa pode ser substituída na Eq. (A2.7) juntamente com uma expressão para obtida da Eq. (A2.2). Dividindo tudo por ρ e rearranjando, a Eq. (A2.7) torna-se:

(A2.8)

Se utilizarmos a derivada material , a Eq. (A2.8) pode ser escrita como:

(A2.9)

O primeiro termo do lado direito da Eq. (A2.9) pode ser relacionado à compressibilidade da água β

(A2.10)

A derivada material no lado direito da Eq. (A2.10) pode ser substituída por uma derivada parcial somente se a desigualdade é satisfeita:

(A2.11)

Assumindo para o lado esquerdo da Eq. (A2.9) que

(A2.12)

Então, substituindo as Eqs. (A2.1) e (A2.10) na Eq. (A2.9) resulta em

(A2.13)

Em um campo de tensões tridimensional, o vetor de velocidade dos grãos = (v_sx, v_sy, v_sz) é relacionado ao vetor de deformação (ou deslocamento do solo) = (u_sx, u_sy, u_sz) por

(A2.14)

Em um campo de tensões unidimensional,

(A2.15)

Se as condições da Eq. (A2.15) são satisfeitas, o termo final da Eq. (A2.13) pode ser expandido (Cooper, 1966; Gambolati & Freeze, 1973; Gambolati, 1973a) como

(A2.16)

em que α é a compressibilidade vertical do meio poroso. A mudança de derivada em torno da igualdade central é válida para um vetor de posição mas não em geral. A derivada material na expressão da direita na Eq. (A2.16) pode ser substituída pela derivada parcial se a Eq. (A2.11) é satisfeita. Nesse caso, a Eq. (A2.13) resulta em

(A2.17)

Como ρ = ρg(h – z) e K é uma constante, a Eq. (A2.17) é simplificada para

(A2.18)

Ou, relembrando que S_s = ρ g(α + nβ) e expandindo a notação vetorial,

(A2.19)

A Equação (A2.19) é idêntica à Eq. (2.75) desenvolvida por Jacob (1940). O desenvolvimento mais rigoroso demonstra que a validade da equação clássica de escoamento é baseada na satisfação das desigualdades das Eqs. (A2.11) e (A2.12) e das condições de tensão da Eq. (A2.15). É improvável que essas condições sejam sempre satisfeitas. Gambolati (1973b) mostra que quando a taxa de consolidação supera a taxa de percolação do fluido , como pode ocorrer em espessas camadas argilosas de baixa permeabilidade e alta compressibilidade, a desigualdade pode não ser satisfeita. Com respeito à condição de tensão, o termo no final da Eq. (A2.13) é realmente a saliente ponta do iceberg que relaciona o campo de escoamento tridimensional ao campo tridimensional de tensões. Biot (1941, 1955) apresentou essa relação e Verruijt (1969) proveu uma clara derivação. Schifman et al. (1969) desenvolveram uma comparação entre o procedimento clássico e de Biot, e Gambolati (1974) analisou a faixa de validade da clássica equação de escoamento.

Considere o problema simples de fluxo de água subterrânea mostrado na Figura 2.25 (a). A equação de fluxo saturado em estado estacionário no plano xy é:

(A3.1)

A definição matemática das condições de contorno é:

em e (A3.2)

em (A3.3)

em (A3.4)

Resolveremos a equação para h(x, y) usando a técnica de separação de variáveis.

Suponha que a solução seja resultado de um produto da seguinte forma:

(A3.5)

A Equação (A3.1) então torna-se

(A3.6)

Dividindo por XY obtém-se

(A3.7)

O lado esquerdo da equação é independente de y. Assim, o lado direito da equação, apesar de sua aparência, deve também ser independente de y, uma vez que é igual identicamente ao lado esquerdo. Similarmente, o lado direito da equação é independente de x, e da mesma forma o lado esquerdo. Se ambos os lados são independentes de x e y, cada lado deve ser igual a uma constante. Dessa forma:

e (A3.8)

A constante G pode ser positiva, negativa ou nula. Todos os três casos resultam em soluções de produtos, mas apenas o caso de G = 0 leva a uma solução que possui significado físico para este problema. Deste modo, as equações (A3.8) resultam em:

e (A3.9)

Essas são equações diferenciais ordinárias, cujas soluções são conhecidas como:

e (A3.10)

A solução produto da Eq. (A3.5) torna-se:

(A3.11)

Pode-se determinar os coeficientes A, B, C e D recorrendo às condições de contorno. Diferenciando a Eq. (A3.11) em relação a y produz

(A3.12)

e recorrendo à Eq. (A3.2) implica que C = 0. A partir da Eq. (A3.11), resulta em

(A3.13)

Valendo-se das condições de contorno das Eqs. (A3.3) e (A3.4) produz F = h₀ e E = –(h₀ – h₁)/x_L. Desse modo, a solução é

(A3.14)

Essa equação é idêntica à Eq. (2.81), apresentada na Seção 2.11 sem demonstração.

É imediatamente claro que a Eq. (A3.14) satisfaz as condições de contorno das Eqs. (A3.3) e (A3.4). Diferenciando com respeito a y produz 0, satisfazendo a Eq. (A3.2). Diferenciando duas vezes com relação a x também produz 0, logo a solução Eq.(A3.14) satisfaz a equação de fluxo [Eq. (A3.1)].

Expressão de Debye-Hückel para atividades iônicas individuais:

å × 108	Íon
2,5
3,0	K+, Cl–,
3,5	OH–, HS–, , F–
4,0	SO42–, PO43–, HPO42–
4,0–4,5	Na+, , ,
4,5	CO32–, SO32–
5	Sr2+, Be2+, S2–
6	Ca2+, Fe2+, Mn2+
8	Mg2+
9	H+, Al3+, Fe3+

Carga	Tamanho,* å	Íons	l = 0,0005	0,001	0,0025	0,005	0,01	0,025	0,05	0,1
1	2.5	Rb+, Cs+, Ag+, , Tl+	0,975	0,964	0,945	0,924	0,898	0,85	0,80	0,75
	3	K+, Cl–, Br–, I–, CN–, , , OH–, F–,	0,975	0,964	0,945	0,925	0,899	0,85	0,805	0,755
	4	Na+, , , , , ,	0,975	0,964	0,947	0,928	0,902	0,86	0,82	0,775
	6	Li+,	0,975	0,965	0,948	0,929	0,907	0,87	0,835	0,80
	9	H+	0,975	0,967	0,950	0,933	0,914	0,88	0,86	0,83
2	4.5	Pb2+, Hg22+, SO42–, CrO42–, CO32–, SO32–, C2O42–, S2O32–, H citrato2–	0,903	0,867	0,805	0,742	0,665	0,55	0,455	0,37
	5	Sr2+, Ba2+, Cd2+, Hg2+, S2–, WO42–	0,903	0,868	0,805	0,744	0,67	0,555	0,465	0,38
	6	Ca2+, Cu2+, Zn2+, Sn2+, Mn2+, Fe2+, Ni2+, Co2+, ftalato2–	0,905	0,870	0,809	0,749	0,675	0,57	0,485	0,405
	8	Mg2+, Be2+	0,906	0,872	0,813	0,755	0,69	0,595	0,52	0,45
3	4	PO43–, Fe(CN)63–, Cr(NH3)63+	0,796	0,725	0,612	0,505	0,395	0,25	0,16	0,095
	9	Al3+, Fe3+, Cr3+, Sc3+, In3+, e terras raras	0,802	0,738	0,632	0,54	0,445	0,325	0,245	0,18

Referências

BERNER, R. A. 1971. Principles of Chemical Sedimentology. McGraw-Hill, New York, 240 pp.

KLOTZ, I. M. 1950. Chemical Thermodynamics. Prentice-Hall, Eaglewood Cliffs, N.J., 369 pp.

MANOV, G. G., R. G. BATES, W. J HAMER, & S. F. ACREE. 1943. Values of the constants in the Debye-Hückel equation for activity coefficients. J. Amer. Chem. Soc., 65, pp. 1765–1767.

=

	erf ()	erfc ()
0	0	1,0
0,05	0,056372	0,943628
0,1	0,112463	0,887537
0,15	0,167996	0,832004
0,2	0,222703	0,777297
0,25	0,276326	0,723674
0,3	0,328627	0,671373
0,35	0,379382	0,620618
0,4	0,428392	0,571608
0,45	0,475482	0,524518
0,5	0,520500	0,479500
0,55	0,563323	0,436677
0,6	0,603856	0,396144
0,65	0,642029	0,357971
0,7	0,677801	0,322199
0,75	0,711156	0,288844
0,8	0,742101	0,257899
0,85	0,770668	0,229332
0,9	0,796908	0,203092
0,95	0,820891	0,179109
1,0	0,842701	0,157299
1,1	0,880205	0,119795
1,2	0,910314	0,089686
1,3	0,934008	0,065992
1,4	0,952285	0,047715
1,5	0,966105	0,033895
1,6	0,976348	0,023652
1,7	0,983790	0,016210
1,8	0,989091	0,010909
1,9	0,992790	0,007210
2	0,995322	0,004678
2,1	0,997021	0,002979
2,2	0,998137	0,001863
2,3	0,998857	0,001143
2,4	0,999311	0,000689
2,5	0,999593	0,000407
2,6	0,999764	0,000236
2,7	0,999866	0,000134
2,8	0,999925	0,000075
2,9	0,999959	0,000041
3	0,999978	0,000022

A equação diferencial parcial que descreve o escoamento estacionário em duas dimensões, em meio homogêneo e isotrópico (Seção 2.11) é a equação de Laplace:

(A6.1)

Para encontrar a equação de diferenças finitas para um nó interior da malha de nós usada para discretizar a região de escoamento, devemos substituir as derivadas parciais de segunda ordem da Eq. (A6.1) por quocientes de diferenças. Vamos considerar o primeiro termo da equação. Lembre-se que a definição de derivada parcial com respeito a x de uma função de duas variáveis h(x, z), é:

(A6.2)

Em um computador, é impossível obter o limite quando , mas é possível aproximar o limite atribuído, arbitrariamente, um pequeno valor para Δx. De fato, podemos fazer isso discretizando uma rede nodal com o espaçamento da malha Δx.

Para qualquer valor de z, por exemplo z₀, podemos expandir h(x, z₀) em uma Série de Taylor em torno do ponto (x₀, z₀), como:

(A6.3)

Se definirmos x = x₀ + Δx (isto é conhecido como diferença progressiva), e desprezarmos os termos de ordem maior que a unidade, podemos aproximar ∂h/∂x por:

(A6.4)

Os termos desprezados da expansão de Taylor representam um erro de truncamento na aproximação por diferenças finitas.

Podemos obter uma expressão similar a (A6.4) substituindo uma diferença regressiva, x = x₀ – Δx, na equação (A6.3). Isto produz:

(A6.5)

Para obter uma aproximação para ∂₂h/∂x₂, escrevemos a equação diferencial em termos de ∂h/∂x, utilizando uma expressão de diferença progressiva.

(A6.6)

e substituímos a expressão de diferença regressiva Eq. (A6.5) na Eq. (A6.6) para obter

(A6.7)

De maneira similar, podemos desenvolver uma equação diferencial para ∂₂h/∂x₂

(A6.8)

Para uma malha quadrada, Δx = Δz, somando as Eqs. (A6.7) e (A6.8), a equação de Laplace, resulta

(A6.9)

Se deixarmos (x₀, z₀) ser o ponto nodal (i, j), a Eq. (A6.9) pode ser rearranjada para resultar:

(A6.10)

que é idêntica à Eq. (5.24).

Toth (1962, 1963) apresentou duas soluções analíticas para o problema de valor de contorno para escoamento permanente vertical, bidimensional, saturado, homogêneo, isotrópico, limitado no topo pela superfície freática e nos outros três lados por fronteiras impermeáveis (área hachurada da Figura 6.1, reproduzida em um sistema de coordenadas xz na Figura A7.1)

Primeiro ele considerou o caso em que a superfície freática é uma linha reta de declividade constante. Nesse caso, a região de escoamento na Figura A7.1 é a região ABCEA. Como não é possível obter uma solução analítica em um domínio trapezoidal, Toth aproximou o campo real de escoamento pela região sombreada ABCDA. Ele projetou os valores de carga hidráulica da superfície real do lençol freático AE sobre a fronteira superior AD da região de solução. A aproximação é satisfatória para α pequeno.

A equação de fluxo é a equação de Laplace:

(7.1)

Para a região limitada por x = s e z = z₀, as condições de contorno são:

em AB e CD

em BC (7.2)

em AD

em que c = tan α.

A solução analítica, obtida por separação de variáveis, é:

(7.3)

Esta equação satisfaz a equação de fluxo (A7.1) e as condições de contorno (A7.2). Quando plotada, leva à rede equipotencial mostrada na Figura A7.1. O escoamento ocorre da fronteira de recarga DF através do domínio até a fronteira de descarga AF.

Toth também considerou o caso em que a superfície freática é representada como uma curva senoidal sobreposta à linha AE (para representar uma topografia ondulada). A condição de contorno torna-se, então:

em AD (A7.4)

em que , e , a’ representa a amplitude da curva senóide e b’ sua frequência.

A solução analítica para este caso assume a forma:

(A7.5)

A rede equipotencial descrita por essa função é similar à Figura 6.2 (b).

Considere um escoamento unidimensional, vertical e homogêneo com o seu limite superior na superfície do solo e o seu limite inferior abaixo da superfície freática. Como descrito na Seção 6.4, a equação do escoamento para esse sistema é:

(A8.1)

e as condições de contorno são:

(A8.2)

no topo, e

(A8.3)

na base.

A solução é obtida em função da carga de pressão, . A posição da superfície freática em qualquer instante de tempo é dada pelo valor de z, no qual . As informações necessárias incluem as condições iniciais , a taxa de precipitação R, a taxa de recarga de água subterrânea Q, e as curvas características não saturadas K(ψ) e C(ψ). Para , K = K₀ e C = 0, sendo , a carga de pressão de entrada de ar.

O esquema numérico de diferenças finitas usado por Rubin & Steinhardt (1963), Liakopoulos (1965b) e Freeze (1969b) é o esquema implícito de Richtmyer (1957). Nesse método, o plano (z, t) é representado por uma malha retangular de pontos j = 1, 2, …, L ao longo do eixo z, e n = 1, 2, … ao longo do eixo t. A distância entre os nós verticais é Δz, e entre os intervalos de tempo, Δz. A solução em qualquer ponto (j, n) da malha é . Sob essa notação, o esquema de diferenças finitas para a Eq. (A8.1) para um nó interior (j = 2 até j = L – 1), depois dos rearranjos pertinentes, é:

(A8.4)

em que os valores de , , e que determinam os valores de K e de C que serão aplicados em um nó em um determinado instante de tempo, são determinados por extrapolação de instantes anteriores, como descrito por Rubin & Steinhardt (1963).

Equações de diferenças finitas que incorporam as condições de contorno (A8.2) e (A8.3) podem ser escritas para os nós nas fronteiras superior e inferior (j = 1, j = L).

Para cada intervalo de tempo, as aproximações por diferenças finitas constituem um sistema de L equações algébricas lineares, com L incógnitas. A fórmula geral é:

(A8.5)

em que os coeficientes A, B, C, e D variam com o nó (j = 1, j = 2 à L – 1, j = L), com o instante de tempo (n = 1, n = 2, n > 2), e com a saturação . As variáveis, das quais A, B, C, e D dependem, são os valores de do instante de tempo anterior, os valores de contorno R e Q, e as relações funcionais K(ψ) e C(ψ).

Os valores de são calculados pela seguinte relação recorrente:

(A8.6)

em que

Os coeficientes E e F são calculados de j = 1 áte j = L, e os ψ’s são calculados forma reversa, de j = L até j = 1, usando a Eq. (A8.6). Para j = 1, E₁ = A₁/B₁ e F₁ = D₁/B₁. Se o nó superior (j = L) está saturado, surgem complicações.

Se K(ψ) e C(ψ) apresentam histerese, essas relações são implementadas no programa computacional na forma de tabelas de valores, representando as curvas de umedecimento e de drenagem. O programa deve localizar a curva correta, considerando se o nó analisado está umedecendo ou secando, e verificar a evolução do nó para determinar o valor de ψ para o qual ocorreu uma mudança de umedecimento para drenagem, ou vice versa.

A equação diferencial parcial que descreve o escoamento transiente através de um meio anisotrópico saturado (Seção 2.11) é:

(A9.1)

Para um aquífero horizontal, confinado de espessura, b, a forma bidimensional da equação (A9.1) se reduz para:

(A9.2)

em que T_x e T_y são os componentes principais do tensor de transmissividade definido por T = Kb, e S é o armazenamento definido por S = S_sb. Para encontrar a equação de diferenças finitas para um nó no interior da malha usada para discretizar a região de escoamento, devemos substituir as derivadas parciais da Eq. (A9.2) por diferenças. Usando as definições desenvolvidas no Apêndice VI, a notação ijk da Seção 8.8 e a Figura 8.26 (c), podemos escrever a expressão de diferenças finitas:

(A9.3)

em que o subscrito  significa que a quantidade entre parênteses é avaliada no ponto médio entre os nós (i, j) e (i + 1, j), e o sobrescrito k significa que a quantidade entre parênteses é avaliada no instante de tempo k. Os termos do lado direito Eq. (A9.3) podem ser aproximados por:

(A9.4A)

(A9.4B)

Se definirmos e  como médias simples na forma

(A9.5)

essas expressões podem ser substituídas nas Eqs. (A9.4), e as Eqs. (A9.4) podem subsequentemente ser substituídas na Eq. (A9.3) para obter:

(A9.6)

Similarmente:

(A9.7)

Finalmente, podemos aproximar o lado direito da equação (A9.2) por:

(A9.8)

Substituindo as Eqs. (A9.6), (A9.7) e (A9.8) para os três termos da Eq. (A9.2) e agrupando os termos, obtemos a equação geral de diferenças finitas para um nó interno em um aquífero heterogêneo e anisotrópico:

(A9.9)

em que

Se o aquífero é homogêneo e isotrópico, T_x = T_y = T para todo (i, j) e S_{i, j} = S para todo (i, j) Sob essas condições, e para uma malha quadrada com Δx = Δy, os coeficientes da Eq. (A9.9) tornam-se:

Se dividirmos por T/Δx², percebe-se que esses coeficientes são os mesmos que aqueles desenvolvidos de uma maneira menos rigorosa na Seção 8.8 e apresentados em conexão com a Eq. (8.57).

Trescott et al. (1976) sugeriram que existem algumas vantagens em utilizar a média harmônica ao invés da média aritmética na Eq. (A9.5). Essa abordagem muda os coeficientes na equação de diferenças finitas, mas não altera os conceitos contidos no desenvolvimento.

A equação (A9.9) é escrita em termos dos valores de carga hidráulica em cinco nós no instante de tempo k, e um nó no instante de tempo k – 1. Essa forma é conhecida como aproximação por diferenças regressivas (ou atrasada). Remson et al. (1971) notam que existem algumas vantagens computacionais em usar a aproximação por diferença central, conhecida como esquema Crank-Nicholson, que utiliza valores de carga hidráulica em cinco nós no instante de tempo k e cinco nós no instante de tempo k – 1. O procedimento implícito de direções alternadas (ADIP) utilizado por Pinder & Bredehoeft (1968) envolve duas equações de diferenças finitas, uma no plano xt e outra no plano yt. Cada uma utiliza valores de carga em três nós no instante de tempo k, e três nós no instante de tempo k – 1.

A equação derivada aqui é uma representação da lei de conservação de massa. A derivação se baseia em Ogata (1970) e Bear (1972). Ela assume que o meio poroso é homogêneo e isotrópico, encontra-se saturado, o fluxo da água em estado estacionário e a lei de Darcy é válida. Sob essa hipótese, o fluxo é descrito pela velocidade linear média, a qual transporta a substância dissolvida por advecção. Se esse fosse o único mecanismo de transporte atuante, os solutos não-reativos transportados pelo fluido se moveriam como um plugue. Na realidade, há um processo de mistura adicional, a dispersão hidrodinâmica (Seção 2.13), que é causada pelas variações na velocidade microscópica dentro de cada canal de poro e de um canal para o outro. Se desejarmos descrever o processo de transporte em uma escala macroscópica, usando parâmetros macroscópicos e levando em consideração o efeito de mistura microscópica, é necessário introduzir um segundo mecanismo de transporte, em adição à advecção, para considerar a dispersão hidrodinâmica.

Para estabelecer uma demonstração matemática da equação de conservação de massa, o fluxo do soluto entrando e saindo de um volume elementar representativo do meio poroso será considerado (Figura A10.1). Em coordenadas cartesianas a descarga específica v possui componentes (v_x, v_y, v_z) e a velocidade linear média possui componentes . A taxa de transporte advectivo é igual a . A concentração de um soluto C é definido como a massa de soluto por unidade de volume da solução. A massa de soluto por unidade de volume do meio poroso é, por consequência, nC. Para um meio homogêneo, a porosidade n é uma constante e ∂(nC)/∂x = ∂nC/∂x

A massa do soluto transportado na direção x pelos dois mecanismos de transporte do soluto pode ser representada como

transporte por advecção = (A10.1)

transporte por dispersão = (A10.2)

onde D_x é o coeficiente de dispersão na direção x e dA é a área da seção elementar do elemento cúbico. O coeficiente de dispersão D_x está relacionado à dispersividade α_x e ao coeficiente de difusão D* pela Eq. (9.4):

(A10.3)

A forma da componente de dispersão incorporada na Eq. (A10.2) é análoga à primeira Lei de Fick.

Se F_x representa a massa total do soluto por unidade de área transportada na direção x por unidade de tempo, então:

(A10.4)

O sinal negativo antes do termo dispersivo indica que o contaminante se move para a zona de menor concentração. Da mesma forma, expressões nas outras duas direções são escritas

(A10.5)

(A10.6)

A quantidade total de soluto entrando no elemento cúbico (Figura A10.1) é

onde os termos parciais indicam a variação espacial da massa de soluto em uma direção específica. A diferença na quantidade que entra e que sai do elemento é, portanto,

Porque a substância dissolvida é assumida como não-reativa, a diferença entre o fluxo que entra e que sai é igual à quantidade da substância acumulada no próprio elemento. A taxa de variação de massa no elemento é

A equação de conservação de massa completa, portanto, torna-se

(A10.7)

Substituindo-se as expressões (A10.4), (A10.5) e (A10.6) na Eq. (A10.7) e cancelando n de ambos lados da equação tem-se que :

(A10.8)

Em meios homogêneos em que a velocidade é estável e uniforme (p.e., se essa não varia no tempo ou espaço), os coeficientes de dispersão D_x, D_y e D_z não variam através do espaço e a Eq. (A10.8) torna-se

(A10.9)

Em uma dimensão

(A10.10)

Em alguma aplicações, a direção unidimensional é tomada como uma coordenada curvilinear na direção do fluxo ao longo da linha de fluxo. A equação de transporte torna-se, então,

(A10.11)

onde l é a direção ao longo da linha de fluxo, D_l o coeficiente de dispersão longitudinal e a velocidade linear média ao longo da linha de fluxo.

Para um problema bidimensional, é possível definir duas direções de coordenadas curvilíneas, S_l e S_t, onde S_l é direcionada ao longo da linha de fluxo e S_t é ortogonal a esta. A equação de transporte resulta em

(A10.12)

onde D_l e D_t são os coeficientes de dispersão longitudinal e transversal, respectivamente.

Se v_l varia ao longo de uma linha de corrente e D_l e D_t variam no espaço, a equação (A10.12) resulta em

(A10.13)

As equações (A10.8) a (A10.13) representam seis formas da equação de advecção-dispersão para o transporte de solutos em meios porosos saturados. A equação (A10.11) é idêntica à Eq. (9.3) na Seção 9.2. A solução de qualquer uma dessas equações fornece a concentração do soluto C em função do espaço e tempo. Ela tomará a forma C(x, y, z, t) para as Eqs. (A10.8) e (A10.9); C(l, t) para a Eq. (A10.11) e C(S_l, S_t, t) para as Eqs. (A10.12) e (A10.13). Existem várias soluções analíticas conhecidas para as formas mais simples da equação de transporte. A Eq. (9.5) da Seção 9.2 é uma solução analítica da Eq. (A10.11), e a Eq. (9.6) é uma solução analítica da Eq. (A10.9). Em muitas situações de campo, análises bi- ou tridimensionais são necessárias. Adicionalmente, as velocidades são raramente uniformes e as dispersões são geralmente variáveis no espaço. Para essas condições, métodos numéricos adaptados para computadores devem ser usados para obter tais soluções.

O coeficiente de dispersão em meio tridimensional em meios porosos, isotrópicos, homogêneos como expresso na Eq. (A10.12) é um tensor simétrico de segunda ordem com nove componentes. D_l e D_t são os termos diagonais da forma bidimensional. As propriedades direcionais do coeficiente de dispersão são causadas pela natureza do processo de dispersão nas direções transversal e longitudinal. Se o meio é anisotrópico, a descrição matemática do processo de dispersão aumenta de complexidade. Não existem soluções numéricas ou analíticas disponíveis para sistemas anisotrópicos. Apenas para meios isotrópicos, o coeficiente de dispersão foi representado com sucesso por métodos experimentais.

É possível estender a equação de transporte para incluir os efeitos de retardação do transporte do soluto por adsorção, reação química, transformação biológica ou decaimento radioativo. Neste caso, o balanço de massa no volume elementar deve incluir um termo “sumidouro”. Para a retardação devido a adsorção, a equação de transporte em um meio homogêneo em um sistema unidimensional ao longo da direção de fluxo toma a forma

(A10.14)

em que ρb é a massa específica (densidade absoluta) do meio poroso, n a porosidade e S a massa do constituinte químico adsorvido a uma unidade de massa da parte sólida do meio poroso. O primeiro termo da Eq. (A10.14) é o termo de dispersão, o segundo é o termo de advecção e o terceiro é o termo de reação. Este termo de reação é tratado na Seção 9.2, em conexão com as Eqs. (9.9) a (9.13). Soluções analíticas para a Eq. (A10.14) são apresentadas por Codell & Schreiber (in press). A Figura 9.12 apresenta exemplos de soluções realizadas por Pickens & Lennox (1976) usando métodos numéricos.