Skip to content

New sorption formulation#2

Open
Manangka wants to merge 16 commits intodevelopfrom
new_sorption_formulation
Open

New sorption formulation#2
Manangka wants to merge 16 commits intodevelopfrom
new_sorption_formulation

Conversation

@Manangka
Copy link
Copy Markdown
Owner

@Manangka Manangka commented Dec 24, 2025

MODFLOW supports modeling of sorption in the transport model using Isotherms.
Isotherms are non-linear functions that establishes a relation between the sorbed concentration and the solute concentration in the groundwater.
$f_{sorption} =-f_m \rho_b \frac{\partial \left(S_w \overline{C}\right)}{\partial t}$

Currently MODFLOW is using the midpoint point rule to discretize the sorption term.
$f_{sorption} = -f_m \rho_b\left(\left(S_w \frac{\partial \overline{C}}{\partial C}\right)^{n+0.5} \frac{C^{n+1} - C^{n}}{\Delta t} +\overline{C}^{n+0.5}\frac{S_w^{n+1} - S_w^{n}}{\Delta t}\right)$

However some investigation showed that this formulation isn't mass conservative.
A test (test_gwt_mst08) has been added to this PR that demonstrates this.

This PR reformulates the sorption term such that it becomes mass conservative. It does this by linearizing the $S_w \overline{C}$ term and freezing the jacobian at the last picard iteration.
$f_{sorption}=-f_m \rho_b \frac{S_w^{\left(k\right)}\left(\frac{\partial \overline{C}}{\partial C}\right)^{\left(k\right)} \left(C^{n+1} - C^{\left(k\right)}\right) + \overline{C}^{\left(k\right)}S_w^{n+1} - \left(S_w \overline{C}\right)^{n}}{\Delta t}$


Further a modification to the Freundlich Isotherm is added to improve stability.
The isotherm becomes problematic for values of a that are smaller than 0.
For those values the derivative goes to infinity causing convergence issues.

The Freundlich Isotherm is given by:
$\overline{C} = K_f C^a$

with its derivative
$\frac{\partial \overline{C}}{\partial C} = a K_f C^{a-1}$

To improve stability a small epsilon is added to keep the derivative finite.
$\overline{C} = K_f \left(C + \epsilon\right)^{a} - K_f \epsilon^{a}$

Note that in the equation above the value of $\overline{C}$ at 0 equals the original Freundlich equation

The corresponding derivative now becomes finite:
$\frac{\partial \overline{C}}{\partial C} = a K_f \left(C + \epsilon\right)^{a - 1}$

Checklist of items for pull request

  • Replaced section above with description of pull request
  • Closed issue #xxxx
  • Referenced issue or pull request #xxxx
  • Added new test or modified an existing test
  • Ran ruff on new and modified python scripts in .doc, autotests, doc, distribution, pymake, and utils subdirectories.
  • Formatted new and modified Fortran source files with fprettify
  • Added doxygen comments to new and modified procedures
  • Updated meson files, makefiles, and Visual Studio project files for new source files
  • Updated definition files
  • Updated develop.toml with a plain-language description of the bug fix, change, feature; required for changes that may affect users
  • Updated input and output guide
  • Removed checklist items not relevant to this pull request

For additional information see instructions for contributing and instructions for developing.

@Manangka Manangka changed the base branch from refactor_isotherm to develop March 23, 2026 13:48
@Manangka Manangka marked this pull request as ready for review March 24, 2026 09:54
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Labels

None yet

Projects

None yet

Development

Successfully merging this pull request may close these issues.

1 participant