なぜ「カップリング」が必要なのか
地下水の流動を計算するだけなら、MODFLOW 6単体で十分です。しかし現実の地下環境では、水が流れると同時に岩石と反応し、鉱物が溶けたり沈殿したりします。熱水系では温度・圧力が高く、この反応がさらに複雑になります。
この熱水系において、流動モデルと反応モデルを切り離して使うと、岩石-水反応を正しく追うことができなくなり（溶解・沈殿現象）、観測データの再現性は保証されません。
そこで本シリーズでは、MODFLOW 6（流動計算）とPHREEQC（地球化学反応計算）をPythonで結合し、1次元カラムスケールでの熱水流動と地球化学反応を同時にシミュレーションする手法を段階的に解説します。
モデルの全体像
本シミュレーションは3つのコアコンポーネントで構成されます。
┌─────────────────┐
│ MODFLOW 6 │ 💧 流動計算
│ (Flopy) │ ・カラム内の圧力・流量
│ │ ・Darcy流速の算出
└────────┬────────┘
│ Darcy流速 (v)
▼
┌─────────────────┐
│ PHREEQC │ 🧪 反応輸送計算
│ (phreeqpy) │ ・移流（TRANSPORT）
│ │ ・鉱物溶解・沈殿（KINETICS）
└────────┬────────┘
│ 濃度・pH・SI・鉱物量
▼
┌─────────────────┐
│ Python │ 📊 統合・可視化
│ (numpy + │ ・時系列プロット
│ matplotlib) │ ・破過曲線の解析
└─────────────────┘3つのコンポーネントの役割を整理すると以下のようになります：
| コンポーネント | ライブラリ | 担当 |
|---|---|---|
| 流動計算 | flopy |
カラム内の流量・水頭・Darcy流速 |
| 反応輸送計算 | phreeqpy |
pH・鉱物量・元素濃度の時間変化 |
| 統合・可視化 | numpy, matplotlib |
データ受け渡し・結果の描画 |
ポイントはMODFLOWが計算したDarcy流速を、そのままPHREEQCの移流速度パラメータに渡すという橋渡しの仕組みです。これがカップリングの核心となります（詳細は第2回で実装します）。
使用ライブラリ
# 流動モデル
import flopy # MODFLOW 6 の Python インターフェース
# 反応計算
from phreeqpy.iphreeqc.phreeqc import IPhreeqc # PHREEQC エンジン
# 数値・可視化
import numpy as np
import matplotlib.pyplot as plt
# ファイル操作
import csv, os, shutilインストールは以下：
pip install flopy
pip install phreeqpy
pip install numpy matplotlib本シリーズでは高温高圧条件（150℃、100atm）に対応した llnl.dat を使用します。これもUSGS/USGS提供のPHREEQCパッケージに含まれているので、同じく作業ディレクトリに配置しておきます。
カラムモデルの設計
物理的なイメージ
円柱形のコアサンプル（あるいは地層の1断面）を想定します。水は一端から注入され、他端から排出されます。内部では岩石鉱物と反応しながら流れます。
Model Discretization
カラム流動モデルの離散化イメージ
MODFLOW 6による離散化
MODFLOWは連続的な空間をグリッドセルに分割して計算します。1次元カラムの場合、DISパッケージで以下のように定義します：
import flopy
# シミュレーション名
sim_name = "column_model"
# --- 基本パラメータ ---
L = 1.0 # カラム長さ [m]
N = 15 # セル分割数
r = 0.05 # カラム半径 [m]
dx = L / N # セル幅 [m]
A = 3.14159 * r**2 # 断面積 [m²]
# --- MODFLOW モデル構築 ---
sim = flopy.mf6.MFSimulation(sim_name=sim_name, exe_name="mf6")
# TDIS：時間設定
tdis = flopy.mf6.ModflowTdis(
sim,
time_units="SECONDS",
nper=1,
perioddata=[(86400, 1, 1.0)]
)
gwf = flopy.mf6.ModflowGwf(sim, modelname=sim_name)
# DISパッケージ：グリッド定義
dis = flopy.mf6.ModflowGwfdis(
gwf,
nlay=1, # 層数：1（1次元）
nrow=1, # 行数：1
ncol=N, # 列数：セル分割数
delr=dx, # セル幅（流れ方向）
delc=1.0, # セル幅（直交方向）：仮想的に1.0
top=1.0,
botm=0.0,
)セル数 N=15 の場合、カラムは15個の計算ブロックに分割されます。各セルの幅は均一で dx = L/N となります。
境界条件の設定
1次元流動モデルには3種類のパッケージを組み合わせます：
| パッケージ | 場所 | 役割 |
|---|---|---|
| WEL（井戸境界） | 入口 Cell 1 | 一定流量 \(Q_{in}\) を注入 |
| CHD（定水頭境界） | 出口 Cell N | 一定水頭 \(h_{out}\) を維持 |
| NPF（節点特性） | 全セル | 透水係数 \(K\) を定義 |
| IC（初期条件） | 全セル | 初期水頭 \(h_{init}\) を設定 |
# --- 物性値 ---
K = 1.0e-5 # 透水係数 [m/s]（均質・等方）
Qin = 1.0e-6 # 注入流量 [m³/s]
h_out = 0.0 # 出口水頭 [m]
h_init = 1.0 # 初期水頭 [m]
# NPF：透水係数
npf = flopy.mf6.ModflowGwfnpf(
gwf,
k=K, # x方向
k22=K, # y方向
k33=K, # z方向
)
# IC：初期条件
ic = flopy.mf6.ModflowGwfic(gwf, strt=h_init)
# WEL：入口境界（定流量）
wel = flopy.mf6.ModflowGwfwel(
gwf,
stress_period_data={0: [[(0, 0, 0), Qin]]}
)
# CHD：出口境界（定水頭）
chd = flopy.mf6.ModflowGwfchd(
gwf,
stress_period_data={0: [[(0, 0, N-1), h_out]]}
)本シリーズでは検証目的のため、透水係数は全セルで一定（均質・等方）とします。実際の地層不均質性を組み込む場合は、k に配列を渡すことで対応できます。
次回予告：実装編（カップリングの核心）
第1回では「何を、なぜ組み合わせるのか」「カラムモデルをどう設計するか」を解説しました。
第2回では以下を実装します：
- MODFLOW 6の実行とCell Budget Fileからの流速抽出
- Courant条件（\(C_r = 1\)）を満たす時間刻みの自動計算
- PHREEQCの初期化（150℃・100atm条件）
- カップリングループ：流速をPHREEQCのTRANSPORTブロックへ動的注入
「なぜCourant数を1に揃えるのか」という数値安定性の話も丁寧に解説する予定です。
まとめ
| 項目 | 設定値 |
|---|---|
| カラム長さ | 1.0 m |
| セル分割数 | 15 |
| 透水係数 | 1.0×10⁻⁵ m/s |
| 注入流量 | 1.0×10⁻⁶ m³/s |
| 温度条件 | 150℃ |
| 圧力条件 | 100 atm（10 MPa） |
| 計算期間 | 24時間 |
本シリーズのコードはすべて GitHubリポジトリ で公開予定です。次回の実装編もお楽しみに。
このシリーズの他の記事：
- #1 モデル設計編（本記事）
- #2 実装編（カップリングの核心）
- #3 解析・可視化編