地下水流動方程式

水文地質学で使用される地下水流動方程式は、帯水層を通る地下水の流れを記述するために使用される数学的関係です。地下水の過渡的流れは、固体内の熱の流れ (熱伝導) を記述する熱伝達で使用される拡散方程式の形式で記述されます。地下水の定常状態の流れは、ポテンシャルフローの形式であるラプラス方程式の形式で記述され、多くの分野で類似しています。

地下水流動方程式は、多くの場合、媒体の特性が実質的に一定であると仮定した小さな代表的要素体積（REV）に対して導出されます。この小さな体積に出入りする水に対して質量バランスが行われ、関係のフラックス項は、層流であることを要求するダルシーの法則と呼ばれる構成方程式を使用して、水頭で表されます。他のアプローチは、カルストや破砕岩（火山性）などの複雑な帯水層の影響を組み込む エージェントベースモデルに基づいています^[1]

過渡的地下水流動方程式を導くには、質量バランスを計算し、ダルシーの法則と併用する必要があります。このバランスは、熱方程式を導くために熱伝達で使用されるエネルギーバランスに似ています。これは、与えられた制御体積に対して、発生源やシンクは別として、質量は生成も破壊もできないという単純な計算です。質量保存則によれば、与えられた時間増分 ( Δt ) において、境界を越えて流入する質量、境界を越えて流出する質量、および体積内の発生源との差が、貯蔵量の変化となります。

${\displaystyle {\frac {\Delta M_{stor}}{\Delta t}}={\frac {M_{in}}{\Delta t}}-{\frac {M_{out}}{\Delta t}}-{\frac {M_{gen}}{\Delta t}}}$

質量は密度×体積として表すことができ、ほとんどの状況下では水は非圧縮性であると考えられます（密度は圧力に依存しません）。境界を越える質量フラックスは体積フラックスになります（ダルシーの法則に見られるように）。テイラー級数を使用して制御体積の境界を越える流入フラックス項と流出フラックス項を表し、発散定理を使用して境界を越えるフラックスを体積全体のフラックスに変換すると、地下水流方程式の最終的な形式（微分形式）は次のようになります。

${\displaystyle S_{s}{\frac {\partial h}{\partial t}}=-\nabla \cdot qG.}$

これは他の分野では拡散方程式または熱方程式として知られており、放物型偏微分方程式（PDE）です。この数学的記述は、水頭の時間変化（左辺）がフラックス（q）とソース項（G ）の負の発散に等しいことを示しています。この方程式には水頭とフラックスの両方が未知数ですが、ダルシーの法則はフラックスを水頭に関連付けているため、フラックス（q）に代入すると次のようになります。

${\displaystyle S_{s}{\frac {\partial h}{\partial t}}=-\nabla \cdot (-K\nabla h)-G.}$

ここで、透水係数（K）が空間的に均一かつ等方性（テンソルではなく）である場合、空間微分からそれを取り出してラプラシアンに簡略化することができ、これにより次の式が得られる。

${\displaystyle S_{s}{\frac {\partial h}{\partial t}}=K\nabla ^{2}hG.}$

これを比貯留量（S _s）で割ると、右辺に水圧拡散係数（ α = K/S _sまたは同等のα = T/S ）が置かれます。水圧拡散係数は、有限の圧力パルスがシステム内を伝播する速度に比例します（ αの値が大きいほど、信号の伝播が速くなります）。地下水流動方程式は次のようになります。

${\displaystyle {\frac {\partial h}{\partial t}}=\alpha \nabla ^{2}hG.}$

ここで、シンク/ソース項Gは同じ単位を持ちますが、適切なストレージ項 (水力拡散係数の置換によって定義される) で除算されます。

特に長方形格子差分モデル（例えば、 USGSが作成したMODFLOW ）を使用する場合、直交座標を扱う。これらの座標では、一般的なラプラシアン演算子は（3次元の流れの場合）具体的には

${\displaystyle {\frac {\partial h}{\partial t}}=\alpha \left[{\frac {\partial ^{2}h}{\partial x^{2}}}+{\frac {\partial ^{2}h}{\partial y^{2}}}+{\frac {\partial ^{2}h}{\partial z^{2}}}\right]-G.}$

MODFLOW コードは、支配的な地下水流方程式の直交3D 形式を離散化してシミュレートします。ただし、ユーザーが希望する場合は、「準 3D」モードで実行するオプションがあります。この場合、モデルはkとS _sではなく、垂直平均TとSを処理します。準 3D モードでは、漏れの概念を使用して 2D 水平層間のフローが計算されます。

もう一つの有用な座標系は3D円筒座標系です（典型的には揚水井が原点に位置する線状源で、Z軸に平行で、収束する放射状の流れを引き起こします）。これらの条件下では、上記の式は次のようになります（rは放射状距離、θは角度）。

${\displaystyle {\frac {\partial h}{\partial t}}=\alpha \left[{\frac {\partial ^{2}h}{\partial r^{2}}}+{\frac {1}{r}}{\frac {\partial h}{\partial r}}+{\frac {1}{r^{2}}}{\frac {\partial ^{2}h}{\partial \theta ^{2}}}+{\frac {\partial ^{2}h}{\partial z^{2}}}\right]-G.}$

この方程式は、原点に位置する揚水井戸 (強度Gのシンク) への流れを表します。この方程式と上記の直交座標バージョンはどちらも地下水の流れに関する基本的な方程式ですが、この点に到達するにはかなりの単純化が必要です。これらの方程式の両方に組み込まれた主な仮定は次のとおりです。

• 帯水層の物質は非圧縮性（圧力の変化によるマトリックスの変化がない、いわゆる沈下）である。
• 水は一定の密度（非圧縮性）である。
• 帯水層にかかる外部荷重（例えば、表土、大気圧）は一定である。
• 1D放射問題では、揚水井戸は漏出のない帯水層を完全に貫通しており、
• 地下水はゆっくりと流れている（レイノルズ数は1未満）、そして
• 透水係数（K）は等方性 スカラーです。

これらの大きな仮定にもかかわらず、地下水流動方程式は、水源とシンクの一時的な分布による帯水層内の水頭の分布をうまく表します。

帯水層に再充填境界条件がある場合、定常状態に達する可能性があり（または多くの場合、近似値として使用できます）、拡散方程式（上記）はラプラス方程式に簡略化されます。

${\displaystyle 0=\alpha \nabla ^{2}h}$

この方程式は、水頭が調和関数であることを述べており、他の分野にも多くの類似点があります。ラプラス方程式は、上記と同様の仮定を使用した手法を使用して解くことができますが、定常流れ場という追加の要件があります。

土木工学や土質力学においてこの方程式を解く一般的な方法は、フローネットを描くグラフィカルな手法を使用することです。フローネットでは、水頭と流れ関数の等高線が曲線グリッドを形成し、複雑な形状を近似的に解くことができます。

揚水井戸への定常流れ（実際には決して発生しませんが、有用な近似値となる場合があります）は、一般にThiem 解と呼ばれます。

上記の地下水流方程式は、3 次元の流れに有効です。非閉じ込め帯水層では、自由表面水位境界条件の存在により、方程式の 3D 形式の解は複雑になります。つまり、水頭の空間分布を解く必要があるだけでなく、この表面の位置も未知数です。支配方程式が線形であっても、これは非線形問題です。

地下水流動方程式の別の定式化は、デュプイ・フォルシュハイマー仮定を適用することによって得られます。この仮定では、水頭は垂直方向に変化しない（つまり、）と仮定します。水平方向の水収支は、帯水層底から不飽和表面まで広がる領域を持つ長い垂直柱に適用されます。この距離は飽和厚さbと呼ばれます。被圧帯水層では、飽和厚さは帯水層の高さHによって決まり、圧力水頭はどこでもゼロではありません。非被圧帯水層では、飽和厚さは地下水面と帯水層底の垂直距離として定義されます。 で、帯水層底がゼロ基準点にある場合、非被圧飽和厚さは水頭に等しくなります（つまり、b=h ）。 ${\displaystyle \partial h/\partial z=0}$${\displaystyle \delta x\delta y}$${\displaystyle \partial h/\partial z=0}$

透水係数と流れの水平成分の両方が帯水層の飽和厚さ全体にわたって均一であると仮定すると（つまり、および）、ダルシーの法則を統合地下水流量、Q _xおよびQ _yで表すことができます。 ${\displaystyle \partial q_{x}/\partial z=0}$${\displaystyle \partial K/\partial z=0}$

${\displaystyle Q_{x}=\int _{0}^{b}q_{x}dz=-Kb{\frac {\partial h}{\partial x}}}$

${\displaystyle Q_{y}=\int _{0}^{b}q_{y}dz=-Kb{\frac {\partial h}{\partial y}}}$

これらを質量バランス式に代入すると、非圧縮飽和地下水流の一般的な 2D 支配方程式が得られます。

${\displaystyle {\frac {\partial nb}{\partial t}}=\nabla \cdot (Kb\nabla h)+N.}$

ここで、n は帯水層の多孔度です。ソース項N (時間あたりの長さ) は、垂直方向への水の追加 (たとえば、再充填) を表します。飽和厚さ、比貯留量、および比収量の正しい定義を組み込むことで、これを閉じ込められた状態と閉じ込められていない状態に対する 2 つの固有の支配方程式に変換できます。

${\displaystyle S{\frac {\partial h}{\partial t}}=\nabla \cdot (Kb\nabla h)+N.}$

（閉じ込め）ここで、S=S _s bは帯水層の貯留係数であり、

${\displaystyle S_{y}{\frac {\partial h}{\partial t}}=\nabla \cdot (Kh\nabla h)+N.}$

（非閉じ込め）、ここでS _{y は}帯水層の比揚水量です。

制限のないケースでは偏微分方程式は非線形ですが、制限のあるケースでは線形であることに注意してください。制限のない定常流れの場合、この非線形性は、PDE をヘッドの 2 乗で表すことによって除去できます。

${\displaystyle \nabla \cdot (K\nabla h^{2})=-2N.}$

あるいは、均質帯水層の場合、

${\displaystyle \nabla ^{2}h^{2}=-{\frac {2N}{K}}.}$

この定式化により、非閉じ込め流の場合に線形 PDE を解くための標準的な方法を適用できます。再充填のない不均質帯水層の場合、閉じ込め/非閉じ込めが混在するケースにポテンシャル フロー法を適用できます。

• 解析要素法
  • 偏微分方程式を解くために使用される数値手法
• デュピュイ＝フォルヒハイマー仮定
  • 鉛直流に関する地下水流動方程式の簡略化
• 地下水エネルギーバランス
  • エネルギーバランスに基づく地下水流動方程式
• リチャーズ方程式

• HF Wang と MP Anderson 地下水モデリング入門: 有限差分法と有限要素法
  • 地下水モデリングの初心者に最適な一冊です。FORTRAN 77の簡単な例とともに、すべての基本概念を網羅しています。
• Freeze, R. Allan; Cherry, John A. (1979). Groundwater Archived 2020-04-06 at the Wayback Machine . Prentice Hall. ISBN 978-0133653120 . 

• USGS 地下水ソフトウェア — MODFLOW のような無料の地下水モデリング ソフトウェア
• 地下水水文学 ( MIT OpenCourseware )