pyUserCalc: A Revised Jupyter Notebook Calculator for Uranium-Series Disequilibria in Basalts

1Introduction¶

Continuous forward melting models are necessary to interpret the origins of empirically-measured U-series isotopic disequilibria in basaltic lavas, but the limited and unreliable availability of reproducible tools for making such calculations remains a persistent problem for geochemists. To date, a number of models have been developed for this task, including classical dynamic melting after McKenzie (1985) and the reactive porous flow model of Spiegelman & Elliott (1993). There have since been numerous approaches to using both the dynamic and porous flow models that range from simplified analytical solutions e.g., Sims et al., 1999Zou & Zindler, 2000 to incremental dynamic melting calculators Stracke et al., 2003, two-porosity calculators Jull et al., 2002Lundstrom et al., 2000Sims et al., 2002, and one-dimensional numerical solutions to reactive porous flow Spiegelman, 2000 and dynamic melting Bourdon et al., 2005Elkins et al., 2019. Unfortunately, some of the approaches published since 1990 lacked publicly available tools that would permit others to directly apply the authors’ methods, and while the more simplified and incremental approaches remain appropriate for asking and approaching some questions, they are insufficient for other applications that require more complex approaches e.g., two-lithology melting Elkins et al., 2019. Other tools like UserCalc that were available to public users Spiegelman, 2000 were limited in application and have since become unavailable.

In light of the need for more broadly accessible and flexible solutions to U-series disequilibrium problems in partial melting, here we present a cloud-server hosted, publicly available numerical calculator for one-dimensional, decompression partial melting. The tool is provided in a Jupyter notebook with importable Python code and can be accessed from a web browser. Users will be able to access and use the tool using a free cloud server account, or on their own computer given any standard Python distribution. As shown below, the notebook is structured to permit the user to select one of two primary model versions, either classical reactive porous flow after Spiegelman & Elliott (1993) and Spiegelman (2000), or a new disequilibrium transport model, developed after the appendix formulas of Spiegelman & Elliott (1993). The new model ranges from pure disequilibrium porous flow transport (i.e., the mass-conserved equivalent of true fractional melting over time) to a “scaled” disequilibrium scenario, where the degree of chemical equilibrium that is reached is determined by the relationship between the rate of chemical reaction and the solid decompression rate (which is, in turn, related to the overall melting rate), in the form of a Damköhler number.

This scaled disequilibrium model resembles the classic dynamic melting model of McKenzie (1985), with the caveat that ours is the first U-series melting model developed for near-fractional, disequilibrium transport where mass is also conserved within a one-dimensional melting regime. That is, rather than controlling the quantity of melt that remains in equilibrium with the solid using a fixed residual porosity, the melt porosity is controlled by Darcy’s Law and mass conservation constraints after Spiegelman & Elliott (1993), and the “near-fractional” scenario is simulated using the reaction rate of the migrating liquid with the upwelling solid matrix.

2Calculating U-series in basalts during mass-conserved, one-dimensional porous flow¶

2.1Solving for equilibrium transport¶

Here we consider several forward melting models that calculate the concentrations and activities of U-series isotopes (²³⁸U, ²³⁰Th, ²²⁶Ra, ²³⁵U, and ²³¹Pa) during partial melting and melt transport due to adiabatic mantle decompression. Following Spiegelman & Elliott (1993), we start with conservation of mass equations for the concentration of a nuclide $i$, assuming chemical equilibrium between melt and solid:

where $t$ is time, $c_i^f$ is the concentration of nuclide $i$ in the melt, $D_i$ is the bulk solid/liquid partition coefficient for nuclide $i$, $\rho_f$ is the density of the fluid and $\rho_s$ is the density of the solid, ϕ is the porosity (local volume fraction of melt), $v$ is the velocity of the melt and $V$ the velocity of the solid in three dimensions, $\lambda_i$ is the decay constant of nuclide $i$, and $(i-1)$ indicates the radioactive parent of nuclide $i$ (see Table 1). Equation (1) states that the change in total mass of nuclide $i$ in both the melt and the solid is controlled by the divergence of the mass flux transported by both phases and by the radioactive decay of both parent and daughter nuclides (i.e., the right hand side of the equation above).

The equilibrium model of Spiegelman & Elliott (1993) assumes that complete chemical equilibrium is maintained between the migrating partial melt and the solid rock matrix along a decompressing one-dimensional column. To close the equations, they assume that melt transport is described by a simplified form of Darcy’s Law for permeable flow through the solid matrix. In one dimension, for a steady-state upwelling column of melting mantle rocks, they defined the one-dimensional melt and solid velocities ($w$ and $W$, respectively), and expressed the melt and solid fluxes as functions of height ($z$) in terms of a constant melting rate $\Gamma_0$:

where $W_0$ is the solid mantle upwelling rate, and $\Gamma_0$ is equivalent to $\rho_s W_0 F_{max}$ divided by the column height $h$ for a maximum degree of melting $F_{max}$.

Assuming an initial condition of secular equilibrium, where the initial activities $\lambda_i c_{i,0}^f D_i$ are equivalent for parent and daughter nuclides, they derived a system of differential equations for the concentration $c_i^f$ in any decay chain, which can be solved numerically using equation (10) from Spiegelman & Elliott (1993):

where $c_i'$ is the scaled melt concentration $(= c_i^f/c_{i,0}^f)$, ζ is the dimensionless fractional height in the scaled column, equal to 0 at the base and 1 at the top, and

is the effective velocity for element $i$.

In their appendix, Spiegelman & Elliott (1993) developed the more general (and, arguably, realistic) form where Γ and $D_i$ are functions of height $z$. The UserCalc model of Spiegelman (2000) then formulated a one-dimensional numerical integration for the concentrations of selected U-series isotopes in continuously produced partial melts with height $z$, after the equilibrium formulas above. The concentration expression derived by Spiegelman (2000) for the equilibrium scenario (formula 6 in Spiegelman (2000)) is:

where $F$ is the degree of melting. Spiegelman (2000) further observed that solving for the natural log of the concentrations normalized to the initial concentration of $i$, $U_i$, rather than the concentrations themselves, is more accurate, particularly for highly incompatible elements (formulas 7-9 in that reference). This is because log concentrations change linearly during melting, rather than exponentially, and are more numerically stable to calculate.

For the formulas above, Spiegelman (2000) defined a series of variables that allow for simpler integration formulas and aid in efficient solution of the model, namely

and substituting from the formulas above

where $D_i^0$ is the initial bulk solid/melt partition coefficient for element $i$, $R_i^{i-1}$ is the ingrowth factor, and $\alpha_i^0$ is the initial degree of secular disequilibrium for nuclide $i$ in the unmelted solid.

$U_i(z) = \ln(c_f(z)/c_f^0)$, the log of the total concentration of nuclide $i$ in the melt, can then be decomposed into

where

is the log concentration of a stable nuclide with the same partition coefficients, and $U^{rad}_i(z)$ is the radiogenic ingrowth component. An alternate way of writing the radiogenic ingrowth component of equation (9) of Spiegelman (2000) is:

where

is the decay constant of nuclide $i$, scaled by the solid transport time ($h/W_0$) across a layer of total height $h$. Note equation (17) is solved over a column of dimensionless height 1 where $\zeta\in[0,1]$.

Using these equations, the UserCalc reactive porous flow calculator accepted user inputs for both $F(z)$ and $D_i(z)$. The method further uses a formula for the melt porosity ($\phi(z)$) based on a Darcy’s Law expression with a relative permeability factor (formula 20 from Spiegelman (2000)):

where $K_r(z)$ is the relative permeability with height $z$, $A_d$ is a permeability calibration function, and $n$ is the permeability exponent. The permeability exponent for a tube-shaped fluid network is expected to be $n = 2$, while for a sheet-shaped network it is $n = 3$; recent measurements of the permeabilities of experimental magmatic melt networks suggest realistic magma migration occurs in a manner intermediate between these two scenarios, with $n = 2.6$ Miller et al., 2014. The relative permeability $K_r$ is calculated with respect to the permeability at the top of the column, i.e. depth $z=z_{final}$:

and allows for locally enhanced flow (e.g., mimicking the effects of a relatively low viscosity fluid).

Our model implementation reproduces and builds on the prior efforts summarized above, using a readily accessible computer language (Python) and web application (Jupyter notebooks).

2.2Solving for complete disequilibrium transport¶

We further present a calculation tool that solves a similar set of equations for pure chemical disequilibrium transport during one-dimensional decompression melting. This model assumes that the solid produces an instantaneous fractional melt in local equilibrium with the solid; however, the melt is not allowed to back-react with the solid during transport, as it would in the equilibrium model above. In the limiting condition defined by stable trace elements (i.e., without radioactive decay), the model reduces to the calculation for an accumulated fractional melt. The model solves for the concentration of each nuclide $i$ in the solid $(s)$ and liquid $(f)$ using equations (26) and (27) of Spiegelman & Elliott (1993):

which maintain conservation of mass for both fluid and solid individually, and do not assume chemical equilibration between the two phases. As above, the equations depend on $F(z)$ and $D_i(z)$, i.e. melt fractions and bulk rock partition coefficients that can vary with depth.

As above, the solid and fluid concentration equations are rewritten in terms of the logs of the concentrations:

and thus

We assume that initial $c_{i,0}^s = D_{i,0} c_{i,0}^f$. As above, the log concentration equations can be broken into stable and radiogenic components, where the stable log concentration equations are:

which are equivalent to a model for a fractionally melted residual solid and an accumulated fractional melt for the liquid.

Reincorporating this with the radiogenic component and scaling all distances by $h$ gives the dimensionless equations:

2.3Solving for transport with chemical reactivity rates¶

The two models described above are end members for complete equilibrium and complete disequilibrium transport. For stable trace elements, these models produce melt compositions that are equivalent to batch melting and accumulated fractional melting e.g., Spiegelman & Elliott, 1993. However, the actual transport of a reactive fluid (like a melt) through a solid matrix can fall anywhere between these end members depending on the rate of transport and re-equilibration between melt and solid, which can be sensitive to the mesoscopic geometry of melt and solid e.g., Spiegelman & Kenyon, 1992. In an intermediate scenario, we envision that some reaction occurs, but chemical equilibration is incomplete due to slow reaction rates relative to the differential transport rates for the fluid and solid. If reaction times are sufficiently rapid to achieve chemical exchange over the lengthscale of interest before the liquid segregates, chemical equilibrium can be achieved; but for reactions that occur more slowly than effective transport rates, only partial chemical equilibrium can occur e.g., Grose & Afonso, 2019Iwamori, 1993Iwamori, 1994Kogiso et al., 2004Liang & Liu, 2016Peate & Hawkesworth, 2005Qin et al., 1992Yang et al., 2000. Such reaction rates can include, for example, the rate of chemical migration over the distance between high porosity veins or channels e.g., Aharonov et al., 1995Jull et al., 2002Spiegelman, 2000Stracke & Bourdon, 2009; or, at the grain scale, the solid chemical diffusivity of elements over the diameter of individual mineral grains e.g., Qin et al., 1992Feineman & DePaolo, 2003Grose & Afonso, 2019Oliveira et al., 2020Van Orman et al., 2002Van Orman et al., 2006.

To model this reactive transport scenario, we start with our equations for disequilibrium transport in a steady-state, one-dimensional conservative system, and add a chemical back-reaction term that permits exchange of elements between the fluid and the solid. The reaction term is scaled by a reactivity rate factor, $\Re$ and expressed in kg/m³/yr. (i.e., the same units as the melting rate). The reactivity rate thus behaves much like the melting rate by governing the rate of exchange between the solid and liquid phases, effectively scaling the degree to which chemical exchange can occur. This new term could simulate a number of plausible scenarios that would physically limit the rate of chemical exchange by transport along a given distance in a linear manner, such as the movement or diffusion of nuclides through the porous solid matrix between melt channels a given distance apart.

First, returning to the conservation of mass equations for a steady-state, one-dimensional, reactive system of stable trace elements, and using $\Gamma(z)$ to represent the melting rate:

where, for an adiabatic upwelling column,

From this, the equations (29) and (30) can be integrated (with appropriate boundary conditions at $z=0$) to give

Next, we expand the concentration equations to include the reactivity factor, and substitute the conservation of total mass determined above:

If we then combine the $\Gamma(z)$ terms and rearrange:

We can now divide the fluid and solid equations by $c_i^f$ and $c_i^s$, respectively, and rearrange the $W_0$ terms:

The first terms on the right-hand side of each of these equations are identical to pure disequilibrium melting, such that if $\Re$ is zero, the equations reduce to the disequilibrium transport case of Spiegelman & Elliott (1993).

To solve, the final terms that involve the reactivity factor can be further rewritten using the definitions for $U_i^f$ and $U_i^s$:

Thus:

and:

Finally, substituting adiabatic upwelling and scaling depth by $h$, and adding radioactive terms gives the full solutions for the dimensionless equations $dU_i/d\zeta$:

where $h$ is the total height of the melting column.

2.3.1The Damköhler number¶

The dimensionless combination

is the Damköhler number, which governs the reaction rate relative to the solid transport time. Damköhler numbers more generally are used to relate the timescales of chemical reactions to the rates of physical transport in a system. If re-equilibration is limited by solid state diffusion, $\Re$ can be estimated using:

where ${\cal D}_i$ is the solid state diffusivity of element $i$, and $d$ is a nominal spacing between melt-channels (this spacing could, for example, be the average grain diameter for grain-scale channels, or 10 cm for closely spaced veins).

In this case (which we will assume for this paper), the Damköhler number can be written

Substituting the definition of $Da$ above yields the final dimensionless ODEs for the disequilbrium transport model:

with initial conditions $U_i^s = U_i^f = 0$.

In the limit where the Damköhler number approaches zero, the above formulas reduce to pure disequilibrium transport, whereas if $Da$ approaches infinity (i.e., infinitely fast reactivity compared to physical transport), the system approaches equilibrium conditions ($c_i^s \rightarrow D_i c_i^f$).

2.3.2Initial conditions¶

Inspection of equation (53) shows that for the initial conditions described above and $F(0)=0$, $\frac{dU_i^f}{d\zeta}$ is ill-defined (at least numerically in a floating-point system). However, taking the limit $\zeta\rightarrow 0$ and applying L’Hôpital’s rule yields

where

The initial radiogenic term also uses the limit from equation (34):

Rearranging equation (55) gives the value for $U_i^{'f}(0)$ for $F=0$ as

3A pyUserCalc Jupyter notebook¶

3.1Code design¶

The UserCalc Python package implements both equilibrium and reactive disequilibrium transport models and provides a set of code classes and utility functions for calculating and visualizing the results of one-dimensional, steady-state, partial melting forward models for both the ²³⁸U and ²³⁵U decay chains. The code package is organized into a set of Python classes and plotting routines, which are documented in the docstrings of the classes and also demonstrated in detail below. Here we briefly describe the overall functionality and design of the code, which is open-source and can be modified to suit an individual researcher’s needs. The code is currently available in a Git repository (https://gitlab.com/ENKI-portal/pyUsercalc), and any future edits or merge requests will be managed through GitLab.

The equilibrium and disequilibrium transport models described above have each been implemented as Python classes with a generic code interface:

'''
Interface:
----------
    model(alpha0,lambdas,D,W0,F,dFdz,phi,rho_f=2800.,rho_s=3300.,method=method,Da=inf)

Parameters:
-----------
    alpha0  :  numpy array of initial activities
    lambdas :  numpy array of decay constants scaled by solid transport time
    D       :  Function D(z) -- returns an array of partition coefficients at scaled height z
    W0      :  float -- Solid mantle upwelling rate
    F       :  Function F(z) -- returns the degree of melting F  
    dFdz    :  Function dFdz(z) --  returns the derivative of F  
    phi     :  Function phi(z) --  returns the porosity  
    rho_f   :  float -- melt density
    rho_s   :  float -- solid density
    method  :  string -- ODE time-stepping scheme to be passed to solve_ivp (one of 'RK45', 'Radau', 'BDF')
    Da      :  float -- Damköhler Number (defaults to \inf, unused in equilibrium model)

Required Method:
----------------
        model.solve(): returns depth and log concentration numpy arrays  z, Us, Uf
'''

which solves the scaled equations (i.e., equation (9) or equations (53) and (54)) for the log concentrations of nuclides $U_i^s$ and $U_i^f$ in a decay chain of arbitrary length, with scaled decay constants $\lambda^\prime_i$ and initial activity ratios $\alpha_i^0$. In the code, we use the variable $z$ for the scaled height in the column (i.e. $z\equiv \zeta$), and the model equations assume a one-dimensional column with scaled height $0\leq z \leq 1$. The bulk partition coefficients $D_i(z)$, degree of melting $F(z)$, melting rate $dF/dz(z)$, and porosity $\phi(z)$ are provided as functions of height in the column. Optional arguments include the melt and solid densities $\rho_f$ and $\rho_s$, the Damköhler number $Da$, and the preferred numerical integration method (see scipy.integrate.solve_ivp). Some of these variables, such as $D_i(z)$ and $F(z)$, are provided by the user as described further below, and are then interpolated using model functions.

UserCalc provides two separate model classes, EquilTransport and DisequilTransport, for the different transport models; the user could add any other model that uses the same interface, if desired. Most users, however, will not access the models directly but rather through the driver class UserCalc.UserCalc, which provides support for solving and visualizing column models for the relevant $^{238}U$ and $^{235}U$ decay chains. The general interface for the UserCalc class is:

''' A class for constructing solutions for 1-D, steady-state, open-system  U-series transport calculations
        as in Spiegelman (2000) and Elkins and Spiegelman (2021).

    Usage:
    ------
     
        us = UserCalc(df,dPdz = 0.32373,n = 2.,tol=1.e-6,phi0 = 0.008,
                  W0 =3.,model=EquilTransport,Da=None,stable=False,method='Radau')

    Parameters:
    -----------
        df     : A pandas dataframe with columns ['P','F', Kr','DU','DTh','DRa','DPa']
        dPdz   : float -- Pressure gradient, to convert pressure P to depth z
        n      : float -- Permeability exponent
        tol    : float -- Tolerance for the ODE solver
        phi0   : float -- Reference melt porosity
        W0     : float -- Upwelling velocity (cm/yr)
        model  : class -- A U-series transport model class (one of EquilTransport or DisequilTransport)
        Da     : float -- Optional Da number for disequilibrium transport model
        stable : bool
             True:  calculates concentrations for non-radiogenic nuclides with same chemical properties 
                    (i.e. sets lambda=0)
             False: calculates the full radiogenic problem
        method: string
                ODE time-stepping method to pass to solve_ivp 
                      (usually one of 'Radau', 'BDF', or 'RK45')
'''

The principal required data input is a spreadsheet containing the degree of melting $F(P)$, relative permeability $K_r(P)$, and bulk partition coefficients for the elements $D_U$, $D_{Th}$, $D_{Ra}$ and $D_{Pa}$ as functions of pressure $P$. The structure of the input data spreadsheet is the same as that described in Spiegelman (2000), and is illustrated in Table 2 below. Because the user provides $F(z)$, $K_r(z)$, and bulk solid $D_i(z)$ input information to the model directly, any considerations such as mineral modes, mineral/melt $D_i$ values, and productivity variations are external to this model and must be developed by the user separately. Once given this input spreadsheet by the user, the code routine initializes the decay constants for the isotopic decay chains and provides functions to interpolate $F(z)$ and $D_i(z)$ and calculate the porosity $\phi(z)$. Once thus initialized, the UserCalc class further provides the following methods:

'''           
    Principal Methods:
    --------
        phi                   :  returns porosity as a function of column height
        set_column_parameters :  resets principal column parameters phi0, n, W0
        solve_1D              :  1D column solution for a single Decay chain 
                                    with arbitrary D, lambda, alpha_0
        solve_all_1D          : Solves a single column model for both 238U and 235U chains.  
                                    returns a pandas dataframe
        solve_grid            : Solves multiple column models for a grid of porosities and upwelling rates
                                    returns a 3-D array of activity ratios    
'''

Of these, the principal user-facing methods are:

• UserCalc.solve_all_1D, which returns a pandas.Dataframe table that contains, at each depth, solutions for the porosity (ϕ), the log concentrations of the specified nuclides in the $^{238}U$ and $^{235}U$ decay chains in both the melt and the solid, and the U-series activity ratios.
• UserCalc.solve_grid, which solves for a grid of one-dimensional solutions for different reference porosities ($phi_0$) and solid upwelling rates ($W_0$) and returns arrays of U-series activity ratios at a specified depth (usually the top of the column), as described in Spiegelman & Elliott (1993).

3.1.1Visualization Functions¶

In addition to the principal classes for calculating U-series activity ratios in partial melts, the UserCalc package also provides functions for visualizing model inputs and outputs. The primary plotting functions include:

• UserCalc.plot_inputs(df): Visualizes the input dataframe to show $F(P)$, $K_r(P)$, and $D_i(P)$.
• UserCalc.plot_1Dcolumn(df): Visualizes the output dataframe for a single one-dimensional melting column.
• UserCalc.plot_contours(phi0,W0,act): Visualizes the output of UserCalc.solve_grid by generating contour plots of activity ratios at a specific depth as functions of the porosity ($\phi_0$) and solid upwelling rate ($W_0$).
• UserCalc.plot_mesh_Ra(Th,Ra,W0,phi0) and UserCalc.plot_mesh_Pa(Th,Pa,W0,phi0): Generates ‘mesh’ plots showing results for different $\phi_0$ and $W_0$ values on (²²⁶Ra/²³⁰Th) vs. (²³⁰Th/²³⁸U) and (²³¹Pa/²³⁵U) vs. (²³⁰Th/²³⁸U) activity diagrams.

Both the primary solver routines and visualization routines will be demonstrated in detail below.

3.1.2Miscellaneous Convenience Functions¶

Finally, the UserCalc module also provides a simple input spreadsheet generator similar to the one provided with the original UserCalc program of Spiegelman (2000). This tool is described more fully in the accompanying Jupyter notebook twolayermodel.ipynb in the Supplemental Materials, and has the interface:

• df = UserCalc.twolayermodel(P, F, Kr, D_lower, D_upper, N=100, P_lambda=1)

3.2An example demonstrating pyUserCalc functionality for a single melting column¶

The Python code cells embedded below provide an example problem that demonstrates the use and behavior of the model for a simple, two-layer upwelling mantle column, with a constant melting rate within each layer and constant $K_r=1$. This example is used to compare the outcomes from the original UserCalc equilibrium model Spiegelman, 2000 to various other implementations of the code, such as pure disequilibrium transport and scaled reactivity rates, as described above.

To run the example code and use this article as a functioning Jupyter notebook, while in a web-enabled browser the user should select each embedded code cell below by mouse-click and then simultaneously type the ‘Shift’ and ‘Enter’ keys to run the cell, after which selection will automatically advance to the following cell. The first cell below imports necessary code libraries to access the Python toolboxes and functions that will be used in the rest of the program:

3.2.1Entering initial input information and viewing input data¶

In the full Jupyter notebook code available in the Git repository and provided here as Supplementary Materials, the user can edit a notebook copy and indicate their initial input data, as has been done for the sample data set below. The name for the user’s input data file should be set in quotes (i.e., replacing the word ‘sample’ in the cell below with the appropriate filename, minus the file extension). This name will be used both to find the input file and to label any output files produced. Our sample file can likewise be downloaded and used as a formatting template for other input files (see Supplementary Materials), and is presented as a useful example below. The desired input file should be saved to a ‘data’ folder in the notebook directory prior to running the code. If desired, a similarly simple two-layer input file can also be generated using the calculator tool provided in the supplementary code.

Once the cell has been edited to contain the correct input file name, the user should run the cell using the technique described above:

The next cell below will read in the input data using the user filename specified above:

The next cell will visualize the input dataframe in Figure 1, using the utility function plot_inputs:

3.2.2Single column equilibrium transport model¶

In its default mode, UserCalc solves the one-dimensional steady-state equilibrium transport model described in Spiegelman (2000). Below we will initialize the model, solve for a single column and plot the results.

First we set the physical parameters for the upwelling column and initial conditions:

Next, we initialize the default equilibrium model:

and run the model for the input code and display the results for the final predicted melt composition in List 1:

The cell below produces Figure 2, which shows the model results with depth:

3.2.3Single column disequilibrium transport model¶

For comparison, we can repeat the calculation using the disequilibrium transport model, and compare the results to the equilibrium model. We first initialize a new model with $Da=0$, which will calculate full disequilibrium transport:

The cells below calculate solutions for this pure disequilibrium scenario, as shown in List 2:

Next we compare the results to our equilibrium calculation above:

The dashed grey curves in Figure 3 illustrate the equilibrium transport solution, which is significantly different from the disequilibrium solution. If we increase the value of $Da$, however, the disequilibrium transport solution should converge towards the equilibrium scenario. To illustrate this, below we calculate the result for $Da=1$:

The outcome of the above calculation (Figure 4, List 3) approaches the equilibrium scenario more closely, as predicted. Below is an additional comparison for $Da=10$:

For $Da=10$ (Figure 5), the activity ratios in the melt are indistinguishable from the equilibrium calculation, suggesting that a Damköhler number of 10 is sufficiently high for a melting system to approach chemical equilibrium, and illustrating that the equilibrium model of Spiegelman & Elliott (1993) and Spiegelman (2000) is the limiting case for the more general disequilibrium model presented here. For this problem, equilibrium transport always provides an upper bound on activity ratios.

3.2.4Stable element concentrations¶

For a stable element, i.e., $\lambda_i = 0$, Spiegelman & Elliott (1993) showed that the equilibrium melting model reduces identically to simple batch melting Shaw, 1970, while the disequilibrium model with $Da=0$ is equivalent to true fractional melting. This presents a useful test of the calculator that verifies the program is correctly calculating stable concentrations. To simulate stable element concentrations for U, Th, Ra, and Pa during equilibrium melting, we can use the same input file example as above and simply test the scenario where $\lambda_i$ values are equal to zero.

First, we impose a “stable” condition that changes all decay constants $\lambda_i = 0$:

For comparison with the results in List 4, we can use the batch melting equation Shaw, 1970 to calculate the concentrations of U, Th, Ra, and Pa using the input values in Table 2 for $F(z)$ and $D_i$, where:

and determine radionuclide activities for the batch melt using the definition of the activity $a$ for a nuclide $i$:

and the initial nuclide activities $a_i^0$, such that:

As the activity ratios in List 5 illustrate, the outcomes of this simple batch melting equation are identical to those produced by the model for equilibrium transport and $\lambda=0$.

Similarly, we can also determine pure disequilibrium melting using the disequilibrium transport model with $\lambda_i=0$. A simple fractional melting problem is easiest to test using constant melt productivity and partitioning behavior, so here we test a simplified, one-layer scenario with constant $dF/dz$ and $D_i$ values:

We note that numerical ordinary differential equation (ODE) solvers may not successfully solve for pure fractional melting with $Da=0$ and stable elements, because the resulting extreme changes in solid concentrations for highly incompatible elements are difficult to resolve using numerical methods. Stable solutions can nonetheless be obtained for very small values of $Da$ that approach $Da=0$, and such solutions still provide a useful test of the disequilibrium transport model. Here we use $Da=10^{-10}$; for such low $Da$ values, the liquid closely approaches the composition of an accumulated fractional melt, and although the liquid and solid outcomes are slightly different from pure fractional melting, the solid is still essentially depleted of all incompatible nuclides.

Similar to our approach for equilibrium and batch melting, we can compare the results of disequilibrium transport for stable elements with pure fractional melting for constant partition coefficients using the definition of aggregated fractional melt concentrations (Figure 7):

or in log units:

3.2.5Considering lithospheric transport scenarios¶

In mantle decompression melting scenarios, melting is expected to cease in the shallow, colder part of the regime where a lithospheric layer is present. The effects of cessation of melting prior to reaching the surface can be envisioned as affecting magma compositions in a number of ways, some of which could be calculated using the models presented here by setting $dF=0$.

There are, however, several limitations when using our transport models to simulate lithospheric melt transport in this way, as the model equations are written to track steady-state decompression and melting. The first limitation is thus the underlying assumption that the solid is migrating and experiences progressive melt depletion in the model, while the solid lithosphere should in fact behave as a rigid matrix that does not experience upwelling. For the disequilibrium transport model with $Da=0$, no chemical reequilibration occurs while $dF=0$, so the lack of solid migration after the cessation of melting does not pose a problem; instead, in the pure disequilibrium transport case, imposing $dF=0$ simply allows for radioactive decay and ingrowth during transport through the lithospheric layer.

The equilibrium transport model, on the other hand, permits full equilibration even if $dF=0$, but the liquid composition does not directly depend on the solid concentration, $c_i^s(z)$, so ongoing chemical reequilibration between the liquid and a modified lithospheric solid could be simulated by modifying the bulk solid/liquid partition coefficients $D_i$. However, the underlying model assumes that the liquid with mass proportion $F_{max}$ reequilibrates with the solid matrix in a steady-state transport regime, at the maximum reference porosity, which may not accurately simulate the transport regime through the fixed lithosphere with no melting. Because it does not directly consider mineral abundances or compositions, the model also does not account for complexities such as low temperature melt-rock reaction or mineral growth.

The case of the scaled disequilibrium transport model with $Da>0$ is the most complex, since the model directly calculates reequilibration of the liquid with a progressively melting solid layer, and thus may not accurately simulate transport through the fixed solid lithosphere. We advise that if the model is used in this way, the results must be interpreted with additional caution.

Finally, calculating a given transport model through the upwelling asthenosphere and into a fixed overlying lithospheric layer neglects an additional, significant limitation: namely that melt-rock interactions, and thus the magma transport style, may be different in the lithosphere than in the melting asthenosphere. As noted above, this could also include low-temperature reactions and the growth of new mineral phases. While it is not possible to change transport models during a single 1D run in the current implementation, one alternative approach is to change the relative permeability, $K_r$, in the lithosphere, in addition to modifying the bulk partition coefficients to reflect lithospheric values. It may also be possible to run a separate, second-stage lithospheric calculation with modified input parameters and revised liquid porosity constraints, but this option is not currently implemented and would require an expansion of the current model.

Despite these caveats, there are some limited scenarios where users may wish to simulate equilibrium or disequilibrium magma transport through a capping layer with constant $dF=0$, constant $\phi = \phi_0$, and revised $D_i$ values for a modified layer mineralogy. The cells below provide options for modifying the existing input data table to impose such a layer. The first cell identifies a final melting pressure $P_{Lithos}$, which is defined by the user in kbar. This value can be set to 0.0 if no lithospheric cap is desired; in the example below, it has been set at 5.0 kbar. There are two overall options for how this final melting pressure could be used to modify the input data. One option (implemented in the Supplementary Materials but not tested here) simply deletes all lines in the input dataframe for depths shallower than $P_{Lithos}$. This is a straightforward option for a one-dimensional column scenario, where melting simply stops at the base of the lithosphere and the composition of the melt product is observed in that position. This is an effective way to limit further chemical interactions after melting has ceased; it fails to account for additional radioactive decay during lithospheric melt transport, but subsequent isotopic decay over a fixed transport time interval could then be calculated using the radioactive decay equations for U-series nuclides.