はじめに：「点」の計算から「旅」の計算へ
#2〜#10 では、ある瞬間の溶液状態を計算してきた。平衡計算は 一枚の写真 と言えよう。
今回の反応経路モデリング（reaction path modeling）は、その写真を何百枚も繋げて作る 映像 に相当する。
すなわち、反応の進行に伴う化学組成や鉱物相の変化を連続的に追跡し、系がどのように進化していくかを記述する手法である。
今回シミュレートするシナリオ
大気 CO₂ を溶かした 酸性雨水（pH 5.6） が地表に降り注ぎ、
土壌 CO₂ ゾーンを通過しながら CO₂ をさらに取り込み、
石灰岩（カルサイト） と接触して少しずつ溶かしながら流れ、
やがて 地下水（pH 7.4、Ca–HCO₃ 型） へと進化する。
この過程で pH・SI・Ca²⁺・HCO₃⁻ がどう変わるかを追う。
土壌 CO₂ ゾーンを通過しながら CO₂ をさらに取り込み、
石灰岩（カルサイト） と接触して少しずつ溶かしながら流れ、
やがて 地下水（pH 7.4、Ca–HCO₃ 型） へと進化する。
この過程で pH・SI・Ca²⁺・HCO₃⁻ がどう変わるかを追う。
ノートこの記事で学ぶこと
REACTIONブロックの基本構文と「少しずつ反応させる」仕組みINCREMENTAL_REACTIONSによるステップ分割の意味- 開放系（CO₂ 補給あり）と閉鎖系（CO₂ 補給なし）の違い
USER_GRAPHとSELECTED_OUTPUTで反応経路を可視化する方法- 石灰岩・花崗岩・黄鉄鉱の 3 シナリオの比較
- Python による反応経路の可視化
理論：反応経路とは何か
PHREEQC の REACTION ブロックの仕組み
REACTION ブロックは、指定した物質を溶液に 少しずつ添加 していく。 各ステップで平衡計算が走り、その都度 pH・SI・化学種が更新される。
図1. REACTION ブロックによる段階的平衡計算の流れ
開放系 vs 閉鎖系
🌬️ 開放系（CO₂ 補給あり）
大気・土壌 CO₂ と常に平衡を保ちながら溶解が進む。
→ pH の上昇が抑えられ、カルサイトがより多く溶ける。
→ 典型例：土壌帯水層の上部
EQUILIBRIUM_PHASES にCO2(g) を追加する。→ pH の上昇が抑えられ、カルサイトがより多く溶ける。
→ 典型例：土壌帯水層の上部
🔒 閉鎖系（CO₂ 補給なし）
溶液内の CO₂ だけで反応が進む。
CO₂ が消費されると H⁺ の供給が止まり
溶解が抑制される。
→ pH が高くなりやすく、溶解量が少ない。
→ 典型例：深部封圧帯水層
CO₂ が消費されると H⁺ の供給が止まり
溶解が抑制される。
→ pH が高くなりやすく、溶解量が少ない。
→ 典型例：深部封圧帯水層
シナリオ 1：石灰岩の溶解（開放系）
地質学的な設定
図2. 石灰岩地帯における雨水→地下水の反応経路
PHREEQCコード（完全版・シナリオ 1）
# ============================================================
# DeepFlow #11 - 反応経路モデリング
# シナリオ1：酸性雨水による石灰岩の溶解（開放系）
# ============================================================
# ---- 初期溶液：土壌通過後の酸性水 ----
SOLUTION 1 "土壌通過雨水"
temp 12
pH 4.5
pe 12
units mol/kgw
-water 1
# ---- 大気 CO₂ との平衡を保つ（開放系） ----
EQUILIBRIUM_PHASES 1
CO2(g) -2.0 10 # log pCO2 = -2.0（土壌 CO₂ ゾーン）
# ---- カルサイトを 100 ステップで少しずつ添加 ----
REACTION 1
CaCO3 1 # 1 mol のカルサイトを...
0.0001 0.0002 0.0003 0.0005 0.0007 0.001 # ステップ量（mol）
0.0015 0.002 0.003 0.004 0.005 0.007
0.010 0.012 0.015 0.018 0.020 0.025
0.030 0.035 0.040 0.045 0.050 0.055
0.060 0.065 0.070 0.075 0.080 0.085
0.090 0.095 0.100 0.110 0.120 0.130
0.140 0.150 0.160 0.170 0.180 0.190
0.200 0.220 0.240 0.260 0.280 0.300
INCREMENTAL_REACTIONS true
# ---- 出力設定 ----
SELECTED_OUTPUT 1
-file limestone_path.txt
-reset false
-step true
-pH true
-temperature true
-totals Ca C(4) Mg
-activities Ca+2 HCO3- CO3-2 H+
-saturation_indices Calcite Aragonite Dolomite CO2(g)
-equilibrium_phases CO2(g)
USER_PUNCH 1
-headings Step pH Ca_mM HCO3_mM SI_Calcite pCO2_log
-start
10 PUNCH STEP_NO, -LA("H+"), TOT("Ca")*1000, TOT("C(4)")*1000, \
SI("Calcite"), SI("CO2(g)")
-end
USER_GRAPH 2
-chart_title "Limestone Dissolution: Reaction Path"
-axis_titles "Ca dissolved (mmol/kgw)" "pH" "SI"
-axis_scale y_axis 4 9.5
-headings Ca_mmol pH SI(Calcite)
-start
10 GRAPH_X TOT("Ca") * 1000
20 GRAPH_Y -LA("H+")
30 GRAPH_SY SI("Calcite")
-end
END
ヒント
INCREMENTAL_REACTIONS true の意味
true にすると、各ステップは「前のステップからの増分」として処理される。 false（デフォルト）では各ステップは初期溶液からの累積量として処理される。 反応経路のシミュレーションでは原則 true を使う。
シナリオ 2：花崗岩の風化（閉鎖系）
花崗岩の主要鉱物である 長石（アルバイト：NaAlSi₃O₈） の風化を追う。
# ============================================================
# シナリオ2：花崗岩の風化（閉鎖系）
# アルバイト（Na長石）の溶解
# ============================================================
SOLUTION 2 "大気平衡雨水"
temp 10
pH 5.65 # 大気 CO₂ との平衡
pe 12
units mol/kgw
-water 1
# ---- アルバイト + カルサイト微量を閉鎖系で溶かす ----
REACTION 2
NaAlSi3O8 1 # アルバイト（長石）
0.0001 0.0002 0.0005 0.001 0.002 0.003
0.005 0.007 0.010 0.012 0.015 0.020
0.025 0.030 0.040 0.050 0.060 0.080
0.100 0.120 0.150 0.200
INCREMENTAL_REACTIONS true
SELECTED_OUTPUT 2
-file granite_path.txt
-reset false
-step true
-pH true
-totals Na Al Si Ca
-activities Al+3 Na+ H4SiO4
-saturation_indices Gibbsite Kaolinite Albite Quartz K-feldspar
USER_PUNCH 2
-headings Step pH Na_mM Al_uM Si_mM SI_Gibbsite SI_Kaolinite SI_Quartz
-start
10 PUNCH STEP_NO, -LA("H+"), TOT("Na")*1000, TOT("Al")*1e6, \
TOT("Si")*1000, SI("Gibbsite"), SI("Kaolinite"), SI("Quartz")
-end
USER_GRAPH 2
-chart_title "花崗岩風化：反応経路"
-axis_titles "溶解ステップ" "log 活量" "SI"
-headings Step Al3+ H4SiO4 SI(Gibbsite) SI(Kaolinite)
-start
10 GRAPH_X STEP_NO
20 GRAPH_Y LA("Al+3"), "Al³⁺"
30 GRAPH_Y LA("H4SiO4"), "H₄SiO₄"
40 GRAPH_SY SI("Gibbsite"), "SI(Gibbsite)"
50 GRAPH_SY SI("Kaolinite"),"SI(Kaolinite)"
-end
END花崗岩風化の反応系列
花崗岩の風化では、溶解が進むにつれて 二次鉱物 が順番に沈殿する：
| 段階 | 主反応 | pH 域 | 生成する二次鉱物 |
|---|---|---|---|
| ① 初期溶解 | 長石 + H₂O + CO₂ → Na⁺ + HCO₃⁻ + Al³⁺ + H₄SiO₄ | 4〜5.5 | なし（溶液中に蓄積） |
| ② Gibbsite 沈殿 | Al³⁺ + 3H₂O → Al(OH)₃ + 3H⁺ | 5.5〜6.5 | Al(OH)₃（Gibbsite） |
| ③ カオリナイト沈殿 | Al³⁺ + Si → Al₂Si₂O₅(OH)₄ | 6.5〜7.5 | カオリナイト（粘土鉱物） |
| ④ 長石との平衡 | 溶解速度 ≈ 沈殿速度 | > 7.5 | スメクタイト・モンモリロナイト |
シナリオ 3：黄鉄鉱の酸化（開放系・酸素存在下）
#6 の AMD テーマを反応経路として追う。
# ============================================================
# シナリオ3：黄鉄鉱の段階的酸化（AMD 形成）
# ============================================================
SOLUTION 3 "中性地下水"
temp 15
pH 7.0
pe 4
units mol/kgw
Na 1.0e-3
Cl 1.0e-3
Alkalinity 1.0e-3 as HCO3
Fe 1.0e-6
S(6) 1.0e-5
# ---- 酸素存在下（大気開放） ----
EQUILIBRIUM_PHASES 3
O2(g) -0.68 10 # 大気 O₂（pO₂ = 0.21 atm, log = -0.68）
# ---- 黄鉄鉱を少しずつ酸化 ----
REACTION 3
FeS2 1
0.0001 0.0002 0.0003 0.0005 0.0008 0.001
0.002 0.003 0.005 0.007 0.010 0.015
0.020 0.025 0.030 0.035 0.040 0.050
0.060 0.080 0.100
INCREMENTAL_REACTIONS true
SELECTED_OUTPUT 3
-file pyrite_path.txt
-reset false
-step true
-pH true
-totals Fe S(6) Al
-activities Fe+2 Fe+3 SO4-2
-saturation_indices Goethite Calcite Gypsum
USER_PUNCH 3
-headings Step pH Fe_mM SO4_mM SI_Goethite
-start
10 PUNCH STEP_NO, -LA("H+"), TOT("Fe")*1000, TOT("S(6)")*1000, \
SI("Goethite")
-end
USER_GRAPH 3
-chart_title "Pyrite Oxidation: Reaction Path (AMD)"
-axis_titles "Reaction step" "pH / log activity" "SI"
-headings Step pH Fe2+ SO4_2- SI(Goethite)
-start
10 GRAPH_X STEP_NO
20 GRAPH_Y -LA("H+")
30 GRAPH_Y LA("Fe+2")
40 GRAPH_Y LA("SO4-2")
50 GRAPH_SY SI("Goethite")
-end
ENDPython で 3 シナリオを比較可視化する
# ============================================================
# reaction_path_plot.py
# 3 シナリオの反応経路を比較プロット
# ============================================================
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
import matplotlib.gridspec as gridspec
# ---- フォント設定（日本語環境） ----
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,
})
# ---- データ読み込み ----
def load_phreeqc_selected(fp, comment="#"):
return pd.read_csv(fp, sep="\t", comment=comment, skipinitialspace=True)
df1 = load_phreeqc_selected("limestone_path.txt")
df2 = load_phreeqc_selected("granite_path.txt")
df3 = load_phreeqc_selected("pyrite_path.txt")
# ======= 3x2 subplot =======
fig = plt.figure(figsize=(14, 10))
gs = gridspec.GridSpec(3, 2, hspace=0.45, wspace=0.35)
# Filter data to exclude initial state (step = -99)
df1_plot = df1[df1["step"] >= 0]
df2_plot = df2[df2["step"] >= 0]
df3_plot = df3[df3["step"] >= 0]
# ======= シナリオ1：石灰岩 =======
ax1a = fig.add_subplot(gs[0, 0])
ax1b = ax1a.twinx()
ax1a.plot(df1_plot["step"], df1_plot["pH"], color="#2563EB", lw=2.2, label="pH")
ax1b.plot(df1_plot["step"], df1_plot["SI_Calcite"], color="#D97706", lw=1.8, ls="--", label="SI(Calcite)")
ax1b.axhline(0, color="#16A34A", lw=1, ls=":", alpha=0.7)
ax1a.set(xlabel="反応ステップ", ylabel="pH", title="シナリオ1：石灰岩溶解（開放系）")
ax1b.set_ylabel("SI", color="#D97706")
lines = ax1a.get_lines() + ax1b.get_lines()
ax1a.legend(lines, [l.get_label() for l in lines], fontsize=9, loc="lower right")
ax1a.grid(True, ls="--", lw=0.5, color="#E5E7EB")
ax1c = fig.add_subplot(gs[0, 1])
ax1c.plot(df1_plot["step"], df1_plot["Ca_mM"], color="#16A34A", lw=2, label="Ca²⁺ (mM)")
ax1c.plot(df1_plot["step"], df1_plot["HCO3_mM"], color="#2563EB", lw=2, ls="--", label="HCO₃⁻ (mM)")
ax1c.set(xlabel="反応ステップ", ylabel="濃度 (mM)",
title="Ca²⁺・HCO₃⁻ の増加")
ax1c.legend(fontsize=9); ax1c.grid(True, ls="--", lw=0.5, color="#E5E7EB")
# ======= シナリオ2：花崗岩 =======
ax2a = fig.add_subplot(gs[1, 0])
ax2a.plot(df2_plot["step"], df2_plot["pH"], color="#92400E", lw=2.2, label="pH")
ax2a.set(xlabel="反応ステップ", ylabel="pH", title="シナリオ2：花崗岩風化（閉鎖系）")
ax2a.grid(True, ls="--", lw=0.5, color="#E5E7EB")
ax2b = fig.add_subplot(gs[1, 1])
ax2b.plot(df2_plot["step"], df2_plot["SI_Gibbsite"], color="#DC2626", lw=1.8, label="SI(Gibbsite)")
ax2b.plot(df2_plot["step"], df2_plot["SI_Kaolinite"], color="#2563EB", lw=1.8, ls="--", label="SI(Kaolinite)")
ax2b.plot(df2_plot["step"], df2_plot["SI_Quartz"], color="#16A34A", lw=1.8, ls="-.", label="SI(Quartz)")
ax2b.axhline(0, color="#9CA3AF", lw=1, ls=":", alpha=0.7)
ax2b.set(xlabel="反応ステップ", ylabel="SI",
title="二次鉱物の飽和指数変化")
ax2b.legend(fontsize=9); ax2b.grid(True, ls="--", lw=0.5, color="#E5E7EB")
# ======= シナリオ3：黄鉄鉱 =======
ax3a = fig.add_subplot(gs[2, 0])
ax3a.plot(df3_plot["step"], df3_plot["pH"], color="#DC2626", lw=2.2, label="pH")
ax3a.axhline(4.0, color="#9CA3AF", lw=1, ls="--", alpha=0.5)
ax3a.text(2, 4.2, "pH 4 (AMD基準)", color="#9CA3AF", fontsize=8)
ax3a.set(xlabel="反応ステップ", ylabel="pH", title="シナリオ3：黄鉄鉱酸化（AMD形成）")
ax3a.grid(True, ls="--", lw=0.5, color="#E5E7EB")
ax3b = fig.add_subplot(gs[2, 1])
ax3b_r = ax3b.twinx()
ax3b.plot(df3_plot["step"], df3_plot["Fe_mM"], color="#EA580C", lw=2, label="Fe (mM)")
ax3b.plot(df3_plot["step"], df3_plot["SO4_mM"], color="#CA8A04", lw=2, ls="--", label="SO₄²⁻ (mM)")
ax3b_r.plot(df3_plot["step"], df3_plot["SI_Goethite"], color="#374151", lw=1.5, ls=":", label="SI(Goethite)")
ax3b.set(xlabel="反応ステップ", ylabel="濃度 (mM)",
title="Fe・SO₄²⁻ 増加と Goethite 過飽和")
lines = ax3b.get_lines() + ax3b_r.get_lines()
ax3b.legend(lines, [l.get_label() for l in lines], fontsize=9)
ax3b_r.set_ylabel("SI(Goethite)", color="#374151")
ax3b.grid(True, ls="--", lw=0.5, color="#E5E7EB")
plt.suptitle("Reaction Path Modeling: 3 Scenario Comparison (PHREEQC)", fontsize=14, y=0.98)
plt.savefig("reaction_paths.svg", bbox_inches="tight")
plt.show()
print("✅ reaction_paths.svg を保存")
3 シナリオの比較結果
図3. 3 シナリオの pH 変化経路 — 石灰岩は上昇、花崗岩は緩やかな上昇、黄鉄鉱は急激な低下
考察：3 シナリオが示すこと
🏔️ 石灰岩（開放系）
pH は 6 → 9.3 まで上昇
SI(Calcite) は 0を超えて増加（≈ +5）
Ca²⁺・HCO₃⁻ は継続的に増加
CO₂ 供給により溶解は持続するが、
pH 上昇によって強い過飽和へ移行
（平衡で止まらない系）
SI(Calcite) は 0を超えて増加（≈ +5）
Ca²⁺・HCO₃⁻ は継続的に増加
CO₂ 供給により溶解は持続するが、
pH 上昇によって強い過飽和へ移行
（平衡で止まらない系）
🪨 花崗岩（閉鎖系）
pH は 7.6〜7.8付近で頭打ち後わずかに低下
Kaolinite・Gibbsite は 強い過飽和
Quartz はゆっくり飽和に接近
二次鉱物は過飽和だが、
実際の沈殿は速度・核形成に制約される
（非平衡状態が持続）
Kaolinite・Gibbsite は 強い過飽和
Quartz はゆっくり飽和に接近
二次鉱物は過飽和だが、
実際の沈殿は速度・核形成に制約される
（非平衡状態が持続）
⚡ 黄鉄鉱（酸化）
pH は 6 → 3 → ~1 まで急落
Fe・SO₄²⁻ は指数的に増加
SI(Goethite) は初期高いが低下
Fe³⁺ による自己触媒で反応加速
酸生成により緩衝能が崩壊し
強酸性条件へ移行
Fe・SO₄²⁻ は指数的に増加
SI(Goethite) は初期高いが低下
Fe³⁺ による自己触媒で反応加速
酸生成により緩衝能が崩壊し
強酸性条件へ移行
「平衡計算」と「反応経路計算」の使い分け
平衡計算（EQUILIBRIUM_PHASES）
系が完全に平衡に達したと仮定したときの状態を求める手法。
現在の水がどの鉱物と平衡にあるか、またはどの鉱物に対して過飽和／不飽和かを評価するために用いる（SI の確認）。
現在の水がどの鉱物と平衡にあるか、またはどの鉱物に対して過飽和／不飽和かを評価するために用いる（SI の確認）。
反応経路計算（REACTION）
有限量の反応物（鉱物・ガス・溶質など）を段階的に加えながら、系の変化を追跡する手法。
例えば、「鉱物をどれだけ溶解させると pH や濃度がどのように変化するか」を評価する際に用いる。
例えば、「鉱物をどれだけ溶解させると pH や濃度がどのように変化するか」を評価する際に用いる。
これらを組み合わせることで、各ステップでの平衡状態を評価しながら反応の進行を追跡でき、地下水の進化過程（風化・溶解・沈殿）を一連のプロセスとして再現することが可能となる。
まとめ：シリーズ全体の統合
ヒント次回 #12「移流分散モデル — 地下水中の物質移動」
TRANSPORT ブロックを使い、帯水層の中を汚染物質や溶存成分が 移流（流れに運ばれる）と分散（拡散）によって広がっていく様子をシミュレートする。 #11 の反応経路にさらに「空間」の次元が加わる、シリーズ最大の統合記事である。
参考文献（References）
Appelo, CAJ, と Dieke Postma. 2005年. Geochemistry, groundwater and pollution. Second. Balkema, Rotterdam, p. 634.
Parkhurst, David L, と CAJ Appelo. 2013年. Description of input and examples for PHREEQC version 3—A computer program for speciation, batch-reaction, one-dimensional transport, and inverse geochemical calculations. US Geological Survey Techniques; Methods, book 6, chap. A43, 497 p.
Yamamoto, S. 1983年. Method of the groundwater survey. Kokon Shoin, Tokyo (in Japanese), 490 p.
Yang, Heejun, T Mishima, S Katazakai, と M Kagabu. 2023年. 「Analytical approach using a chemical equilibrium formula and geochemical modeling for alkalinity measurements of small natural water samples」. Applied Geochemistry 148: 105535.
- #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 の応用）
DeepFlow | 地球科学シミュレーションの深みへ