Pular para o conteúdo principal

Modelo de águas subterrâneas

1. Visão geral

O componente de águas subterrâneas de HydroPol2D representa a dinâmica acoplada de:

  • armazenamento na zona vadose,
  • recarregar da infiltração,
  • fluxo lateral saturado de águas subterrâneas,
  • exfiltração de águas subterrâneas para a superfície,
  • e interação das águas subterrâneas com as células dos rios.

A implementação atual combina três componentes complementares:

  1. um módulo vadose/recarga em camadas, que move a água através do armazenamento próximo à superfície, zona raiz e zona de transmissão usando um fechamento Darcy reduzido consistente com Richards com base nas relações constitutivas van Genuchten-Mualem,

  2. um módulo de troca de ascensão capilar, que pode mover a água para cima, de águas subterrâneas rasas para o armazenamento vadoso disponível, e

  3. um solucionador bidimensional de fluxo de água subterrânea Boussinesq, que direciona o fluxo saturado lateralmente através do domínio e retorna a exfiltração para a superfície da terra quando o lençol freático sobe acima da elevação da superfície local.

Esta estrutura permite que HydroPol2D represente o acoplamento dinâmico entre:

  • águas superficiais,
  • armazenamento vadoso,
  • recarga de águas subterrâneas,
  • fluxo saturado de águas subterrâneas,
  • e fluxo de retorno com excesso de saturação.

Figura 1. Representação conceitual do componente aquífero em HydroPol2D.


2. Estrutura Conceitual Geral

O modelo de águas subterrâneas é acoplado a dois relógios. A troca da zona vadosa é avaliada a cada intervalo de tempo da água superficial porque a infiltração e a evapotranspiração precisam de feedback rápido da umidade do solo. As cabeceiras das águas subterrâneas saturadas são atualizadas com menos frequência, geralmente em um intervalo de tempo alvo diário, a menos que a estabilidade ou os limites de mudança de carga exijam uma atualização mais curta. A lógica é:

  1. calcular a posição do lençol freático a partir da atual altura da água subterrânea;
  2. determinar o armazenamento não saturado máximo permitido pela profundidade atual do lençol freático;
  3. calcular a drenagem/recarga das camadas vadosas para as águas subterrâneas durante a atual etapa de superfície;
  4. calcular a ascensão capilar das águas subterrâneas para as camadas vadosas quando ativado,
  5. acumular o volume líquido de troca no programador de águas subterrâneas,
  6. atualizar as cabeças saturadas de água subterrânea quando o cronograma for devido,
  7. resolver o fluxo lateral de águas subterrâneas usando uma equação 2D Boussinesq quando o fluxo de base está ativo;
  8. calcular a exfiltração de águas subterrâneas para a superfície;
  9. armazenamento vadoso correto se o aumento do lençol freático reduzir o armazenamento insaturado disponível,
  10. atualizar a profundidade da água superficial em conformidade.

As principais variáveis ​​são:

  • hh = carga hidráulica da água subterrânea [L][\mathrm{L}]
  • z0z_0 = elevação da base do aquífero ou elevação do leito rochoso [L][\mathrm{L}]
  • H=hz0H = h - z_0 = espessura saturada [L][\mathrm{L}]
  • RR = taxa de recarga para águas subterrâneas [LT1][\mathrm{L}\,\mathrm{T}^{-1}]
  • qexfq_{\mathrm{exf}} = fluxo de exfiltração das águas subterrâneas para a superfície [LT1][\mathrm{L}\,\mathrm{T}^{-1}]
  • SyS_y = rendimento específico [][-]
  • KK = condutividade hidráulica saturada [LT1][\mathrm{L}\,\mathrm{T}^{-1}]

2.1 Agendador assíncrono

O agendador assíncrono fica ativo por padrão quando a modelagem de águas subterrâneas está habilitada. Os principais controles são:

ControlarPadrãoSignificado
flag_groundwater_async1Usa atualizações mais lentas de águas subterrâneas saturadas em vez de atualizar cabeças a cada passo da superfície.
groundwater_target_dt_min1440Intervalo alvo de atualização das águas subterrâneas, em minutos.
groundwater_min_dt_min1Intervalo mínimo de atualização das águas subterrâneas, em minutos.
groundwater_max_head_change_m0,25Força uma atualização quando o câmbio acumulado mudaria a cabeça nesse valor.
groundwater_courant0,25Fator de estabilidade para a estimativa do intervalo de tempo das águas subterrâneas Boussinesq.
flag_capillary_rise1Permite a ascensão capilar em camadas simples quando o armazenamento em camadas do solo está presente.

A recarga pendente, a ascensão capilar e a troca líquida são armazenadas como variáveis ​​de estado do escalonador. Quando uma atualização das águas subterrâneas é executada, a profundidade líquida de troca acumulada é dividida pelo intervalo decorrido das águas subterrâneas para formar o termo de recarga passado para o solucionador saturado.

3. Posição do lençol freático e capacidade da zona não saturada

HydroPol2D determina a posição do lençol freático em relação à superfície da terra usando a altura atual da água subterrânea. Deixar:

  • zsurfz_{\mathrm{surf}} = elevação da superfície terrestre [L][\mathrm{L}]
  • DsoilD_{\mathrm{soil}} = profundidade do solo [L][\mathrm{L}]
  • hh = carga hidráulica da água subterrânea [L][\mathrm{L}]

A espessura saturada acima do leito rochoso é calculada como

zgw=h(zsurfDsoil)z_{\mathrm{gw}} = h - \left(z_{\mathrm{surf}} - D_{\mathrm{soil}}\right)

A profundidade da superfície da terra até o lençol freático é então

zwt=Dsoilzgwz_{\mathrm{wt}} = D_{\mathrm{soil}} - z_{\mathrm{gw}}

com a restrição

zwt0z_{\mathrm{wt}} \ge 0

onde:

  • zgwz_{\mathrm{gw}} = espessura saturada medida para cima a partir da base do aquífero [L][\mathrm{L}]
  • zwtz_{\mathrm{wt}} = profundidade até o lençol freático abaixo da superfície do solo [L][\mathrm{L}]

Esta variável é central para o acoplamento entre o armazenamento vadoso e as águas subterrâneas.

3.1 Armazenamento máximo não saturado

A água máxima que pode ser armazenada na zona não saturada é definida usando a convenção atual armazenamento ativo acima do conteúdo de água residual:

SUZ,max=zwt(θsatθr)S_{\mathrm{UZ,max}} = z_{\mathrm{wt}} \left( \theta_{\mathrm{sat}} -\theta_r \right)

onde:

  • SUZ,maxS_{\mathrm{UZ,max}} = armazenamento não saturado máximo [L][\mathrm{L}]
  • θsat\theta_{\mathrm{sat}} = conteúdo volumétrico saturado de água [][-]
  • θr\theta_r = teor volumétrico de água residual [][-]

Assim, à medida que o lençol freático sobe, a zona vadosa diminui e a capacidade de armazenamento disponível para infiltração diminui. Quando o armazenamento em camadas do solo está ativo, o código avalia este limite através das capacidades explícitas dos armazenamentos próximos à superfície, zona radicular e zona de transmissão após truncamento pela posição atual do lençol freático.

4. Troca de Vadose em Camadas Padrão para Água Subterrânea

A recarga é calculada a partir do armazenamento vadoso atual após o módulo de infiltração já ter movido a água da superfície para o solo. Isso evita infiltração de contagem dupla. O caminho padrão atual é dividido em camadas:

  • armazenamento próximo à superfície,
  • armazenamento na zona raiz,
  • armazenamento na zona de transmissão,
  • recarga de águas subterrâneas.

4.1 Acoplamento em degrau superficial

A troca vadose é avaliada em cada etapa da água superficial porque a infiltração e a evapotranspiração precisam de feedback imediato da umidade do solo. Na implementação atual, a infiltração é primeiro aceita pelo módulo de infiltração, armazenada nas camadas vadosas e só então drenada para baixo em direção às águas subterrâneas.

4.2 Estado constitutivo em camadas

Para uma camada vadosa genérica \ell com armazenamento ativo SS_{\ell} e espessura zz_{\ell}, HydroPol2D avalia

θ=θr,+Sz\theta_{\ell} = \theta_{\mathrm{r},\ell}+ \frac{S_{\ell}}{z_{\ell}} Se,=θθr,θsat,θr,S_{e,\ell} = \frac{ \theta_{\ell} - \theta_{\mathrm{r},\ell} }{ \theta_{\mathrm{sat},\ell} - \theta_{\mathrm{r},\ell} } K=Ksat,Se,v[1(1Se,1/m)m]2K_{\ell} = K_{\mathrm{sat},\ell} S_{e,\ell}^{\ell_v} \left[ 1 - \left(1 - S_{e,\ell}^{1/m_{\ell}}\right)^{m_{\ell}} \right]^2

com

m=11nm_{\ell} = 1 -\frac{1}{n_{\ell}}

onde:

  • θ\theta_{\ell} = teor de água da camada [][-]
  • Se,S_{e,\ell} = saturação efetiva da camada [][-]
  • KK_{\ell} = condutividade insaturada da camada [LT1][\mathrm{L}\,\mathrm{T}^{-1}]
  • v\ell_v = Parâmetro de conectividade de poros Mualem [][-]

Esta é a mesma família constitutiva van Genuchten-Mualem usada pelo módulo de infiltração.

4.3 Cascata de drenagem descendente

O caminho de recarga em camadas é calculado como uma cascata descendente:

  1. a drenagem próxima à superfície vai para a zona radicular quando existe armazenamento radicular,
  2. a drenagem da zona raiz vai para a zona de transmissão quando existe armazenamento de transmissão,
  3. a drenagem da zona de transmissão torna-se recarga de água subterrânea.

Conceitualmente, a profundidade de drenagem de uma camada ao longo de um intervalo de tempo pode ser resumida como

ΔD\aproxKΔt\Delta D_{\ell} \aprox K_{\ell}\,\Delta t

sujeito a:

  • água disponível na camada doadora,
  • capacidade restante na camada receptora,
  • e a espessura atual da camada e os parâmetros do solo.

A recarga líquida descendente que atinge as águas subterrâneas é então

Rdown=ΔSgwΔtR_{\mathrm{down}} = \frac{\Delta S_{\rightarrow \mathrm{gw}}}{\Delta t}

onde ΔSgw\Delta S_{\rightarrow \mathrm{gw}} é a transmissão total de água para as águas subterrâneas durante a atual etapa de superfície. A drenagem direta da zona próxima à superfície ou radicular para as águas subterrâneas ocorre apenas quando as camadas vadosas inferiores estão ausentes.

4.4 Formulação substituta legada

HydroPol2D ainda contém o antigo fechamento de recarga de balde representativo para casos em que os dados explícitos de armazenamento em camadas não estão disponíveis. Nesse modo de fallback, um único estado de armazenamento vadoso é convertido em um conteúdo representativo de água e drenado com um fechamento semelhante ao Darcy. Essa formulação deve ser tratada como legado/substituto, e não como o caminho teórico principal para o modelo atual.

5. Ascensão Capilar Opcional

Quando flag_capillary_rise = 1, HydroPol2D permite que águas subterrâneas rasas subam para o armazenamento vadoso disponível. A implementação atual é uma simples regra de ascensão capilar em camadas, e não uma solução completa da equação de Richards. A troca ascendente aumenta com:

  • profundidade mais rasa do lençol freático,
  • maior déficit de armazenamento na camada receptora,
  • e maior condutividade saturada nessa camada.

Conceitualmente, a transferência capilar ascendente para a camada \ell é limitada por

Cmin ⁣(S,def,Ksat,fprox(zwt)fdry(S))C_{\ell} \le \min\!\left( S_{\ell,\mathrm{def}},\, K_{\mathrm{sat},\ell} f_{\mathrm{prox}}(z_{\mathrm{wt}}) f_{\mathrm{dry}}(S_{\ell}) \right)

onde:

  • S,defS_{\ell,\mathrm{def}} é o déficit de armazenamento da camada restante
  • fproxf_{\mathrm{prox}} é um fator monotônico de proximidade do lençol freático
  • fdryf_{\mathrm{dry}} é um fator de secura monotônico

No código atual, a ascensão capilar é avaliada primeiro na zona de transmissão, depois na zona raiz e apenas na camada próxima à superfície quando não existe armazenamento vadoso inferior. A troca líquida vadose-água subterrânea acumulada pelo programador é, portanto,

Rnet=RdownCupR_{\mathrm{net}} = R_{\mathrm{down}} - C_{\mathrm{up}}

onde CupC_{\mathrm{up}} é o fluxo de ascensão capilar ascendente total.

6. Modelo bidimensional de fluxo de águas subterrâneas Boussinesq

Uma vez calculada a recarga, HydroPol2D direciona as águas subterrâneas lateralmente usando a função PHTOKEN1XYZ_2D_explicit. Este solucionador representa o sistema de águas subterrâneas saturadas através de uma aproximação de fluxo não confinado integrada em profundidade.

6.1 Equação governante

O roteamento das águas subterrâneas segue a equação 2D Boussinesq de forma conservadora:

Syht=F+RS_y \frac{\partial h}{\partial t} =- \nabla\cdot\mathbf{F} + R

onde:

  • SyS_y = rendimento específico [][-]
  • hh = carga hidráulica da água subterrânea [L][\mathrm{L}]
  • F\mathbf{F} = vetor de fluxo de águas subterrâneas [L2T1][\mathrm{L}^2\,\mathrm{T}^{-1}]
  • RR = taxa de recarga [LT1][\mathrm{L}\,\mathrm{T}^{-1}]

HydroPol2D avalia esta equação numericamente calculando os fluxos nas interfaces das células e, em seguida, atualizando explicitamente a altura da água subterrânea.

6.2 Espessura saturada

A espessura saturada ativa é definida como

H=max(hz0,0)H = \max(h - z_0, 0)

onde:

  • HH = espessura saturada [L][\mathrm{L}]
  • z0z_0 = elevação da base do aquífero ou elevação do leito rochoso [L][\mathrm{L}]

Isto reforça a condição física de que não existe armazenamento saturado abaixo da base do aquífero.

6.3 Fluxos de interface

HydroPol2D calcula os fluxos de águas subterrâneas nas interfaces das células usando a lei de Darcy com ponderação de espessura saturada. Para a direção xx:

Fx=KxHxhxF_x = -K_x H_x \frac{\partial h}{\partial x}

e para a direção yy:

Fy=KyHyhyF_y = -K_y H_y \frac{\partial h}{\partial y}

onde:

  • Fx,FyF_x, F_y = fluxos por unidade de largura [L2T1][\mathrm{L}^2\,\mathrm{T}^{-1}]
  • Kx,KyK_x, K_y = condutividade hidráulica nas interfaces das células [LT1][\mathrm{L}\,\mathrm{T}^{-1}]
  • Hx,HyH_x, H_y = espessuras saturadas de interface [L][\mathrm{L}]

Na implementação, a condutividade hidráulica da interface é tomada como a média aritmética das condutividades nas células adjacentes:

Kx=Ki+Ki+12Ky=Kj+Kj+12K_x = \frac{K_i + K_{i+1}}{2} \qquad K_y = \frac{K_j + K_{j+1}}{2}

e a espessura saturada da interface também é calculada:

Hx=Hi+Hi+12Hy=Hj+Hj+12H_x = \frac{H_i + H_{i+1}}{2} \qquad H_y = \frac{H_j + H_{j+1}}{2}

6.4 Limitação conservadora de fluxo

Para evitar instabilidade numérica e esgotamento não físico do armazenamento de água subterrânea, HydroPol2D aplica um limitador de fluxo para que nenhuma célula possa perder mais água do que contém durante uma subetapa. Por exemplo, na direção xx, o fluxo máximo permitido é

Fx,max=SyHdonorΔtF_{x,\max} = \frac{S_y H_{\mathrm{donor}}}{\Delta t}

e o fluxo final da interface é limitado como

Fxsign(Fx)min(Fx,Fx,max)F_x \leftarrow \operatorname{sign}(F_x)\,\min\left(|F_x|,\,F_{x,\max}\right)

Uma expressão semelhante é usada na direção yy. Esta é uma característica importante do solucionador de águas subterrâneas HydroPol2D porque garante a admissibilidade da massa local durante a integração explícita.

6.5 Divergência dos fluxos de águas subterrâneas

A divergência do campo de fluxo é então calculada a partir dos fluxos de interface como

F=Fxx+Fyy\nabla\cdot\mathbf{F} = \frac{\partial F_x}{\partial x}+ \frac{\partial F_y}{\partial y}

usando diferenças finitas conservadoras na grade raster.

6.6 Atualização explícita do preditor-corretor

O solucionador Boussinesq usa uma estratégia explícita de preditor-corretor.

Etapa do preditor

Uma primeira estimativa da carga é calculada usando os fluxos atuais:

h=ht+ΔtSy(Ft)h^{\ast} = h^t + \frac{\Delta t}{S_y}\left(-\nabla \cdot \mathbf{F}^t\right)

Etapa do corretor

Os fluxos são recalculados a partir do estado do preditor, calculados em média e usados ​​na atualização final:

ht+Δt=ht+ΔtSy(Favg+R)h^{t+\Delta t} = h^t+ \frac{\Delta t}{S_y} \left( -\nabla \cdot \mathbf{F}^{\mathrm{avg}} + R \right)

onde:

  • hh^{\ast} = carga hidráulica preditora [L][\mathrm{L}]
  • ht+Δth^{t+\Delta t} = carga hidráulica final após o intervalo de tempo [L][\mathrm{L}]
  • Favg\mathbf{F}^{\mathrm{avg}} = média dos fluxos da interface do preditor e do corretor

Este procedimento preditor-corretor melhora a estabilidade e reduz a difusão numérica em relação a uma simples atualização explícita de Euler.

7. Tempo de atualização das águas subterrâneas e etapa adaptativa do tempo

HydroPol2D usa uma lógica de passo de tempo de águas subterrâneas de dois níveis.

7.1 Agendador assíncrono externo

A troca vadose é acumulada a cada etapa da água superficial, mas as alturas saturadas das águas subterrâneas são atualizadas apenas quando o cronograma das águas subterrâneas é devido. O intervalo de atualização das águas subterrâneas é o mínimo de:

  • o intervalo alvo do usuário,
  • o limite explícito de estabilidade das águas subterrâneas,
  • o intervalo implícito na mudança de carga acumulada máxima permitida,
  • e o intervalo parcial final no final da corrida.

Se flag_groundwater_async = 0, a atualização saturada da água subterrânea é forçada a cada passo da superfície.

7.2 Limite explícito de estabilidade das águas subterrâneas

O limite externo do escalonador é baseado em uma estimativa de difusividade explícita para a equação Boussinesq não confinada:

D=KHSyD = \frac{K H}{S_y} Δtgw,stable=CgwLgrid24Dmax\Delta t_{\mathrm{gw,stable}} = \frac{ C_{\mathrm{gw}} \, L_{\mathrm{grid}}^2 }{ 4D_{\max} }

onde:

  • DD = difusividade das águas subterrâneas [L2T1][\mathrm{L}^2\,\mathrm{T}^{-1}]
  • CgwC_{\mathrm{gw}} = Courant de águas subterrâneas/fator de segurança [][-]
  • Lgrid=min(Δx,Δy)L_{\mathrm{grid}} = \min(\Delta x,\Delta y)
  • DmaxD_{\max} = difusividade máxima do domínio sobre as células ativas das águas subterrâneas

Este é o limite usado pela lógica estimate_groundwater_stable_timestep atual.

7.3 Subpasso interno Boussinesq

Quando uma atualização das águas subterrâneas é acionada, o programador converte a profundidade líquida de troca acumulada em uma taxa média de recarga ao longo do intervalo decorrido das águas subterrâneas e passa-a para o solucionador saturado. O solucionador Boussinesq então avança nesse intervalo usando subetapas internas do preditor-corretor. Essas subetapas internas são adicionalmente limitadas por uma condição baseada na velocidade da forma

Δtinnermin ⁣(CvΔxux,CvΔyuy)\Delta t_{\mathrm{inner}} \le \min\!\left(\frac{C_v \Delta x}{|u_x|},\, \frac{C_v \Delta y}{|u_y|} \right)

onde:

  • ux,uyu_x, u_y são velocidades das águas subterrâneas de Darcy
  • CvC_v é o fator de segurança interno baseado na velocidade

O resultado é uma estratégia de dois níveis:

  1. um limite externo do programador baseado na estabilidade explícita das águas subterrâneas e na mudança de carga permitida,
  2. um loop de subetapa preditor-corretor interno dentro do solucionador Boussinesq.

8. Condições Limites

O solucionador de águas subterrâneas suporta:

  • Condições de contorno de Dirichlet,
  • condições de contorno sem fluxo,
  • mascaramento de domínio.

8.1 Limites de Dirichlet

Se uma máscara de Dirichlet for fornecida, a cabeça é prescrita diretamente:

h=hDh = h_{\mathrm{D}}

nas células de limite selecionadas, onde:

  • hDh_{\mathrm{D}} = limite prescrito [L][\mathrm{L}]

8.2 Limites sem fluxo

Os limites sem fluxo são impostos por condições de gradiente zero no perímetro do domínio. Na prática, HydroPol2D copia a cabeça do vizinho interior válido mais próximo para a célula do perímetro, aplicando assim uma condição discreta de Neumann.

8.3 Máscara de captação

Todos os cálculos de águas subterrâneas estão restritos ao domínio computacional válido. As células fora da máscara de captação são definidas como NaN.

9. Exfiltração de águas subterrâneas para a superfície

Depois de atualizar a altura do lençol freático, HydroPol2D calcula a exfiltração onde o lençol freático se eleva acima da superfície efetiva do terreno. A elevação da superfície usada para exfiltração é

hsurf=z0+Dsoilh_{\mathrm{surf}} = z_0 + D_{\mathrm{soil}}

onde:

  • hsurfh_{\mathrm{surf}} = elevação efetiva da superfície terrestre para emergência de água subterrânea [L][\mathrm{L}]

A exfiltração é então calculada como

qexf=max(0,hhsurf)SyΔtq_{\mathrm{exf}} = \max\left(0,\,h - h_{\mathrm{surf}}\right) \frac{S_y}{\Delta t}

onde:

  • qexfq_{\mathrm{exf}} = fluxo de exfiltração [LT1][\mathrm{L}\,\mathrm{T}^{-1}]

Esta expressão significa que se a altura da água subterrânea exceder a superfície da terra, o excesso de armazenamento saturado é expelido para cima e transferido para o sistema de águas superficiais. A cabeceira da água subterrânea é então cortada de modo que

hhsurfh\le h_{\mathrm{surf}}

o que evita que as águas subterrâneas permaneçam acima da superfície da terra após a exfiltração ter sido contabilizada.

10. Intercâmbio Rio-Aquífero

A interface do solucionador atual inclui:

  • uma máscara de rio,
  • um termo de condutância do rio,
  • estágio do rio,
  • elevação do leito do rio,

que são projetados para apoiar a troca rio-aquífero. Na implementação atual, a variável de espaço reservado

qriverq_{\mathrm{river}}

é inicializado e transportado pela interface de acoplamento como o termo de troca de rio. Isto significa que a interface está presente, mas a troca rio-aquífero deve atualmente ser descrita como um espaço reservado/gancho estrutural, e não uma reivindicação de teoria autônoma madura no mesmo nível que recarga, fluxo lateral de águas subterrâneas, ascensão capilar ou exfiltração.

11. Correção do excesso de saturação após elevação do lençol freático

Após a atualização do Boussinesq, o HydroPol2D recalcula a profundidade do lençol freático e, portanto, o novo armazenamento máximo não saturado:

SUZ,maxnew=zwtnew(θsatθr)S_{\mathrm{UZ,max}}^{\,\mathrm{new}} = z_{\mathrm{wt}}^{\,\mathrm{new}} \left( \theta_{\mathrm{sat}} -\theta_r \right)

Se o armazenamento vadoso atualizado exceder esta nova capacidade, o excesso é devolvido à superfície como excesso de água de saturação:

Sexcess=max(SSUZ,maxnew,0)S_{\mathrm{excess}} = \max\left( S - S_{\mathrm{UZ,max}}^{\,\mathrm{new}}, 0 \right)

A correção é então:

SSSexcessS \leftarrow S - S_{\mathrm{excess}}

e a profundidade da água superficial aumenta em

hsurfnew=hsurfold+Sexcessh_{\mathrm{surf}}^{\,\mathrm{new}} = h_{\mathrm{surf}}^{\,\mathrm{old}}+ S_{\mathrm{excess}}

Na implementação padrão em camadas, a mesma correção é aplicada após recalcular a capacidade somada admissível da camada vadose. Esta é uma parte muito importante do HydroPol2D porque garante a consistência de massa entre as zonas vadosa e saturada quando o lençol freático sobe durante um intervalo de tempo.

12. Acoplamento de volta às águas superficiais

O modelo de águas subterrâneas retroalimenta o sistema de águas superficiais de duas maneiras:

  1. exfiltração, através de qexfq_{\mathrm{exf}},
  2. excesso de saturação, através de SexcessS_{\mathrm{excess}}.

A profundidade da água superficial é atualizada conforme

dt=dt+Δtqexfd_t = d_t + \Delta t\,q_{\mathrm{exf}}

com a conversão de unidade apropriada no código de m\mathrm{m} para mm\mathrm{mm}. Isto significa que HydroPol2D permite que as águas subterrâneas influenciem ativamente a geração de inundações e a humidade próxima da superfície, em vez de tratá-las como um processo desconectado de limites inferiores.

13. Balanço de Massa

O solucionador de águas subterrâneas realiza uma verificação explícita do balanço de massa durante cada atualização. O volume de entrada de recarga é

Vrecharge=\somaRΔtAcellV_{\mathrm{recharge}} = \soma R\,\Delta t\,A_{\mathrm{cell}}

O volume de perda por exfiltração é

Vexf=\somaqexfΔtAcellV_{\mathrm{exf}} = \soma q_{\mathrm{exf}}\,\Delta t\,A_{\mathrm{cell}}

A mudança no armazenamento de água subterrânea é

ΔVstorage=\somaSymax(ht+Δtz0,0)Acell\somaSymax(htz0,0)Acell\Delta V_{\mathrm{storage}} = \soma S_y\,\max(h^{t+\Delta t}-z_0,0)\,A_{\mathrm{cell}}- \soma S_y\,\max(h^t-z_0,0)\,A_{\mathrm{cell}}

O resíduo do balanço de massa é então

εgw=VrechargeVexfΔVstorage\varepsilon_{\mathrm{gw}} = V_{\mathrm{recharge}}- V_{\mathrm{exf}}- \Delta V_{\mathrm{storage}}

onde:

  • εgw\varepsilon_{\mathrm{gw}} = erro de balanço de massa da água subterrânea [L3][\mathrm{L}^3]

Esta quantidade é armazenada como erro de diagnóstico de águas subterrâneas em HydroPol2D. Quando a água subterrânea assíncrona está ativa, a recarga acumulada e a troca capilar pendentes não são tratadas como perda de modelo. Eles são armazenados como troca de água subterrânea em trânsito até que a próxima atualização saturada consuma o volume do escalonador. Esta distinção é importante para interpretar o diagnóstico do balanço hídrico de todo o modelo.

14. Resumo

O modelo de águas subterrâneas HydroPol2D combina:

  • um módulo de recarga de zona vadose em camadas baseado em um fechamento Darcy – van Genuchten – Mualem,
  • simples ascensão capilar das águas subterrâneas rasas para o armazenamento vadoso disponível,
  • um programador assíncrono de águas subterrâneas para atualizações mais lentas de carga saturada,
  • um modelo 2D Boussinesq de fluxo de águas subterrâneas para fluxo lateral saturado,
  • passo de tempo explícito adaptativo,
  • exfiltração superficial quando a água subterrânea sobe acima da superfície do solo,
  • e uma correção do excesso de saturação que preserva a massa quando o lençol freático reduz a capacidade de armazenamento vadoso.

Esta estrutura permite que HydroPol2D represente um sistema de subsuperfície acoplado dinamicamente no qual:

  • a infiltração não desaparece simplesmente em um balde estático,
  • recarga afeta a altura das águas subterrâneas,
  • o fluxo de água subterrânea redistribui a água lateralmente,
  • e as águas subterrâneas podem retornar à superfície terrestre e influenciar a geração de fluxo superficial.

O resultado é um modelo de águas subterrâneas fisicamente interpretável, conservador em massa e totalmente integrado à arquitetura hidrológica e hidrodinâmica do HydroPol2D.