pyUserCalc: A revised Jupyter notebook calculator for uranium-series disequilibria in basalts

Lynne J. Elkins; Marc Spiegelman

doi:https://doi.org/10.1029/2020EA001619

Original

AGU Notebooks Now

Download PDF Download TEX Learn More

Back to Article

pyUserCalc: A revised Jupyter notebook calculator for uranium-series disequilibria in basalts

Download Notebook

Contents

0.1Lynne J. Elkins¹, Marc Spiegelman²¶

¹ University of Nebraska-Lincoln, Lincoln, NE, USA, lelkins@unl.edu

² Lamont-Doherty Earth Observatory of Columbia University, Palisades, NY, USA

1Summary¶

This Jupyter notebook calculates U-series disequilibria in partial melts using the equations after Elkins and Spiegelman (2021). Users should cite that reference when publishing outcomes of this calculator tool using the following citation:

Elkins, L.J., & Spiegelman, M. (2021). pyUserCalc: A revised Jupyter notebook calculator for uranium-series disequilibria in basalts. Earth and Space Science, 8, e2020EA001619. Elkins & Spiegelman (2021).

The notebook determines (²³⁰Th/²³⁸U), (²²⁶Ra/²³⁰Th), and (²³¹Pa/²³⁵U) in partial melts for pure equilibrium or disequilibrium porous flow transport in a 1D decompressing mantle melting column. The disequilibrium transport model considers reactivity rates relative to the solid decompression rate using a Damköhler number. The equations for these models are briefly summarized below, but the user is referred to the full Elkins and Spiegelman (2021) manuscript for more detailed definitions and derivations.

1.1Equilibrium porous flow¶

As described in the Elkins and Spiegelman (2021) main text, this equilibrium transport model recreates the solutions of Spiegelman and Elliott (1993) for U-series disequilibrium calculations during partial melting. As a brief summary, the concentration expression from Spiegelman (2000) for the equilibrium scenario (formula 6 in that reference) is:

\frac{dc_i^f}{dz} = \frac{-c_i^f(z)}{F(z) + (1 - F(z)) D_i(z)} \frac{d}{dz} [F(z) + (1 - F(z)) D_i(z)] + \frac{\lambda_{i-1}\overline{\rho D_{i-1}} c_{i-1}^f(z) - \lambda_{i}\overline{\rho D_{i}} c_{i}^f(z)}{\rho_s W_0 [F(z) + (1 - F(z)) D_i(z)]}

(1)

The model solves for the log of the concentration, $U_i^f$ :

U_i^f = \ln\left(\frac{c_i^f}{c_{i,0}^f}\right)

(2)

where

\frac{dU_i^f}{dz} = \frac{1}{c_i^f(z)} \frac{dc_i^f}{dz}

(3)

and thus

\frac{dU_i^f}{dz} = \frac{-1}{F(z) + (1 - F(z)) D_i(z)} \frac{d}{dz} [F(z) + (1-F(z))D_i(z)] + \frac{\lambda_i}{w_{eff}^i} [R_i^{i-1} \exp[U_{i-1}^f(z)-U_i^f(z)]-1]

(4)

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

U_i(z) = U^{stable}_i(z) + U^{rad}_i(z)

(5)

1.2Disequilibrium porous flow¶

The disequilibrium transport model solves the set of equations for disequilibrium transport after Spiegelman and Elliott (1993), with an additional reactivity term for a scaled reactivity rate $\Re$ . The log concentrations are solved similar to the expressions for equilibrium transport above, but with modified mass conservation that tracks pure disequilibrium flow. These modified expressions, expressed in terms of the scaled column height ζ, are:

\begin{align} \frac{dU_i^f}{d\zeta} &= \frac{1}{F(\zeta)}\left[\frac{dF}{d\zeta}(\frac{D_i^0}{D_i(\zeta)}\exp(U_i^s(\zeta) - U_i^f(\zeta)) - 1)\right] - \frac{\Re h}{\rho_s W_0 F(\zeta)}\left[1 - \frac{D_i^0}{D_i(\zeta)}\exp(U_i^s(\zeta) - U_i^f(\zeta))\right] + \frac{\rho_f \phi}{\rho_s F(\zeta)}\lambda'_i \left[\frac{D_i^0 \alpha_{i-1}}{D_{i-1}^0 \alpha_i} exp[U_{i-1}^f(\zeta) - U_i^f(\zeta)]-1\right] \\ \frac{dU_i^s}{d\zeta} &= \frac{1}{(1 - F(\zeta))}\left[\frac{dF}{d\zeta} (1 - \frac{1}{D_i(\zeta)})\right]+ \frac{\Re h}{\rho_s W_0 D_i(\zeta)(1 - F(\zeta))}\left[\frac{D_i(\zeta)}{D_i^0} \exp(U_i^f(\zeta) - U_i^s(\zeta)) - 1 \right] + \frac{1-\phi}{1-F(\zeta)}\lambda'_i \left[\frac{\alpha_{i-1}}{\alpha_i} exp[U_{i-1}^s(\zeta) - U_i^s(\zeta)]-1\right]\\ \end{align}

(6)

and the reactivity rate scales with the Damköhler number as a function of column height, upwelling rate, and solid density:

Da = \frac{\Re h}{\rho_s W_0}

(7)

2The calculator tool¶

The python code cells embedded below implement the models described above and in the main text of Elkins and Spiegelman (2021). A copy of this .ipynb file should be saved to a user directory, along with the “UserCalc.py” driver file and one folder named “data”. The input file should be a comma-delimited text file written in the same format as the “sample” file provided, and should be saved to the “data” directory. This Jupyter notebook can then be run either native on the user’s computer using a python software package, such as Jupyter Notebook in the Anaconda package, or from a cloud account in a shared online JupyterLab or JupyterHub environment like the ENKI platform.

Once this notebook has been opened, select each embedded code cell 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. Cells may be edited prior to running to specify the model calculations desired. Note that when modifying and running the model repeatedly, it may be necessary to restart the kernel for each fresh start.

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:

# Select this cell with by mouseclick, and run the code by simultaneously typing the 'Shift' + 'Enter' keys.
# If the browser is able to run the Jupyter notebook, a number [1] will appear to the left of the cell.

import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
%matplotlib inline
import os
import shutil

# Import UserCalc:
import UserCalc

2.1Enter initial input information and view input data¶

Edit the cell below with the name of the input data file set in quotes, as has been done for the “sample” file below, and then run the cell. The two subsequent cells will display the data in table and figure formats to check that the input data are correct. Note that if an input file has previously been run, this code will overwrite the old output directory with a fresh one, so any files from previous runs that need to be saved should be copied or relocated first.

runname='sample'
output = '{}'.format(runname)
if os.path.exists(output):
    shutil.rmtree(output)
os.mkdir(output)

input_file = 'data/{}.csv'.format(runname)
df = pd.read_csv(input_file,skiprows=1,dtype=float)
df

UserCalc.plot_inputs(df)

2.2Modifying inputs for lithospheric transport¶

For scenarios where a modified transport scheme through a lithospheric layer is desired, the cells below will truncate the 1D melting regime at a final melting pressure, $P_{Lithos}$ . Edit the cell to impose a $P_{Lithos}$ value; set this value to 0.0 if no lithospheric cap is desired. It is also possible to define the pressure of the base of the crust ( $P_{crust}$ ) in case an additional compositional layer is desired.

There are two options for how the final melting pressure is imposed. The first option deletes all rows from the input table for depths shallower than $P_{Lithos}$ , completely truncating the melting calculation.

The second (default) option changes the degree of melting increments ( $dF$ ) to a value of 0 for all depths shallower than $P_{Lithos}$ , but continues to track transport to the surface. For full disequilibrium transport with $Da=0$ , this scenario will simulate ongoing radioactive decay and ingrowth during lithospheric transport after melting has ceased. For equilibrium or partial disequilibrium transport models, the partial melt and the solid will continue to interact and reequilibrate, but no new melt will be produced; this may or may not be a sensible choice depending on the problems being explored (see main text for further discussion).

Both of the described options generate alternative input dataframes, without overwriting the original data input file. Any of the options can then be selected to run using transport models below.

Plithos = 5.0
Pcrust = 2.2

df_nolith = df[df.P > Plithos]

df_lith = df.copy()
Pfinal = df_lith.iloc[(df_lith['P']-Plithos).abs().idxmin()]
F_max = Pfinal[1].tolist()
df_lith.loc[(df_lith['P'] < Plithos),['F']] = F_max

The two cells below offer additional options for modifying melt equilibration during lithospheric transport by imposing new, fixed solid-melt partition coefficients in the lithospheric layer (and, optionally, the crust). This is only relevant if reactive lithospheric transport will be calculated below; otherwise, these cell will have no effect on model results. Select the desired option in the first cell below.

The option called “new” imposes new partition coefficient values at depths shallower than $P_{Lithos}$ , and, if desired, again at depths shallower than $P_{crust}$ . (Note, however, that very abrupt and large changes in $D_i$ may be difficult for the model to accurately extrapolate and can produce errors, so for a calculation with large partitioning changes at the asthenosphere-lithosphere interface, it can work better to generate a separate input file with a more gradual adjustment in $D_i$ values at layer boundaries.) The alternative option, called “old,” sets all lithosphere partition coefficients equal to the partition coefficient values in the input file at $P_{Lithos}$ . After editing the partition coefficients to use desired values, run both cells to implement the changes.

# Indicate how you wish to impose lithospheric partition coefficients by indicating "new" or "old":
define_lith_Ds = 'old'

if define_lith_Ds == 'new':
    
    # Edit here to impose new melt-rock partition coefficients in the lithosphere and/or crust
    D_U_lith = 1.0405e-3
    D_Th_lith = 1.3795e-3
    D_Ra_lith = 0.00001
    D_Pa_lith = 0.00001

    D_U_crust = 1.0405e-3
    D_Th_crust = 1.3795e-3
    D_Ra_crust = 0.00001
    D_Pa_crust = 0.00001

elif define_lith_Ds == 'old':
    
    # These rows will define lithosphere partition coefficients based on the final melting pressure
    D_U_lith = Pfinal[3].tolist()
    D_Th_lith = Pfinal[4].tolist()
    D_Ra_lith = Pfinal[5].tolist()
    D_Pa_lith = Pfinal[6].tolist()

    D_U_crust = Pfinal[3].tolist()
    D_Th_crust = Pfinal[4].tolist()
    D_Ra_crust = Pfinal[5].tolist()
    D_Pa_crust = Pfinal[6].tolist()
    
# Implement the changes:
df_lith.loc[(df_lith['P'] < Plithos),['DU']] = D_U_lith
df_lith.loc[(df_lith['P'] < Plithos),['DTh']] = D_Th_lith
df_lith.loc[(df_lith['P'] < Plithos),['DRa']] = D_Ra_lith
df_lith.loc[(df_lith['P'] < Plithos),['DPa']] = D_Pa_lith

df_lith.loc[(df_lith['P'] < Pcrust),['DU']] = D_U_crust
df_lith.loc[(df_lith['P'] < Pcrust),['DTh']] = D_Th_crust
df_lith.loc[(df_lith['P'] < Pcrust),['DRa']] = D_Ra_crust
df_lith.loc[(df_lith['P'] < Pcrust),['DPa']] = D_Pa_crust

2.3Set up 1D model runs¶

It is necessary to tell the model which data input scenario will be run below. The options available are to calculate transport models using the original data input file (“original”), for input values modified to truncate the run at the base of the lithosphere (“no_lith”), or for input values modified to consider lithospheric transport (“lith”).

# Specify the data input scenario that should be run below. Options are 'original', 'no_lith', and 'lith'.
data_input_option = 'no_lith'

# Run this cell to implement the data input scenario for the rest of the model run.
if data_input_option == 'original':
    df_input = df

elif data_input_option == 'lith':
    df_input = df_lith
    
elif data_input_option == 'no_lith':
    df_input = df_nolith

# Generate a figure for the input scenario, and save it to the output folder:
UserCalc.plot_inputs(df_input)
plt.savefig("{}/{}_inputs.pdf".format(runname,runname),transparent=True,dpi=300,bbox_inches='tight')

Before running the model, check and edit the cell below to change the default physical input parameters to those desired. Run the cell to save the parameter values.

# Maximum melt porosity:
phi0 = 0.008

# Solid upwelling rate in cm/yr. (to be converted to km/yr. in the driver function):
W0 = 3.

# Permeability exponent:
n = 2.

# Solid and liquid densities in kg/m3:
rho_s = 3300.
rho_f = 2800.

# Initial activity values (default is 1.0):
alpha0_238U = 1.
alpha0_235U = 1.
alpha0_230Th = 1.
alpha0_226Ra = 1.
alpha0_231Pa = 1.
alpha0_all = np.array([alpha0_238U, alpha0_230Th, alpha0_226Ra, alpha0_235U, alpha0_231Pa])

The next cell specifies which of the three possible model scenarios will be run. The Damköhler number for the scaled reactivity rate transport model can also be specified.

# Indicate below which transport models should be run by typing "yes" or "no".
equilibrium = 'yes'
disequilibrium = 'yes'
scaled = 'yes'

# If using the scaled reactivity model, provide the desired Da number value here:
Da_number = 0.1

The next cell initializes and runs 1D calculations for the specified transport models and input data.

Note that for input data files that generate very rapid concentration changes, the disequilibrium model may fail because the ODE solver cannot solve the problem. To try to prevent this, we recommend first testing $Da=10^{-10}$ instead of $Da=0$ by default. This will often avoid the error and should produce a comparable result, as explored in Elkins and Spiegelman (2021). If errors persist, revising the input data to create more gradual changes in $F(P)$ or $D_i(P)$ may be necessary.

The default model run scenario considers radioactive decay during transport. To instead determine instantaneous batch or fractional melting without radioactive decay, change the “stable” value to “True” in the cells below.

# Run this cell to calculate results for the transport models specified above.
if equilibrium == 'yes':
    us_eq = UserCalc.UserCalc(df_input,stable=False)
    df_out_eq = us_eq.solve_all_1D(phi0,n,W0,alpha0_all)

if disequilibrium == 'yes':
    us_diseq = UserCalc.UserCalc(df_input,model=UserCalc.DisequilTransport,Da=1.0e-10,stable=False)
    df_out_diseq = us_diseq.solve_all_1D(phi0,n,W0,alpha0_all)

if scaled == 'yes':
    us_diseqda = UserCalc.UserCalc(df_input,model=UserCalc.DisequilTransport,Da=Da_number,stable=False)
    df_out_diseqda = us_diseqda.solve_all_1D(phi0,n,W0,alpha0_all)

The cells below will generate figures showing output results for all transport models at the top of the melting column. Run these cells to verify a logical outcome:

# Equilibrium model result:
if equilibrium == 'yes':
    df_out_eq.tail(n=1)

# Pure disequilibrium model result:
if disequilibrium == 'yes':
    df_out_diseq.tail(n=1)

# Scaled disequilibrium model result:
if scaled == 'yes':
    df_out_diseqda.tail(n=1)

Run the next set of cells to view the model output results as a series of depth figures, and to save them to the output folder.

# Equilibrium model results figure:
if equilibrium == 'yes':
    fig = UserCalc.plot_1Dcolumn(df_out_eq)
    plt.savefig("{}/{}_1D_depths_eq.pdf".format(runname,runname),transparent=True,dpi=300,bbox_inches='tight')

# Pure disequilibrium model results figure:
if disequilibrium == 'yes':
    fig = UserCalc.plot_1Dcolumn(df_out_diseq)
    plt.savefig("{}/{}_1D_depths_diseq.pdf".format(runname,runname),transparent=True,dpi=300,bbox_inches='tight')

# Scaled disequilibrium model results figure:
if scaled == 'yes':
    fig = UserCalc.plot_1Dcolumn(df_out_diseqda)
    plt.savefig("{}/{}_1D_depths_diseq_Da={}.pdf".format(runname,runname,us_diseqda.Da),transparent=True,dpi=300,bbox_inches='tight')

Finally, run the appropriate cells below to save the results as .csv files in the output folder:

# Save equilibrium results:
if equilibrium == 'yes':
    df_out_eq.to_csv("{}/{}_1D_solution_eq.csv".format(runname,runname))

# Save pure disequilibrium results:
if disequilibrium == 'yes':
    df_out_diseq.to_csv("{}/{}_1D_solution_diseq.csv".format(runname,runname))

# Save scaled disequilibrium results:
if scaled == 'yes':
    df_out_diseqda.to_csv("{}/{}_1D_solution_diseq_Da={}.csv".format(runname,runname,us_diseqda.Da))

3Batch operations¶

If batch operations are being calculated, edit and run the cells below to calculate outcomes over a range of maximum residual melt porosity (ϕ) and the solid mantle upwelling rate ( $W_0$ ) values. Edit the first cell to select whether to define the specific ϕ and $W_0$ values as evenly spaced log grid intervals or using manually specified values (default). All upwelling rates are entered in units of cm/yr. Other input parameters will match those set above, and the results are automatically saved to the output folder. This set of batch operations automatically calculates outputs for 1) the equilibrium, 2) the pure disequilibrium, and 3) the scaled disequilibrium model with a Damköhler number (using the value chosen for the single column calculations above).

# Edit the rows below to specify reference porosity and solid upwelling rates for batch runs.

# phi0 = np.logspace(-3,-2,11)
phi0 = np.array([0.001, 0.002, 0.005, 0.01, 0.02])

# W0 = np.logspace(-1,1,11)
W0 = np.array([0.5, 1., 2., 5., 10., 20., 50.])

# Enable time counter to visually track batch calculation progress:
import time
tic = time.perf_counter()
toc = time.perf_counter()

# Calculate the gridded values for the equilibrium model:
if equilibrium == 'yes':
    act_eq = us_eq.solve_grid(phi0,n,W0,us_eq.D_238,us_eq.lambdas_238,us_eq.alphas_238)
    Th_eq = act_eq[0]
    Ra_eq = act_eq[1]
    df_eq = pd.DataFrame(Th_eq)
    df_eq.to_csv("{}/{}_batch_Th_eq.csv".format(runname,runname))
    df_eq = pd.DataFrame(Ra_eq)
    df_eq.to_csv("{}/{}_batch_Ra_eq.csv".format(runname,runname))
    
    act_eq_235 = us_eq.solve_grid(phi0,n,W0,us_eq.D_235,us_eq.lambdas_235,us_eq.alphas_235)
    Pa_eq = act_eq_235[0]
    df_eq = pd.DataFrame(Pa_eq)
    df_eq.to_csv("{}/{}_batch_Pa_eq.csv".format(runname,runname))

# Calculate the gridded values for the pure disequilibrium model:
if disequilibrium == 'yes':
    act_diseq = us_diseq.solve_grid(phi0,n,W0,us_diseq.D_238,us_diseq.lambdas_238,us_diseq.alphas_238)
    Th_diseq = act_diseq[0]
    Ra_diseq = act_diseq[1]
    df_diseq = pd.DataFrame(Th_diseq)
    df_diseq.to_csv("{}/{}_batch_Th_diseq.csv".format(runname,runname))
    df_diseq = pd.DataFrame(Ra_diseq)
    df_diseq.to_csv("{}/{}_batch_Ra_diseq.csv".format(runname,runname))
    
    act_diseq_235 = us_diseq.solve_grid(phi0,n,W0,us_diseq.D_235,us_diseq.lambdas_235,us_diseq.alphas_235)
    Pa_diseq = act_diseq_235[0]
    df_diseq = pd.DataFrame(Pa_diseq)
    df_diseq.to_csv("{}/{}_batch_Pa_diseq.csv".format(runname,runname))

# Calculate the gridded values for the scaled disequilibrium model:
if scaled == 'yes':
    act_diseqda = us_diseqda.solve_grid(phi0,n,W0,us_diseqda.D_238,us_diseqda.lambdas_238,us_diseqda.alphas_238)
    Th_diseqda = act_diseqda[0]
    Ra_diseqda = act_diseqda[1]
    df_diseqda = pd.DataFrame(Th_diseqda)
    df_diseqda.to_csv("{}/{}_batch_Th_diseq_Da={}.csv".format(runname,runname,us_diseqda.Da))
    df_diseqda = pd.DataFrame(Ra_diseqda)
    df_diseqda.to_csv("{}/{}_batch_Ra_diseq_Da={}.csv".format(runname,runname,us_diseqda.Da))
    
    act_diseqda_235 = us_diseqda.solve_grid(phi0,n,W0,us_diseqda.D_235,us_diseqda.lambdas_235,us_diseqda.alphas_235)
    Pa_diseqda = act_diseqda_235[0]
    df_diseqda = pd.DataFrame(Pa_diseqda)
    df_diseqda.to_csv("{}/{}_batch_Pa_diseq_Da={}.csv".format(runname,runname,us_diseqda.Da))

The cells below create and export contour figures as described in the Elkins and Spiegelman (2021) main text, first for equilibrium and then for disequilibrium and scaled reactivity transport.

# Equilibrium contour plots for U-Th-Ra chain:
if equilibrium == 'yes':
    UserCalc.plot_contours(phi0,W0,act_eq, figsize=(12,12))
    plt.suptitle('Equilibrium transport',y=0.92,fontsize=22)
    plt.savefig("{}/{}_contour_Th_Ra_eq.pdf".format(runname,runname),transparent=True,dpi=300,bbox_inches='tight')

# Equilibrium contour plots for U-Pa chain:
if equilibrium == 'yes':
    UserCalc.plot_contours(phi0,W0,act_eq_235)
    plt.suptitle('Equilibrium transport',fontsize=22)
    plt.savefig("{}/{}_contour_Pa_eq.pdf".format(runname,runname),transparent=True,dpi=300,bbox_inches='tight')

# Pure disequilibrium contour plots for U-Th-Ra chain:
if disequilibrium == 'yes':
    UserCalc.plot_contours(phi0,W0,act_diseq, figsize=(12,12))
    plt.suptitle('Disequilibrium transport, Da = {}'.format(us_diseq.Da),y=0.92,fontsize=22)
    plt.savefig("{}/{}_contour_Th_Ra_diseq.pdf".format(runname,runname),transparent=True,dpi=300,bbox_inches='tight')

# Pure disequilibrium contour plots for U-Pa chain:
if disequilibrium == 'yes':
    UserCalc.plot_contours(phi0,W0,act_diseq_235)
    plt.suptitle('Disequilibrium transport, Da = {}'.format(us_diseq.Da),fontsize=22)
    plt.savefig("{}/{}_contour_Pa_diseq.pdf".format(runname,runname),transparent=True,dpi=300,bbox_inches='tight')

# Scaled disequilibrium contour plots for U-Th-Ra chain:
if scaled == 'yes':
    UserCalc.plot_contours(phi0,W0,act_diseqda, figsize=(12,12))
    plt.suptitle('Disequilibrium transport, Da = {}'.format(us_diseqda.Da),y=0.92,fontsize=22)
    plt.savefig("{}/{}_contour_Th_Ra_diseq_Da={}.pdf".format(runname,runname,us_diseqda.Da),transparent=True,dpi=300,bbox_inches='tight')

# Scaled disequilibrium contour plots for U-Pa chain:
if scaled == 'yes':
    UserCalc.plot_contours(phi0,W0,act_diseqda_235)
    plt.suptitle('Disequilibrium transport, Da = {}'.format(us_diseqda.Da),fontsize=22)
    plt.savefig("{}/{}_contour_Pa_diseq_Da={}.pdf".format(runname,runname,us_diseqda.Da),transparent=True,dpi=300,bbox_inches='tight')

Finally, the remaining cells display and export grid or mesh figures as described in the Elkins and Spiegelman (2021) main text, again showing equilibrium results followed by disequilibrium and scaled reactivity results.

# Equilibrium grid plots for U-Th-Ra chain:
if equilibrium == 'yes':
    UserCalc.plot_mesh_Ra(Th_eq,Ra_eq,W0,phi0)
    plt.title('Equilibrium transport',y=1.1)
    plt.savefig("{}/{}_grid_Ra_Th_eq.pdf".format(runname,runname),transparent=True,dpi=300,bbox_inches='tight')

# Equilibrium grid plots for U-Pa chain:
if equilibrium == 'yes':
    UserCalc.plot_mesh_Pa(Th_eq,Pa_eq,W0,phi0)
    plt.title('Equilibrium transport',y=1.1)
    plt.savefig("{}/{}_grid_Pa_Th_eq.pdf".format(runname,runname),transparent=True,dpi=300,bbox_inches='tight')

# Pure disequilibrium grid plots for U-Th-Ra chain:
if disequilibrium == 'yes':
    UserCalc.plot_mesh_Ra(Th_diseq,Ra_diseq,W0,phi0)
    plt.title('Disequilibrium transport, Da = {}'.format(us_diseq.Da),y=1.1)
    plt.savefig("{}/{}_grid_Ra_Th_diseq.pdf".format(runname,runname),transparent=True,dpi=300,bbox_inches='tight')

# Pure disequilibrium grid plots for U-Pa chain:
if disequilibrium == 'yes':
    UserCalc.plot_mesh_Pa(Th_diseq,Pa_diseq,W0,phi0)
    plt.title('Disequilibrium transport, Da = {}'.format(us_diseq.Da),y=1.1)
    plt.savefig("{}/{}_grid_Pa_Th_diseq.pdf".format(runname,runname),transparent=True,dpi=300,bbox_inches='tight')

# Scaled disequilibrium grid plots for U-Th-Ra chain:
if scaled == 'yes':
    UserCalc.plot_mesh_Ra(Th_diseqda,Ra_diseqda,W0,phi0)
    plt.title('Disequilibrium transport, Da = {}'.format(us_diseqda.Da),y=1.1)
    plt.savefig("{}/{}_grid_Ra_Th_diseq_Da={}.pdf".format(runname,runname,us_diseqda.Da),transparent=True,dpi=300,bbox_inches='tight')

# Scaled disequilibrium grid plots for U-Pa chain:
if scaled == 'yes':
    UserCalc.plot_mesh_Pa(Th_diseqda,Pa_diseqda,W0,phi0)
    plt.title('Disequilibrium transport, Da = {}'.format(us_diseqda.Da),y=1.1)
    plt.savefig("{}/{}_grid_Pa_Th_diseq_Da={}.pdf".format(runname,runname,us_diseqda.Da),transparent=True,dpi=300,bbox_inches='tight')

References¶

Elkins, L. J., & Spiegelman, M. (2021). pyUserCalc: A Revised Jupyter Notebook Calculator for Uranium‐Series Disequilibria in Basalts. Earth and Space Science, 8(12). 10.1029/2020ea001619

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

Example: Two Layer Model