PHREEQC from Scratch #12: Advection-Dispersion Model — Solute Transport in Groundwater

Simulating solute transport in groundwater using the TRANSPORT block. Combining advection, dispersion, and reaction to track everything from contaminant diffusion to carbonate dissolution fronts in 2D (space × time).
Geochemistry
PHREEQC
Author

DeepFlow

Published

May 18, 2026

Introduction: The Grand Integration of the Series

From #2 to #11, we calculated the chemistry at “a single point and a single moment”.

  • Speciation (#2): Solution composition at one point
  • EQUILIBRIUM_PHASES (#3〜#5): Mineral equilibrium at one point
  • REACTION (#11): Changes over time

This time, the TRANSPORT block adds the dimension of space to this. As groundwater flows through an aquifer, it reacts with minerals, transports contaminants, and forms concentration fronts—we will calculate this entire process all at once.

Everything built up in this series is utilized here
Activity coefficients (#9) × Saturation indices (#10) × Reaction paths (#11)
　　　　　　↓
TRANSPORT: Space × Time × Chemical Reaction
NoteWhat you will learn in this article
  • The three concepts of advection, dispersion, and reaction
  • Basic syntax of the TRANSPORT block
  • Scenario 1: Movement of a saline front (pure advection-dispersion)
  • Scenario 2: Calcite dissolution front (advection + chemical reaction)
  • Scenario 3: Heavy metal contaminant diffusion and soil adsorption
  • Visualization of concentration profiles using Python

Theory: Advection-Dispersion Equation

Solute transport in groundwater is described by the Advection-Dispersion Equation (ADE):

\[\frac{\partial C}{\partial t} = D \frac{\partial^2 C}{\partial x^2} - v \frac{\partial C}{\partial x} + R\]

Term Name Meaning PHREEQC Parameter
∂C/∂t Time change Rate of change in concentration over time -time_step
D ∂²C/∂x² Dispersion term Diffusion/dispersion due to concentration gradient -dispersivity
v ∂C/∂x Advection term Mass transport by groundwater flow -shifts (velocity × time)
R Reaction term Mineral dissolution/precipitation, adsorption, etc. EQUILIBRIUM_PHASES etc.

PHREEQC’s Approach: Operator Splitting

PHREEQC does not solve the ADE directly; instead, it uses a method where advection, dispersion, and reaction are calculated alternately:

Operator Splitting Concept Diagram ① Advection Shift cell forward (Repeat -shifts times) ② Dispersion Mix with adjacent cells (Controlled by dispersivity) ③ Reaction Equilibrate in each cell (EQUILIBRIUM_PHASES etc.) Next time step Repeat ①-③ -shifts times Record results for each cell 1 time step = 1 shift = Time for groundwater to move one cell

Basic Syntax of the TRANSPORT Block

TRANSPORT
    -cells            20        # Number of cells (spatial divisions)
    -length           0.1       # Length of each cell (m)
    -shifts           40        # Number of time steps (advection steps)
    -time_step        864000    # Time per step (seconds) = 10 days
    -flow_direction   forward   # Flow direction (forward/back/diffusion_only)
    -boundary_conditions  flux  flux  # Boundary conditions (flux/constant)
    -dispersivity     0.015     # Dispersivity (m)
    -diffusion_coefficient 1e-9 # Diffusion coefficient (m²/s)
    -punch_cells      1-20      # Cells to output
    -punch_frequency  10        # Output frequency (every 10 steps)
    -print_cells      1 5 10 20 # Cells to display on screen
TipBalancing Cell Number and Computational Accuracy

Increasing the number of cells improves accuracy but also increases computation time. Determine the cell length so that the Peclet number (Pe = v·Δx/D) is 2 or less. For Pe > 2, numerical dispersion occurs, overestimating physical dispersion.

Scenario 1: Movement of a Saline Front (Pure Advection-Dispersion)

Setup

A NaCl solution (0.1 mol/L) is injected into an aquifer initially filled with pure water (20 cells, 0.5 m each = 10 m total length). We will track how the saline front moves and disperses.

# ============================================================
#  DeepFlow #12 - Scenario 1: Saline Front Movement
#  Injecting NaCl solution into a pure water aquifer
# ============================================================

# ---- Initial Solution: Pure Water (fills all 20 cells) ----
SOLUTION 0  "Injection Water (NaCl 0.1 mol/L)"
    temp      25
    pH        7
    units     mol/kgw
    Na        0.1
    Cl        0.1

SOLUTION 1-20  "Initial Aquifer (Pure Water)"
    temp      25
    pH        7
    -water    1

# ---- TRANSPORT Settings ----
TRANSPORT
    -cells            20
    -lengths          20*0.5    # Repeat 0.5 m for 20 cells (PHREEQC array notation, not multiplication)
    -shifts           40        # 40 time steps
    -time_step        86400     # 1 day/step. Apparent velocity: 0.5 m / 86400 s = 5.8e-6 m/s
    -flow_direction   forward
    -boundary_conditions  flux  flux
    -dispersivities   20*0.05   # Repeat 0.05 m for 20 cells (same as above)
    -diffusion_coefficient 1e-9
    -punch_cells      1-20
    -punch_frequency  10        # Output at steps 10, 20, 30, 40

# ---- Output Settings ----
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

# ---- Graph Settings ----
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

Scenario 2: Calcite Dissolution Front (Advection + Chemical Reaction)

Setup

Acidic water (pH 4.5) is injected into an aquifer containing calcite (initial state: Ca–HCO₃ water in equilibrium with calcite). We track how the dissolution front moves, observing pH, Ca²⁺, and SI(Calcite). What happens when acidic water enters a limestone aquifer?

# ============================================================
#  Scenario 2: Calcite Dissolution Front
#  Injecting acidic water into a limestone aquifer
# ============================================================

# ---- Injection Water (Acidic + CO2) ----
SOLUTION 0  "Injected Acidic Water"
    temp      12
    pH        4.5
    pe        12
    units     mol/kgw
    -water    1

EQUILIBRIUM_PHASES 0
    CO2(g)   -2.0   10   # Soil CO2

# ---- Initial Groundwater (Calcite + CO2 Equilibrium) ----
SOLUTION 1-20  "Initial Groundwater"
    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    # Repeat 0.5 m for 20 cells (PHREEQC array notation, not multiplication)
    -shifts           60
    -time_step        86400
    -flow_direction   forward
    -boundary_conditions  flux flux
    -dispersivities   20*0.05   # Repeat 0.05 m for 20 cells
    -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

# ---- Graph Settings (1: Water Chemistry) ----
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

# ---- Graph Settings (2: Mineral) ----
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

Interpreting the Results

Observation Geochemical Meaning
pH front moves downstream Propagation of the low pH zone due to acidic water advection. Dispersion makes the front gradual.
Ca²⁺ rises ahead of the front Ca²⁺ is supplied by calcite dissolution. It dissolves just enough to satisfy equilibrium conditions.
SI(Calcite) = 0 zone exists A zone where acid is consumed and equilibrium with calcite is re-established. Downstream of here maintains the original groundwater composition.
Zone where calcite mass decreases Cells where minerals have been depleted by dissolution. Long-term, porosity increases and hydraulic conductivity changes.

Scenario 3: Heavy Metal (Cd²⁺) Contamination Migration and Soil Adsorption

Setup

We assume a case where Cd²⁺ contaminated water from mine drainage enters an aquifer. We incorporate the retardation effect due to cation exchange in the soil.

# ============================================================
#  Scenario 3: Heavy Metal Contamination and Adsorption
#  Cd²⁺ Advection + Retardation by Cation Exchange
# ============================================================

# ---- Injection Water: Cd²⁺ Contaminated Water ----
SOLUTION 0  "Cd Contaminated Water"
    temp      15
    pH        6.5
    units     mol/kgw
    Na        1.0e-3
    Cl        1.0e-3
    Cd        1.0e-5    # 10 μmol/L of Cd²⁺
    -water    1

# ---- Initial Aquifer Water: Clean Groundwater ----
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

# ---- Soil Cation Exchanger (All Cells) ----
# Using exchanger definitions from phreeqc.dat:
#   NaX, KX, CaX2, MgX2, CdX2 etc.
# Since X- is not defined natively, we initialize via -equilibrate
EXCHANGE 1-20
    NaX      0.0018   # Na occupies exchange sites initially
    CaX2     0.0001   # Ca also slightly occupies
    -equilibrate  1   # Equilibrate with SOLUTION 1 to set initial state

# ---- TRANSPORT ----
TRANSPORT
    -cells            20
    -lengths          20*0.5    # Repeat 0.5 m for 20 cells (PHREEQC array notation, not multiplication)
    -shifts           80
    -time_step        86400
    -flow_direction   forward
    -boundary_conditions  flux  flux
    -dispersivities   20*0.05   # Repeat 0.05 m for 20 cells
    -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
ImportantWhat is the Retardation Factor?

When adsorption occurs, the migration speed of the contaminant becomes slower than the groundwater velocity. This ratio is the retardation factor R, expressed as:

\[R = 1 + \frac{\rho_b \cdot K_d}{\theta}\]

\(\rho_b\): Bulk density, \(K_d\): Distribution coefficient, \(\theta\): Porosity.

Physical Meaning

The retardation factor R means:

Groundwater velocity ÷ Contaminant migration velocity. For example:

  • R=1: Same speed as water (no adsorption)

  • R=5: 1/5 the speed of water

  • R=10: 1/10 the speed of water

Heavy metals like Cd²⁺ strongly adsorb to soil, resulting in a high R, often migrating at only 1/5 to 1/10 the speed of groundwater. The EXCHANGE block in PHREEQC calculates this effect thermodynamically and automatically.

Visualizing Concentration Profiles with Python

# ============================================================
#  transport_plot.py
#  Visualize spatial concentration profiles for 3 scenarios
# ============================================================
import numpy as np
import os
import matplotlib
import matplotlib.font_manager as fm
import matplotlib.pyplot as plt
import matplotlib.cm as cm

# ---- Font Settings ----
plt.rcParams.update({
    "font.family": "sans-serif",
    "font.sans-serif": ["Arial", "Helvetica", "DejaVu Sans"],
    "axes.unicode_minus": False,
    "figure.dpi":      150,
})

#  Numerical Core: Reproducing PHREEQC's Operator Splitting
#  ① Advection: Shift cell one to the right
#  ② Dispersion: Mixing between adjacent cells (Mixing ratio = α/(dx+α))

def adv_disp(nx, dx, alpha, n_shifts, C0, C_init,
             punch_freq, R=1):
    """
    1D Advection-Dispersion Simulator (Operator Splitting)

    Parameters
    ----------
    nx         : Number of cells
    dx         : Cell length (m)
    alpha      : Dispersivity (m)
    n_shifts   : Number of time steps (shifts)
    C0         : Injected concentration
    C_init     : Initial aquifer concentration
    punch_freq : Record every N shifts
    R          : Retardation factor (Adsorption present = R>1)

    Returns
    -------
    list of (shift, concentration_array)
    """
    C = np.full(nx, float(C_init))
    f_mix = alpha / (dx + alpha)   # Dispersion mixing ratio
    results = []

    for s in range(1, n_shifts + 1):
        # ① Advection (considering Retardation R: if R=5, shift once every 5 steps)
        if s % R == 0:
            C_adv = np.empty(nx)
            C_adv[0] = C0          # Inlet boundary: Injected concentration
            C_adv[1:] = C[:-1]     # Shift one cell right
        else:
            C_adv = C.copy()       # No shift (adsorption delay)

        # ② Dispersion (Mixing with adjacent cells)
        C_new = C_adv.copy()
        for i in range(nx):
            left  = C_adv[i - 1] if i > 0      else C0          # Left boundary
            right = C_adv[i + 1] if i < nx - 1 else C_adv[i]   # Right boundary (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):
    """
    pH Front Simulator (Advection-Dispersion + Calcite Buffer Approximation)

    Calcite Buffer Approximation:
      Calculate pH based on explicit tracking of acid and calcite amounts
    """
    acid = np.zeros(nx)                        # Acid concentration in cell (mol/L)
    calcite = np.full(nx, float(calcite_init)) # Remaining calcite (mol/L)    
    f_mix = alpha / (dx + alpha)
    results = []

    for s in range(1, n_shifts + 1):
        # ① Advection (Acid movement)
        acid_adv = np.empty(nx)
        acid_adv[0] = acid_in      # Injected water contains a constant amount of acid
        acid_adv[1:] = acid[:-1]

        # ② Dispersion (Acid mixing)
        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

        # ③ Chemical Reaction (Acid neutralization by calcite and depletion)
        acid_new = acid_d.copy()
        for i in range(nx):
            if acid_new[i] > 0:
                # Amount neutralized is the minimum of present acid and remaining calcite
                reacted = min(acid_new[i], calcite[i])
                acid_new[i] -= reacted
                calcite[i] -= reacted

        acid = acid_new.copy()

        # ④ pH Calculation
        pH_out = np.empty(nx)
        for i in range(nx):
            if calcite[i] > 0:
                # If calcite remains, it's buffered and maintains initial pH
                pH_out[i] = pH_init
            else:
                # In depleted zones, pH drops according to acid concentration
                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

#  Scenario Parameter Setup

x = np.arange(0.5, 20.5) * 0.5   # Cell center coordinates (m)

# Scenario 1: Saline Front (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
)

# Scenario 2: Calcite Dissolution Front (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
)

# Scenario 3: Cd²⁺ Contamination (Conservative Tracer + Adsorption Retardation 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   # Retardation with R=5
)

#  Plotting

fig, axes = plt.subplots(1, 3, figsize=(15, 5))
fig.patch.set_facecolor("white")

# ---- Scenario 1: Saline Front ----
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} d")
ax.set(
    xlabel="Distance (m)",
    ylabel="Na (mM)",
    title="Scenario 1: Saline Front Movement\nNaCl 0.1 mol/L Injected",
)
ax.legend(title="Days Elapsed", 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")

# ---- Scenario 2: Calcite Dissolution Front ----
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} d")
ax.axhline(7.4, color="#9CA3AF", lw=1, ls=":", alpha=0.8)
ax.text(0.2, 7.55, "Initial 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, "Injected pH 4.5", color="#93C5FD", fontsize=8.5)
ax.set(
    xlabel="Distance (m)",
    ylabel="pH",
    title="Scenario 2: Calcite Dissolution Front\nAcidic Water pH 4.5 Injected",
)
ax.legend(title="Days Elapsed", 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")

# ---- Scenario 3: Cd Contamination & Adsorption Retardation ----
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="Conservative Tracer (R=1)")
ax.plot([], [], color="#DC2626", lw=2,            label="Cd (R=5, Adsorption Retardation)")
ax.set(
    xlabel="Distance (m)",
    ylabel="Concentration (umol/L)",
    title="Scenario 3: Cd Contamination Migration\nRetardation by Cation Exchange",
)
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: Concentration Profiles for 3 Scenarios",
             fontsize=13, y=0.99)
plt.tight_layout()
plt.savefig("transport_profiles.png", dpi=150,
            bbox_inches="tight", facecolor="white")
plt.show()
print("Saved transport_profiles.png")

Comparing the 3 Scenarios

The figure above compares the concentration profiles (front shape and movement) across the three scenarios. It visually demonstrates how the addition of chemical reactions and adsorption alters the transport mechanism compared to pure advection and dispersion.

1. Saline Front (Pure Advection & Dispersion)

The left panel (Scenario 1) shows the transport of a conservative (non-reactive) solute, NaCl. The defining characteristic is the smooth “S-shaped curve” of the concentration front. Advection pushes the center of the front downstream with the groundwater flow, while dispersion gradually smears the concentration gradient, widening the front over time. This illustrates the most fundamental solute transport behavior.

2. pH Front (Advection-Dispersion + Dissolution & Buffering)

The middle panel (Scenario 2) shows the pH front resulting from acidic water injection. Unlike the smooth S-curve in Scenario 1, the front here forms an extremely sharp, “vertical step (step function)”. This occurs because the pH is completely buffered and maintained at its initial value (7.4) until the calcite in the aquifer is fully dissolved. Once the calcite is depleted, the pH drops abruptly to the injected value (4.5). This sharp deformation of the front is a classic hallmark of a reactive transport front.

3. Cd²⁺ Contamination (Advection-Dispersion + Adsorption Retardation)

The right panel (Scenario 3) shows the migration of the heavy metal Cd²⁺, which undergoes cation exchange (adsorption) with the soil. Compared to a conservative tracer moving perfectly with the water (dashed line, R=1), the Cd²⁺ front (solid lines) exhibits a similar S-shape but remains much further behind. This clear “retardation” physically illustrates how adsorption reactions effectively slow down the migration velocity of contaminants.

TRANSPORT Block Essential Parameter Guide

Parameter Unit Meaning Typical Value
-cells — Spatial divisions. More = higher accuracy & time 10~100
-length m Length of each cell 0.01~10
-shifts — Time steps (= number of times fluid passes a cell) 20~200
-time_step s Time per step. Calculated as length / velocity 86400 (1 day)
-dispersivity m Dispersivity α. Empirical rule: ~1/10 of migration distance 0.01~1
-stagnant — Dual-porosity model (mobile/immobile zones) Important in fractured rock

Summary: Series Completion

PHREEQC Series Completion Map #1~#6 Basics & Speciation Equilibrium, Mixing, Redox "Chemistry at 1 point" #7~#10 Solubility, Activity, SI Python Visualization "Deepening Thermodynamics" #11 Reaction Path Modeling REACTION Block "Temporal Changes" #12 ← Current TRANSPORT Block Advection, Dispersion, Reaction "Space × Time × Chemistry" Series Completion Now capable of fully simulating the "groundwater journey" in PHREEQC
TipTo those who have made it this far

Starting from installation in #1, walking through speciation, mixing, equilibrium, redox, solubility, activity, saturation index, reaction paths, and now advection-dispersion—you have traversed the core knowledge necessary to perform practical geochemical simulations in PHREEQC.

Your next steps could open doors to inverse modeling (the INVERSE_MODELING block) using actual field data (water quality analysis values), or moving forward to FloPy × PHREEQC coupling (like the hydrothermal simulation showcased at the beginning of the series).

References

Appelo, CAJ, and Dieke Postma. 2005. Geochemistry, Groundwater and Pollution. Second. Balkema, Rotterdam, p. 634.
Parkhurst, David L, and 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, and 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.

DeepFlow | Science beneath the surface