Skip to content

mibiremo API reference

phreeqc

Python interface to PhreeqcRM for biogeochemical reactive transport modeling.

PhreeqcRM is a reaction module developed by the U.S. Geological Survey (USGS) for coupling biogeochemical calculations with transport models. It provides a high-performance interface to PHREEQC modeling capabilities for reactive transport simulations in environmental and hydrological applications.

PhreeqcRM enables:

  • Multi-threaded biogeochemical calculations for large-scale transport models
  • Equilibrium and kinetic biogeochemical reactions in porous media
  • Parallel processing for computationally intensive reactive transport

This interface wraps the phreeqcrm pip package, exposing four high-level convenience methods (create, initialize_phreeqc, run_initial_from_file, get_selected_output_df) and a public self.rm attribute (the underlying phreeqcrm.PhreeqcRM instance) for advanced users who need direct access to the full PhreeqcRM API.

PhreeqcRM documentation and source code can be found at:

Last revision: 10/04/2026

PhreeqcRM

Python interface to PhreeqcRM for geochemical reactive transport modeling.

This class facilitates coupling between transport codes and geochemical reaction calculations by managing multiple reaction cells, each representing a grid cell in the transport model. The PhreeqcRM approach allows efficient parallel processing of geochemical calculations across large spatial domains.

The class handles
  • Creation and initialization of PhreeqcRM instances
  • Loading thermodynamic databases (PHREEQC format)
  • Setting up initial chemical conditions from input files
  • Running equilibrium and kinetic geochemical reactions
  • Transferring concentrations between transport and reaction modules
  • Managing porosity, saturation, temperature, and pressure fields
  • Retrieving calculated properties and concentrations
Typical workflow
  1. Create instance and call create() method to initialize with grid size
  2. Load thermodynamic database with initialize_phreeqc()
  3. Set initial conditions with run_initial_from_file()
  4. In transport time loop:
  5. Transfer concentrations to reaction module with rm.SetConcentrations()
  6. Advance time with rm.SetTime() and rm.SetTimeStep()
  7. Run reactions with rm.RunCells()
  8. Retrieve updated concentrations with rm.GetConcentrations()

Attributes:

Name Type Description
nxyz int

Number of grid cells in the reactive transport model.
n_threads int

Number of threads for parallel geochemical processing.
rm PhreeqcRM

The underlying PhreeqcRM instance. Use this attribute to access the full PhreeqcRM API directly (e.g., rm.SetTime(t), rm.GetConcentrations()).
components ndarray

Array of component names for transport.
species ndarray

Array of aqueous species names in the system.

Examples:

See page Examples for usage examples.

Source code in mibiremo/phreeqc.py 
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
class PhreeqcRM:
    """Python interface to PhreeqcRM for geochemical reactive transport modeling.

    This class facilitates coupling between transport codes and geochemical
    reaction calculations by managing multiple reaction cells, each representing
    a grid cell in the transport model. The PhreeqcRM approach allows efficient
    parallel processing of geochemical calculations across large spatial domains.

    The class handles:
        - Creation and initialization of PhreeqcRM instances
        - Loading thermodynamic databases (PHREEQC format)
        - Setting up initial chemical conditions from input files
        - Running equilibrium and kinetic geochemical reactions
        - Transferring concentrations between transport and reaction modules
        - Managing porosity, saturation, temperature, and pressure fields
        - Retrieving calculated properties and concentrations

    Typical workflow:
        1. Create instance and call create() method to initialize with grid size
        2. Load thermodynamic database with initialize_phreeqc()
        3. Set initial conditions with run_initial_from_file()
        4. In transport time loop:
           - Transfer concentrations to reaction module with rm.SetConcentrations()
           - Advance time with rm.SetTime() and rm.SetTimeStep()
           - Run reactions with rm.RunCells()
           - Retrieve updated concentrations with rm.GetConcentrations()

    Attributes:
        nxyz (int): Number of grid cells in the reactive transport model.
        n_threads (int): Number of threads for parallel geochemical processing.
        rm (phreeqcrm.PhreeqcRM): The underlying PhreeqcRM instance. Use this
            attribute to access the full PhreeqcRM API directly (e.g.,
            ``rm.SetTime(t)``, ``rm.GetConcentrations()``).
        components (numpy.ndarray): Array of component names for transport.
        species (numpy.ndarray): Array of aqueous species names in the system.

    Examples:
        See page [Examples](examples.md) for usage examples.
    """

    def __init__(self):
        """Initialize PhreeqcRM instance.

        Creates a new PhreeqcRM object with default values. The instance must be
        created using the create() method before it can be used for calculations.
        """
        self._initialized = False
        self.nxyz = 1
        self.n_threads = 1
        self.rm = None
        self.components = None
        self.species = None

    def create(self, nxyz=1, n_threads=1) -> None:
        """Creates a PhreeqcRM reaction module instance.

        Initializes the PhreeqcRM library and creates a reaction module with
        the specified number of grid cells and threads.
        This method must be called before any other PhreeqcRM operations.

        Args:
            nxyz (int, optional): Number of grid cells in the model. Must be
                positive. Defaults to 1.
            n_threads (int, optional): Number of threads for parallel processing.
                Use -1 for automatic detection of CPU count. Defaults to 1.

        Raises:
            RuntimeError: If PhreeqcRM instance creation fails.

        Examples:
            >>> rm = PhreeqcRM()
            >>> rm.create(nxyz=100, n_threads=4)
        """
        if n_threads == -1:
            n_threads = os.cpu_count()

        self.n_threads = n_threads
        self.nxyz = nxyz
        try:
            self.rm = phreeqcrm.PhreeqcRM(nxyz, n_threads)
            self._initialized = True
        except Exception as e:
            raise RuntimeError(f"Failed to create PhreeqcRM instance: {e}")

    def initialize_phreeqc(
        self,
        database_path,
        units_solution=2,
        units=1,
        porosity=1.0,
        saturation=1.0,
        multicomponent=True,
    ) -> None:
        """Initialize PhreeqcRM with database and default parameters.

        Loads a thermodynamic database and sets up the PhreeqcRM instance with
        standard parameters for geochemical calculations. This is a convenience
        method that handles common initialization tasks.

        Args:
            database_path (str): Path to the PHREEQC database file (.dat format).
                Common databases include phreeqc.dat, Amm.dat, pitzer.dat.
            units_solution (int, optional): Units for solution concentrations.
                1 = mol/L, 2 = mmol/L, 3 = μmol/L. Defaults to 2.
            units (int, optional): Units for other phases (Exchange, Surface,
                Gas, Solid solutions, Kinetics). Defaults to 1.
            porosity (float, optional): Porosity value assigned to all cells.
                Must be between 0 and 1. Defaults to 1.0.
            saturation (float, optional): Saturation value assigned to all cells.
                Must be between 0 and 1. Defaults to 1.0.
            multicomponent (bool, optional): Enable multicomponent diffusion
                by saving species concentrations. Defaults to True.

        Raises:
            RuntimeError: If the PhreeqcRM instance is not initialized or if
                the database fails to load.

        Examples:
            >>> rm = PhreeqcRM()
            >>> rm.create(nxyz=100)
            >>> rm.initialize_phreeqc("phreeqc.dat", units_solution=1)
        """
        if not self._initialized:
            raise RuntimeError("PhreeqcRM instance not initialized. Call create() first.")

        status = self.rm.LoadDatabase(database_path)
        if status < 0:
            raise RuntimeError(f"Failed to load Phreeqc database (error code: {status})")

        self.rm.SetComponentH2O(False)
        self.rm.SetRebalanceFraction(0.5)
        self.rm.SetUnitsSolution(units_solution)
        self.rm.SetUnitsPPassemblage(units)
        self.rm.SetUnitsExchange(units)
        self.rm.SetUnitsSurface(units)
        self.rm.SetUnitsGasPhase(units)
        self.rm.SetUnitsSSassemblage(units)
        self.rm.SetUnitsKinetics(units)
        self.rm.SetPorosity(porosity * np.ones(self.nxyz))
        self.rm.SetSaturationUser(saturation * np.ones(self.nxyz))
        self.rm.SetFilePrefix("phr")
        self.rm.OpenFiles()
        if multicomponent:
            self.rm.SetSpeciesSaveOn(True)

    def run_initial_from_file(self, pqi_file, ic):
        """Set up initial conditions from PHREEQC input file and initial conditions array.

        Loads initial geochemical conditions by running a PHREEQC input file and
        mapping the defined solutions, phases, and other components to the grid cells.
        This method also retrieves component and species information for later use.

        Args:
            pqi_file (str): Path to the PHREEQC input file (.pqi format) containing
                definitions for solutions, equilibrium phases, exchange, surface,
                gas phases, solid solutions, and kinetic reactions.
            ic (numpy.ndarray): Initial conditions array with shape (nxyz, 7) where
                each row corresponds to a grid cell and columns represent:
                - Column 0: Solution ID
                - Column 1: Equilibrium phase ID
                - Column 2: Exchange ID
                - Column 3: Surface ID
                - Column 4: Gas phase ID
                - Column 5: Solid solution ID
                - Column 6: Kinetic reaction ID
                Use -1 for unused components.

        Raises:
            RuntimeError: If the PHREEQC input file fails to run.
            ValueError: If initial conditions array has incorrect shape or cannot
                be converted to integer array.

        Examples:
            >>> import numpy as np
            >>> ic = np.array([[1, -1, -1, -1, -1, -1, -1]])  # Only solution 1
            >>> rm.run_initial_from_file("initial.pqi", ic)
        """
        status = self.rm.RunFile(True, True, True, pqi_file)
        if status < 0:
            raise RuntimeError(f"Failed to run Phreeqc input file (error code: {status})")

        if ic.shape != (self.nxyz, 7):
            raise ValueError(f"Initial conditions array must have shape ({self.nxyz}, 7), got {ic.shape}")

        if not isinstance(ic, np.ndarray):
            try:
                ic = np.array(ic).astype(np.int32)
            except Exception as e:
                raise ValueError("Initial conditions must be convertible to a numpy array of integers") from e

        ic1 = ic.flatten("F").astype(np.int32)
        self.rm.InitialPhreeqc2Module(ic1)

        self.rm.FindComponents()
        self.components = np.array(list(self.rm.GetComponents()))
        self.species = np.array(list(self.rm.GetSpeciesNames()))

        self.rm.SetTime(0.0)
        self.rm.SetTimeStep(0.0)
        self.rm.RunCells()

    def get_selected_output_df(self) -> pd.DataFrame:
        """Retrieve selected output data as a pandas DataFrame.

        Extracts the current selected output data from PhreeqcRM and formats it
        as a pandas DataFrame with appropriate column headers. Selected output
        typically includes calculated properties like pH, pe, ionic strength,
        activities, saturation indices, and user-defined calculations.

        Returns:
            pandas.DataFrame: DataFrame containing selected output data with
                rows representing grid cells and columns representing the
                selected output variables defined in the PHREEQC input.

        Examples:
            >>> df = rm.get_selected_output_df()
            >>> print(df.columns)  # Show available output variables
            >>> print(df['pH'])     # Access pH values for all cells
        """
        headings = list(self.rm.GetSelectedOutputHeadings())
        ncolsel = len(headings)
        so = np.array(self.rm.GetSelectedOutput())
        return pd.DataFrame(so.reshape(ncolsel, self.nxyz).T, columns=headings)

__init__()

Initialize PhreeqcRM instance.

Creates a new PhreeqcRM object with default values. The instance must be created using the create() method before it can be used for calculations.

Source code in mibiremo/phreeqc.py 
74
75
76
77
78
79
80
81
82
83
84
85
def __init__(self):
    """Initialize PhreeqcRM instance.

    Creates a new PhreeqcRM object with default values. The instance must be
    created using the create() method before it can be used for calculations.
    """
    self._initialized = False
    self.nxyz = 1
    self.n_threads = 1
    self.rm = None
    self.components = None
    self.species = None

create(nxyz=1, n_threads=1)

Creates a PhreeqcRM reaction module instance.

Initializes the PhreeqcRM library and creates a reaction module with the specified number of grid cells and threads. This method must be called before any other PhreeqcRM operations.

Parameters:

Name Type Description Default
nxyz int

Number of grid cells in the model. Must be positive. Defaults to 1.

 1
n_threads int

Number of threads for parallel processing. Use -1 for automatic detection of CPU count. Defaults to 1.

 1

Raises:

Type Description
RuntimeError

If PhreeqcRM instance creation fails.

Examples:

>>> rm = PhreeqcRM()
>>> rm.create(nxyz=100, n_threads=4)
Source code in mibiremo/phreeqc.py 
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
def create(self, nxyz=1, n_threads=1) -> None:
    """Creates a PhreeqcRM reaction module instance.

    Initializes the PhreeqcRM library and creates a reaction module with
    the specified number of grid cells and threads.
    This method must be called before any other PhreeqcRM operations.

    Args:
        nxyz (int, optional): Number of grid cells in the model. Must be
            positive. Defaults to 1.
        n_threads (int, optional): Number of threads for parallel processing.
            Use -1 for automatic detection of CPU count. Defaults to 1.

    Raises:
        RuntimeError: If PhreeqcRM instance creation fails.

    Examples:
        >>> rm = PhreeqcRM()
        >>> rm.create(nxyz=100, n_threads=4)
    """
    if n_threads == -1:
        n_threads = os.cpu_count()

    self.n_threads = n_threads
    self.nxyz = nxyz
    try:
        self.rm = phreeqcrm.PhreeqcRM(nxyz, n_threads)
        self._initialized = True
    except Exception as e:
        raise RuntimeError(f"Failed to create PhreeqcRM instance: {e}")

get_selected_output_df()

Retrieve selected output data as a pandas DataFrame.

Extracts the current selected output data from PhreeqcRM and formats it as a pandas DataFrame with appropriate column headers. Selected output typically includes calculated properties like pH, pe, ionic strength, activities, saturation indices, and user-defined calculations.

Returns:

Type Description
DataFrame

pandas.DataFrame: DataFrame containing selected output data with rows representing grid cells and columns representing the selected output variables defined in the PHREEQC input.

Examples:

>>> df = rm.get_selected_output_df()
>>> print(df.columns)  # Show available output variables
>>> print(df['pH'])     # Access pH values for all cells
Source code in mibiremo/phreeqc.py 
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
def get_selected_output_df(self) -> pd.DataFrame:
    """Retrieve selected output data as a pandas DataFrame.

    Extracts the current selected output data from PhreeqcRM and formats it
    as a pandas DataFrame with appropriate column headers. Selected output
    typically includes calculated properties like pH, pe, ionic strength,
    activities, saturation indices, and user-defined calculations.

    Returns:
        pandas.DataFrame: DataFrame containing selected output data with
            rows representing grid cells and columns representing the
            selected output variables defined in the PHREEQC input.

    Examples:
        >>> df = rm.get_selected_output_df()
        >>> print(df.columns)  # Show available output variables
        >>> print(df['pH'])     # Access pH values for all cells
    """
    headings = list(self.rm.GetSelectedOutputHeadings())
    ncolsel = len(headings)
    so = np.array(self.rm.GetSelectedOutput())
    return pd.DataFrame(so.reshape(ncolsel, self.nxyz).T, columns=headings)

initialize_phreeqc(database_path, units_solution=2, units=1, porosity=1.0, saturation=1.0, multicomponent=True)

Initialize PhreeqcRM with database and default parameters.

Loads a thermodynamic database and sets up the PhreeqcRM instance with standard parameters for geochemical calculations. This is a convenience method that handles common initialization tasks.

Parameters:

Name Type Description Default
database_path str

Path to the PHREEQC database file (.dat format). Common databases include phreeqc.dat, Amm.dat, pitzer.dat.

 required
units_solution int

Units for solution concentrations. 1 = mol/L, 2 = mmol/L, 3 = μmol/L. Defaults to 2.

 2
units int

Units for other phases (Exchange, Surface, Gas, Solid solutions, Kinetics). Defaults to 1.

 1
porosity float

Porosity value assigned to all cells. Must be between 0 and 1. Defaults to 1.0.

 1.0
saturation float

Saturation value assigned to all cells. Must be between 0 and 1. Defaults to 1.0.

 1.0
multicomponent bool

Enable multicomponent diffusion by saving species concentrations. Defaults to True.

 True

Raises:

Type Description
RuntimeError

If the PhreeqcRM instance is not initialized or if the database fails to load.

Examples:

>>> rm = PhreeqcRM()
>>> rm.create(nxyz=100)
>>> rm.initialize_phreeqc("phreeqc.dat", units_solution=1)
Source code in mibiremo/phreeqc.py 
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
def initialize_phreeqc(
    self,
    database_path,
    units_solution=2,
    units=1,
    porosity=1.0,
    saturation=1.0,
    multicomponent=True,
) -> None:
    """Initialize PhreeqcRM with database and default parameters.

    Loads a thermodynamic database and sets up the PhreeqcRM instance with
    standard parameters for geochemical calculations. This is a convenience
    method that handles common initialization tasks.

    Args:
        database_path (str): Path to the PHREEQC database file (.dat format).
            Common databases include phreeqc.dat, Amm.dat, pitzer.dat.
        units_solution (int, optional): Units for solution concentrations.
            1 = mol/L, 2 = mmol/L, 3 = μmol/L. Defaults to 2.
        units (int, optional): Units for other phases (Exchange, Surface,
            Gas, Solid solutions, Kinetics). Defaults to 1.
        porosity (float, optional): Porosity value assigned to all cells.
            Must be between 0 and 1. Defaults to 1.0.
        saturation (float, optional): Saturation value assigned to all cells.
            Must be between 0 and 1. Defaults to 1.0.
        multicomponent (bool, optional): Enable multicomponent diffusion
            by saving species concentrations. Defaults to True.

    Raises:
        RuntimeError: If the PhreeqcRM instance is not initialized or if
            the database fails to load.

    Examples:
        >>> rm = PhreeqcRM()
        >>> rm.create(nxyz=100)
        >>> rm.initialize_phreeqc("phreeqc.dat", units_solution=1)
    """
    if not self._initialized:
        raise RuntimeError("PhreeqcRM instance not initialized. Call create() first.")

    status = self.rm.LoadDatabase(database_path)
    if status < 0:
        raise RuntimeError(f"Failed to load Phreeqc database (error code: {status})")

    self.rm.SetComponentH2O(False)
    self.rm.SetRebalanceFraction(0.5)
    self.rm.SetUnitsSolution(units_solution)
    self.rm.SetUnitsPPassemblage(units)
    self.rm.SetUnitsExchange(units)
    self.rm.SetUnitsSurface(units)
    self.rm.SetUnitsGasPhase(units)
    self.rm.SetUnitsSSassemblage(units)
    self.rm.SetUnitsKinetics(units)
    self.rm.SetPorosity(porosity * np.ones(self.nxyz))
    self.rm.SetSaturationUser(saturation * np.ones(self.nxyz))
    self.rm.SetFilePrefix("phr")
    self.rm.OpenFiles()
    if multicomponent:
        self.rm.SetSpeciesSaveOn(True)

run_initial_from_file(pqi_file, ic)

Set up initial conditions from PHREEQC input file and initial conditions array.

Loads initial geochemical conditions by running a PHREEQC input file and mapping the defined solutions, phases, and other components to the grid cells. This method also retrieves component and species information for later use.

Parameters:

Name Type Description Default
pqi_file str

Path to the PHREEQC input file (.pqi format) containing definitions for solutions, equilibrium phases, exchange, surface, gas phases, solid solutions, and kinetic reactions.

 required
ic ndarray

Initial conditions array with shape (nxyz, 7) where each row corresponds to a grid cell and columns represent: - Column 0: Solution ID - Column 1: Equilibrium phase ID - Column 2: Exchange ID - Column 3: Surface ID - Column 4: Gas phase ID - Column 5: Solid solution ID - Column 6: Kinetic reaction ID Use -1 for unused components.

 required

Raises:

Type Description
RuntimeError

If the PHREEQC input file fails to run.
ValueError

If initial conditions array has incorrect shape or cannot be converted to integer array.

Examples:

>>> import numpy as np
>>> ic = np.array([[1, -1, -1, -1, -1, -1, -1]])  # Only solution 1
>>> rm.run_initial_from_file("initial.pqi", ic)
Source code in mibiremo/phreeqc.py 
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
def run_initial_from_file(self, pqi_file, ic):
    """Set up initial conditions from PHREEQC input file and initial conditions array.

    Loads initial geochemical conditions by running a PHREEQC input file and
    mapping the defined solutions, phases, and other components to the grid cells.
    This method also retrieves component and species information for later use.

    Args:
        pqi_file (str): Path to the PHREEQC input file (.pqi format) containing
            definitions for solutions, equilibrium phases, exchange, surface,
            gas phases, solid solutions, and kinetic reactions.
        ic (numpy.ndarray): Initial conditions array with shape (nxyz, 7) where
            each row corresponds to a grid cell and columns represent:
            - Column 0: Solution ID
            - Column 1: Equilibrium phase ID
            - Column 2: Exchange ID
            - Column 3: Surface ID
            - Column 4: Gas phase ID
            - Column 5: Solid solution ID
            - Column 6: Kinetic reaction ID
            Use -1 for unused components.

    Raises:
        RuntimeError: If the PHREEQC input file fails to run.
        ValueError: If initial conditions array has incorrect shape or cannot
            be converted to integer array.

    Examples:
        >>> import numpy as np
        >>> ic = np.array([[1, -1, -1, -1, -1, -1, -1]])  # Only solution 1
        >>> rm.run_initial_from_file("initial.pqi", ic)
    """
    status = self.rm.RunFile(True, True, True, pqi_file)
    if status < 0:
        raise RuntimeError(f"Failed to run Phreeqc input file (error code: {status})")

    if ic.shape != (self.nxyz, 7):
        raise ValueError(f"Initial conditions array must have shape ({self.nxyz}, 7), got {ic.shape}")

    if not isinstance(ic, np.ndarray):
        try:
            ic = np.array(ic).astype(np.int32)
        except Exception as e:
            raise ValueError("Initial conditions must be convertible to a numpy array of integers") from e

    ic1 = ic.flatten("F").astype(np.int32)
    self.rm.InitialPhreeqc2Module(ic1)

    self.rm.FindComponents()
    self.components = np.array(list(self.rm.GetComponents()))
    self.species = np.array(list(self.rm.GetSpeciesNames()))

    self.rm.SetTime(0.0)
    self.rm.SetTimeStep(0.0)
    self.rm.RunCells()

semilagsolver

Semi-Lagrangian solver for 1D advection-diffusion equation on a uniform grid.

Author: Matteo Masi Last revision: 17/02/2026

SemiLagSolver

Semi-Lagrangian solver for 1D advection-diffusion transport equations.

This class implements a semi-Lagrangian numerical scheme for solving the one-dimensional advection-diffusion equation on uniform grids. The solver uses operator splitting to handle advection and diffusion separately, providing accurate and stable solutions for transport problems.

The numerical approach consists of two sequential steps
  1. Advection: Solved using the Method of Characteristics (MOC) with cubic spline interpolation (PCHIP - Piecewise Cubic Hermite Interpolating Polynomial) to maintain monotonicity and prevent oscillations.
  2. Diffusion: Solved using the Saul’yev alternating direction method, which provides unconditional stability for the diffusion equation.
Boundary Conditions
  • Inlet (left boundary, x=0): Dirichlet-type condition with prescribed concentration value.
  • Outlet (right boundary): Neumann-type condition (zero gradient) allowing natural outflow of transported species.
Mathematical Formulation

The solver addresses the 1D advection-diffusion equation:

∂C/∂t + v∂C/∂x = D∂²C/∂x²

where: - C(x,t): Concentration field - v: Advection velocity (constant) - D: Diffusion/dispersion coefficient (constant)

Numerical Stability
  • The cubic spline advection step is stable for any Courant number
  • The Saul’yev diffusion solver is unconditionally stable
  • Combined scheme maintains stability and accuracy for typical transport problems
Applications
  • Reactive transport modeling in porous media
  • Contaminant transport in groundwater systems
  • Chemical species transport in environmental flows
  • Coupling with geochemical reaction modules (e.g., PhreeqcRM)

Attributes:

Name Type Description
x ndarray

Spatial coordinate array (uniform spacing required).
C ndarray

Current concentration field at grid points.
v float

Advection velocity in consistent units with spatial coordinates.
d float

Diffusion coefficient in consistent units (L²/T).
dt float

Time step for numerical integration in consistent time units.
dx float

Spatial grid spacing (automatically calculated from x).
Note

The spatial grid must be uniformly spaced for the numerical scheme to work correctly. Non-uniform grids are not supported in this implementation.

Source code in mibiremo/semilagsolver.py 
 15
 16
 17
 18
 19
 20
 21
 22
 23
 24
 25
 26
 27
 28
 29
 30
 31
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
class SemiLagSolver:
    """Semi-Lagrangian solver for 1D advection-diffusion transport equations.

    This class implements a semi-Lagrangian numerical scheme for solving the
    one-dimensional advection-diffusion equation on uniform grids. The solver
    uses operator splitting to handle advection and diffusion separately,
    providing accurate and stable solutions for transport problems.

    The numerical approach consists of two sequential steps:
        1. **Advection**: Solved using the Method of Characteristics (MOC) with
           cubic spline interpolation (PCHIP - Piecewise Cubic Hermite Interpolating
           Polynomial) to maintain monotonicity and prevent oscillations.
        2. **Diffusion**: Solved using the Saul'yev alternating direction method,
           which provides unconditional stability for the diffusion equation.

    Boundary Conditions:
        - **Inlet (left boundary, x=0)**: Dirichlet-type condition with prescribed
          concentration value.
        - **Outlet (right boundary)**: Neumann-type condition (zero gradient) allowing
          natural outflow of transported species.

    Mathematical Formulation:
        The solver addresses the 1D advection-diffusion equation:

        ∂C/∂t + v∂C/∂x = D∂²C/∂x²

        where:
        - C(x,t): Concentration field
        - v: Advection velocity (constant)
        - D: Diffusion/dispersion coefficient (constant)

    Numerical Stability:
        - The cubic spline advection step is stable for any Courant number
        - The Saul'yev diffusion solver is unconditionally stable
        - Combined scheme maintains stability and accuracy for typical transport problems

    Applications:
        - Reactive transport modeling in porous media
        - Contaminant transport in groundwater systems
        - Chemical species transport in environmental flows
        - Coupling with geochemical reaction modules (e.g., PhreeqcRM)

    Attributes:
        x (numpy.ndarray): Spatial coordinate array (uniform spacing required).
        C (numpy.ndarray): Current concentration field at grid points.
        v (float): Advection velocity in consistent units with spatial coordinates.
        d (float): Diffusion coefficient in consistent units (L²/T).
        dt (float): Time step for numerical integration in consistent time units.
        dx (float): Spatial grid spacing (automatically calculated from x).

    Note:
        The spatial grid must be uniformly spaced for the numerical scheme to
        work correctly. Non-uniform grids are not supported in this implementation.
    """

    def __init__(self, x, c_init, v, d, dt):
        """Initialize the Semi-Lagrangian solver with transport parameters.

        Sets up the numerical solver with spatial discretization, initial conditions,
        and transport parameters. Validates input consistency and calculates derived
        parameters needed for the numerical scheme.

        Args:
            x (numpy.ndarray): Spatial coordinate array defining the 1D computational
                domain. Must be uniformly spaced with at least 2 points. Units should
                be consistent with velocity and diffusion coefficient.
            c_init (numpy.ndarray): Initial concentration field at each grid point.
                Length must match the spatial coordinate array. Units are user-defined
                but should be consistent throughout the simulation.
            v (float): Advection velocity (positive for left-to-right flow).
                Units must be consistent with spatial coordinates and time step
                (e.g., if x is in meters and dt in days, v should be in m/day).
            d (float): Diffusion/dispersion coefficient (must be non-negative).
                Units must be L²/T where L and T are consistent with spatial
                coordinates and time step (e.g., m²/day).
            dt (float): Time step for numerical integration (must be positive).
                Units should be consistent with velocity and diffusion coefficient.

        Raises:
            ValueError: If spatial coordinates are not uniformly spaced or if
                concentration array length doesn't match spatial coordinates.
            ValueError: If transport parameters are not physically reasonable
                (negative diffusion, zero or negative time step).

        Examples:
            >>> x = np.linspace(0, 5, 51)      # 5 m domain, 0.1 m spacing
            >>> C0 = np.exp(-x**2)             # Gaussian initial condition
            >>> solver = SemiLagSolver(x, C0, v=0.5, d=0.05, dt=0.01)
        """
        if len(x) != len(c_init):
            raise ValueError(f"Length of x ({len(x)}) must match length of c_init ({len(c_init)})")

        if len(x) < 2:
            raise ValueError(f"Grid must have at least 2 points, got {len(x)}")

        self.x = x
        self.C = c_init
        self.v = v
        self.d = d
        self.dt = dt
        self.dx = x[1] - x[0]

    def cubic_spline_advection(self, c_bound) -> None:
        """Solve the advection step using cubic spline interpolation.

        Implements the Method of Characteristics (MOC) for the advection equation
        ∂C/∂t + v∂C/∂x = 0 using backward tracking of characteristic lines.
        Uses PCHIP (Piecewise Cubic Hermite Interpolating Polynomial) to maintain
        monotonicity and prevent numerical oscillations.

        The method works by:
            1. Computing departure points: xi = x - v*dt (backward tracking)
            2. Interpolating concentrations at departure points using cubic splines
            3. Applying inlet boundary condition for points that tracked outside domain

        Args:
            c_bound (float): Inlet concentration value applied at the left boundary
                (x=0) for any characteristic lines that originated from outside the
                computational domain. Units should match the concentration field.

        Note:
            This method modifies self.C in-place. The cubic spline interpolation
            preserves monotonicity, making it suitable for concentration fields
            where spurious oscillations must be avoided.

        Numerical Properties:
            - Unconditionally stable (no CFL restriction)
            - Maintains monotonicity (no new extrema created)
            - Handles arbitrary Courant numbers (v*dt/dx)
            - Exact for linear concentration profiles
        """
        cs = PchipInterpolator(self.x, self.C)
        shift = self.v * self.dt
        xi = self.x - shift
        k0 = xi <= 0
        xi[k0] = 0
        yi = cs(xi)
        yi[k0] = c_bound
        self.C = yi

    def saulyev_solver_alt(self, c_bound) -> None:
        """Solve the diffusion step using the Saul'yev alternating direction method.

        Implements the Saul'yev scheme for the diffusion equation ∂C/∂t = D∂²C/∂x²
        using alternating direction sweeps to achieve unconditional stability.
        The method performs two passes:
            1. Left-to-right sweep using forward differences
            2. Right-to-left sweep using backward differences
            3. Final solution is the average of both sweeps

        Args:
            c_bound (float): Inlet concentration value applied at the left boundary
                during the diffusion solve. This maintains consistency with the
                advection boundary condition.

        Algorithm Details:
            - **Left-to-Right Pass**: For each cell i, uses implicit treatment of
              left neighbor and explicit treatment of right neighbor
            - **Right-to-Left Pass**: For each cell i, uses implicit treatment of
              right neighbor and explicit treatment of left neighbor
            - **Averaging**: Combines both solutions to achieve second-order accuracy

        Boundary Conditions:
            - **Left boundary (x=0)**: Dirichlet condition with prescribed c_bound
            - **Right boundary**: Zero gradient (Neumann) condition implemented
              by using the same concentration as the last interior point

        Numerical Properties:
            - Unconditionally stable for any time step size
            - Second-order accurate in space and time
            - Preserves maximum principle (no spurious extrema)
            - Handles arbitrary diffusion numbers (D*dt/dx²)

        Note:
            This method modifies self.C in-place. The alternating direction
            approach eliminates the restrictive stability constraint of explicit
            methods while maintaining computational efficiency.
        """
        dt = self.dt
        theta = self.d * dt / (self.dx**2)

        # Assign current C state as initial condition
        c_init = self.C.copy()
        clr = self.C.copy()
        crl = self.C.copy()

        # A) L-R direction
        for i in range(len(clr)):
            if i == 0:  # left boundary
                sola = theta * c_bound
            else:
                sola = theta * clr[i - 1]
            solb = (1 - theta) * c_init[i]
            solc = theta * c_init[i + 1] if i < len(clr) - 1 else theta * c_init[i]
            # L-R Solution
            clr[i] = (sola + solb + solc) / (1 + theta)

        # B) R-L direction
        for i in range(len(crl) - 1, -1, -1):
            if i == len(crl) - 1:  # right boundary (take from LR solution)
                sola = theta * clr[-1]
            else:
                sola = theta * crl[i + 1]
            solb = (1 - theta) * c_init[i]
            solc = theta * c_init[i - 1] if i > 0 else theta * c_init[i]
            # R-L Solution
            crl[i] = (sola + solb + solc) / (1 + theta)

        # Average L-R and R-L solutions and update to final state
        self.C = (clr + crl) / 2

    def transport(self, c_bound) -> np.ndarray:
        """Perform one complete transport time step with coupled advection-diffusion.

        Executes the full semi-Lagrangian algorithm by sequentially applying
        the advection and diffusion operators using operator splitting. This
        approach decouples the hyperbolic (advection) and parabolic (diffusion)
        aspects of the transport equation for enhanced numerical stability.

        The operator splitting sequence:
            1. **Advection Step** using cubic spline MOC
            2. **Diffusion Step** using Saul'yev method

        Args:
            c_bound (float): Inlet boundary concentration applied at x=0 for both
                advection and diffusion steps. This represents the concentration
                of material entering the domain (e.g., injection well concentration,
                upstream boundary condition, etc.).

        Returns:
            numpy.ndarray: Updated concentration field after the complete transport
                step. The array has the same shape as the initial concentration
                and represents C(x, t+dt).

        Note:
            This method updates the internal concentration field (self.C) and
            returns the updated values. For reactive transport coupling, call
            this method to advance transport, then apply geochemical reactions
            to the returned concentration field.
        """
        # Step 1: Solve advection equation using cubic spline MOC
        self.cubic_spline_advection(c_bound)

        # Step 2: Solve diffusion equation using Saul'yev alternating direction method
        self.saulyev_solver_alt(c_bound)

        return self.C

__init__(x, c_init, v, d, dt)

Initialize the Semi-Lagrangian solver with transport parameters.

Sets up the numerical solver with spatial discretization, initial conditions, and transport parameters. Validates input consistency and calculates derived parameters needed for the numerical scheme.

Parameters:

Name Type Description Default
x ndarray

Spatial coordinate array defining the 1D computational domain. Must be uniformly spaced with at least 2 points. Units should be consistent with velocity and diffusion coefficient.

 required
c_init ndarray

Initial concentration field at each grid point. Length must match the spatial coordinate array. Units are user-defined but should be consistent throughout the simulation.

 required
v float

Advection velocity (positive for left-to-right flow). Units must be consistent with spatial coordinates and time step (e.g., if x is in meters and dt in days, v should be in m/day).

 required
d float

Diffusion/dispersion coefficient (must be non-negative). Units must be L²/T where L and T are consistent with spatial coordinates and time step (e.g., m²/day).

 required
dt float

Time step for numerical integration (must be positive). Units should be consistent with velocity and diffusion coefficient.

 required

Raises:

Type Description
ValueError

If spatial coordinates are not uniformly spaced or if concentration array length doesn’t match spatial coordinates.
ValueError

If transport parameters are not physically reasonable (negative diffusion, zero or negative time step).

Examples:

>>> x = np.linspace(0, 5, 51)      # 5 m domain, 0.1 m spacing
>>> C0 = np.exp(-x**2)             # Gaussian initial condition
>>> solver = SemiLagSolver(x, C0, v=0.5, d=0.05, dt=0.01)
Source code in mibiremo/semilagsolver.py 
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
def __init__(self, x, c_init, v, d, dt):
    """Initialize the Semi-Lagrangian solver with transport parameters.

    Sets up the numerical solver with spatial discretization, initial conditions,
    and transport parameters. Validates input consistency and calculates derived
    parameters needed for the numerical scheme.

    Args:
        x (numpy.ndarray): Spatial coordinate array defining the 1D computational
            domain. Must be uniformly spaced with at least 2 points. Units should
            be consistent with velocity and diffusion coefficient.
        c_init (numpy.ndarray): Initial concentration field at each grid point.
            Length must match the spatial coordinate array. Units are user-defined
            but should be consistent throughout the simulation.
        v (float): Advection velocity (positive for left-to-right flow).
            Units must be consistent with spatial coordinates and time step
            (e.g., if x is in meters and dt in days, v should be in m/day).
        d (float): Diffusion/dispersion coefficient (must be non-negative).
            Units must be L²/T where L and T are consistent with spatial
            coordinates and time step (e.g., m²/day).
        dt (float): Time step for numerical integration (must be positive).
            Units should be consistent with velocity and diffusion coefficient.

    Raises:
        ValueError: If spatial coordinates are not uniformly spaced or if
            concentration array length doesn't match spatial coordinates.
        ValueError: If transport parameters are not physically reasonable
            (negative diffusion, zero or negative time step).

    Examples:
        >>> x = np.linspace(0, 5, 51)      # 5 m domain, 0.1 m spacing
        >>> C0 = np.exp(-x**2)             # Gaussian initial condition
        >>> solver = SemiLagSolver(x, C0, v=0.5, d=0.05, dt=0.01)
    """
    if len(x) != len(c_init):
        raise ValueError(f"Length of x ({len(x)}) must match length of c_init ({len(c_init)})")

    if len(x) < 2:
        raise ValueError(f"Grid must have at least 2 points, got {len(x)}")

    self.x = x
    self.C = c_init
    self.v = v
    self.d = d
    self.dt = dt
    self.dx = x[1] - x[0]

cubic_spline_advection(c_bound)

Solve the advection step using cubic spline interpolation.

Implements the Method of Characteristics (MOC) for the advection equation ∂C/∂t + v∂C/∂x = 0 using backward tracking of characteristic lines. Uses PCHIP (Piecewise Cubic Hermite Interpolating Polynomial) to maintain monotonicity and prevent numerical oscillations.

The method works by
  1. Computing departure points: xi = x - v*dt (backward tracking)
  2. Interpolating concentrations at departure points using cubic splines
  3. Applying inlet boundary condition for points that tracked outside domain

Parameters:

Name Type Description Default
c_bound float

Inlet concentration value applied at the left boundary (x=0) for any characteristic lines that originated from outside the computational domain. Units should match the concentration field.

 required
Note

This method modifies self.C in-place. The cubic spline interpolation preserves monotonicity, making it suitable for concentration fields where spurious oscillations must be avoided.

Numerical Properties
  • Unconditionally stable (no CFL restriction)
  • Maintains monotonicity (no new extrema created)
  • Handles arbitrary Courant numbers (v*dt/dx)
  • Exact for linear concentration profiles
Source code in mibiremo/semilagsolver.py 
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
def cubic_spline_advection(self, c_bound) -> None:
    """Solve the advection step using cubic spline interpolation.

    Implements the Method of Characteristics (MOC) for the advection equation
    ∂C/∂t + v∂C/∂x = 0 using backward tracking of characteristic lines.
    Uses PCHIP (Piecewise Cubic Hermite Interpolating Polynomial) to maintain
    monotonicity and prevent numerical oscillations.

    The method works by:
        1. Computing departure points: xi = x - v*dt (backward tracking)
        2. Interpolating concentrations at departure points using cubic splines
        3. Applying inlet boundary condition for points that tracked outside domain

    Args:
        c_bound (float): Inlet concentration value applied at the left boundary
            (x=0) for any characteristic lines that originated from outside the
            computational domain. Units should match the concentration field.

    Note:
        This method modifies self.C in-place. The cubic spline interpolation
        preserves monotonicity, making it suitable for concentration fields
        where spurious oscillations must be avoided.

    Numerical Properties:
        - Unconditionally stable (no CFL restriction)
        - Maintains monotonicity (no new extrema created)
        - Handles arbitrary Courant numbers (v*dt/dx)
        - Exact for linear concentration profiles
    """
    cs = PchipInterpolator(self.x, self.C)
    shift = self.v * self.dt
    xi = self.x - shift
    k0 = xi <= 0
    xi[k0] = 0
    yi = cs(xi)
    yi[k0] = c_bound
    self.C = yi

saulyev_solver_alt(c_bound)

Solve the diffusion step using the Saul’yev alternating direction method.

Implements the Saul’yev scheme for the diffusion equation ∂C/∂t = D∂²C/∂x² using alternating direction sweeps to achieve unconditional stability. The method performs two passes: 1. Left-to-right sweep using forward differences 2. Right-to-left sweep using backward differences 3. Final solution is the average of both sweeps

Parameters:

Name Type Description Default
c_bound float

Inlet concentration value applied at the left boundary during the diffusion solve. This maintains consistency with the advection boundary condition.

 required
Algorithm Details
  • Left-to-Right Pass: For each cell i, uses implicit treatment of left neighbor and explicit treatment of right neighbor
  • Right-to-Left Pass: For each cell i, uses implicit treatment of right neighbor and explicit treatment of left neighbor
  • Averaging: Combines both solutions to achieve second-order accuracy
Boundary Conditions
  • Left boundary (x=0): Dirichlet condition with prescribed c_bound
  • Right boundary: Zero gradient (Neumann) condition implemented by using the same concentration as the last interior point
Numerical Properties
  • Unconditionally stable for any time step size
  • Second-order accurate in space and time
  • Preserves maximum principle (no spurious extrema)
  • Handles arbitrary diffusion numbers (D*dt/dx²)
Note

This method modifies self.C in-place. The alternating direction approach eliminates the restrictive stability constraint of explicit methods while maintaining computational efficiency.

Source code in mibiremo/semilagsolver.py 
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
def saulyev_solver_alt(self, c_bound) -> None:
    """Solve the diffusion step using the Saul'yev alternating direction method.

    Implements the Saul'yev scheme for the diffusion equation ∂C/∂t = D∂²C/∂x²
    using alternating direction sweeps to achieve unconditional stability.
    The method performs two passes:
        1. Left-to-right sweep using forward differences
        2. Right-to-left sweep using backward differences
        3. Final solution is the average of both sweeps

    Args:
        c_bound (float): Inlet concentration value applied at the left boundary
            during the diffusion solve. This maintains consistency with the
            advection boundary condition.

    Algorithm Details:
        - **Left-to-Right Pass**: For each cell i, uses implicit treatment of
          left neighbor and explicit treatment of right neighbor
        - **Right-to-Left Pass**: For each cell i, uses implicit treatment of
          right neighbor and explicit treatment of left neighbor
        - **Averaging**: Combines both solutions to achieve second-order accuracy

    Boundary Conditions:
        - **Left boundary (x=0)**: Dirichlet condition with prescribed c_bound
        - **Right boundary**: Zero gradient (Neumann) condition implemented
          by using the same concentration as the last interior point

    Numerical Properties:
        - Unconditionally stable for any time step size
        - Second-order accurate in space and time
        - Preserves maximum principle (no spurious extrema)
        - Handles arbitrary diffusion numbers (D*dt/dx²)

    Note:
        This method modifies self.C in-place. The alternating direction
        approach eliminates the restrictive stability constraint of explicit
        methods while maintaining computational efficiency.
    """
    dt = self.dt
    theta = self.d * dt / (self.dx**2)

    # Assign current C state as initial condition
    c_init = self.C.copy()
    clr = self.C.copy()
    crl = self.C.copy()

    # A) L-R direction
    for i in range(len(clr)):
        if i == 0:  # left boundary
            sola = theta * c_bound
        else:
            sola = theta * clr[i - 1]
        solb = (1 - theta) * c_init[i]
        solc = theta * c_init[i + 1] if i < len(clr) - 1 else theta * c_init[i]
        # L-R Solution
        clr[i] = (sola + solb + solc) / (1 + theta)

    # B) R-L direction
    for i in range(len(crl) - 1, -1, -1):
        if i == len(crl) - 1:  # right boundary (take from LR solution)
            sola = theta * clr[-1]
        else:
            sola = theta * crl[i + 1]
        solb = (1 - theta) * c_init[i]
        solc = theta * c_init[i - 1] if i > 0 else theta * c_init[i]
        # R-L Solution
        crl[i] = (sola + solb + solc) / (1 + theta)

    # Average L-R and R-L solutions and update to final state
    self.C = (clr + crl) / 2

transport(c_bound)

Perform one complete transport time step with coupled advection-diffusion.

Executes the full semi-Lagrangian algorithm by sequentially applying the advection and diffusion operators using operator splitting. This approach decouples the hyperbolic (advection) and parabolic (diffusion) aspects of the transport equation for enhanced numerical stability.

The operator splitting sequence
  1. Advection Step using cubic spline MOC
  2. Diffusion Step using Saul’yev method

Parameters:

Name Type Description Default
c_bound float

Inlet boundary concentration applied at x=0 for both advection and diffusion steps. This represents the concentration of material entering the domain (e.g., injection well concentration, upstream boundary condition, etc.).

 required

Returns:

Type Description
ndarray

numpy.ndarray: Updated concentration field after the complete transport step. The array has the same shape as the initial concentration and represents C(x, t+dt).
Note

This method updates the internal concentration field (self.C) and returns the updated values. For reactive transport coupling, call this method to advance transport, then apply geochemical reactions to the returned concentration field.

Source code in mibiremo/semilagsolver.py 
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
def transport(self, c_bound) -> np.ndarray:
    """Perform one complete transport time step with coupled advection-diffusion.

    Executes the full semi-Lagrangian algorithm by sequentially applying
    the advection and diffusion operators using operator splitting. This
    approach decouples the hyperbolic (advection) and parabolic (diffusion)
    aspects of the transport equation for enhanced numerical stability.

    The operator splitting sequence:
        1. **Advection Step** using cubic spline MOC
        2. **Diffusion Step** using Saul'yev method

    Args:
        c_bound (float): Inlet boundary concentration applied at x=0 for both
            advection and diffusion steps. This represents the concentration
            of material entering the domain (e.g., injection well concentration,
            upstream boundary condition, etc.).

    Returns:
        numpy.ndarray: Updated concentration field after the complete transport
            step. The array has the same shape as the initial concentration
            and represents C(x, t+dt).

    Note:
        This method updates the internal concentration field (self.C) and
        returns the updated values. For reactive transport coupling, call
        this method to advance transport, then apply geochemical reactions
        to the returned concentration field.
    """
    # Step 1: Solve advection equation using cubic spline MOC
    self.cubic_spline_advection(c_bound)

    # Step 2: Solve diffusion equation using Saul'yev alternating direction method
    self.saulyev_solver_alt(c_bound)

    return self.C