Introduction: What is Inverse Modeling?
So far in this series, we have performed calculations such as “What kind of water (C) do we get if we mix water (A) and water (B)?” or “What happens if we dissolve limestone in this water?” These are calculations that predict the “result” from the “cause,” which is known as Forward Modeling.
However, in actual groundwater investigations and hot spring field sites, the situation is the exact opposite. All we have in our hands is the “resulting water quality data (C)” that has already mixed and sprung forth, and we cannot directly observe what happened deep underground.
“What kind of water (A) and water (B) originally mixed, at what ratio, and what minerals dissolved (or precipitated) during that time to create this hot spring (C)?”
The calculation that deduces the “cause (recipe)” from this “result” is called Inverse Modeling.
Difference Between Forward and Inverse Modeling
Look at the figure below. The fundamental difference between the two approaches is evident at a glance.
Inverse analysis retrospectively estimates underground reaction processes from the observed final water quality.
As the figure shows, in forward modeling, the arrows flow from left to right (cause → result), while in inverse modeling, the arrows go from right to left (result → cause). PHREEQC automatically performs this reverse deduction by solving simultaneous mass conservation equations.
🍹 Analogy: Reverse-Engineering a Mixed Juice Recipe
Let’s compare the concept of inverse modeling to analyzing a mixed juice.
You have a “finished mixed juice (final water quality)” in front of you. Upon analysis, it contains 10g of sugar and 50mg of vitamin C. In your kitchen, you have the following ingredients (initial water quality/minerals): 1. Apple juice (2g sugar, 5mg vitamin C / 100ml) 2. Orange juice (1g sugar, 20mg vitamin C / 100ml) 3. Gum syrup (pure sugar)
Now, if you ask PHREEQC, “How was this mixed juice made?” (by running inverse modeling), PHREEQC will solve a system of equations based on the law of mass conservation and provide an answer (model) like this:
PHREEQC’s Answer (Computational Model) “If you mix 200ml of apple juice and 200ml of orange juice, and finally dissolve 4g of gum syrup, it will match the ingredients perfectly!”
In this way, inverse modeling is a tool like a master detective in geochemistry that automatically calculates the “mixing ratio of waters” and the “amount of mineral change (dissolved or precipitated)” required to bridge the gap between the “starting water” and the “arriving water.”
This time, let’s actually use inverse modeling to unravel an underground drama, taking as an example the process of “rainwater infiltrating into the ground and becoming groundwater in a carbonate aquifer,” which is common in groundwater research.
Practical Case Study: Evolution from Rainwater to Carbonate Aquifer Groundwater
To experience the power of inverse modeling, we will analyze a scenario where “rainwater infiltrates underground, passes through a carbonate rock (such as limestone) aquifer, and its water quality evolves.”
Let’s assume the following reactions occurred during the groundwater flow process:
- Takes in CO₂ from the soil.
- Calcite and Dolomite from the limestone dissolve.
- Trace amounts of Gypsum and Halite in the strata dissolve.
PHREEQC Code and Water Quality Data
First, let’s look at the input file (PHREEQC code) to calculate these. We input the initial rainwater as the starting solution (SOLUTION 1) and the observed groundwater as the final solution (SOLUTION 2).
TITLE Carbonate aquifer inverse model
# 1. Starting Solution: Rainwater
SOLUTION 1 Rainwater
temp 15
pH 5.6
units mg/L
Ca 2
Mg 0.5
Na 2
Cl 3.08
S(6) 2
Alkalinity 5 as HCO3
# 2. Final Solution: Aquifer Groundwater
SOLUTION 2 Groundwater
temp 15
pH 7.3
units mg/L
Ca 55.2
Mg 4.87
Na 10
Cl 15.42
S(6) 21.2
Alkalinity 165 as HCO3
# Inverse Modeling Settings
INVERSE_MODELING 1
-solutions 1 2
-phases
Calcite
Dolomite
Gypsum
Halite
CO2(g)
-uncertainty 0.10
ENDRegarding the Error (uncertainty) Setting To prevent errors (such as No models found) due to analytical errors in practice, the -uncertainty (tolerance) is broadened from the default 0.05 (5%) to 0.10 (10%). If the charge balance is off in actual analytical data, fine-tuning may be necessary, such as widening this further or adding the charge keyword to specific ions.
How to Read INVERSE_MODELING Results and Geochemical Interpretation
When you run the code above, PHREEQC performs mass balance calculations and returns several output blocks. Let’s focus on “Phase mole transfers,” the core part of the inverse analysis, to interpret the results.
Extracting the Execution Results
Phase mole transfers: Minimum Maximum Formula (Approximate SI in solution 1, 2 at 288 K, 1 atm)
Calcite 9.481e-04 0.000e+00 0.000e+00 CaCO3 ( -4.83, -0.31)
Dolomite 1.798e-04 0.000e+00 0.000e+00 CaMg(CO3)2 ( -9.95, -1.37)
Gypsum 1.999e-04 0.000e+00 0.000e+00 CaSO4:2H2O ( -4.47, -2.26)
Halite 3.482e-04 0.000e+00 0.000e+00 NaCl ( -9.70, -8.35)
CO2(g) 1.089e-03 0.000e+00 0.000e+00 CO2 ( -1.92, -2.14)Note: The unit of the output value is mol/kgw (moles per 1 kg of water). For example, 9.481e-04 means \(9.481 \times 10^{-4}\), or 0.000948 mol (approx. 0.95 mmol).
Meaning of the Sign (Plus/Minus)
The most important thing when reading Phase mole transfers is the sign (plus/minus).
| Sign | Physical Meaning | Impact on Water Quality |
|---|---|---|
| + (Positive) | Minerals dissolve, or gases are absorbed | The component increases in the water |
| - (Negative) | Minerals precipitate, or gases degas (release) | The component decreases in the water |
In this output, all values are positive (+). This means that all specified minerals and gases resulted in being “dissolved/absorbed into the water.”
Detailed Interpretation of Each Mineral
Let’s explain what kind of geochemical story each result is telling.
① Calcite: Dissolved approx. 0.95 mmol - Reaction: \(\text{CaCO}_3 \rightarrow \text{Ca}^{2+} + \text{CO}_3^{2-}\) - Interpretation: This is the main factor for the increase in Ca concentration and alkalinity (HCO3) in the groundwater. Looking at the SI (Saturation Index) on the far right, it approaches 0 from Solution 1 (-4.83) to Solution 2 (-0.31), which is consistent with the natural phenomenon of “Calcite dissolving and reaching a saturated state.”
② Dolomite: Dissolved approx. 0.18 mmol - Reaction: \(\text{CaMg}(\text{CO}_3)_2 \rightarrow \text{Ca}^{2+} + \text{Mg}^{2+} + 2\text{CO}_3^{2-}\) - Interpretation: This explains the observed increase in Mg concentration. Since Dolomite is also undersaturated (SI: -9.95 → -1.37), the result that a dissolution reaction occurred is thermodynamically reasonable.
③ Gypsum: Dissolved approx. 0.20 mmol - Reaction: \(\text{CaSO}_4\cdot 2\text{H}_2\text{O} \rightarrow \text{Ca}^{2+} + \text{SO}_4^{2-} + 2\text{H}_2\text{O}\) - Interpretation: This compensates for the increase in sulfate ions (SO4) contained in the groundwater. It suggests that gypsum lenses interbedded in the strata or trace sulfate minerals have dissolved.
④ Halite: Dissolved approx. 0.35 mmol - Interpretation: Explains the increase in Na and Cl.
⑤ CO2(g): Absorbed approx. 1.1 mmol - Reaction: \(\text{CO}_2 + \text{H}_2\text{O} \rightarrow \text{H}_2\text{CO}_3\) - Interpretation: This indicates that when rainwater passed through the soil, it absorbed a large amount of soil CO2 from plant root respiration and microbial activity. This CO2 dissolved in the water to become carbonic acid, playing the role of an “acid” that dissolves the aforementioned Calcite and Dolomite.
What to Do When the Model Fails (No models found)
In this example, the data is intentionally adjusted so that the mass balance matches, but when actual observation data is put in, it very often results in Number of models found: 0 (no solution found) at first.
The main causes and workarounds are as follows:
- Missing Candidate Minerals
- For example, if Na and Cl are increasing but you forget to put Halite in
-phases, the calculations won’t balance and it will result in an error. Look at the changes in water quality and add minerals that supplement the missing components.
- For example, if Na and Cl are increasing but you forget to put Halite in
- Analytical Data Errors (Charge Imbalance)
- If the charge balance is broken due to analytical errors of ions, PHREEQC cannot find a solution. In this case, measures such as widening
-uncertaintyto 0.10–0.20 or having PHREEQC fine-tune the charge balance by attaching thechargekeyword to a specific ion (e.g., Cl or Alkalinity) are necessary.
- If the charge balance is broken due to analytical errors of ions, PHREEQC cannot find a solution. In this case, measures such as widening
Case Study 2: Seawater Intrusion and Reverse Ion Exchange in Coastal Groundwater
A model very frequently used in coastal groundwater research is one that elucidates “water-rock reactions and ion exchange occurring after mixing” in addition to “mixing of freshwater and seawater.”
Here, using observation data from a hypothetical coastal well, let’s analyze water quality evolution that cannot be explained by simple mixing alone through inverse modeling.
1. Observation Data and Estimation of “Seawater Mixing Ratio by Cl⁻”
We have prepared the following three water quality data sets (freshwater, seawater, and observed coastal groundwater).
| Parameter (mg/L) | Freshwater | Seawater | Coastal_Well |
|---|---|---|---|
| pH | 7.2 | 8.1 | 7.4 |
| Na | 20 | 10500 | 500 |
| Ca | 40 | 400 | 180 |
| Mg | 10 | 1300 | 45 |
| Cl | 30 | 19000 | 1000 |
| SO4 | 15 | 2700 | 120 |
| Alkalinity | 180 | 140 | 240 |
In practice, we first roughly estimate the mixing ratio of seawater \(f\) using chloride ion (Cl⁻), which is a conservative tracer.
\[ f = \frac{Cl_{sample} - Cl_{fresh}}{Cl_{sea} - Cl_{fresh}} = \frac{1000 - 30}{19000 - 30} \approx 0.051 \]
Where: - \(f\): Mixing ratio of seawater (fraction) - \(Cl_{sample}\): Chloride concentration in the observed groundwater sample (mg/L) - \(Cl_{fresh}\): Chloride concentration in the freshwater end-member (mg/L) - \(Cl_{sea}\): Chloride concentration in the seawater end-member (mg/L)
In other words, it can be estimated that about 5.1% is seawater and about 94.9% is freshwater.
2. Limits of Simple Mixing Models and Suspicion of “Reverse Ion Exchange”
If they were simply mixed, what should the concentrations of Na and Ca be?
- Theoretical value of Na = (0.949 × 20) + (0.051 × 10500) = 554 mg/L
- Theoretical value of Ca = (0.949 × 40) + (0.051 × 400) = 58 mg/L
However, the observed values are Na=500 mg/L (less than the theoretical value) and Ca=180 mg/L (abnormally higher than the theoretical value). This is strong evidence that “Reverse Ion Exchange”, famous in coastal aquifers, is occurring. It is thought that large amounts of seawater-derived Na⁺ adsorbed onto clay minerals in the strata, and instead, Ca²⁺ from the strata was pushed out into the water.
3. Inverse Modeling in PHREEQC
To quantitatively prove this assumption, we input the following code into PHREEQC. The key is to define virtual exchange minerals (NaX, CaX2) in PHASES to account for cation exchange.
TITLE Coastal groundwater
# Definition of virtual Cation Exchange Phases
PHASES
NaX
NaX = Na+ + X-
log_k 0.0
CaX2
CaX2 = Ca+2 + 2X-
log_k 0.0
# 1. Freshwater End-Member
SOLUTION 1 Freshwater
units mg/L
pH 7.2
Na 20; Ca 40; Mg 10; Cl 30; S(6) 15; Alkalinity 180 as HCO3
# 2. Seawater End-Member
SOLUTION 2 Seawater
units mg/L
pH 8.1
Na 10500; Ca 400; Mg 1300; Cl 19000; S(6) 2700; Alkalinity 140 as HCO3
# 3. Observed Groundwater
SOLUTION 3 Coastal_Well
units mg/L
pH 7.4
Na 500; Ca 180; Mg 45; Cl 1000; S(6) 120; Alkalinity 240 as HCO3
# Inverse Modeling
INVERSE_MODELING 1
-solutions 1 2 3
-phases
Calcite
Dolomite
Gypsum
CO2(g)
NaX
CaX2
-mineral_water true
-uncertainty 0.10
-balances
Cl
ENDTwo Important Options
-mineral_water true: Normal inverse analysis solves “Solution 1 → Final point,” but including this option explicitly directs it to solve a reaction model involving mixing: “Solution 1 + Solution 2 + Mineral Reactions = Final point.” This is essential for coastal groundwater analysis.-balances Cl: Instructs that for the specified component (here, chloride ion Cl), the mass balance is closed “solely by mass conservation due to simple mixing,” regardless of mineral dissolution/precipitation. Since Cl is treated as a conservative tracer that does not react with rocks, including this setting makes the mixing ratio calculation more stable.
4. Interpretation of Execution Results
Let’s decipher the Phase mole transfers of the applicable model from the calculation output.
Solution fractions: Fraction
Solution 1 (Freshwater) 9.499e-01
Solution 2 (Seawater) 5.014e-02
Phase mole transfers: Moles
Calcite 3.273e-03
Dolomite -1.311e-03
Gypsum -3.554e-04
CO2(g) 2.748e-04
NaX -2.738e-03
CaX2 1.369e-03① Proof of Seawater Mixing Ratio The proportion of Solution 2 (seawater) was 5.014e-02 (about 5.0%), perfectly matching the 5.1% calculated manually from Cl⁻ at the beginning.
② Quantification of Reverse Ion Exchange NaX is negative (adsorption), and CaX2 is positive (desorption). Specifically, it was proven that approximately 2.74 mmol of Na⁺ adsorbed from the water to the strata, and instead, approximately 1.37 mmol of Ca²⁺, equivalent to half that charge, was released into the water.
③ Complex Reactions of Carbonate Minerals and the Truth About Ca Increase It became clear that the observed “abnormal increase in Ca concentration (theoretical 58 mg/L → observed 180 mg/L)” was significantly contributed to not only by the release of Ca²⁺ from reverse ion exchange but also by the dissolution of Calcite (+3.27 mmol). Furthermore, as a large amount of Ca²⁺ was supplied to the water, the water tended to become supersaturated with respect to carbonates and sulfates containing Mg. Consequently, while Calcite dissolved, Dolomite and Gypsum shifted to precipitation (-). This complex chain reaction, which could never be explained by simple mixing or “reverse ion exchange alone,” has been brilliantly reconciled and clarified by PHREEQC’s mass balance calculations.
Conclusion
Let’s summarize the key points of inverse modeling introduced this time.
| Item | Details |
|---|---|
| Purpose | To reverse-calculate mixing ratios and underground mineral reactions from observed final water quality. |
| Core of Output | The sign of Phase mole transfers (+ for dissolution, - for precipitation) and its amount. |
| Selecting from Multiple Models | Prioritize geological basis + minimum residuals + minimal model. |
| Setting Uncertainty | Start from 0.10 to 0.20 for problems with large errors. |
| Unique Strength | Can thermodynamically verify and prove “how water quality evolved.” |
Inverse modeling is a detective game where you hypothesize, “Let’s suspect Calcite precipitation” or “Gypsum must have dissolved here,” and proceed through trial and error while feeding models. However, when a solution is finally found, it is not merely a calculation result; it becomes a solid proof that you have “logically deciphered the invisible underground reaction processes.”
Next time, we will step into the world of KINETICS and REACTION (Reaction Path Modeling), which simulate how these chemical reactions, such as “precipitation” and “dissolution,” progress over time. Stay tuned!
References
- #1 Installation and Initial Calculation
- #2 Analyzing Seawater with Speciation
- #3 Mixing and EQUILIBRIUM_PHASES
- #4 Considering Reaction Rates with KINETICS
- #5 Mixing Groundwater and Seawater
- #6 Pyrite Oxidation and AMD
- #7 Solubility Diagrams (Gibbsite)
- #8 Analyzing and Visualizing PHREEQC Results with Python
- #9 Ionic Strength and Activity Coefficients
- #10 Mastering the Saturation Index (SI)
- #11 Reaction Path Modeling
- #12 Advection-Dispersion Model
- #13 Redox Sequences
- #14 Groundwater Diagnosis of Nitrate Pollution
- [#15 Inverse Modeling — Deducing “Invisible Underground Processes” from Water Quality Data]
DeepFlow | Science beneath the surface