はじめに:シリーズ最大の統合
#2 から #11 まで、私たちは「ある場所・ある瞬間」の化学を計算してきた。
今回の TRANSPORT ブロックはこれに空間の次元を加える。 地下水が帯水層の中を流れながら、鉱物と反応し、汚染物質を運び、濃度前線を形成する—— その全プロセスを一度に計算する。
↓
TRANSPORT:空間 × 時間 × 化学反応
- 移流(advection)・分散(dispersion)・反応(reaction)の3つの概念
TRANSPORTブロックの基本構文- シナリオ1:塩水フロントの移動(純粋な移流分散)
- シナリオ2:カルサイト溶解前線(移流+化学反応)
- シナリオ3:重金属汚染の拡散と土壌吸着
- Python による濃度プロファイルの可視化
理論:移流分散方程式
地下水中の溶質移動は 移流分散方程式(ADE: Advection-Dispersion Equation) で記述される:
\[\frac{\partial C}{\partial t} = D \frac{\partial^2 C}{\partial x^2} - v \frac{\partial C}{\partial x} + R\]
| 項 | 名称 | 意味 | PHREEQCパラメータ |
|---|---|---|---|
| ∂C/∂t | 時間変化 | 濃度の時間的変化率 | -time_step |
| D ∂²C/∂x² | 分散項 | 濃度勾配による拡散・分散 | -dispersivity |
| v ∂C/∂x | 移流項 | 地下水の流れによる物質輸送 | -shifts(流速×時間) |
| R | 反応項 | 鉱物溶解・沈殿・吸着など | EQUILIBRIUM_PHASES等 |
PHREEQCのアプローチ:オペレータスプリッティング
PHREEQCは ADE を直接解くのではなく、移流・分散・反応を交互に計算する手法を取る:
TRANSPORT ブロックの基本構文
TRANSPORT
-cells 20 # セル数(空間方向の分割数)
-length 0.1 # 各セルの長さ(m)
-shifts 40 # タイムステップ数(移流ステップ数)
-time_step 864000 # 1ステップの時間(秒)= 10日
-flow_direction forward # 流れの方向(forward/back/diffusion_only)
-boundary_conditions flux flux # 境界条件(flux/constant)
-dispersivity 0.015 # 分散長(m)
-diffusion_coefficient 1e-9 # 拡散係数(m²/s)
-punch_cells 1-20 # 出力するセル番号
-punch_frequency 10 # 何ステップごとに出力するか
-print_cells 1 5 10 20 # 画面表示するセル
セル数を増やすと精度が上がるが計算時間も増える。 ペクレ数(Pe = v·Δx/D) が 2 以下になるようにセル長を決めよう。 Pe > 2 では数値分散が生じて物理的な分散を過大評価する。
シナリオ1:塩水フロントの移動(純粋な移流分散)
設定
純水に満たされている帯水層(20セル、各0.5m = 総延長10m)に、 NaCl 溶液(0.1 mol/L)を注入する。 塩水フロントがどのように移動・拡散するかを追う。
# ============================================================
# DeepFlow #12 - シナリオ1:塩水フロントの移動
# 純水の帯水層に NaCl 溶液を注入
# ============================================================
# ---- 初期溶液:純水(全20セルを満たす) ----
SOLUTION 0 "注入水(NaCl 0.1 mol/L)"
temp 25
pH 7
units mol/kgw
Na 0.1
Cl 0.1
SOLUTION 1-20 "初期帯水層(純水)"
temp 25
pH 7
-water 1
# ---- TRANSPORT 設定 ----
TRANSPORT
-cells 20
-lengths 20*0.5 # 0.5 m を 20セル分繰り返す(※掛け算ではなく、0.5を20回並べるPHREEQC特有の記法)
-shifts 40 # 40 タイムステップ
-time_step 86400 # 1 日/ステップ 単純計算(見かけの移動速度)で0.5 m / 86400 s = 5.8 * 10^-6 m/s
-flow_direction forward
-boundary_conditions flux flux
-dispersivities 20*0.05 # 分散長 0.05 m を 20セル分繰り返す(同上)
-diffusion_coefficient 1e-9
-punch_cells 1-20
-punch_frequency 10 # 10・20・30・40 ステップで出力
# ---- 出力設定 ----
SELECTED_OUTPUT 1
-file saltfront.txt
-distance true
-time true
-totals Na Cl
USER_PUNCH 1
-headings dist_m time_d Na_mM Cl_mM
-start
10 PUNCH DIST, TIME/86400, TOT("Na")*1000, TOT("Cl")*1000
-end
# ---- グラフ描画設定 ----
USER_GRAPH 1
-headings Distance Na_mM Cl_mM
-chart_title "Salt Front Movement"
-axis_titles "Distance (m)" "Concentration (mM)"
-axis_scale x_axis 0 10
-axis_scale y_axis 0 120
-initial_solutions false
-start
10 GRAPH_X DIST
20 GRAPH_Y TOT("Na")*1000, TOT("Cl")*1000
-end
END
シナリオ2:カルサイト溶解前線(移流+化学反応)
設定
カルサイトを含む帯水層(初期:カルサイトと平衡した Ca–HCO₃ 水)に 酸性水(pH 4.5) を注入する。 溶解前線がどのように移動するか、pH・Ca²⁺・SI(Calcite) を追う。簡単に石灰岩地下水に酸性の水が入ったらどうなるか?
# ============================================================
# シナリオ2:カルサイト溶解前線
# 酸性水を石灰岩帯水層に注入
# ============================================================
# ---- 注入水(酸性+CO2) ----
SOLUTION 0 "注入酸性水"
temp 12
pH 4.5
pe 12
units mol/kgw
-water 1
EQUILIBRIUM_PHASES 0
CO2(g) -2.0 10 # 土壌CO2
# ---- 初期地下水(カルサイト+CO2平衡) ----
SOLUTION 1-20 "初期地下水"
temp 12
pH 7.2
units mol/kgw
Ca 2.0e-3
Alkalinity 4.0e-3 as HCO3
-water 1
EQUILIBRIUM_PHASES 1-20
Calcite 0.0 0.005
CO2(g) -2.0 10
# ---- Advection-Dispersion ----
TRANSPORT
-cells 20
-lengths 20*0.5 # 0.5 m を 20セル分繰り返す(※掛け算ではなく、PHREEQC特有の配列記法)
-shifts 60
-time_step 86400
-flow_direction forward
-boundary_conditions flux flux
-dispersivities 20*0.05 # 分散長 0.05 m を 20セル分繰り返す
-diffusion_coefficient 1e-9
-punch_cells 1-20
-punch_frequency 20
# ---- Output ----
SELECTED_OUTPUT 1
-file calcite_front.txt
-distance true
-time true
-pH true
-totals Ca C(4)
-saturation_indices Calcite
-equilibrium_phases Calcite
USER_PUNCH 1
-headings dist_m time_d pH Ca_mM SI_Calcite Calcite_mol
-start
10 PUNCH DIST, TIME/86400, -LA("H+"), TOT("Ca")*1000, SI("Calcite"), EQUI("Calcite")
-end
# ---- グラフ描画設定 (1: 水質) ----
USER_GRAPH 1
-headings Distance pH Ca_mM
-chart_title "Calcite Dissolution Front (Water Chemistry)"
-axis_titles "Distance (m)" "pH" "Ca (mM)"
-axis_scale x_axis 0 10
-initial_solutions false
-start
10 GRAPH_X DIST
20 GRAPH_Y -LA("H+")
30 GRAPH_SY TOT("Ca")*1000
-end
# ---- グラフ描画設定 (2: 鉱物) ----
USER_GRAPH 2
-headings Distance SI_Calcite Calcite_mol
-chart_title "Calcite Dissolution Front (Mineral)"
-axis_titles "Distance (m)" "SI (Calcite)" "Calcite (mol)"
-axis_scale x_axis 0 10
-initial_solutions false
-start
10 GRAPH_X DIST
20 GRAPH_Y SI("Calcite")
30 GRAPH_SY EQUI("Calcite")
-end
END
結果の読み方
| 観察されること | 地球化学的意味 |
|---|---|
| pH 前線が下流へ移動 | 酸性水の移流による pH 低下域の伝播。分散によって前線がなだらかになる |
| Ca²⁺ が前線前方で上昇 | カルサイトの溶解によって Ca²⁺ が供給される。平衡条件を満たす量だけ溶解 |
| SI(Calcite) = 0 ゾーンが存在 | 酸が消費されカルサイトと再び平衡になった帯域。ここより下流は元の地下水組成を保つ |
| カルサイト量が減少するゾーン | 溶解によって鉱物が消耗されたセル。長期的には空隙率が増加して透水係数が変化する |
シナリオ3:重金属(Cd²⁺)汚染の移動と土壌吸着
設定
鉱山廃水由来の Cd²⁺ 汚染水が帯水層に侵入した場合を想定する。 土壌の陽イオン交換(ion exchange) による吸着遅延効果を組み込む。
# ============================================================
# シナリオ3:重金属汚染の移動と土壌吸着
# Cd²⁺ の移流 + 陽イオン交換による遅延
# ============================================================
# ---- 注入水:Cd²⁺ 汚染水 ----
SOLUTION 0 "Cd汚染水"
temp 15
pH 6.5
units mol/kgw
Na 1.0e-3
Cl 1.0e-3
Cd 1.0e-5 # 10 μmol/L の Cd²⁺
-water 1
# ---- 初期帯水層水:清浄地下水 ----
SOLUTION 1-20
temp 15
pH 7.2
units mol/kgw
Na 1.5e-3
Ca 0.5e-3
Cl 1.0e-3
Alkalinity 2.0e-3 as HCO3
-water 1
# ---- 土壌の陽イオン交換体(全セル) ----
# phreeqc.dat の交換体定義を使用:
# NaX, KX, CaX2, MgX2, CdX2 など
# X- は phreeqc.dat では定義されないため -equilibrate で初期化
EXCHANGE 1-20
NaX 0.0018 # Na が初期条件で交換サイトを占める
CaX2 0.0001 # Ca も微量占有
-equilibrate 1 # SOLUTION 1 と平衡化して初期状態を決定
# ---- TRANSPORT ----
TRANSPORT
-cells 20
-lengths 20*0.5 # 0.5 m を 20セル分繰り返す(※掛け算ではなく、PHREEQC特有の配列記法)
-shifts 80
-time_step 86400
-flow_direction forward
-boundary_conditions flux flux
-dispersivities 20*0.05 # 分散長 0.05 m を 20セル分繰り返す
-diffusion_coefficient 1e-9
-punch_cells 1-20
-punch_frequency 20
SELECTED_OUTPUT 3
-file cd_transport.txt
-distance true
-time true
-totals Cd Na Ca
-molalities Cd+2
USER_PUNCH 3
-headings dist_m time_d Cd_umol Cd2_activity Na_mM
-start
10 PUNCH DIST, TIME/86400, TOT("Cd")*1e6, ACT("Cd+2"), TOT("Na")*1000
-end
USER_GRAPH 1
-headings Distance Cd_umol Na_mM
-chart_title "Heavy Metal Transport with Ion Exchange"
-axis_titles "Distance (m)" "Cd (umol/L)" "Na (mM)"
-axis_scale x_axis 0 10
-initial_solutions false
-start
10 GRAPH_X DIST
20 GRAPH_Y TOT("Cd")*1e6
30 GRAPH_SY TOT("Na")*1000
-end
END
吸着がある場合、汚染物質の移動速度は地下水の流速より遅くなる。 この比率が遅延係数 R であり、次式で表される:
\[R = 1 + \frac{\rho_b \cdot K_d}{\theta}\]
\(\rho_b\):嵩密度、\(K_d\):分配係数、\(\theta\):空隙率。
物理的な意味
遅延係数 R は
地下水の流速 ÷ 汚染物質の移動速度を意味する。例えば:
R=1:水と同じ速度(吸着なし)
R=5:水の1/5の速度
R=10:水の1/10の速度
Cd²⁺ のような重金属は土壌に強く吸着するため R が大きく、 地下水の1/5〜1/10の速度でしか移動しないことが多い。 PHREEQC の EXCHANGE ブロックはこの効果を熱力学的に自動計算する。
Python で濃度プロファイルを可視化する
# ============================================================
# transport_plot.py
# 3シナリオの空間濃度プロファイルを可視化
# ============================================================
import numpy as np
import os
import matplotlib
import matplotlib.font_manager as fm
import matplotlib.pyplot as plt
import matplotlib.cm as cm
# ---- フォント設定(日本語環境) ----
plt.rcParams.update({
"font.family": "sans-serif",
"font.sans-serif": ["MS Gothic", "Noto Sans CJK JP", "IPAexGothic", "DejaVu Sans"],
"axes.unicode_minus": False,
"figure.dpi": 150,
})
# 数値計算コア:PHREEQCのオペレータスプリッティングを再現
# ① 移流:セルを1つ右にシフト
# ② 分散:隣接セル間のミキシング(混合率 = α/(dx+α))
def adv_disp(nx, dx, alpha, n_shifts, C0, C_init,
punch_freq, R=1):
"""
1次元移流分散シミュレータ(オペレータスプリッティング)
Parameters
----------
nx : セル数
dx : セル長(m)
alpha : 分散長(m)
n_shifts : タイムステップ数(シフト数)
C0 : 注入水濃度
C_init : 初期帯水層濃度
punch_freq : 何シフトごとに記録するか
R : 遅延係数(吸着あり=R>1)
Returns
-------
list of (shift, concentration_array)
"""
C = np.full(nx, float(C_init))
f_mix = alpha / (dx + alpha) # 分散混合率
results = []
for s in range(1, n_shifts + 1):
# ① 移流(遅延係数 R を考慮:R=5 なら5ステップに1回シフト)
if s % R == 0:
C_adv = np.empty(nx)
C_adv[0] = C0 # 入口境界:注入水濃度
C_adv[1:] = C[:-1] # 1セル右へシフト
else:
C_adv = C.copy() # シフトなし(吸着遅延)
# ② 分散(隣接セルとのミキシング)
C_new = C_adv.copy()
for i in range(nx):
left = C_adv[i - 1] if i > 0 else C0 # 左境界
right = C_adv[i + 1] if i < nx - 1 else C_adv[i] # 右境界(flux BC)
C_new[i] = C_adv[i] + f_mix * (left - 2 * C_adv[i] + right) / 2
C = C_new.clip(min(C0, C_init), max(C0, C_init))
if s % punch_freq == 0:
results.append((s, C.copy()))
return results
def ph_front(nx, dx, alpha, n_shifts, pH_in, pH_init,
punch_freq, calcite_init=0.005, acid_in=0.002):
acid = np.zeros(nx) # セル内の酸の濃度 (mol/L)
calcite = np.full(nx, float(calcite_init)) # カルサイト残量 (mol/L)
f_mix = alpha / (dx + alpha)
results = []
for s in range(1, n_shifts + 1):
# ① 移流(酸の移動)
acid_adv = np.empty(nx)
acid_adv[0] = acid_in # 注入水は一定の酸を含む
acid_adv[1:] = acid[:-1]
# ② 分散(酸の混合)
acid_d = acid_adv.copy()
for i in range(nx):
L = acid_adv[i - 1] if i > 0 else acid_in
R = acid_adv[i + 1] if i < nx - 1 else acid_adv[i]
acid_d[i] = acid_adv[i] + f_mix * (L - 2 * acid_adv[i] + R) / 2
# ③ 化学反応(カルサイトによる酸の中和と枯渇)
acid_new = acid_d.copy()
for i in range(nx):
if acid_new[i] > 0:
# 中和できる量は、存在する酸とカルサイト残量の少ない方
reacted = min(acid_new[i], calcite[i])
acid_new[i] -= reacted
calcite[i] -= reacted
acid = acid_new.copy()
# ④ pH の計算
pH_out = np.empty(nx)
for i in range(nx):
if calcite[i] > 0:
# カルサイトが残っていれば緩衝されて初期pHを維持
pH_out[i] = pH_init
else:
# 枯渇ゾーンでは酸の濃度に応じてpHが低下
ratio = min(acid[i] / acid_in, 1.0)
pH_out[i] = pH_init - ratio * (pH_init - pH_in)
if s % punch_freq == 0:
results.append((s, pH_out.copy()))
return results
# 各シナリオのパラメータ設定
x = np.arange(0.5, 20.5) * 0.5 # セル中心座標(m)
# シナリオ1:塩水フロント(NaCl 0.1 mol/L = Na 100 mM)
res1 = adv_disp(
nx=20, dx=0.5, alpha=0.5, n_shifts=40,
C0=100.0, C_init=0.0, punch_freq=10
)
# シナリオ2:カルサイト溶解前線(pH)
res2 = ph_front(
nx=20, dx=0.5, alpha=0.3, n_shifts=60,
pH_in=4.5, pH_init=7.4, punch_freq=20,
calcite_init=0.005, acid_in=0.002
)
# シナリオ3:Cd²⁺ 汚染(保守的トレーサー + 吸着遅延 R=5)
res3_cons = adv_disp(
nx=20, dx=0.5, alpha=0.3, n_shifts=80,
C0=10.0, C_init=0.0, punch_freq=20, R=1
)
res3_cd = adv_disp(
nx=20, dx=0.5, alpha=0.3, n_shifts=80,
C0=10.0, C_init=0.0, punch_freq=20, R=5 # R=5 で遅延
)
# 作図
fig, axes = plt.subplots(1, 3, figsize=(15, 5))
fig.patch.set_facecolor("white")
# ---- シナリオ1:塩水フロント ----
ax = axes[0]
blues = cm.Blues(np.linspace(0.35, 1.0, len(res1)))
for (step, C), col in zip(res1, blues):
ax.plot(x, C, color=col, lw=2, label=f"{step} 日")
ax.set(
xlabel="距離 (m)",
ylabel="Na (mM)",
title="シナリオ1:塩水フロント移動\nNaCl 0.1 mol/L 注入",
)
ax.legend(title="経過日数", fontsize=8.5, loc="upper left",
framealpha=0.9, edgecolor="#E5E7EB")
ax.set_ylim(-2, 108)
ax.grid(True, ls="--", lw=0.5, color="#E5E7EB")
ax.set_facecolor("white")
# ---- シナリオ2:カルサイト溶解前線 ----
ax = axes[1]
oranges = cm.Oranges(np.linspace(0.35, 1.0, len(res2)))
for (step, pH), col in zip(res2, oranges):
ax.plot(x, pH, color=col, lw=2, label=f"{step} 日")
ax.axhline(7.4, color="#9CA3AF", lw=1, ls=":", alpha=0.8)
ax.text(0.2, 7.55, "初期 pH 7.4", color="#9CA3AF", fontsize=8.5)
ax.axhline(4.5, color="#93C5FD", lw=1, ls=":", alpha=0.8)
ax.text(0.2, 4.65, "注入水 pH 4.5", color="#93C5FD", fontsize=8.5)
ax.set(
xlabel="距離 (m)",
ylabel="pH",
title="シナリオ2:カルサイト溶解前線\n酸性水 pH 4.5 注入",
)
ax.legend(title="経過日数", fontsize=8.5, loc="lower right",
framealpha=0.9, edgecolor="#E5E7EB")
ax.set_ylim(3.8, 7.9)
ax.grid(True, ls="--", lw=0.5, color="#E5E7EB")
ax.set_facecolor("white")
# ---- シナリオ3:Cd 汚染と吸着遅延 ----
ax = axes[2]
grays = cm.Greys(np.linspace(0.35, 0.75, len(res3_cons)))
reds = cm.Reds(np.linspace(0.35, 1.0, len(res3_cd)))
for i, ((step, Cc), (_, Cd)) in enumerate(zip(res3_cons, res3_cd)):
ax.plot(x, Cc, color=grays[i], lw=1.5, ls="--")
ax.plot(x, Cd, color=reds[i], lw=2)
ax.plot([], [], color="gray", lw=1.5, ls="--", label="保守的トレーサー (R=1)")
ax.plot([], [], color="#DC2626", lw=2, label="Cd (R=5, 吸着遅延)")
ax.set(
xlabel="距離 (m)",
ylabel="濃度 (umol/L)",
title="シナリオ3:Cd 汚染移動\n陽イオン交換による遅延",
)
ax.legend(fontsize=8.5, loc="upper right",
framealpha=0.9, edgecolor="#E5E7EB")
ax.grid(True, ls="--", lw=0.5, color="#E5E7EB")
ax.set_facecolor("white")
fig.suptitle("PHREEQC TRANSPORT:3シナリオの濃度プロファイル",
fontsize=13, y=0.99)
plt.tight_layout()
plt.savefig("transport_profiles.png", dpi=150,
bbox_inches="tight", facecolor="white")
plt.show()
print("transport_profiles.png を保存しました")シナリオ1, 2, 3の結果比較(Pythonコードの結果)

上図は、各シナリオにおける物質の濃度プロファイル(前線の形状と移動)を比較したものです。単なる移流・分散に加え、化学反応や吸着が作用することで、物質の移動メカニズムがどのように変化するかを視覚的に確認できます。
1. 塩水フロント(純粋な移流・分散)
左図(シナリオ1)は、保守的(非反応性)な溶質である NaCl の移動を示しています。 前線の濃度プロファイルがなだらかな 「S字曲線」 を描いているのが特徴です。地下水流による「移流」が前線の中心位置を下流へ押し流しつつ、「分散」が時間経過とともに濃度勾配を均して前線の幅を広げていくという、最も基本的な物質移動の挙動が現れています。
2. pH 前線(移流・分散 + 溶解・緩衝反応)
中央図(シナリオ2)は、酸性水の侵入に伴う pH 前線です。シナリオ1のようななだらかなS字曲線ではなく、極めてシャープな 「垂直の段差(ステップ関数)」 になっていることが分かります。 これは、侵入した酸が帯水層内のカルサイトを溶解し尽くすまでは pH が初期値(7.4)に完全に緩衝・維持され、カルサイトが枯渇した瞬間に pH が注入水の値(4.5)まで急激に低下するためです。化学反応によって前線が鋭く変形する典型的な反応前線の例です。
3. Cd²⁺ 汚染(移流・分散 + 吸着遅延)
右図(シナリオ3)は、土壌への陽イオン交換(吸着)を伴う重金属 Cd²⁺ の移動です。 水と一緒に移動する保守的トレーサー(破線、R=1)と比較すると、Cd²⁺ のプロファイル(実線)は前線の形状自体はS字曲線で似ているものの、到達位置が大きく手前に留まっています。これが 「遅延(Retardation)」 であり、吸着反応が汚染物質の移動速度を物理的に遅らせる効果を明確に示しています。
TRANSPORT の重要パラメータ早見表
| パラメータ | 単位 | 意味 | 典型値 |
|---|---|---|---|
| -cells | — | 空間分割数。多いほど精度↑・計算時間↑ | 10〜100 |
| -length | m | 各セルの長さ | 0.01〜10 |
| -shifts | — | タイムステップ数(= セルを流体が通過する回数) | 20〜200 |
| -time_step | 秒 | 1ステップの時間。length/流速 で計算 | 86400(1日) |
| -dispersivity | m | 分散長α。移動距離の約1/10が経験則 | 0.01〜1 |
| -stagnant | — | デュアルポロシティモデル(移動層/非移動層) | 亀裂岩盤で重要 |
まとめ:シリーズ全体の完成
#1 のインストールから始まり、Speciation・混合・平衡・酸化還元・溶解度・活量・飽和指数・反応経路、そして今回の移流分散まで——PHREEQCで地球化学シミュレーションを実践するために必要な知識の核心を一通り歩いた。
次のステップとして、実際の野外データ(水質分析値)を使った逆モデリング(INVERSE_MODELING ブロック)や、FloPy × PHREEQC の結合(シリーズ冒頭の熱水シミュレーション)へ進む道が開けている。
参考文献(References)
- #1 インストールと最初の計算
- #2 Speciationで海水を解析する
- #3 MixingとEQUILIBRIUM_PHASES
- #4 カルサイト−CO₂水反応
- #5 炭酸地下水と海水の混合
- #6 黄鉄鉱の酸化(AMD)
- #7 溶解度ダイアグラム(Gibbsite)
- #8 Pythonでの可視化
- #9 イオン強度と活量係数
- #10 飽和指数(SI)の使いこなし
- #11 反応経路モデリング(REACTION block の応用)
- #12 移流分散モデル(TRANSPORT block の応用) ← 今回
- #13 酸化還元シークエンス — 地下水が「還元」されていく順番
DeepFlow | 地球科学シミュレーションの深みへ