### Define system parameters and InterfaceSystem setup

Source: https://pyequion2.readthedocs.io/en/latest/tutorials.html

Sets up the InterfaceSystem with elements, defines index maps, and specifies pipe parameters and transport parameters for the pipe flow.

```python
elements = ['Ca', 'C', 'Na', 'Cl', 'Mg']
intsys = InterfaceSystem(elements, from_elements=True)
intsys.set_interface_phases()
index_map = {el: i for i, el in enumerate(elements)} #For creating the solution vector
reverse_index_map = {i: el for i, el in enumerate(elements)} #For creating the solution vector

TK = 298.15
pipe_diameter = 0.05 #m
flow_velocity = 1.0
pipe_length = 80.0 #m
pipe_time = pipe_length/flow_velocity

transport_params = {'type': 'pipe',
                    'shear_velocity': shear_velocity(flow_velocity, pipe_diameter, TK)}
```

--------------------------------

### Install PyEquion2 with pip

Source: https://pyequion2.readthedocs.io/en/latest/install.html

Use this command to install the PyEquion2 package from PyPI using pip.

```sh
pip install pyequion2
```

--------------------------------

### Retrieve total carbon concentration

Source: https://pyequion2.readthedocs.io/en/latest/tutorials.html

Gets the total molal concentration of carbon in the solution.

```python
>>> solution.elements_molals['C']
0.12677211939028027
```

--------------------------------

### Set up EquilibriumSystem with elements Na and C

Source: https://pyequion2.readthedocs.io/en/latest/tutorials.html

Initializes the equilibrium system with elements 'Na' and 'C' using from_elements=True.

```python
>>> from pyequion2 import EquilibriumSystem
>>> EquilibriumSystem(['Na', 'C'], from_elements=True)
```

--------------------------------

### Initialize EquilibriumSystem with extended Debye activity model

Source: https://pyequion2.readthedocs.io/en/latest/tutorials.html

Sets up the system with the EXTENDED_DEBYE activity model at initialization.

```python
>>> eqsys = EquilibriumSystem(['CaCl2', 'NaHCO3'], activity_model="EXTENDED_DEBYE")
```

--------------------------------

### Set up CO2 degassing system

Source: https://pyequion2.readthedocs.io/en/latest/tutorials.html

Imports EquilibriumSystem, sets up a CO2 feed, defines a molal balance, and solves the equilibrium at 298.15 K. The solution shows significant degassing.

```python
>>> from pyequion2 import EquilibriumSystem #Import the necessary module
>>> eqsys = EquilibriumSystem(['CO2']) #We set up the feed components of our system
>>> molal_balance = {'C':0.5} #Set up the balances
>>> TK = 298.15
>>> solution, _ = eqsys.solve_equilibrium_elements_balance_phases(TK, molal_balance)
```

--------------------------------

### Initialize EquilibriumSystem from elements

Source: https://pyequion2.readthedocs.io/en/latest/tutorials.html

Sets up the equilibrium system using element names instead of feed components.

```python
>>> eqsys = EquilibriumSystem(['Ca', 'Cl', 'Na', 'C'], from_elements=True)
```

--------------------------------

### Set up equilibrium system with precipitation

Source: https://pyequion2.readthedocs.io/en/latest/tutorials.html

Imports the EquilibriumSystem module and defines the system with elements Ca, Na, C, Cl, along with molal balance and temperature/pressure conditions.

```python
>>> from pyequion2 import EquilibriumSystem #Import the necessary module
>>> eqsys = EquilibriumSystem(['Ca', 'Na', 'C', 'Cl'], from_elements=True)
>>> molal_balance = {'Ca':0.028, 'C':0.065, 'Na':0.065, 'Cl':0.056}
>>> TK = 298.15 #Temperature in Kelvin
>>> PATM = 1.0 #Pressure in atm
```

--------------------------------

### Set up and solve ODE system

Source: https://pyequion2.readthedocs.io/en/latest/tutorials.html

Defines the initial element balance vector and uses scipy.integrate.solve_ivp to solve the differential equation over the pipe time.

```python
initial_elements_balance = {'Ca':0.028, 'C':0.065, 'Na':0.075, 'Cl':0.056, 'Mg':0.02}
initial_elements_vector = np.hstack([initial_elements_balance[reverse_index_map[i]]
                                     for i in range(len(initial_elements_balance))])

start_time = time.time()
sol = scipy.integrate.solve_ivp(f, (0.0, pipe_time), initial_elements_vector)
elapsed_time = time.time() - start_time
```

--------------------------------

### Set up interface phases and transport parameters

Source: https://pyequion2.readthedocs.io/en/latest/tutorials.html

Sets up solid phases at the interface, retrieves bulk solute molals, and defines transport parameters for a pipe flow with shear velocity 0.05 m/s.

```python
>>> intsys.set_interface_phases()
>>> molals_bulk = solution.solute_molals
>>> transport_params = {'type': 'pipe',
                        'shear_velocity': 0.05}
```

--------------------------------

### Initialize EquilibriumSystem with feed components

Source: https://pyequion2.readthedocs.io/en/latest/tutorials.html

Creates an equilibrium system from feed components, automatically including all corresponding aqueous species.

```python
>>> eqsys = EquilibriumSystem(['CaCl2', 'NaHCO3'])
```

--------------------------------

### Generate and read solution log file

Source: https://pyequion2.readthedocs.io/en/latest/tutorials.html

Save a log file of the solution using savelog and then read it back. The output shows the full log content including conditions, component concentrations, properties, and phase information.

```python
>>> solution.savelog("solutionlog.txt")
>>> with open("solutionlog.txt", "r") as f:
...     print(f.read())
CONDITIONS
aqueous
T - 298.15 K
P - 1.00 atm
[Ca]=2.80e-02 mol/kg H2O
[C]=6.50e-02 mol/kg H2O
[Na]=6.50e-02 mol/kg H2O
[Cl]=5.60e-02 mol/kg H2O
----------------------------------------
[COMPONENT]    [CONCENTRATION (mol/kg H2O)]    [ACTIVITY (mol/kg H2O)]    [MOLE FRACTION]
H2O    5.551e+01    1.000e+00    9.963014e-01
Na+    6.374e-02    4.881e-02    1.144102e-03
Cl-    5.597e-02    4.157e-02    1.004663e-03
HCO3-    5.505e-02    4.243e-02    9.880605e-04
Ca++    2.134e-02    7.677e-03    3.829556e-04
CaHCO3+    5.332e-03    4.154e-03    9.569780e-05
CO2    1.745e-03    1.800e-03    3.132396e-05
CaCO3    1.319e-03    1.360e-03    2.367630e-05
NaHCO3    1.129e-03    1.165e-03    2.027042e-05
CO3--    2.989e-04    1.055e-04    5.364689e-06
NaCO3-    1.252e-04    9.586e-05    2.247923e-06
CaCl2    1.286e-05    1.326e-05    2.308680e-07
Na2CO3    1.145e-06    1.181e-06    2.055222e-08
OH-    7.258e-07    5.365e-07    1.302778e-08
CaOH+    8.823e-08    6.754e-08    1.583637e-09
H+    2.319e-08    1.887e-08    4.161924e-10
NaOH    1.658e-08    1.709e-08    2.975389e-10
C    6.500e-02
Ca    2.800e-02
Cl    5.600e-02
H    1.111e+02
Na    6.500e-02
O    5.570e+01
----------------------------------------
PROPERTIES
pH = 7.724
I = 0.133 mol/kg H2O
conductivity = 0.827 S/m
----------------------------------------
[PHASE]    [AMOUNT mol/kg H2O]
Aragonite    0.0
Calcite    0.0
Halite    0.0
Vaterite    0.0
CO2(g)    0.0
H2O(g)    0.0
----------------------------------------
[PHASE]    [SUPERSATURATION]    [SI]
Aragonite    175.54843195012336   2.2443969545735136
Calcite    244.43691952856932   2.388166801965749
Halite    5.4607121083262775e-05   -4.262750719099147
Vaterite    66.33386579002101   1.8217353077514495
```

--------------------------------

### Check interface phases considered

Source: https://pyequion2.readthedocs.io/en/latest/tutorials.html

Lists the solid phases that will be considered at the interface, showing Calcite, Dolomite, and Halite.

```python
>>> intsys.interface_phases
['Calcite', 'Dolomite', 'Halite']
```

--------------------------------

### Import EquilibriumSystem from pyequion2

Source: https://pyequion2.readthedocs.io/en/latest/tutorials.html

Imports the main API interface for setting up and calculating electrolytical balance systems.

```python
>>> from pyequion2 import EquilibriumSystem
```

--------------------------------

### Initialize solution stats for initial guess holding

Source: https://pyequion2.readthedocs.io/en/latest/tutorials.html

Initializes dictionaries to hold the solution of each step as the initial guess for the next step, improving solver efficiency.

```python
>>> solution_stats = {'res': None, 'x': 'default'}
>>> solution_stats_int = {'res': None, 'x': 'default'}
```

--------------------------------

### Import InterfaceSystem for transport equilibrium

Source: https://pyequion2.readthedocs.io/en/latest/tutorials.html

Imports InterfaceSystem and initializes it with element components. InterfaceSystem inherits from EquilibriumSystem.

```python
>>> from pyequion2 import InterfaceSystem
>>> intsys = InterfaceSystem(['Ca', 'C', 'Na', 'Cl', 'Mg'], from_elements=True)
```

--------------------------------

### Access solution properties

Source: https://pyequion2.readthedocs.io/en/latest/tutorials.html

Access pH, ionic strength (in mol/kg H2O), and electrical conductivity (in S/m) of the solution. The output shows the values for each property.

```python
>>> solution.ph
7.724341513364517
>>> solution.ionic_strength #mol/kg H2O
0.13338239393747503
>>> solution.electrical_conductivity #S/m
0.8269708424123452
```

--------------------------------

### Inspect species in the equilibrium system

Source: https://pyequion2.readthedocs.io/en/latest/tutorials.html

Lists all aqueous species included in the system after initialization.

```python
>>> eqsys.species
['H2O', 'CO2', 'CO3--', 'Ca++', 'CaCO3', 'CaCl2', 'CaHCO3+',
 'CaOH+', 'Cl-', 'H+', 'HCO3-', 'Na+', 'Na2CO3', 'NaCO3-',
 'NaHCO3', 'NaOH', 'OH-']
```

--------------------------------

### Calculate electrolytical balance with solve_equilibrium_mixed_balance

Source: https://pyequion2.readthedocs.io/en/latest/tutorials.html

Sets up an EquilibriumSystem with feed components CaCl2 and NaHCO3, defines a molal balance, and solves for the equilibrium at 298.15 K and 1 atm. Returns a solution object and solution statistics.

```python
from pyequion2 import EquilibriumSystem #Import the necessary module
eqsys = EquilibriumSystem(['CaCl2', 'NaHCO3']) #We set up the feed components of our system
molal_balance = {'Ca':0.028, 'C':0.065, 'Na':0.065, 'Cl':0.056} #Set up the balances
TK = 298.15 #Temperature in Kelvin
PATM = 1.0 #Pressure in atm
#Returns the solution class and tuple of solution statistics
solution, solution_stats = eqsys.solve_equilibrium_mixed_balance(TK, molal_balance=molal_balance, PATM=PATM)
```

--------------------------------

### Solve bulk equilibrium for interface system

Source: https://pyequion2.readthedocs.io/en/latest/tutorials.html

Defines element balances and solves the bulk equilibrium using solve_equilibrium_mixed_balance.

```python
>>> elements_balance = {'Ca':0.028, 'C':0.065, 'Na':0.075, 'Cl':0.056, 'Mg':0.02}
>>> TK = 298.15
>>> solution, solution_stats = intsys.solve_equilibrium_mixed_balance(TK, molal_balance=elements_balance)
```

--------------------------------

### Change activity model after initialization

Source: https://pyequion2.readthedocs.io/en/latest/tutorials.html

Switches the activity model to PITZER after the system has been created.

```python
>>> eqsys.set_activity_functions(activity_model="PITZER")
```

--------------------------------

### Specify temperature and pressure

Source: https://pyequion2.readthedocs.io/en/latest/tutorials.html

Sets the temperature in Kelvin and pressure in atm for the calculation. Pressure is optional and defaults to 1 atm.

```python
>>> TK = 298.15 #Temperature in Kelvin
>>> PATM = 1.0 #Pressure in atm
```

--------------------------------

### Inspect reaction fluxes at interface

Source: https://pyequion2.readthedocs.io/en/latest/tutorials.html

Shows the flux of solid formation at the surface for each mineral phase.

```python
>>> solution_int.reaction_fluxes
{'Aragonite': 0.0,
 'Calcite': 0.0011847758387262064,
 'Dolomite': 0.002514465853801235,
 'Halite': 1e-200,
 'Vaterite': 0.0}
```

--------------------------------

### Import pyequion2 modules

Source: https://pyequion2.readthedocs.io/en/latest/tutorials.html

Imports the InterfaceSystem class and the water_properties module from pyequion2.

```python
>>> from pyequion2 import InterfaceSystem
>>> from pyequion2 import water_properties #Collection of a bunch of water properties
```

--------------------------------

### Solve equilibrium with mixed balance

Source: https://pyequion2.readthedocs.io/en/latest/tutorials.html

Calls solve_equilibrium_mixed_balance with the temperature, molal_balance, and activities_balance_log to compute the solution.

```python
>>> solution, solution_stats = \
...           eqsys.solve_equilibrium_mixed_balance(TK,
...                                                 molal_balance=molal_balance,
...                                                 activities_balance_log=activities_balance_log)
```

--------------------------------

### Solve equilibrium with phases

Source: https://pyequion2.readthedocs.io/en/latest/tutorials.html

Calls solve_equilibrium_elements_balance_phases to compute equilibrium assuming the most stable phase of each possible solid precipitates if possible.

```python
>>> solution, solution_stats = eqsys.solve_equilibrium_elements_balance_phases(TK, molal_balance)
```

--------------------------------

### Check pH of solution

Source: https://pyequion2.readthedocs.io/en/latest/tutorials.html

Prints the pH of the solution to confirm it matches the fixed value.

```python
>>> solution.ph
9.000000000000147
```

--------------------------------

### Print species molals and activities

Source: https://pyequion2.readthedocs.io/en/latest/tutorials.html

Iterates over the solution's species and prints their molal concentrations and activities.

```python
>>> for specie in solution.species:
...     print(solution.molals[specie], solution.activities[specie])
H2O 55.508435061791985 1.0
CO2 0.0017451987860942807 0.0017996296574881963
CO3-- 0.00029889095285012476 0.00010547453992485908
Ca++ 0.021336180473077816 0.007676960165455775
CaCO3 0.0013191134282240777 0.0013602551559961564
CaCl2 1.286269802134016e-05 1.3263871725652214e-05
CaHCO3+ 0.005331755169106101 0.004154442718456503
CaOH+ 8.823157066306216e-08 6.753529488353837e-08
Cl- 0.05597427460395732 0.041566218194101646
H+ 2.318795121973022e-08 1.8865072850393035e-08
HCO3- 0.05504929827427368 0.04242877615220498
Na+ 0.06374309497676348 0.04880999552873938
Na2CO3 1.1450565517797437e-06 1.1807696329515286e-06
NaCO3- 0.00012524189134245758 9.586419010815542e-05
NaHCO3 0.0011293564415574976 0.0011645798531929594
NaOH 1.6577233005177517e-08 1.7094259056857954e-08
OH- 7.258362733865141e-07 5.365250097984246e-07
```

--------------------------------

### Import numpy and scipy.integrate

Source: https://pyequion2.readthedocs.io/en/latest/tutorials.html

Imports the numpy library for numerical operations and scipy.integrate for ODE solving.

```python
>>> import numpy as np
>>> import scipy.integrate #For ODE solving
```

--------------------------------

### Solve interface equilibrium

Source: https://pyequion2.readthedocs.io/en/latest/tutorials.html

Solves the interface equilibrium using the bulk molals and transport parameters.

```python
>>> solution_int, res = intsys.solve_interface_equilibrium(TK,
...                                                        molals_bulk,
...                                                        transport_params)
```

--------------------------------

### Inspect element reaction fluxes

Source: https://pyequion2.readthedocs.io/en/latest/tutorials.html

Shows the flux of each element due to reactions at the interface.

```python
>>> solution_int.elements_reaction_fluxes
{'H': 0.0,
 'O': 0.01864112263898603,
 'Cl': 1e-200,
 'C': 0.0062137075463286765,
 'Ca': 0.0036992416925274415,
 'Na': 1e-200,
 'Mg': 0.002514465853801235}
```

--------------------------------

### Solve equilibrium with solve_equilibrium_mixed_balance

Source: https://pyequion2.readthedocs.io/en/latest/tutorials.html

Calls the method to calculate the electrolytical balance, returning the solution and statistics.

```python
>>> solution, solution_stats = eqsys.solve_equilibrium_mixed_balance(TK, molal_balance=molal_balance, PATM=PATM)
```

--------------------------------

### Dissolve CO2 gas by increasing pressure

Source: https://pyequion2.readthedocs.io/en/latest/tutorials.html

Solves the equilibrium at 100 atm and shows that the CO2 gas molality becomes negligible (1e-200), indicating full dissolution.

```python
>>> PATM = 100.0
>>> solution_high_pressure, _ = eqsys.solve_equilibrium_elements_balance_phases(TK,
...                                                                             molal_balance,
...                                                                             PATM=PATM)
>>> solution_high_pressure.gas_molals['CO2(g)']
1e-200
```

--------------------------------

### Convert temperature to Kelvin

Source: https://pyequion2.readthedocs.io/en/latest/tutorials.html

Converts 35°C to Kelvin by adding 273.15.

```python
>>> TK = 35 + 273.15 #35 ºC to Kelvin
```

--------------------------------

### Inspect saturation indexes

Source: https://pyequion2.readthedocs.io/en/latest/tutorials.html

Access the saturation_indexes attribute to see the saturation index for each possible solid precipitate. The output is a dictionary with mineral names and their SI values.

```python
>>> solution.saturation_indexes
{'Aragonite': 2.2443969545735136,
 'Calcite': 2.388166801965749,
 'Halite': -4.262750719099147,
 'Vaterite': 1.8217353077514495}
```

--------------------------------

### Define non-stoichiometric molal balance

Source: https://pyequion2.readthedocs.io/en/latest/tutorials.html

Shows a molal balance where the carbon amount differs from sodium, demonstrating flexibility.

```python
>>> molal_balance = {'Ca':0.028, 'C':0.075, 'Na':0.065, 'Cl':0.056} #molal_balance['C'] != molal_balance['Na']
```

--------------------------------

### Verify conservation including precipitation

Source: https://pyequion2.readthedocs.io/en/latest/tutorials.html

Adds the precipitated calcite to the dissolved elements to show conservation of Ca and C.

```python
>>> solution.solid_molals['Calcite'] + solution.elements_molals['Ca']
0.028
>>> solution.solid_molals['Calcite'] + solution.elements_molals['C']
0.06499999999999997
```

--------------------------------

### Inspect residual of solution

Source: https://pyequion2.readthedocs.io/en/latest/tutorials.html

Computes the maximum absolute residual from solution_stats to verify convergence.

```python
>>> import numpy as np
>>> np.max(np.abs(solution_stats['res']))
8.505418591653324e-13
```

--------------------------------

### Check gas molals

Source: https://pyequion2.readthedocs.io/en/latest/tutorials.html

Shows the gas molals, indicating no degassing in the system.

```python
>>> solution.gas_molals
{'CO2(g)': 1.0714181266404203e-26, 'H2O(g)': 0.0}
```

--------------------------------

### Check dissolved elements after precipitation

Source: https://pyequion2.readthedocs.io/en/latest/tutorials.html

Displays the dissolved molals of Ca and C after precipitation, which are not conserved alone.

```python
>>> solution.elements_molals['Ca'], solution.elements_molals['C']
(0.0370571521476301, 5.7152147630119065e-05)
```

--------------------------------

### Inspect transport fluxes at interface

Source: https://pyequion2.readthedocs.io/en/latest/tutorials.html

Shows the flux of each species due to transport, with some reactive species having negative fluxes.

```python
>>> solution_int.transport_fluxes
{'CO2': -2.195876510859801e-06,
 'CO3--': 0.0017140194883697559,
 'Ca++': 0.0010649319732573606,
 'CaCO3': 0.002586944230083052,
 'CaHCO3+': 4.6679707750352626e-05,
 'CaOH+': 6.857814366765381e-07,
 'Cl-': 0.0,
 'H+': -3.963201178383444e-10,
 'HCO3-': -4.7309531127778806e-05,
 'Mg++': 0.0009919544743964446,
 'MgCO3': 0.0014404475305928642,
 'MgHCO3+': 6.879935653506196e-05,
 'MgOH+': 1.3264492276864372e-05,
 'Na+': -0.0004119691355775593,
 'Na2CO3': 5.351388380819752e-06,
 'NaCO3-': 0.00041290482611060334,
 'NaHCO3': -1.1933573855193434e-05,
 'NaOH': 2.951065605094775e-07,
 'OH-': 3.759842968655567e-05}
```

--------------------------------

### Inspect element molalities

Source: https://pyequion2.readthedocs.io/en/latest/tutorials.html

Access the elements_molals attribute to see the molality of each element in the solution. The output is a dictionary with element symbols as keys and molalities as values.

```python
>>> solution.elements_molals
{'H': 111.07838138730193,
 'O': 55.70169069365096,
 'Ca': 0.028,
 'Na': 0.06499999999999999,
 'Cl': 0.056,
 'C': 0.06499999999999999}
```

--------------------------------

### Inspect CO2 gas molality after degassing

Source: https://pyequion2.readthedocs.io/en/latest/tutorials.html

Displays the gas molality of CO2(g) from the solution, showing a high value indicating degassing.

```python
>>> solution.gas_molals['CO2(g)']
0.49972244867277366
```

--------------------------------

### Check precipitating phases

Source: https://pyequion2.readthedocs.io/en/latest/tutorials.html

Displays the solid molals for each possible precipitate. The value for calcite is the mols of calcite per unit of liquid H2O.

```python
>>> solution.solid_molals
{'Aragonite': 0.0,
 'Calcite': 0.023947112655582508,
 'Halite': 3.507133834338337e-116,
 'Vaterite': 0.0}
```

--------------------------------

### Define molal balance for elements

Source: https://pyequion2.readthedocs.io/en/latest/tutorials.html

Specifies the amount of each element in molals, allowing non-stoichiometric balances.

```python
>>> molal_balance = {'Ca':0.028, 'C':0.065, 'Na':0.065, 'Cl':0.056}
```

--------------------------------

### Fix pH by setting H+ log-activity

Source: https://pyequion2.readthedocs.io/en/latest/tutorials.html

Fixes the pH by setting the log-activity of H+ to -9.0, as pH is the negative base-10 logarithm of H+ activity.

```python
>>> activities_balance_log = {'H+': -9.0} #Fix our pH by fixing H+ log-activity
```

--------------------------------

### Define molal balance for Na

Source: https://pyequion2.readthedocs.io/en/latest/tutorials.html

Sets the molal balance for Na to 0.150 molal, as Na is non-volatile and its amount remains as added.

```python
>>> molal_balance = {'Na': 0.150}
```

--------------------------------

### Inspect solution statistics residual

Source: https://pyequion2.readthedocs.io/en/latest/tutorials.html

Checks the residual of the solution as an error metric, using numpy to find the maximum absolute residual.

```python
>>> import numpy as np
>>> residual = solution_stats['res']
>>> np.max(np.abs(solution_stats['res']))
7.513989430663059e-13
```

--------------------------------

### Define functions for shear velocity calculation

Source: https://pyequion2.readthedocs.io/en/latest/tutorials.html

Defines reynolds_number, darcy_friction_factor, and shear_velocity functions to compute flow characteristics. The darcy_friction_factor uses laminar flow formula for Reynolds < 2300 and Blasius correlation otherwise.

```python
def reynolds_number(flow_velocity, pipe_diameter, TK=298.15): #Dimensionless
    kinematic_viscosity = water_properties.water_kinematic_viscosity(TK)
    return flow_velocity*pipe_diameter/kinematic_viscosity


def darcy_friction_factor(flow_velocity, pipe_diameter, TK=298.15):
    reynolds = reynolds_number(flow_velocity, pipe_diameter, TK)
    if reynolds < 2300:
        return 64/reynolds
    else: #Blasius
        return 0.316*reynolds**(-1./4)

def shear_velocity(flow_velocity, pipe_diameter, TK=298.15):
    f = darcy_friction_factor(flow_velocity, pipe_diameter, TK)
    return np.sqrt(f/8.0)*flow_velocity
```

--------------------------------

### Define right-hand side of differential equation

Source: https://pyequion2.readthedocs.io/en/latest/tutorials.html

Defines the function f(t, y) that computes the derivative of element concentrations based on equilibrium and interface equilibrium solutions. It uses global solution_stats and solution_stats_int to hold initial guesses.

```python
def f(t, y):
    global solution_stats
    global solution_stats_int
    elements_balance = {el: y[index_map[el]] for el in elements}
    solution, solution_stats = intsys.solve_equilibrium_mixed_balance(TK,
                                                                      molal_balance=elements_balance,
                                                                      tol=1e-6,
                                                                      initial_guess=solution_stats['x'])
    molals_bulk = solution.solute_molals
    solution_int, solution_stats_int = intsys.solve_interface_equilibrium(TK,
                                                                          molals_bulk,
                                                                          transport_params,
                                                                          tol=1e-6,
                                                                          initial_guess=solution_stats_int['x'])
    elements_reaction_fluxes = solution_int.elements_reaction_fluxes
    wall_scale = 4/(pipe_diameter*water_properties.water_density(TK))
    dy = -wall_scale*np.hstack(
        [elements_reaction_fluxes[reverse_index_map[i]]
         for i in range(y.shape[0])])
    return dy
```

--------------------------------

### Define molal balance with specific ion concentration

Source: https://pyequion2.readthedocs.io/en/latest/tutorials.html

Specifies the amount of a specific ion (Ca++) in the molal balance, fixing its concentration.

```python
>>> molal_balance = {'Ca++':0.028, 'C':0.075, 'Na':0.065, 'Cl':0.056} #Our solution will necessary have 0.028 molals of Ca++ ions.
```

=== COMPLETE CONTENT === This response contains all available snippets from this library. No additional content exists. Do not make further requests.