Thermolab: A Thermodynamics Laboratory for Nonlinear Transport Processes in Open Systems

We developed a numerical thermodynamics laboratory called “Thermolab” to study the effects of the thermodynamic behavior of nonideal solution models on reactive transport processes in open systems. The equations of the state of internally consistent thermodynamic data sets are implemented in MATLAB functions and form the basis for calculating Gibbs energy. A linear algebraic approach is used in Thermolab to compute Gibbs energy of mixing for multicomponent phases to study the impact of the nonideality of solution models on transport processes. The Gibbs energies are benchmarked with experimental data, phase diagrams, and other thermodynamic software. Constrained Gibbs minimization is exemplified with MATLAB codes and iterative refinement of composition of mixtures may be used to increase precision and accuracy. All needed transport variables such as densities, phase compositions, and chemical potentials are obtained from Gibbs energy of the stable phases after the minimization in Thermolab. We demonstrate the use of precomputed local equilibrium data obtained with Thermolab in reactive transport models. In reactive fluid flow the shape and the velocity of the reaction front vary depending on the nonlinearity of the partitioning of a component in fluid and solid. We argue that nonideality of solution models has to be taken into account and further explored in reactive transport models. Thermolab Gibbs energies can be used in Cahn‐Hilliard models for nonlinear diffusion and phase growth. This presents a transient process toward equilibrium and avoids computational problems arising during precomputing of equilibrium data.

3 of 44 documenting the algorithms. Matlab also provides a prototype coding development platform. Short MATLAB codes serve as "flow charts" used in the past to document the algorithms. The advantage of the MATLAB/ OCTAVE codes compared to the "flow charts" is that those are actually working computer programs, and can be copy-pasted from the figure and executed. Once working these codes can also be easily translated to any preferred programming language for optimization and supercomputing. The focus in this contribution is entirely on the technicalities with limited example applications. We document the nonideality of the models and the impact on nonlinear transport processes in open systems by showing examples of how this can be achieved with Thermolab. Future studies will be investigating the impact of the nonideality on mass transport.

Background and Motivation
A typical set of equations that governs the physics of the reactive transport based on Beinlich et al. (2020) is shown to motivate the use of equilibrium thermodynamics in transport processes. We consider the one-dimensional transport of a mobile species described by the total mass concentration balance, assuming nonmoving, nondeforming solid (v s (x,t) = 0): with total mass concentration in the time derivative above defined as: where ρ f and ρ s are density of fluid and solid, respectively, C f and C s are concentration in mass fraction of a species in fluid and solid, respectively, and ϕ is volume fraction of fluid, all varying as a function of time and space.
In Equation 1, the Darcy flux, q D and q Cf diffusion flux in fluid are given by: where, P f is fluid pressure, and is chemical potential of the species in the fluid. Gravitational effects in the Darcy flux have been neglected. Furthermore, the transient permeability k and effective diffusivity ef f of the mobile species are defined as: where k 0 is background permeability, η f is fluid viscosity, is the diffusion coefficient of the species in the fluid, and the latter three parameters are assumed constant.
For this simplifying case, conservation of total mass of the system becomes: with total mass density defined as: Porosity evolution is governed by the time integrated mass conservation equation of an immobile species in the solid, ρ im (e.g., Beinlich et al., 2020;Malvoisin et al., 2015;Plümper et al., 2017): 4 of 44 where ϕ(x,t 0 ) and ρ im (x,t 0 ) are the initial values of porosity and immobile species density at t 0 and may vary in space.
In Equation 7, the expression for ρ n (x,t) in Malvoisin et al. (2015) and Plümper et al. (2017), where fluid was pure water, is replaced with ρ im (x,t), which denotes the mass density of any immobile species in the solid: The full set of equations above presents a system of three equations for eight unknowns: fluid and solid density (ρ f ), and (ρ s ), weight fraction of a mobile component in fluid and solid (C f ), and (C s ), chemical potential of the component in the fluid (μ Cf ), and the weight fraction of an immobile component in the solid ( im ), porosity (ϕ), and fluid pressure (P f ). All unknowns vary in time and space. Background permeability (k 0 ), fluid viscosity (η f ), and diffusion coefficient of mobile component C f in the fluid ( ) are given and kept constant. Porosity, fluid pressure, and either the fluid or the solid composition of each component in the system are not thermodynamically constrained and the three equations are solved for ϕ, P f , and C s . The remaining unknowns, solid and fluid compositions and densities, are found from thermodynamic relationships and in general vary as a function of temperature, pressure, and composition. The equations above imply the assumption of constant temperature, which may be justified in the case of a small modeling domain where temperature is nearly instantly homogenized. Furthermore, we neglect the effect of varying pressure on the thermodynamic properties arguably justified in the case pressure variations are kept small.
In complex mineral, melt, and fluid solutions, the nonideal mixing behavior leads to nonlinear thermodynamic relationships. Parameterization of the various variables to obtain analytical expressions to close the system of transport equations may be possible in some cases , but may be difficult or impossible due to the complex topologies and discontinuities of phase assemblage boundaries (Stixrude & Lithgow-Bertelloni, 2022). In that case these thermodynamic properties either have to be directly calculated (i.e., on-the fly) or a precomputed lookup table may be prepared to interpolate the needed values. Here, we compute these relationships with Thermolab, benchmarked with phase diagrams for rocks, fluids, and melts. Using MATLAB codes, the first part of this contribution documents the technical details of the equilibrium thermodynamic calculations and the treatment of solid solutions. This is considered essential to understand the effects on nonlinear transport processes. We include a description of methods to sufficiently resolve the equilibrium compositions and densities of the involved phases. This is needed to avoid numerical artifacts and instability in the transport codes. In the last part, we get back to a detailed example of the implementation and solution of the system of transport equations in which all quantities evolve through time. For a full list of symbols and notation see Table 1.

Gibbs Energy Calculation
The starting point of equilibrium thermodynamic calculations in Thermolab is the Gibbs energy of a mineral, melt, or fluid. All other required properties are derived from the Gibbs energy. The calculation of Gibbs energies requires in the first place an internally consistent thermodynamic data set for pure solids, fluids, gases, or aqueous species, often referred to as endmembers. Second, it needs solution models that describe the energy of mixing between endmembers dissolved in a phase, also referred to as mixing models or activity-composition relationships (e.g., Ganguly, 2020;Holland & Powell, 2003). Thermolab currently has built in several internally consistent thermodynamic data sets for minerals, melts, and fluids, the most extensive and up-to-date are the Holland and Powell endmember databases , 2011. The SUPCRT (dslop98) mineral database Johnson et al., 1992) is implemented to allow calculations with additional minerals not included in the Holland and Powell data sets. For water (and/or CO 2 ), several equations of state (EOS) are implemented, including the EOS of the International Association for the Properties of Water and Steam (IAPWS) (Wagner & Pruss, 2002, revised version, 2018, the EOS of Johnson and Norton (1991), Z. G. Zhang and Duan (2005), C. Zhang and Duan (2009), Pitzer and Sterner (1994), and the CORK EOS, from Holland and Powell (1991). For aqueous species, we have implemented the  formulation, and the Tanger and Helgeson (1988) formulation as implemented in SUPCRT92 (Johnson et al., 1992) using the 1998 database including more than 1,300 aqueous species. Furthermore, it includes the Deep Earth Water (DEW) model (Sverjensky et al., 2014) with the recent updates and additional species from Huang and Sverjensky (2019), Aranovich et al. (2020), and the Miron data set for aqueous species (Miron et al., 2016 the endmember basis for use in thermodynamics of mixtures and for the local equilibrium calculations. In the future, this can be extended with additional thermodynamic endmember data sets. The following documents the main code as shown in an example in Figure 1. The mathematical notation follows the convention that bold font variables present vectors, bold font capital variables represent matrices, and all others are scalars variables.
Meaning of symbols and variables are summarized in Table 1.

Figure 1.
Thermolab code example to calculate the Gibbs energies of a phase for a given composition at fixed temperature (T) and pressure (P). Here is an example for the amphibole model of Green et al. (2016). This code is a possible MATLAB translation of the equations and methods documented in the Gibbs energy calculation section of the main text. The first line in the code clears MATLAB memory and figure. In Lines 3 and 4, T and P are defined. Lines 5 and 6 specify, respectively, for which phase to calculate the Gibbs energy and the name of the Excel file in which the solution model data is stored. Line 7 specifies the name of the elements in the phase, and line 8 the corresponding composition in moles for which to calculate the Gibbs energy. Line 10 is a call to read the needed thermodynamic data. Lines 12 and 13 are needed for aqueous species and compute the density and dielectric constant of water. Lines 15-24 compute the Gibbs energy of the endmembers in the solution (e.g., Equations 9, 12 and 17). Proportions and site fractions are calculated in lines 26-33. Numerical errors are removed in lines 35 and 36. In Lines 29 and 30, a grid of all possible states of order-disorder is generated for the particular composition of the phase specified in line 8. For all these compositions the mechanical, ideal, and nonideal mixing energy is calculated in lines 38, 39, and 40, respectively. The nonideal Gibbs energy, line 40, is a call to an external function that comprises a collection of nonideal mixing functions that can be opted. Line 43 finds the Gibbs energy that is minimum and thus finds the state of order-disorder. Line 44 displays the corresponding proportions of the endmembers.

Endmembers
An internally consistent thermodynamic database contains the Gibbs energy of an endmember, for example, a mineral, gas, fluid, or melt species with a fixed composition at reference conditions ( 0 ). Alternatively, it holds enthalpy and entropy data from which the Gibbs energy of an endmember at reference conditions can be obtained: where 0 is enthalpy and 0 is entropy at reference temperature T r = 298.15°K and pressure P r = 10 5 Pa.
To calculate the Gibbs energy at elevated temperature the entropy is obtained by integration with respect to temperature at reference pressure (P = P r ): where c P is the heat capacity at constant reference pressure (P r ).
Then, the integral of entropy at reference pressure with respect to temperature is subtracted from the Gibbs energy at reference conditions to get the Gibbs energy at reference pressure and elevated temperature: To obtain the Gibbs energy at elevated pressure the volume, v 0 , at fixed temperature is integrated with respect to pressure and added to the Gibbs energy at reference pressure for each temperature: An additional temperature and pressure dependent excess Gibbs energy term ( 0 exc ) is added to account for Landau phase transitions, ordering reactions, or in case of aqueous species for Gibbs energy of solvation.

Solids
An example of a heat capacity expression to be used in Equation 10 consists of a four-parameter (a, b, c, and d) polynomial fit to experimental heat capacity data as used in the data set of Powell (1998, 2011): For the pressure dependence of Gibbs energy an equation of state (EOS) is used to relate volume to pressure. For example, the Murnaghan equation of state used by  (with K′ set to 4): where is bulk modulus at reference pressure and 0 is volume at reference pressure and given temperature.
Note that the reference pressure P r is added to ensure the limit of standard reference conditions. The EOS can be rearranged for volume: Substituting Equations 9, 10, 11, and 15 into Equation 12, and using a suitable heat capacity formula such as Equation 13, the integrals can be evaluated analytically or numerically. Lines 18-21 in Figure 1 show a MATLAB example of using Equations 9-12. In the example, the integrals are evaluated in a separate function. To use a different thermodynamic database, the heat capacity equation and EOS can be replaced with the appropriate expression. A number of important endmembers such as quartz are treated with additional volume and entropy terms to account for first or second order phase transitions such as heat capacity anomalies or order-disorder in the crystal lattice (e.g., sillimanite Holland and Powell (1996)). This then requires an excess Gibbs energy contribution, which in principle can also be calculated according to Equation 12 using different volume and heat capacities and integration limits (e.g., Berman & Brown, 1985). Alternatively, an expression for this additional energy is given . In the example in Figure 1, this is represented by a call to a MATLAB function (line 20 in Figure 1). For the specific details, we refer to the original papers documenting the internally consistent databases (Berman, 1988;, 2011Johnson et al., 1992).

Fluids
For molecular fluids such as H 2 O and CO 2 , the EOS, for example, Equation 14, is replaced by an appropriate fluid EOS. As water is one of the most important fluids on our planet, extensive work has been done on the thermodynamic formulation. The main difference between solids is that the fluid EOS, especially for water, usually cannot easily be rearranged for volume as there are multiple volumes possible for a single pressure in the region of coexisting fluid and gas. An example of this is given by the Modified Redlich Kwong (MRK) EOS, on which Holland and Powell (1991) base the fluid Gibbs energies to be compatible with the extensively used thermodynamically consistent data set of : Because this is a multivalued function a suitable algorithm must be used to calculate volume as function of P-T in the two-phase region (e.g., where gas and liquid coexist). In the single-phase region, the volume for the stable phase (gas or liquid) should be determined. Holland and Powell (1991) rearranged the equation as a cubic in volume after which the correct root must be found in each phase region. Rather than solving which phase (e.g., gas, liquid, or both) is stable, the different regions were predefined and the correct volume was selected based on the P-T conditions. The MRK formulation in Equation 16 is extended with a virial contribution and Holland and Powell (1991) give detailed instructions on how to calculate the volume as function of P-T for H 2 O, CO 2 , and several other COH species. The advantage of the CORK formulation is that it can be integrated analytically and details in the original paper of Holland and Powell (1991) are sufficient for reproducing the values for use in Thermolab. The more updated thermodynamically consistent data set of Holland and Powell (2011), which replaces the  data set, uses the EOS of Pitzer and Sterner (1994).
Calculating the Gibbs energy of fluids then follows Equations 9-12 as for solids, with an entropy integral that is consistent with the data set for solid, melt, and/or gas endmembers. The Powell (1998, 2011) data sets have parameters for heat capacity, fitted simultaneously with the minerals to ensure the Gibbs energy at room pressure and elevated temperature can be calculated consistently.
In principle, the CORK or Pitzer and Sterner (1994) EOS can be used in conjunction with the SUPCRT mineral database, however, for the entropy integral a thermodynamically consistent heat capacity formulation is needed to obtain valid thermodynamic calculations. The SUPCRT database does not contain heat capacity values for H 2 O. The specific details needed for this calculation in the original papers of Johnson and Norton (1991) and Johnson et al. (1992) could not be found. Similarly, the DEW spreadsheet does not contain H 2 O entries and combines Delany and Helgeson (1978) and an internal routine for use only above at least 0.1 GPa, likely as it is intended to be used for Deep Earth applications. Moreover, at low P and T, for example, for shallow processes, the CORK EOS should not be used (Holland & Powell, 1991). We found that by using the NIST Shomate heat capacity equations (Shomate, 1954), using parameters from the NIST website for liquid water and gas, the Gibbs energies retrieved from SUPCRT can be reproduced accurately. The simplicity of the formulation and up-to-date online documentation of the parameters is of advantage. For calculations at elevated P, the numerical integration of the IAPWS or the EOS of Johnson and Norton (1991) give satisfactory results at low T, whereas the CORK EOS can be used above 100°C (Holland & Powell, 1991).

Aqueous Species
Gibbs energies of aqueous species are calculated in Thermolab following the formulation of Tanger and Helgeson (1988) as also outlined in Johnson et al. (1992). The resulting Gibbs energies have been compared to the output from SUPCRT and are reproduced within ∼1 J precision. The DEW model (Sverjensky et al., 2014) and the Miron data set (Miron et al., 2016(Miron et al., , 2017 use the same formulation, with refitted parameters for some of the endmember species. Although fundamentally the aqueous species are also calculated according to Equation 12, there is a difference to the solid, melt, and fluid endmembers because the aqueous species need the density and dielectric constant of the solvent. The contribution to Gibbs energy of solvation therefore needs the density and dielectric constant of water for the Born equation (Figure 1, Lines 12-13). This contribution is included in the excess Gibbs energy term (Figure 1, Line 20). Furthermore, it should be noted that these Gibbs energies are usually on a molality based scale. Note that in the example for amphibole in Figure 1, the properties of water are irrelevant, but to maintain some generality they have been left in the example. When computing the Gibbs energy of aqueous species, we change the name of the phase from "Amphibole" into any desired aqueous species or fluid mixture defined in a solution model database file. In Thermolab, the user can choose from a number of EOS and dielectric constants to be used in the aqueous species Gibbs energy (Fernández et al., 1997;Johnson & Norton, 1991;Sverjensky et al., 2014). The aqueous species endmembers of  offer a restricted set of species compared to the SUPCRT data, but it has the advantage that it is fitted in the internally consistent data set of the mineral, melt, and gas endmembers of Powell (1998, 2011). Their formulation is based on the Anderson density equation (Anderson et al., 1991) and uses the CORK EOS for the water density. These aqueous species endmember data can also be used in fluid mixtures (Evans & Powell, 2007).

Dependent Endmembers
Endmembers can also be formed from a linear combination of several other endmembers ( Figure 1, Line 23): where the ν holds the stoichiometric reaction coefficients. A formation energy may be associated with such a reaction and this can be captured in any functional form in the excess Gibbs energy term of Equation 12.

Gibbs Energy of Mixtures
The Gibbs energy of a phase (mineral, gas, fluid, or melt species) that can form a mixture between several endmembers is represented by the sum of mechanical, ideal, and nonideal mixing energies ( Figure 1, Line 41): The mechanical part of the mixing is the sum of the Gibbs energy of all endmembers in the mixture (g 0 ), weighted by their proportions (p) (Figure 1, Line 38).
The ideal mixing Gibbs energy is given by a linear combination of the site fractions (z) in the mixture multiplied with their logarithm and the site multiplicity (m) (Figure 1, Line 39): where the symbol indicates elementwise multiplication of vectors and matrices.
Using the site fractions of the endmembers (Zt), the last term in Equation 20 ensures that the ideal Gibbs energy is zero in the limit of the pure endmember. This for example occurs when a solution model is defined with certain endmembers having site fractions between 0 and 1 (e.g., Appendix in Tajčmanová et al., 2009). How to obtain m, p, z, and Zt is explained in the following section.
Nonideal Gibbs energy can be calculated in the simplest way with a binary mixing formula: This is essentially a sum of multiplications of binary pairs of endmembers in the mixture multiplied with their interaction parameters. These interaction parameters can depend on composition to get an asymmetric mixing formulation, referred to as subregular mixing model (Ganguly, 2020). When p is replaced with a size parameter adjusted proportion and using the appropriate W the asymmetric formulism of Holland and Powell (2003) can be used. If the molar volumes of the endmembers are used as size parameter, it is essentially a Van Laar mixing model (Van Laar, 1906), as for example, used by Aranovich and Newton (1999). The W and any size parameters are in general also temperature and pressure dependent such as also the molar volume used in Van Laar mixing models. Any other regular or subregular solution models can be used in place of Equation 21 and ternary interaction terms may be added. In Thermolab, it is possible to expand the codes with a variety of mixing formulas for the nonideality by adding the appropriate formulation to the function called in line 40 in Figure 1.

Solution Models
The data for m and Z t are retrieved from a site occupancy table, stored for example, in an Excel spreadsheet.
In this spreadsheet, also data for W and any size parameters for asymmetric formalism models (Holland & Powell, 2003) and the model type are stored. In Figure 1, the data is loaded in the beginning of the code (Line 10, Figure 1) from a function that reads the Excel data. Site fractions z are obtained by multiplying the transpose of the site fraction speciation matrix with a column vector of proportions: After specifying proportions of the endmembers in a solution, the Gibbs energy of that particular phase and composition at a given P and T can be calculated from the above.

Site Occupancy
The site occupancy of the mixture can be represented by a table listing the occupancy of crystallographic sites for each endmember. For example, a binary olivine solution with only one crystallographic site M, and a fixed silica tetrahedral site (which has composition SiO 4 ), is shown in Table 2. It can be represented by a matrix: The columns represent occupancy of each species on a crystallographic site for the endmember in the rows of the table. Site fractions z can be retrieved from this table by dividing each site occupancy over the sum of moles of the species on that site. Unless defined otherwise by the authors of a solution model, the sum of moles of all species on each site gives the site multiplicity, and can be represented by a vector: As the site multiplicity data is not always determined from S t because of decisions by the authors of the solution model, we store the multiplicity values in the spreadsheet. Examples where the theoretical site multiplicity is replaced by a different value is for example the T1 site in Amphibole (Diener et al., 2007), where the authors have taken 1 instead of 4 (Appendix A).  Elementwise division of each row in Table 2 or in matrix Equation 23 with  this multiplicity vector leads to a table of site fractions represented by the  matrix (Table 3): Using Equation 22 and spelling this out in matrix-vector notation, it is essentially a set of equations where the coefficients in front of the proportions are directly read from the site fraction table (Table 3) and then transposed: The site fractions become:

Ordering
In fact, in olivine, Mg and Fe can interchange on two molecular sites M1 and M2 (Cemič, 2005). This exchange can take place while total chemistry of the olivine is fixed. To account for variations in distribution of Mg and Fe over the M1 and M2 sites a third endmember to account for the ordering on the sites is introduced. Introducing the endmember "Ordered Olivine" which has a composition composed of 0.5 mol forsterite and 0.5 mol fayalite and having an "ordered" configuration with all Mg on M1 and all Fe on M2. We may write the site speciation as shown in Table 4.
The Gibbs energy of this endmember can be made out of two endmembers, 0.5 mol of forsterite and 0.5 mol of fayalite and an energy of reaction may be added.
After this, we get three proportions that define the Gibbs energy. To reach the configuration with all Fe on M1 and all Mg on M2, a negative proportion of the ordered olivine must be used (Powell & Holland, 1999).
Using the speciation matrix and Equation 22 we get:

Worked Example
Using the multiplicities and site fractions for the olivine example without order-disorder we get using Equation 20: (Note that the T-site fraction cancels due to the logarithms because the site fraction is always 1).

From Equation 21
, the nonideal mixing energy becomes: In case we use a  model for olivine, the nonideal parameters are equal to each other and reduce to the simple binary symmetric mixing parabola: For the example with order-disorder the ideal mixing Gibbs energy becomes: Using Equation 21, with three endmembers, the nonideal Gibbs energy becomes:

Proportions From Site Fractions
The proportions can be written as function of site fractions by solving the system of equations in Equation 29 using the independent equations and the fact that proportions sum up to 1. Starting from the equation that proportions sum up to one, and adding sequentially only equations from Equation 29 that increase the rank gives a closed system of equations. This leads to: Which after solving for proportions gives: Hence, there are two independent site fractions together completely defining the composition and ordering of the olivine in this example.

Proportions From Composition
From the olivine Table 4, we can derive a list of compositions of each endmember by summing up the moles of each atom over the sites.
Again, this table can serve to set up a system of equations and may be added by the constraint that the sum of the proportions equals 1 (first equation): However, this table only has two independent equations (rank = 2). We need another equation to solve for all three proportions. A simple solution may be to treat the ordered endmember as a known variable on the right hand side.
Our extended system can be written: Starting from the first required equation and adding equations until we have reached the rank of the matrix gives us directly the first 2 and the last equation. This means the set of independent equations becomes: The matrix in Equation 39 can be inverted and used to find the expressions for the endmember proportions as function of composition: Instead of varying the ordered endmember proportion, it is also possible to vary one of the site fractions. For example, we can have the site fraction of Mg on M1 as variable in addition to the compositional variables. This way, we can fix the bulk composition of the mineral while varying the site occupancies due to ordering. To this end, we add the equation for site fraction of Mg on M1 (Equation 1 in Equation 29), to the original system of equations: Here only the first, second, and last equation are independent and so we have: Solved for the three unknown proportions: These equations guarantee a fixed olivine composition, while changing the distribution of Fe and Mg on M1 and M2. The site fraction equation that will describe this distribution, or the order-disorder, can be obtained by ), is convenient as it guarantees the minimum amount of compositionally independent variables. A MATLAB code that automates this procedure is used in Thermolab and detailed in Appendix A. This function gives the independent components for each solution model as well as the site fraction to vary as order parameter and the matrix of converting composition to proportion. A worked example for a more complex mineral, Amphibole with 11 endmembers is given in Appendix A. It is also possible to calculate site fractions from bulk composition as shown in the Appendix in Vrijmoed and Podladchikov (2015).

Aqueous Fluids
For mixtures of aqueous species, an approach similar to the ideal olivine example can be taken. Treating the fluid mixture with a one-site model in which all species mix in a similar fashion the mechanical and ideal mixture is calculated as in Equations 19 and 20, and for the nonideality, the Helgeson Kirkham Flowers (HKF) extended Debye-Hückel activity is used (Dolejš & Wagner, 2008;Helgeson et al., 1981). The Gibbs energy of the fluid mixture is on a molality scale when the formulation of Johnson et al. (1992) and the HKF activity model (Helgeson et al., 1981) is used and must be converted to a molar based scale to be compatible with the solids. This is taken care of by adding a correction term to the molality based Gibbs energy of the pure aqueous endmember as in Dolejš (2013): where 0 is the mole fraction and 0 is the molality based Gibbs energy, M w is the molar mass of H 2 O. We note that in Equation 12, it is unspecified what concentration scale the Gibbs energy is based on. This information must be known a priori and the conversion above must be applied appropriately.

Equilibrium Calculations Using Reactions
Gibbs energies calculated with the Thermolab codes are benchmarked with THERMOCALC (Holland & Powell, 2011), SUPCRT92 (Johnson et al., 1992), the DEW spreadsheet (Sverjensky et al., 2014), or data and phase diagrams from the relevant publications (Miron et al., 2016) and Perple_X (Connolly, 2005). Additionally, comparing with experimental data is essential to increase confidence in the method.
With the Gibbs energies of endmembers, equilibrium conditions for simple reactions can be calculated. Results for several solid state and fluid-solid reactions are shown in Figure 2. The curves are produced by making use of the equilibrium condition: Solubility of aqueous silica, by plotting the equilibrium constant of the dissolution of quartz into water using the  database. Comparison between data from Manning (1994) shows good agreement between model and experiment. (c) Comparison between experimental data and the SUPCRT (dashed lines) and Miron (solid lines) data set for the dissociation of KCl in water. This reproduces the figure shown by Miron et al. (2016), to demonstrate their improved fit to the data. (d) Activity diagram at room pressure, temperature, for basic weathering reactions using the equilibrium constants of the reaction calculated with Thermolab using the (dslop98) SUPCRT aqueous species and minerals data and Gibbs energy of water from NIST. Here a procedure to remove metastable extensions of reactions is performed in an automated way (i.e., automated Schreinemakers analysis). Gbs = Gibbsite, Ms = Muscovite, Mic = Microcline, Kln = Kaolinite, Prl = Pyrophyllite, SiO2 (am.) = amorphous silica.
with the change of Gibbs energy of the reaction given by: where the ν hold the stoichiometric coefficients in the reaction and g 0 are the endmember Gibbs energies of each phase in the reaction.
For pure phases (e.g., endmembers), the equilibrium constant K eq will be 1 and the logarithmic term disappears leaving only Δ 0 = 0. Using the endmember Gibbs energies and stoichiometric coefficient of reactions the contour line of Δ 0 = 0 can be plotted to visualize the reaction in P-T space (Figure 2a). In a simple dissolution reaction K eq will be equal to solubility and so we can use Δ 0 to calculate the solubility of quartz (Figure 2b). An example of a single dissociation reaction is shown in Figure 2c. By varying activities in K eq and using them as axis on a diagram the contour lines where the right hand side is 0 are used to obtain an activity-activity plot ( Figure 2d).

Gibbs Minimization in MATLAB
Using the ΔG of reaction, the metastable part of individual reactions will also be plotted (e.g., Figure 2a). For an equilibrium phase diagram, a Schreinemakers analysis or similar procedure may follow to draw the stable reactions (e.g., Figure 2d). Additionally, the phases and reactions to be calculated must be known beforehand. However, in many cases it can be desired to predict which reactions may take place. Instead of predicting the equilibrium line of a chemical reaction we calculate the stable phase assemblage at a point in P-T by determining which mineral has the minimum Gibbs energy. In the example code (Figure 3a), the calculation of Gibbs energy of the endmembers is done in line 9, a call to the main Thermolab function, which includes the code in Figure 1 to focus on the minimization algorithm. If all phases have the same composition (i.e., they are polymorphs), then the "min" function can be used (Figures 3a and 3b). Instead of locating the coexisting phases in a reaction in P-T-X space, the Gibbs minimization delivers the stable phases everywhere except for on the reaction line. Therefore, the reaction lines reflect the resolution of the P-T grid for which the Gibbs energies of each phase are calculated. At sufficiently high resolution, the reaction lines become smooth.
Natural chemical systems, for example, rocks, generally, are multicomponent systems. In this case, mass balance must be considered while finding the minimum of Gibbs energy. For a more thorough thermodynamic analysis of this, see for example, Connolly (2017). Constrained Gibbs energy minimization is employed in the majority of phase diagram calculation software (Connolly, 2009;Gordon & McBride, 1994). In MATLAB, the simplest way this can be achieved is using linear programming (function "linprog," e.g., Dantzig et al., 1955). With this function, minimization can be done under the constraints of mass balance equalities, which for completeness is shortly outlined in the following.
The function "linprog" is used to find the minimum Gibbs energy of the system by solving the following optimization problem, where α and g are vectors of the molar amount and the Gibbs energies of all phases, respectively. The "linprog" algorithm then searches for the components in α that make the minimum Gibbs energy of the system, while respecting the following mass balance equality: where n sys is a vector of total moles of each component in the system and N phs holds the molar composition of the phases as detailed below.
Equation 48 states that the sum of the amount of each phase multiplied by the composition of a component should equal the sum of that component in the system. A final requirement is that no negative amount α is allowed as that is physically meaningless: As a worked example we may consider the binary system SiO 2 -MgO. Possible phases that can be built from these components are quartz, periclase, enstatite, and forsterite. Spelled out for this case the equations read in matrix form: Figure 3. Two examples of Gibbs minimization to calculate phase equilibrium with Thermolab. The Gibbs energy calculation shown in Figure 1 is used here as a MATLAB function to improve clarity of the algorithm. (a) Code showing unconstrained Gibbs minimization used to calculate the aluminosilicate phase diagram. The Gibbs energies of the phases are obtained in line 9, using a unique identifier to distinguish them from the same phases in another database. Andalusite from the tc-ds55 THERMOCALC  data set (and, tc-ds55), kyanite (ky,tc-ds55), and sillimanite (sill,tc-ds55) in the order in which they are specified in line 5. The result is a 3D array with T in first dimension, P in second dimension, and phase in third dimension. Line 11 then finds the minimum in the third dimension, which results in a value (val) and index (ind) for each P-T. The index corresponds to the list of phases in line 5. (b) Resulting plot from code in panel (a). A color coded visualization of the phase index that gives the lowest Gibbs energy is shown using the pcolor function in MATLAB (code: line 13). The color bar indicates that the blue area has index 1, which then means it is the first phase in the list (line 5) that is stable (i.e., 1 = andalusite, 2 = kyanite, and 3 = sillimanite). The resolution of phase diagram can be increased by decreasing the step size of temperature and pressure (line 3 and 4). (c) MATLAB example to calculate chemical equilibrium using constrained Gibbs energy minimization. For the first lines see panel (a). In lines 7-9, the system composition is specified using total moles of SiO 2 and MgO. Then the phase compositions in moles is expressed in the matrix N phs with in each column the phase corresponding to the list in line 5. Each row holds the components (SiO 2 and MgO). This information is used as equality constraints in the function linprog which is called in line 14. It then finds the minimum Gibbs energy of the system out of the four possible phases, while satisfying the equality constraints. The second and third input arguments for "linprog" are empty as they are reserved for inequalities and the fourth and fifth arguments are for the matrix of coefficients and right-hand side vector, respectively. To restrict the linear programming search for the amount of each phase to positive values we input also the lower bound vector of zeros (LB), for example, Equation 49. Then alph will hold the amount of each phase that will make minimum system Gibbs energy while obeying system composition. The phases that have an alph above zero represent the equilibrium assemblage for the given P-T-X conditions. (d) G-X diagram in which the Gibbs energy of the phases are plotted against system composition (for plotting the Gibbs energy is normalized over the sum of the moles of component in each phase: i.e., the sum over the rows of N phs ). . Note that in Figure 3c single elements have been used as components instead of the oxides, which leads to the same result. After defining the lower bounds of alpha's as zero, linprog finds the alpha's that make the minimum of Gibbs energy. The results are plotted in Figure 3d and shows that for a system composition between enstatite and forsterite, alpha's of enstatite and forsterite are nonzero, which means that those are the stable phases and the magnitude of the alpha's give the molar amount of each stable phase in order of definition of the list of phases in Figure 3c, Line 5.
Performing a minimization at each point in a grid of pressure and temperature and evaluating at each P-T point the stable phase assemblage leads to a phase diagram in which also the amount of each phase is obtained from the mass balance constraint. Figure 4a shows that this produces the same result in case of the pure Al 2 SiO 5 system, as in the unconstrained minimization, but that it can now be used to do multicomponent systems such as dehydration of antigorite in SiO 2 , MgO, and H 2 O (Figure 4b).
Most minerals, fluids, and melts are not just pure phases, but can form mixtures of endmembers. The basic example above is shown because essentially the Gibbs minimization approach used in Thermolab is the same for systems in which phases occur that have a variable composition. In the following we shortly outline how we approximate the equilibrium in Thermolab to a reasonable degree by using the same linear programming method.

Linearization of Mixtures
We compute the Gibbs energy of a set of discrete compositions of the mixture and add them to the list of endmembers as discrete entities, being treated as a fixed composition phase. In the following, all these entities that are the consequence of discretizing a real phase that can form a thermodynamic mixture will be arbitrarily called phase-compounds. With this approach the minimization algorithm remains unchanged from the above case for pure phases in Figure 3c. Figure 5 shows a complete code example to do a calculation with mineral solid solutions using the approach described here.
An important step in preparing the mixtures for Gibbs minimization is to generate a grid holding a set of discrete phase compositions representing the mineral, fluid, melt, or gas. The compositions need to cover the full range of possible site fractions or proportions that can exist for a mixture. Fixing the pressure and temperature, Gibbs energy can be calculated either by specifying the site fractions or the proportions. The simplest is to use the site fractions for grid generation because they range from zero to one. Using proportions is more complex because, although proportions sum up to one for each mixture, they can also be negative. Thus, it requires knowledge of the range of values for proportions in each different solution. The grid generation is done in the function call in Line 12, Figure 5, by varying the site fractions by default from 0 to 1, and generating a multidimensional Cartesian grid over the correct number of independent site fractions, using the MATLAB function "ndgrid" inside. With this brutal way, many nonphysical site fractions are generated, but these are removed inside the grid generation function. For example, on a crystallographic site with three site fractions, always one of them is dependent because the sum should equal one. So, there are two independent site fractions that both can vary between 0 and 1. When creating a Cartesian grid with "ndgrid," it is possible to have a combi- Figure 5. Complete example code for the calculation of chemical equilibrium using Gibbs energy minimization including solid solutions, aqueous species, and pure phases in Thermolab for the case of serpentinite with addition of carbon (Beinlich et al., 2020). After defining a P and T of interest, the solution model Excel spreadsheet is selected (line 5). The chemical components in the system are defined in line 6 and the composition of the system is specified in line 7. The names of the phases to be considered in the calculation are specified in lines 8-10, where the tc-ds55 identifies the thermodynamic data set to be used (here THERMOCALC data set 55). The static thermodynamic data (td) for all phases are loaded in line 11 and stored in one structure array (this replaces the command in Figure 1 nation of the two independent site fractions being for example, 0.9 and 0.5. This is not a feasible composition because the sum of those fractions exceeds 1. So, all combinations in the grid that produced values with a sum of the site fractions above one are excluded. In summary, we just generate Gibbs energies for the full range of possible site fractions in the solution including all states of ordering and internal speciation. For amphibole or clinopyroxene from for example, Green et al. (2016), only maximal 3-6 discrete compositions per dimension can be used, when executed on a single PC. The resulting set of Gibbs energies and compositions is used in place of the example in Figure 3c. The call to the main code "tl_gibbs_energy," in Figure 5 Line 17, now includes the complete Gibbs energy code example Figure 1.
As an example of including solutions, we show a variety of phase diagrams, all produced with the same code for a discrete set of P-T values. The first example in Figure 6a, applies the olivine solution model as described above together with a binary melt mixing model, replacing the solid endmembers with their liquid equivalents (liquid forsterite and fayalite). For olivine, a grid is created where 1 and 2 vary independently and for the melt, the site fraction of forsterite liquid is varied. The Margules parameters in the melt are adjusted until they fit the experimental data (Bowen & Schairer, 1935). The T-X diagram including CO 2 fluid using the mixing model of Aranovich and Newton (1999) is produced by changing the bulk CO 2 in the system. When the amount of fluid components (H 2 O) in the calculation is orders of magnitude higher than the solid components, the diagram approximates the T-X CO2 diagram as the system is largely dominated by fluid. The results show good agreement with the experimental data from Aranovich and Newton (1999). A more complex calculation for a metapelite in KFMASH is done as a benchmark with the Perple_X and THERMOCALC diagrams and shows good agreement, improving the confidence of the method. Note that in complex systems, it becomes increasingly more demanding to get high-resolution phase diagrams. For example, a computation with the melt model of Holland et al. (2018) yields a low resolution diagram. Additionally, Cr 2 O 3 , Fe 2 O 3 , and K 2 O were omitted. However, the topology of the diagram is overall similar. In particular, the melt and pyroxene models used in the reproduction of phase diagram for basalt melting are computationally challenging with the method of linearizing the solutions with discrete hypothetical phase compounds.

Postprocessing the Minimization Results
To precompute a lookup table for the local equilibrium the unknown variables need to be retrieved from the minimization results. The results of multiple computations can be saved in one data set to be postprocessed separately. An example for one "linprog" result is given in Figure 7.
Using the discrete phase compound linearization approach, mixtures (solid solutions, fluids, melts, and gas mixtures) are split into individual phases having a fixed composition. The minimization algorithm will determine which of those discrete phases are stable, and for mixtures this often results in multiple discrete phase compounds with different composition of that mixture being stable. This is a mathematical consequence of solving the optimization program and the thermodynamic meaning is that we have a divariant field in which the composition of mixtures may change. Thermodynamically, it can happen that two or three phases are stable as a result of a miscibility gap. However, this needs to be determined by an algorithm that distinguishes the discrete phase compound compositions from each other. If they are significantly different for a given resolution of the discretization, then they are true separate phases. In the other case, the properties are obtained from a weighted average into one composition for the true stable phase. Afterwards, only the stable phase amount and the composition are retained. A clustering algorithm is used to determine if we have multiple or single phases stable for each particular mixture (Figure 7, Line 10).
Considering a given equilibrium calculation with linprog, the main result is a vector alpha, which holds the total amount of mole of each phase that has been considered in the system. The clustering algorithm first takes care of all distinct phases and removes any phases with zero alpha (i.e., those phases are found not stable by the linprog algorithm).
Mole fraction amount of the ith stable phase in the system is found by normalizing the molar amount of stable phase by the total: where np is the number of stable phases in this case.
Volume fraction ϕ of phase i is then found from molar volume and mole fraction amount: where the molar volume of each phase can be found using the numerical derivative of the Gibbs energy: The weight fraction of each phase is found from: in which the molar mass of all phases can be found by: where molm holds the molar masses of the elements.
The density of phase i can then be computed: Composition in weight is found from: Chemical potential is calculated following: where μ nc is the chemical potential of the last dependent proportions (as proportions sum up to 1) and is calculated from:

Improving the Calculations
Increasing the resolution of the P-T grid improves the smoothness of the lines in the diagram, but it does not ultimately lead to a better calculation. This is illustrated in Figure 8, where the resolution of the T-X grid for the olivine diagram is kept constant, while the number of discrete compositions for which the olivine and melt are calculated is increased. This means that the mixture is better approximated, and the result is that the compositions and modal abundances become better resolved and the step-like behavior that is present in the calculation with lower compositional grid gradually disappears. Although this is feasible for calculations that involve one or two simple mixtures, like binary olivine and melt, the computational demand required for more complex chemical systems and minerals makes such high-resolution compositional grids unpractical for computation. Even if  computers with abundant memory are employed, we found that linprog will not find a solution for systems with more than 2 million discrete phases in the minimization.
Another option to improve the calculation is therefore to refine compositions in an iterative approach as also suggested and employed since decades (Connolly, 2009;Rossi et al., 2009;W. B. White et al., 1958). In Thermolab, this is open for the user to implement or improve. Here, we describe a simple approach that currently is employed as a working example. First, an initial minimization is done with a suitable compositional starting grid. This initial starting grid is also subject to possible modification and improvement; however, it is possible to simply set up a coarse resolution initial grid so that each minimization is fast. Each discrete phase which has an alpha above zero is recognized as a stable phase. And, the index in the list of nonzero alpha is used to find the composition of the mixture in terms of the endmember proportions and site fractions. In the next step, the proportions are used to generate a denser (i.e., higher resolution) compositional grid around the stable proportions found in the first minimization. For a binary, a simple graphical explanation of this procedure in Thermolab is illustrated in Figure 9. The new higher resolution compositional grid is added to the already existing grid and taken to the next iteration. We checked that the Gibbs energy of the system is always decreasing during the iterations. The iterations are stopped when the Gibbs energy of the system is not changing anymore within limits of machine-precision, or when the compositional spacing in the refined grid reaches this tolerance. A number of parameters are manually set by the user to control the initial refinement window, the number of refined discrete compounds generated, and the factor with which this window is decreased during each iteration. This leads to improved calculation of compositions, modal abundances, as shown in Figure 10 for the metapelite (KFMASH) P-T diagram from Figure 6c.

Investigating Nonlinear Transport Processes in Open Systems
Coupling local equilibrium thermodynamics to transport processes such as the system of equations introduced above has been discussed previously (e.g., Malvoisin et al., 2015). The method described here follows the approach employed in recent studies of Malvoisin et al. (2015), Plümper et al. (2017), and Beinlich et al. (2020).
For the transport models a complete set of equilibria for all possible external conditions, T, P, X, may be precomputed and stored in a lookup table. To this end, loops over P, T, and X can be programmed around the linprog minimization and stored in a database. As an example, we show the soapstone formation in serpentinite from Beinlich et al. (2020). To the initial bulk serpentinite composition, carbon (C) is added to produce a lookup table of 150 different bulk compositions at fixed T and P. Employing the refinement method the Beinlich et al. (2020) thermodynamic calculations can be improved in smoothness and computation speed (Figures 12b, 12d and 12f). However, the main results are the same. The important aspect is that the Beinlich et al. (2020) calculation was more robust, because no iterative refinement was used and hence any chance of reaching a local-minimum rather than the true minimum was avoided. Thus, the results of Beinlich et al. (2020) served as benchmark for the refinement algorithm described above.

Reactive Transport Example
Regardless of the assumptions or the numerical method used to solve the equations, the thermodynamic properties need to be calculated using either a computation on the fly, a precomputed lookup table (Malvoisin et al., 2015), or by parameterization . After retrieving the density of solid, fluid, and concentration of fluid and solid, and the chemical potential, the system of equations described above can be solved. Here, a possible numerical implementation in MATLAB is shown (Figure 11), using explicit finite differences and a lookup table approach. The purpose is to focus on the concept of coupling the transport processes to the precomputed local equilibrium table to show how the effect of the nonideal solution models can be taken into account in reactive transport. The precomputed thermodynamic relationships are loaded in the beginning of the code. After setting up the physical and numerical parameters and defining a Cartesian grid the initial conditions. The simulation starts with a constant initial porosity distribution throughout the domain. Fluid pressure is initialized such that flow would start from the left boundary toward the right into the model. Initial solid concentration is chosen such that everywhere in the domain a serpentinite is stable and only the boundary on the left has a different composition. Such that a fluid with a composition out of equilibrium with the rock on the right will enter the model. The fluid reacts with the initial rock to form new mineral assemblages progressively transforming the rock from the left to the right in the model.
The unknowns, ρ f , ρ s , im , C f, and μ Cf are found from interpolation in the lookup table at each time step (Figure 11, Lines 36-40) at a given P-T (0.3 GPa and 300°C). Here, the fluid only contains CO 2 and for example, Mg can be used as immobile species in the solid im to find porosity according to Equation 7. From this, transient permeability and effective diffusion coefficient can be computed and used in the Darcy and diffusion fluxes. These are used in the total mass balance (5) and total mass concentration balance (1) to update concentrations and fluid pressure.
In the example, the fluid pressure is found from finding the steady state solution of the total mass balance using an iterative pseudotransient method. The total mass concentration C tot is found from an explicit finite different formulation. From this new C tot , Equation 2 is rearranged to find the updated C s . In the postprocessing stage the thermodynamic lookup table can also be used to find the modal abundances of the minerals to plot the evolution of the mineralogy through time.
Details of the coding example are found in the caption of Figure 11 and results of the transport code are shown in Figure 12. The nonlinearity of the thermodynamic relationship between solid and fluid composition is demonstrated in Figures 12c, 12d and 12f. The shape of the curves is a result of both the nonideality of solution models as well as changes in the phase assemblage. Depending on the function of solid composition versus fluid composition, the reaction fronts behave very differently in terms of shape and velocity (Figures 12a, 12c and 12e). In the case that the incoming fluid composition is low in CO 2 in equilibrium with antigorite, talc and magnesite no sharp fronts develop (Figure 12a) but modal abundance may change. When incoming CO 2 fluid composition is higher and antigorite disappears from the system, a sharp reaction front propagates, and when the CO 2 in the fluid increases further, talc also disappears and a second reaction front appears. The shape of the reaction fronts is controlled by the curvature of the C f -C s relation as well as the advection velocity versus speed of diffusion.

Cahn-Hilliard Exsolution
Natural processes usually are not in equilibrium and do not strictly follow the path that can be followed from a phase diagram. Instead, there will be a process that develops toward the state of global equilibrium depicted on a phase diagram although it may never reach it (e.g., due to low temperature effectively arresting diffusion, or other kinetic effects). With Thermolab, the Gibbs energies can be directly used instead of first calculating thermodynamic equilibrium and the process toward equilibrium can be modeled. A demonstration of this is given by solving the Cahn-Hilliard equations in which driving force for diffusion is chemical potential and an uphill diffusion process causes an initial random homogeneous system to develop into an equilibrium phase assemblage (Cahn & 10.1029/2021GC010303 29 of 44 Figure 11. Hilliard, 1958). A simplified version of the Cahn-Hilliard equation consists of a balance of mass concentration and a flux equation for the diffusion. Assuming a system without fluid, and assuming constant densities a basic mass concentration balance can be written (Figure 13, Line 36): where is the concentration of species A in the solid.
Flux of concentration in solid of a component can be defined as a function of gradients in chemical potential differences or a single chemical potential in conjunction with the Gibbs-Duhem relation (e.g., p. 80 in Lebon et al., 2008;Nauman & He, 2001), here the former is used (Figure 13, Line 35): The chemical potential difference between species A and B in the solid, − , can be conveniently expressed by: where the last term in Equation 64 is the interfacial free energy contribution as introduced by Cahn and Hilliard (1958). Similar to Figure 7, the chemical potential can be calculated with numerical differentiation of the Gibbs energy (see Figure 13, Lines 30-32).
These equations describe the spinodal decomposition of a phase with strong nonideal behavior, such as feldspar, due to nonlinear diffusion driven by chemical potential gradients. A useful review is found in Nauman and He (2001). See also caption Figure 14 for details.

Discussion
The motivation behind the development of Thermolab is to study the effects of the nonideality of solution models in reactive transport processes in open systems. Most natural materials including minerals, rocks, melt, fluids, and gases display some degree of nonideal mixing behavior (e.g., Ganguly, 2020). We showed that for the case study of the soapstone formation studied by Beinlich et al. (2020) both the shape and the velocity of the reaction front vary strongly depending on the nonlinearity of the partitioning of carbon between fluid and solid ( Figure 12). For evaluating potential risks during transport of nuclear waste material understanding this nonlinear behavior is crucial and a study of Shao et al. (2009) is line with this conclusion.
For mathematical analysis, the formulation of complex solution models is important to show in a transparent manner such that nonlinearities in transport processes can be studied. Therefore we presented the linear algebraic approach that is used in Thermolab to compute Gibbs energy of mixing for arbitrary multicomponent phases (Figure 1). The starting point was a crystallographic or structural model of a solution and a predefined set of endmembers as developed by previous workers (e.g., Green et al., 2016;Palin et al., 2016). More in-depth  discussion on how to define new crystallographic and speciation models for solid solutions is given by Myhill and Connolly (2021), after which the model should then be fitted to experimental data.
To ensure the reliability of the Gibbs energy calculation and equilibrium calculation methods, Thermolab is benchmarked versus phase diagrams and exemplified in use for reactive transport projects. Constrained Gibbs energy minimization is discussed following existing approaches (Connolly, 2005;W. B. White et al., 1958) and the results are then used to calculate equilibrium compositions and retrieve the so-called isotherm that can be used in a transport code.
The isotherm, introduced in chromatographic studies of metasomatism, is a relation between fluid and solid composition at fixed P-T (Hofmann, 1972). An approach in which the system is divided in solid and fluid and for which the mass conservations equations have been summed up to eliminate reaction source terms has proven very useful for reaction front propagation studies (Orr, 2005). Shape and velocity of propagating reaction fronts strongly depend on the nonlinearity of solutions models (Guy, 1993). The approach of precomputing the equilibrium compositions and using the results as an isotherm in the transport codes is useful to study the effects of the nonlinearity of solution models on reactive transport. As the studies of Hofmann (1972) and Guy (1993) focus on single phases, our results show this is functioning similarly on multiphase systems such as rocks. The steep reaction fronts act similar to isotherms with miscibility gaps (Fletcher & Hofmann, 1974), but the solid solutions are responsible for the continuity of the isotherm across the "jump-like" curves. Previous studies usually introduced simple isotherms to study the behavior during transport. Thermolab has been motivated by the need to generate realistic flow functions using complex solution models.
Solid solutions with strong nonideal behavior result in phase separation due to nonlinear diffusion when combined with transport models. The Cahn-Hilliard model (Cahn & Hilliard, 1958) is used to stabilize the numerical solutions by introducing an energy penalty for generating surfaces between phases introduced by a surface energy parameter. This method is widely applied in material science (Nauman & He, 2001) and also some geological applications have been studied ). For geological materials other than feldspar, the behavior of such nonlinear diffusion systems can be further investigated with Thermolab.
Mixing of databases is possible in a flexible framework like Thermolab, however, it may not be recommended since internally consistent databases are not necessarily consistent among each other. Nevertheless, to use aqueous species in fluids together with most up-to-date solution models of minerals, the only possibility to date is to combine for example SUPCRT databases with Powell (1998, 2011). It has been argued that such a combination may serve as a good approximation (Dolejš, 2013), when appropriate EOS and dielectric constant for water are used. However, refitting the aqueous species databases from SUPCRT in conjunction with mineral database of  is likely the more reliable approach to combine data sets (Miron et al., 2016(Miron et al., , 2017. The equilibrium relation between C f and C s (i.e., the isotherm). The steep slope at C s between ∼0.045 and 0.052 is enlarged in the inset. Here antigorite disappears (the soapstone on the left in panel (c)), and the steepness of the slope on the isotherm causes a sharp soapstone front. The shallow slopes below C s = 0.045 correspond to the transitional front where serpentinite is partially transformed to soapstone. (e) Development of an additional reaction front, forming listvenite (quartz-magnesite rock), as incoming fluid composition is now higher, ∼0.086 weight fraction dissolved carbon. (f) As the incoming fluid composition lies on another steep part of the isotherm, the listvenite front is also sharp. Inset shows the curvature of the isotherm corresponding to the listvenite stable assemblage. Endmember data from the tc-ds55 data set , solution models: CO 2 -H 2 O Fluid (Aranovich & Newton, 1999), Antigorite (Padrόon-Navarta et al., 2013), Chlorite , Talc , Magnesite, and Dolomite (R. W. White et al., 2003). t* = dimensionless time.
Figure 13. Cahn-Hilliard example code for 1D binary exsolution. The length of the domain and diffusion coefficient, assumed constant for simplicity here, are defined in lines 3-4, after which the thermodynamic data is loaded for a given P-T. The phase for which Gibbs energy is calculated is input in line 9, and the components in the system can be given in line 10. Line 11 initializes static thermodynamic data as in Figure 5. Line 12 calculates the Gibbs energy of the endmembers in the solution (this is done to increase performance as they do not vary here). Numerical parameters are setup in lines 14-18, followed by preprocessing for generating the grid and the time step. An initial setup consists of a homogeneous concentration with small random perturbations. The physical process is modeled in a time loop from lines 27 to 38. Lines 28-31 show how Thermolab is used without precomputed lookup tables. A small increment of proportions of the endmembers is used to make a numerical derivative to calculate the chemical potential in line 32, with Cahn-Hilliard addition of energy regularization added in line 33, corresponding to Equation 64. The chemical potential gradients drive diffusion in the system via the diffusion flux in line 35, Equation 63 and line 36 is the mass concentration balance in simplified form, Equation 62. Lines 36-37 are boundary conditions that ensure total concentration in the system remains constant. Here, making the concentrations equal in a mass conservative way (using the mean of the three grid points on the boundary), the chemical potentials will be the same at these nodes and hence it will lead to no flux boundary condition. Lines 41-45 are plotting the results as the model runs. Line 48 visualizes the coarsening of the separate exsolved phases through time.

Conclusions
With Thermolab it is possible to reproduce published phase diagrams involving complex solution models and these solution models can then be used in transport codes to investigate the effects of nonlinearity on the open system processes. Complex flow functions can be retrieved from Thermolab and used in mathematical analysis and numerical models of reactive transport. Gibbs energy can be directly used to construct chemical potentials as a driving force for nonlinear diffusion leading to phase separation.
The modeling framework provided by Thermolab allows users to add custom functions, include new databases, design solution models and also improve procedures to calculate equilibrium. Limits of the discretization of solutions approach are due to computer memory restrictions, which for complex solution models increases computation time to unpractical durations. Computers with increased memory would be a solution, however, we found that the linprog algorithm will stop to converge when about 2 million discrete phases or more in systems with more than six components are used. Refinement strategies are a potential solution to this problem; however, these compromise the robustness of the result by possibly missing the global minimum of the Gibbs energy. Thus, Thermolab leaves the door open for the development of faster and more robust future algorithms by a transparent open source coding environment, which due to the compactness of MATLAB example code can be translated to other programming languages with minimum effort.

Appendix A
An example of the automated selection of compositional or site fraction variables to calculate proportions for a mineral having 11 endmembers is given for amphibole: From the site occupancy table, we obtain the following site fraction table by dividing each occupancy over the sum of moles on the site. For example, the M2 site reaches a sum of 2 mol when the Mg, Fe Al Fe 3 , and Ti are summed. Then the moles of each occupancy are normalized by dividing by this sum to get fractions on the site. See Table A2 for the result. Figure 14. Example of using Thermolab without precomputed thermodynamic data. Physical processes develop toward thermodynamic equilibrium using the Cahn-Hilliard model at 500°C, 0.1 GPa. The 1D example code for a binary system is given in Figure 13. An initially homogenous concentration distribution in a single phase develops into a two phase system due to uphill diffusion as a result of the nonlinearity of solution models (the nonideality that leads to exsolution of phases).      The system of equations describing the relation between proportion and site fractions that is obtained automatically using the MATLAB script from Equation A1 is given in matrix form by: 1 1 1 1 1 1 1 1 1   1 1 0 1 1 1 1 1 1 Figure A1 shows the MATLAB script used to retrieve the matrix in Equation A2 above.
And from the inverse of the matrix in Equation A2 above the proportions can be found as function of site fractions:   Figure A1. Example of automated method for determining independent compositional variables, site fractions, and order-disorder variables for an amphibole with 11 endmembers .
Endmember proportions can also be obtained from the bulk composition using the table of moles of elements of the endmembers: The automatic selection of independent compositional variables and independent site fractions to describe ordering in the mineral starts from taking the first equation and adding from the top down only the equations that increase the rank of the system until the rank of the selected system of equations equals the number of endmember proportions. For Equation A4 above, this results in the following system of equations: