Pythonコード

import flopy

import numpy as np

import matplotlib.pyplot as plt

import os

import shutil

import sys

import time

import csv

from phreeqpy.iphreeqc.phreeqc_dll import IPhreeqc

# ==========================================

# 0. 準備

# ==========================================

def get_exe_path():

    try:

        base_dir = os.path.dirname(os.path.abspath(__file__))

    except NameError:

        base_dir = os.getcwd()

    exe_name = "mf6.exe"

    exe_path = os.path.join(base_dir, exe_name)

    if not os.path.exists(exe_path):

        print(f"エラー: {exe_path} が見つかりません。")

        sys.exit(1)

    db_name = "llnl.dat"

    db_path = os.path.join(base_dir, db_name)

    if not os.path.exists(db_path):

        print(f"エラー: {db_name} が見つかりません。")

        sys.exit(1)

    print(f"Using database: {db_path}")

    return exe_path, db_path, base_dir

# ==========================================

# 1. Flopyによる流動モデル

# ==========================================

def setup_and_run_flow(exe_path, length, radius, n_cells, flow_rate_m3s, k_perms):

    name = "column_flow"

    ws = "workspace_column"

    base_dir = os.path.dirname(exe_path)

    ws_path = os.path.join(base_dir, ws)

    if os.path.exists(ws_path): shutil.rmtree(ws_path)

    os.makedirs(ws_path, exist_ok=True)

    area = np.pi * radius**2

    delr = length / n_cells

    sim = flopy.mf6.MFSimulation(sim_name=name, exe_name=exe_path, sim_ws=ws_path)

    tdis = flopy.mf6.ModflowTdis(sim, time_units='SECONDS', nper=1, perioddata=[(1.0, 1, 1.0)])

    ims = flopy.mf6.ModflowIms(sim)

    gwf = flopy.mf6.ModflowGwf(sim, modelname=name, save_flows=True)

    sim.register_ims_package(ims, [gwf.name])

    dis = flopy.mf6.ModflowGwfdis(gwf, nlay=1, nrow=1, ncol=n_cells,

                                  delr=delr, delc=np.sqrt(area), top=length, botm=0.0)

    npf = flopy.mf6.ModflowGwfnpf(gwf, k=k_perms, save_specific_discharge=True)

    ic = flopy.mf6.ModflowGwfic(gwf, strt=1.0)

    wel_spd = {0: [[(0, 0, 0), flow_rate_m3s]]}

    wel = flopy.mf6.ModflowGwfwel(gwf, stress_period_data=wel_spd)

    chd_spd = {0: [[(0, 0, n_cells-1), 1.0]]}

    chd = flopy.mf6.ModflowGwfchd(gwf, stress_period_data=chd_spd)

    budget_file = f"{name}.cbc"

    oc = flopy.mf6.ModflowGwfoc(gwf,

                                budget_filerecord=budget_file,

                                head_filerecord=f"{name}.hds",

                                saverecord=[('HEAD', 'ALL'), ('BUDGET', 'ALL')])

    sim.write_simulation()

    sim.run_simulation(silent=True)

    cbc_path = os.path.join(ws_path, budget_file)

    try:

        cbb = flopy.utils.CellBudgetFile(cbc_path, precision='double')

        spdis = cbb.get_data(text='DATA-SPDIS')[0]

        qx = spdis['qx']

        avg_q = np.mean(np.abs(qx))

        return avg_q

    except Exception as e:

        print(f"Error reading budget: {e}")

        return None

# ==========================================

# 2. PhreeqPyによる反応輸送モデル

# ==========================================

def run_reactive_transport_kinetics(darcy_velocity, length, n_cells, porosity, flow_rate_ml_min, db_path):

    print("--- Step 2: Running PHREEQC Transport ---")

    velocity = darcy_velocity / porosity

    delr = length / n_cells

    dt = delr / velocity

    # ---------------------------------------------------------

    # 時間設定の変更

    # ---------------------------------------------------------

    target_hours = 24.0  # 目標時間 [時間]

    target_seconds = target_hours * 3600.0

    total_steps = int(target_seconds / dt)

    print(f"Target Time: {target_hours} hours")

    print(f"Flow Rate  : {flow_rate_ml_min} mL/min")

    print(f"Delta t    : {dt:.2f} sec (Time per step)")

    print(f"Total Steps: {total_steps}")

    print(f"Conditions : 150 C, 10 MPa")

    phreeqc = IPhreeqc()

    phreeqc.load_database(db_path)

    # PHASES

    phases_def = """

    PHASES

        Amor.Si

            SiO2 = SiO2(aq)

            log_k -2.71

        K-Feldspar

            KAlSi3O8 + 4H+ = K+ + Al+++ + 3SiO2(aq) + 2H2O

            log_k -0.275

            delta_h -53.53 kJ

        Magnetite

            Fe2SiO4 + 4H+ = 2Fe++ + SiO2(aq) + 2H2O

            log_k 3.0

        Fe(OH)3

            Fe(OH)3 + 3H+ = Fe+3 + 3H2O

            log_k 3.0

        Boehmite

            AlOOH + 3H+ = Al+3 + 2H2O

            log_k 8.5

    """

    phreeqc.run_string(phases_def)

    # RATES

    rates_def = """

    RATES

        Albite

        -start

        10  k_acid = 10^-10.2 * EXP(-65.0 / 0.008314 * (1/TK - 1/298.15))

        20  k_neut = 10^-12.6 * EXP(-69.8 / 0.008314 * (1/TK - 1/298.15))

        30  k_base = 10^-15.6 * EXP(-71.0 / 0.008314 * (1/TK - 1/298.15))

        40  rate = (k_acid * ACT("H+")^0.46 + k_neut + k_base * ACT("H+")^-0.57) * (1 - SR("Albite"))

        50  moles = rate * TIME * PARM(1)

        60  SAVE moles

        -end

        Anorthite

        -start

        40  rate = (10^-10.2 * ACT("H+")^0.46 + 10^-12.6) * (1 - SR("Anorthite"))

        50  moles = rate * TIME * PARM(1)

        60  SAVE moles

        -end

        Enstatite

        -start

        40  rate = (10^-10.2 * ACT("H+")^0.5) * (1 - SR("Enstatite"))

        50  moles = rate * TIME * PARM(1)

        60  SAVE moles

        -end

        Magnetite

        -start

        40  rate = 10^-11.0 * (1 - SR("Magnetite"))

        50  moles = rate * TIME * PARM(1)

        60  SAVE moles

        -end

        Diopside

        -start

        40  rate = (10^-10.2 * ACT("H+")^0.5) * (1 - SR("Diopside"))

        50  moles = rate * TIME * PARM(1)

        60  SAVE moles

        -end

        K-Feldspar

        -start

        40  rate = (10^-10.2 * ACT("H+")^0.46) * (1 - SR("K-Feldspar"))

        50  moles = rate * TIME * PARM(1)

        60  SAVE moles

        -end

        Amor.Si

        -start

        40  rate = 10^-13.0 * (1 - SR("Amor.Si"))

        50  moles = rate * TIME * PARM(1)

        60  SAVE moles

        -end

    """

    phreeqc.run_string(rates_def)

    # ---------------------------------------------------------

    # 初期状態 & KINETICS

    # ---------------------------------------------------------

    # 【修正】 parms (表面積) を大幅に大きくしました (1e2 ~ 1e4)

    # これにより反応速度が上がり、濃度が数桁上昇します。

    # 計算負荷を下げるため -step_divide をさらに強化しています。

    initial_sol_def = f"""

    SOLUTION 1-{n_cells} Initial Pore Water

        units mmol/kgw

        pH 7.0

        temp 150

        pressure 100

    EQUILIBRIUM_PHASES 1-{n_cells}

        Calcite   0 0

        Fe(OH)3   0 0

        Boehmite  0 0

        Dolomite  0 0

        Magnesite 0 0

    KINETICS 1-{n_cells}

    -step_divide 1e7     # 分割数を増やして収束しやすくする

    -bad_step_max 1000

    Albite

        -formula  NaAlSi3O8

        -m        0.1285

        -parms    1e3    # 0.1 -> 1000

    Anorthite

        -formula  CaAl2Si2O8

        -m        0.1285

        -parms    1e2    # 0.1 -> 100

    Enstatite

        -formula  MgSiO3

        -m        0.072

        -parms    1e2    # 0.1 -> 100

    Magnetite

        -formula  Fe2SiO4

        -m        0.054

        -parms    1e2    # 0.1 -> 100

    Diopside

        -formula  CaMg(SiO3)2

        -m        0.052

        -parms    1e3    # 0.1 -> 1000

    K-Feldspar

        -formula  KAlSi3O8

        -m        0.03

        -parms    1e2    # 0.1 -> 100

    Amor.Si

        -formula  SiO2

        -m        0.417

        -parms    1e4    # 0.1 -> 10000

    SAVE solution 1-{n_cells}

    END

    """

    phreeqc.run_string(initial_sol_def)

    # 流入水

    injection_sol_def = """

    SOLUTION 0 Inflow Water

        units mmol/kgw

        pH 7.0 charge

        temp 150

        pressure 100

    EQUILIBRIUM_PHASES 0

        CO2(g) 1.3 10.0

    SAVE solution 0

    END

    """

    phreeqc.run_string(injection_sol_def)

    # データ格納用辞書

    results = {

        'time': [], 'pH': [],

        'Na': [], 'Ca': [], 'K': [], 'Mg': [], 'Fe': [], 'Al': [], 'Si': [],

        'Calcite': [], 'Fe(OH)3': [], 'Boehmite': [], 'Magnetite': []

    }

    phreeqc.run_string("""

    SELECTED_OUTPUT

        -reset false

        -pH

        -totals Na Ca K Mg Fe Al Si

        -equilibrium_phases Calcite Fe(OH)3 Boehmite

        -kinetic_reactants Magnetite

    END

    """)

    start_time = time.time()

    col_map = {}

    # ループ開始

    for step in range(1, total_steps + 1):

        transport_step = f"""

        TRANSPORT

            -cells {n_cells}

            -shifts 1

            -time_step {dt}

            -punch_cells {n_cells}

            -punch_frequency 1

        END

        """

        try:

            phreeqc.run_string(transport_step)

        except Exception as e:

            print(f"\nError at step {step}: {e}")

            break

        elapsed = time.time() - start_time

        # 進捗表示（ステップ数が多いので見やすく）

        sys.stdout.write(f"\rProgress: [{step}/{total_steps}] ({elapsed:.1f}s) - {(step/total_steps)*100:.1f}%")

        sys.stdout.flush()

        res = phreeqc.get_selected_output_array()

        if res and len(res) > 1:

            header = res[0]

            data = res[-1]

            if not col_map:

                h_up = [str(x).upper() for x in header]

                def find_idx(key):

                    for i, h in enumerate(h_up):

                        if key.upper() == h: return i

                    for i, h in enumerate(h_up):

                        if key.upper() in h: return i

                    return None

                col_map['pH'] = find_idx('pH')

                for elem in ['Na', 'Ca', 'K', 'Mg', 'Fe', 'Al', 'Si']:

                    col_map[elem] = find_idx(elem)

                for min_name in ['Calcite', 'Fe(OH)3', 'Boehmite', 'Magnetite']:

                    col_map[min_name] = find_idx(min_name)

            # 時間 (時間単位)

            current_time_hr = step * dt / 3600.0

            results['time'].append(current_time_hr)

            results['pH'].append(data[col_map['pH']])

            for elem in ['Na', 'Ca', 'K', 'Mg', 'Fe', 'Al', 'Si']:

                idx = col_map[elem]

                val = data[idx] if idx is not None else 0.0

                results[elem].append(val * 1000) # mmol/kgw

            for min_name in ['Calcite', 'Fe(OH)3', 'Boehmite', 'Magnetite']:

                idx = col_map[min_name]

                val = data[idx] if idx is not None else 0.0

                results[min_name].append(val)

    print("\nCalculation finished.")

    return results

# ==========================================

# メイン実行部

# ==========================================

if __name__ == "__main__":

    L = 0.15          

    r = 0.01          

    N = 15            

    Q_ml_min = 1.0    

    phi = 0.35        

    K = 1e-4          

    Q_m3s = (Q_ml_min * 1e-6) / 60.0

    mf6_path, db_path, base_dir = get_exe_path()

    q_val = setup_and_run_flow(mf6_path, L, r, N, Q_m3s, K)

    if q_val is not None:

        print(f"Darcy Velocity: {q_val:.2e} m/s")

        res = run_reactive_transport_kinetics(q_val, L, N, phi, Q_ml_min, db_path)

        times = res['time']

        if len(times) > 0:

            # 保存

            output_file = os.path.join(base_dir, "simulation_results_24h.dat")

            print(f"Saving results to: {output_file}")

            try:

                with open(output_file, 'w', newline='', encoding='utf-8') as f:

                    writer = csv.writer(f)

                    header_row = ['Time(hr)', 'pH'] + \

                                 [f"{e}(mmol/kgw)" for e in ['Na', 'Ca', 'K', 'Mg', 'Fe', 'Al', 'Si']] + \

                                 [f"{m}(mol)" for m in ['Calcite', 'Fe(OH)3', 'Boehmite', 'Magnetite']]

                    writer.writerow(header_row)

                    for i in range(len(times)):

                        row = [times[i], res['pH'][i]]

                        row += [res[e][i] for e in ['Na', 'Ca', 'K', 'Mg', 'Fe', 'Al', 'Si']]

                        row += [res[m][i] for m in ['Calcite', 'Fe(OH)3', 'Boehmite', 'Magnetite']]

                        writer.writerow(row)

                print("File saved successfully.")

            except Exception as e:

                print(f"Error saving file: {e}")

            # グラフ描画

            fig, (ax1, ax2, ax3) = plt.subplots(3, 1, figsize=(10, 15), sharex=True)

            # 1. 溶存濃度 (対数軸)

            elements = ['Na', 'Ca', 'K', 'Mg', 'Fe', 'Al', 'Si']

            colors = plt.cm.tab10(np.linspace(0, 1, len(elements)))

            for i, elem in enumerate(elements):

                ax1.plot(times, res[elem], label=elem, color=colors[i], linewidth=2)

            ax1.set_yscale('log')

            ax1.set_ylabel('Concentration (mmol/kgw)')

            ax1.set_title('1. Aqueous Species Concentrations (Log Scale)')

            ax1.legend(loc='upper right', bbox_to_anchor=(1.15, 1))

            ax1.grid(True, which="both", alpha=0.3)

            # 2. pH

            ax2.plot(times, res['pH'], 'k-', linewidth=2, label='pH')

            ax2.set_ylabel('pH')

            ax2.set_title('2. pH Change')

            ax2.grid(True, alpha=0.3)

            ax2.legend()

            # 3. 鉱物量

            minerals = ['Calcite', 'Fe(OH)3', 'Boehmite', 'Magnetite']

            colors_min = ['blue', 'brown', 'cyan', 'black']

            for i, mn in enumerate(minerals):

                ax3.plot(times, res[mn], label=mn, color=colors_min[i], linewidth=2)

            ax3.set_ylabel('Mineral Amount (mol)')

            ax3.set_xlabel('Time (hours)')

            ax3.set_title('3. Mineral Phase Changes')

            ax3.legend()

            ax3.grid(True, alpha=0.3)

            plt.tight_layout()

            plt.show()

        else:

            print("No results to plot.")