ARTEMIS reactive transport/enhanced weathering model in soils

ARTEMIS Version 1.0 A Reactive Transport Enhanced-weathering Model In Soils User Guide Lyla L. Taylor 22 November 2025 Leverhulme Centre for Climate Change Mitigation* Director: David J. Beerling School of Biosciences University of Sheffield *Funded by the Leverhulme Trust (Leverhulme Research Centre grant number RC-2015-029 to David J. Beerling) Contents 1 Introduction 1 1.1 Installation on Linux systems ............................... 1 1.2 Development history .................................... 2 1.3 Two geochemical models .................................. 4 1.4 Disclaimer .......................................... 4 2 Input files 5 2.1 Rock-MIP parameter files ................................. 5 2.2 Rock-MIP initialChemistry file ............................... 5 2.3 Phase data file ....................................... 7 2.4 Initial soil layer data file .................................. 9 2.5 Feedstock PSD file ..................................... 13 2.6 Timetable of events ..................................... 14 2.7 Timeseries input files .................................... 14 2.8 Depthwise time series input files .............................. 15 2.9 Observations files ...................................... 16 2.10 Files in subdirectory “phreeqcdbincludefiles” ....................... 17 3 Model parameters and settings 19 3.1 The optional defaultsettings.m script ........................... 19 3.1.1 Surface area method ................................ 19 3.1.2 Evapotranspiration method ............................ 23 3.1.3 EXCHANGE method ................................ 23 3.1.4 MIX direction .................................... 24 3.1.5 Soil gas phase options ............................... 24 4 Energy Farm example 27 4.1 Input files and scripts .................................... 27 i ii 4.2 How to do an example Energy Farm run ......................... 28 4.3 Loading and displaying the Energy Farm results ..................... 32 4.4 Energy Farm calibration runs ............................... 32 4.4.1 Phosphorus dynamics during calibration runs .................. 34 5 Navigating the software 37 5.1 The preprocessor: setupdatastruc.m ............................ 37 5.1.1 Files to be written prior to the first model run .................. 38 5.1.2 checkphasetype.m .................................. 39 5.1.3 implicitcalc.m .................................... 40 5.2 PHREEQC wrapper: getphreeqcrun.m .......................... 41 5.3 Data visualisation utilities ................................. 42 Bibliography 45 List of Figures 1.1 Overall steps for running this model ............................ 3 4.1 Energy Farm calibration runs ............................... 34 4.2 Topsoil P pools during calibrating runs .......................... 35 4.3 P uptake during calibration runs ............................. 36 5.1 Preprocessor: functions called ............................... 39 5.2 PHREEQC wrapper: functions called ........................... 42 5.3 PHREEQC input file for the first timestep ........................ 43 5.4 PHREEQC input file for subsequent timesteps ...................... 44 iii List of Tables 1.1 Software required to run ARTEMIS ............................ 3 2.1 Input .csv files for the preprocessor ............................ 6 2.2 Time variables allowed for timeseries and timetable data ................ 6 2.3 Rock-MIP parameter files: example data ......................... 7 2.4 Input data for all solid inorganic phases to be modelled ................. 8 2.5 Allowed designations for kinetic reaction orders in phase table ............. 10 2.6 Soil layer file column headings ............................... 10 2.7 Soil layer input parameters ................................. 11 2.8 Timetable structure arrays ................................. 15 2.9 Time series files of forcing data .............................. 16 2.10 The phreeqcdbincludefiles directory ............................ 18 3.1 Feedstock applications and extended runs ......................... 20 3.2 Parameters controlling reactive surface area (RSA). ................... 20 3.3 Preprocessor parameters related to sorption ....................... 21 3.4 Additional parameters used in this model ......................... 22 3.5 Options for setting the SAmethod keyword ........................ 23 3.6 Options for setting the ETmethod keyword ........................ 24 3.7 Options for the soilco2opt keyword ............................ 25 3.8 Other keywords for soil gases ............................... 26 4.1 Energy Farm calibration scripts and settings ....................... 27 4.2 Five-year Energy Farm runs ................................ 28 4.3 MATLAB scripts for Energy Farm runs through 2070 .................. 28 4.4 Example output files .................................... 30 4.5 Data files for Energy Farm runs .............................. 31 4.6 Mineralogy of the blueridge metabasalt .......................... 33 iv Chapter 1 Introduction ARTEMIS is a deterministic, process-based reactive transport model (RTM) for calculating carbon dioxide removal (CDR) following enhanced rock weathering (ERW) in soils. It aims to account for as many soil processes as possible, as realistically as possible, allowing users to investigate process effects by turning individual processes on or off. These aims are, of course, only partially fulfilled, and there are pitfalls associated with this level of complexity. For example, users should be wary of trying to calibrate every process at a particular field site, given the considerable heterogeneity of managed soils. Nevertheless, this code may prove useful to the CDR, ERW and geochemical communities. This user guide provides instructions for installing ARTEMIS (Section 1.1). It also explains the structure and parameterizations for ARTEMIS Version 1.0, along with some examples showing how to run the model. It assumes some experience with both MATLAB and with the underlying geochemical speciation code (PHREEQC, Parkhurst and Appelo, 2013). A preprocessor and PHREEQC wrapper comprise the bulk of the MATLAB code described here (Figure 1.1). The steps shown in Figure 1.1 will be described in subsequent chapters. Several publicly-available MATLAB software packages are also required to run ARTEMIS (Table 1.1). The code runs under the Linux operating system and has been successfully tested on two highperformance computers (HPCs) with different recent versions of MATLAB (Section 1.2). There is no guarantee whatsoever that ARTEMIS will run on any other system; I am not a Windows or Mac user, so have never tried to install or run it on those operating systems. 1.1 Installation on Linux systems Note: these instructions assume basic knowledge of both Linux and MATLAB, and that MATLAB is already installed. To install ARTEMIS, download the compressed tarball and unpack it with tar xzvf: tar xzvf ARTEMIS_GMD_code_v1.0_Sep-2025.tar.gz A directory “soilmodel” will be created containing all the m-files (MATLAB code). There are several subdirectories: 1 2 phreeqcdbincludefiles contains the ARTEMIS PHREEQC databases along with some additional .csv files (Section 2.10). examples contains the data and scripts for the Energy Farm example runs (see Chapter 4). There is no Makefile, but there are two environment variables for the home directory and a scratch directory in this code which should be defined: $HOME and $SCRATCH (the default directory for model output files). For example, I have defined $SCRATCH in my .bashrc file: export SCRATCH=/mnt/parscratch/users/$USER Finally, the MATLAB startup.m file needs to add the main model directory to the MATLAB path, along with the directories containing the MATLAB code for StoichTools and the SkillMetrics toolbox (Table 1.1). 1.2 Development history Development of ARTEMIS commenced in 2020 with the aim of modelling the effect of rock dust treatments in a soy-maize-maize field experiment run by the Leverhulme Centre for Climate Change Mitigation (LC3M) (Kantola et al., 2017; Kantola et al., 2023; Beerling et al., 2024). The interface with PHREEQC (Parkhurst and Appelo, 2013) was loosely based on MATLAB software developed to model CO2consumption for the Hubbard Brook wollastonite treatment (Taylor et al., 2021) but was also informed by the PHREEQC script written by Peter Wade for the mesocosm study of Kelland et al. (2020). ARTEMIS has always included more complex biologically-mediated processes such as nutrient uptake, decomposition, nitrogen cycling and soil pCO2based on input soil respiration. It is partly based on established land and soil models (Neitsch et al., 2011; Lawrence et al., 23 March 2020) and on earlier geochemical models of the long-term global carbon cycle that I developed at the University of Sheffield (Taylor et al., 2016; Taylor et al., 2012). Subsequent beta-test or prototype versions of ARTEMIS were adapted with the aim of modelling several other LC3M experiments at both laboratory and field scale, and to produce runs for a model intercomparison project for enhanced weathering (Rock-MIP, Taylor et al., 2023; Taylor et al., 2024). Version 1.0 distils code from several of these experiment-specific prototypes, so there is no guarantee that users will find ARTEMIS easy to follow or modify. ARTEMIS was developed on a succession of high-performance computers (HPCs) at the University of Sheffield. The Energy Farm example simulations (see Chapter 4) were run in batch mode on a Dell PowerEdge C6420 running SUSE Liberty Linux 7 with the slurm job scheduling system and MATLAB version R2022a. That HPC was decommissioned almost immediately after completing those example runs, but before ARTEMIS was to be published. The code was transferred to its current home, a Dell PowerEdge R650 also running SUSE Liberty Linux 7 with slurm and MATLAB version R2023b. The “best” calibrated example run for the Energy Farm (see Table 4.1) was tested on the new HPC to check that the outputs had not changed. ARTEMIS does not require parallel processing or multiple cores. I have occasionally done short runs on the MATLAB command line on my desktop computer running various versions of Ubuntu and 3 MATLAB, but this is not the most convenient way to undertake a larger study. In batch mode on the HPC, the long Energy Farm runs covering the period 2016–2070 (Section 4.1 and Table 4.3) and treating apatite as an “implicit” phase (see Section 4.2 and the GMD manuscript) each took ≃7hours to run, while the five-year run with kinetic apatite and the cvode implicit solver required 3–4 hours. The run with kinetic apatite and the Runge-Kutta solver ran for over 2.5 hours before crashing after 1270 out of 1829 days. The quickest run times were 21–23 minutes for the five-year runs without sorption; most of remaining five-year runs required 33–60 minutes. Therefore, the ability to run simultaneous jobs in the background is useful for ARTEMIS. Figure 1.1: Overall steps for running this model. See Chapter 4for examples of how to run the software. Table 1.1: Software required to run ARTEMIS Software Required? Language Purpose MATLAB required Running the preprocessor and PHREEQC wrapper https://uk.mathworks.com/products/matlab.html PHREEQC required C and C++ Geochemical calculations https://www.usgs.gov/software/phreeqc-version-3 StoichTools required MATLAB chemical notation parsing https://uk.mathworks.com/matlabcentral/fileexchange/29774-stoichiometry-tools LC3M PSD software optional MATLAB Particle size distribution tracking, developed by Mark Lomas (Beerling et al., 2020) Not yet publicly available. Contact David Beerling for more information. 4 1.3 Two geochemical models ARTEMIS is the second RTM to be published by the LC3M. Several larger-scale LC3M publications (Beerling et al., 2020; Kantzas et al., 2022; Beerling et al., 2025) employed a different RTM developed independently by mathematician Mark Lomas under the auspices of the late geochemist Steven A. Banwart. Mark’s code (Beerling et al., 2025, https://doi.org/10.5281/zenodo.10940280) calculates potential CDR based on feedstock weathering, and is far more amenable to large-scale studies and emulation than ARTEMIS as it requires fewer inputs and shorter runtimes. Although I did not write the code for that model, I provided advice about how to code the alkalinity calculations (Beerling et al., 2020), and developed equations describing the overall effect of nitrogen cycling on pH (Kantzas et al., 2022) and, more recently, cation exchange (currently under development). 1.4 Disclaimer This documentation, along with the code and output files for the example model runs, is intended to accompany a model description paper intended for submission to Geoscientific Model Development (GMD). It is freely available on a “warts and all” basis, with no guarantees of any kind and no technical support. I will however try to answer questions relating to the ARTEMIS code and the GMD manuscript. I would like to thank my coauthors on the GMD manuscript describing ARTEMIS Version 1.0 and indeed other past and present members of the LC3M team who provided many helpful suggestions, data and support over the years. The name ARTEMIS was chosen by the LC3M director and manuscript coauthor, David Beerling; it is a better name than any I had come up with. However, neither he nor any other coauthors or colleagues are responsible for the errors, omissions or other problems with the Version 1.0 code or documentation that will undoubtedly come to light. I am entirely to blame for those. 11 Table 2.7 (see page 11) shows the required parameters that ought to appear in this table, along with some optional ways of entering some of the data. For example, the model requires volumetric porosity, field capacity and wilting point, as well as the initial water content. These may be specified for each layer, or they may be calculated from the soil texture. As these parameters may not all have been measured, the preprocessor can calculate them using the required soil texture data and one of the following pedotransfer functions (from the Community Land Model Version 5 as represented by Jinyun Tang in main/FuncPedotransferMod.F90): Cosby5 Cosby et al. (1984, their Table 5), default for this preprocessor and for CLM5 Cosby4 Cosby et al. (1984, their Table 4) NoilhanLacarrere1995 Noilhan and Lacarrère (1995) The desired pedotransfer function may be specified in defaultsettings.m (Section 3.1); otherwise Cosby5 will be used. By default, if no initial soil moisture is given, the model starts at field capacity for all cells. Data related to the cation exchange capacity are required, but there are several different options: Provide CEC and exchangeable cation data If the exchangeable cation data do not add to 100% of the given CEC, exchangeable acidity will be added (default: as protons). Provide CEC only Without exchangeable cation data, ARTEMIS will equilibrate the PHREEQC EXCHANGE block with the solution. Provide exchangeable cation data only CEC will be estimated as the sum of the data provided. Solid mineral or amorphous soil phase data may be provided in this file. These should be given in weight percent of the mineral soil (as sand/silt/clay percentages are) for each layer. The function setupsoilphases.m will look for fields corresponding to the phases in ds.phasetable in ds.cell; these phases may optionally have e.g. _native appended. Alternatively, their treatment as kinetic or equilibrium phases may be indicated by appending _kinetic or _eq to the phase name, e.g., Enstatite_kinetic, Calcite_eq. Kinetic treatment of phases also depends on the presence of both thermodynamic and kinetic parameters in the phases file; equilibrium phases require only thermodynamic data. Table 2.7: Soil layer input parameters topdepth Top depth cm Required bottomdepth Bottom depth cm Required bulkden Bulk density g cm−3Required parameter definition units options default purpose Continued on next page 12 Table 2.7: Soil layer input parameters (Continued) sandpct Sand wgt% Required for pedotransfer function and native minerals siltpct Silt wgt% Required for pedotransfer function claypct Clay wgt% Required for pedotransfer function TOC Total organic carbon wgt% Required Sorption to organic matter TON Total organic nitrogen wgt% Required Decomposition Moisture parameters used internally that can be given explicitly (default: use pedotransfer function) por Porosity volume fraction pedotransf. Maximum water content fc Field capacity volume fraction pedotransf. Maximum water content without flow wp Wilting point volume fraction pedotransf. Minimum water content (hygroscopic water) H2O Initial water volume fraction field capacity Initial water content at start of model run Cation exchange CEC Cation exchange capacity cmolc kg−1 soil Together with the bulk density and layer soil volume, CEC sets the size of the (clay) exchanger in the PHREEQC EXCHANGE block. Ca_exch Exchangeable Ca cmolc kg−1 soil or %CEC if CEC given Provides initial amount of an exchangeable element on the exchanger (clays as represented in the PHREEQC EXCHANGE block). Generally, <X>_exch for any exchangeable element <X>. Initial soil solutions (default: flood soil with rainwater) pH_soln pH Initial soil solution pH Ca_soln Ca mol L−1Initial soil solution elemental concentration. Generally, <X>_soln for element <X>, which may also include Amm (ammonium). Preferably all major ions should be given. parameter definition units options default purpose Continued on next page 13 Table 2.7: Soil layer input parameters (Continued) Solid kinetically-dissolving native soil phases (RSA est. from sand/silt/clay) not feedstocks Kfeldspar_kinetic Soil mineral wgt% Initial amount of a native soil phase represented kinetially. Generally, <PHASE>_kinetic where <PHASE> matches the name of a phase in the phase input file (Table 2.4). Secondary phases to include as EQUILIBRIUM_PHASES Calcite_eq Soil mineral mol g−1 soil Initial moles for equilibrium phase <PHASE> corresponding to a phase in the phase input file (Table 2.4). parameter definition units options default purpose 2.5 Feedstock PSD file An optional file containing measured particle size distribution data can be provided. This should be a .csv file with row names in column 1. These row headings should match the string describing each line in the description below. Subsequent columns contain contents indicated by column 1. NB strings are NOT quoted and must NOT contain commas. field There must be a column heading called radius or diameter providing the upper bin edge radius or diameter. Subsequent column headings give the name of each feedstock, of which there must be at least one. label This row contains a text label describing each column. units This row contains the units for the data found in each column. Units should be microns for radius or diameter columns and weight percent in each bin for feedstock columns. BETN2_m2g This metadata row contains the BET (Brunauer et al., 1938) surface area for each feedstock, in m2rock g−1rock. Enter 0 for the radius and diameter columns, or where the data are missing. geo_m2g This metadata row contains the geometric surface area for each feedstock, in m2rock g−1 rock. Enter 0 for the radius and diameter columns, or where the data are missing. lambda This metadata row contains BET/geo surface area ratios. Enter 0 for the radius and diameter columns, or where the data are missing. format This row contains the string %f in every column, indicating that MATLAB should treat the data as floating-point numbers. data These rows contain the measured data as described by the label and units rows. For feedstock columns, rows contain the weight percentage of feedstock in each bin. 14 The format of this file is a legacy from earlier model versions which has largely been deprecated elsewhere as it does not take advantage of MATLAB’s built-in routines to read .csv data. It was designed to make the units etc clear and allow inclusion of metadata (such as the surface area data above). An example of such a file is examples/EF_data/Lewis_feedstocks_PSD.csv 2.6 Timetable of events This file schedules feedstock and fertiliser applications, if any, as well as providing observed peak harvest data. It must include sufficient timestamp information to assign experimental day number and decimal year to each event (see Table 2.2). The resulting structure should have the fields shown in Table 2.8, some of which will be populated by the preprocessor (e.g., runindex) as noted, but others should be columns in the timetable file with headings exactly matching the variable names shown. The peaktot, peakroot and harvest data are the biomasses of the crop at “peak” (when biomass is largest) and harvest time. If these data are provided, and the crop is included in the file phreeqcdbincludefiles/SWAT_AppendixA_crop_params4.csv, the preprocessor can generate biomass and leaf area index (LAI) curves. This file contains all the data in Appendix A of the SWAT model (Neitsch et al., 2011) and some additional columns such as an alternative crop name and a flag indicating whether the crop is a legume or not. The daily change in biomass and LAI curves allow the preprocessor to calculate evapotranspiration and nutrient uptake from the soil layers. A limited set of common fertilisers will be recognised and formulas and molar masses will be generated for them, allowing their inclusion in the PHREEQC KINETICS block. They are assumed to be pelletised, but a RATES file will need to be written for them. There is a parameter “writetrRATESfiles” which, if set (ds.parm.writetrRATESfiles=1), will allow the creation of a RATES file for these treatments. Each fertiliser will have a PHREEQC RATE that is like the one in phreeqcdbincludefiles/Pellets.dat; these will all be written to Pellets.RATES. Run the preprocessor once to write this file (see Section 4.2 for examples). This keyword should be switched off for subsequent model runs (ds.parm.writetrRATESfiles=0) to prevent simultaneous batch jobs from trying to rewrite this file. This keyword appears in the examples/EF_CLM5/defaultsettings.m file (Section 3.1). Tillage and sowing depths need to be positive numbers; they should be zero for other timetable entries. The amount of feedstock and fertilisers in each model layer will be determined by the fraction of the layer which is above the given tillage depth. 2.7 Timeseries input files These .csv files should contain at least one column fully specifying the day for which the data are valid (see Table 2.2). If monthly data are provided they will be assumed to be valid for calendar day 15 and will be resampled. The sampling of the data in such a file must be the same for all variables in it; separate timeseries files can be provided for data with different sampling. Variables (other than time variables) that will be recognised by the preprocessor are shown in Table 2.9. The first line of the files gives the units. There may be several acceptable options for the units (see Table 2.9). For e.g. hydrological variables, units may be given per year, per month, or per 15 Table 2.8: Timetable structure arrays Variable name Value or units Description expID ID string for desired experiment (default: use all entries) treatment Name of feedstock, fertiliser or other treatment applied doseform feedstock, pellets Form of the applied treatment. Feedstocks are assumed to have powdered/granular form. At present, only pellets are recognised for fertilisers. dosage Amount of applied treatment (blank if none) doseunits t/ha, t/m2, kg/ha, kg/m2, mol/ha, mol/m2 Units for dosage applied (t, kg or mol, per ha or m2) molmass g/mol Molar mass of the treatment (1 if moles were the given units) formula Chemical formula for treatment (blank if none) feedstock Name of feedstock applied (blank if none, will be created from ds.timetable.treatment and ds.feedstocktable) crop Name of crop planted, harvested or sampled (fallow if none) runindex Index linking timetable events to arrays in ds.run (will be created when timetable file is read) molpm2 mol/m2 Calculated moles of treatment added per square meter of land (blank if none) dosekg kg/m2 sowdepth cm Sowing depth tillagedepth cm Depth of tillage peaktot g m2 land Total (dry) biomass peakabove g m2 land Aboveground (dry) biomass peakroot g m2 land Root (dry) biomass harvest g m2 land Yield (dry) at harvest day. The second line gives the name of the variable provided and these must be valid MATLAB variable names as described in the beginning of Chapter 2. Variable names that are meaningful to the preprocessor and model are shown in Table 2.9. 2.8 Depthwise time series input files The name of a depthwise time series input file should provide the name of the variable for which the data are intended. Examples include: depthtimeseries_soiltmp.csv Soil temperature data (units: Celsius) These .csv files should contain at least one column fully specifying the day for which the data are valid (see Table 2.2). If monthly data are provided they will be assumed to be valid for calendar 16 Table 2.9: Time series files of forcing data Line 2 heading Allowed units Description <timestamp> see Table 2.2 airtemp Celsius Air temperature precip mm/unit time Precipitation ET mm/unit time Total evapotranspiration from soil transpiration mm/unit time Total transpiration (default: calculate from ET) evaporation mm/unit time Total evaporation (default: calculate from ET) infiltration mm/unit time Infiltration (default: use precip) NPP gC m−1 landunit time−1Net primary productivity (or use peak biomass in timetable 2.8 soilresp gC m−1 landunit time−1Soil respiration (or estimate using biomass or NPP) day 15. The data will be resampled to provide daily values. The sampling of the data in such a file must be the same for all variables in it; separate timeseries files can be provided for data with different sampling. Variables (other than time variables) that will be recognised by the preprocessor are shown in Table 2.9. The first line of the files gives the units in the first column, assumed to be valid for all (nontimestamp) columns. There may be several acceptable options for the units (see Table 2.9). For (e.g.), hydrological variables, units may be given per year, per month, or per day. Columns that contain layer data may give the depth for which they are valid and the depth units, e.g. “15cm”. The second line gives the layer the name of the variable provided and these must be valid MATLAB variable names as described in the beginning of Chapter 2. Aside from the usual time-related variables (Table 2.2), these fields should contain the string “lay” followed by an integer giving the layer number, e.g., “lay1”. Any number of layers can be defined. These do not need to correspond to either the layers defined in the init_layers.csv file or the PHREEQC cells that will be defined in the model (Section 2.4). For example, the data can come from another model where different standard layers are defined. These data will be resampled to provide values for each PHREEQC cell in the model. 2.9 Observations files The names of files containing observations should include the string “obs_”. The observations files used for the Energy Farm example runs include: EF_obs_soil_drains_B0_SOTON.csv Drain solution chemistry (ppm) for the large “B0” plot 17 under soy-maize-maize rotation. EF_obs_soil_drains_insitu_SOTON.csv Drain solution data including temperature (◦C), pH, DOC (mg C/L), PO4 (µmol/L) for the large plots under soy-maize-maize rotation. EF_obs_soil_lysimeters_SOTON.csv Solution chemistry (ppm) for the both large “B0” plot (25 cm deep) and the small plots (50 cm deep) under soy-maize-maize rotation. EF_obs_soil_pH_largeplots_fromIlsaKantola_10Feb2021.csv Soil pH for depth ranges 0– 10cm and 10–30cm, for large plots under soy-maize-maize rotation. EF_obs_total_cCaMgrelease_means_stdev_acrossallblocks_TiCat.csv Ca and Mg feedstock release rates (tCO2/ha) from Beerling et al. (2024, Dataset S06). These are .csv files with headings on line 1 and units on line 2. Observations occupy the remaining lines of the file, generally one per sample. They include a datecollected column with a collection date in a MATLAB datetime format (e.g., yyyy-mm-dd, dd/mm/yyyy). The collection date is necessary because the preprocessor will assign each observation to a particular index in the model run. Likewise, files with a “soil” designation include depth (cm), and the preprocessor will assign a cell (layer) index to each observation. Data in files with a “total” designation are for the entire soil column and do not require a depth. For geochemical data, elements can be included. The strings “NO3”, “NH4” and “SO4” are also recognised, although NH4 will be replaced with redox-decoupled “AmmH”. Acceptable units for solutes include ppm, mg/L, mol/L, mmol/L, mM, umol/L and {\mu}mol/L. It is possible to include standard deviations for display purposes. For example, cCarelease and cCarelease_stdev along with “navg” giving the number of samples averaged are columns in the file EF_obs_total_cCaMgrelease_means_stdev_acrossallblocks_TiCat.csv which will allow calculation of the standard error for cumulative Ca release from feedstock. 2.10 Files in subdirectory “phreeqcdbincludefiles” These files include some bespoke PHREEQC databases with redox-decoupled nitrogen and some additional master species such as Ti(OH)4, Urea, and some organic acids. They are based on the Tipping_Hurley.dat database which ships with recent versions of PHREEQC and includes sorption to organic matter. The basic database is THAmmOrg.dat; THAmmOrgAl.dat is a version which includes organic matter complexation with Al from Erlandsson et al. (2016). THAmmOrgAlP.dat is an untested version which includes some P sorption data. There are several phase data files available in the phreeqcdbincludefiles/ directory (Table 2.10). Most of the kinetic data in these files comes from Palandri and Kharaka (2004) while most of the thermodynamic data are from the THERMODDEM (Blanc et al., 2012) Mineral mass and density data are from webmineral.com or were calculated using the StoichTools package (Table 1.1). The origins of the data for each phase are indicated in columns headed “gnotes”, “thermnotes”, “ratenotes”. 18 Table 2.10: The phreeqcdbincludefiles directory Filename default? Description PHREEQC databases and data-blocks THAmmOrg.dat default See text. Essentially the PHREEQC database Tipping_Hurley.dat with redox-decoupled nitrogen and a few additional species. THAmmOrgAl.dat As THAmmOrg.dat, but including Al complexation with organic matter (Erlandsson et al., 2016). THAmmOrgAlP.dat UNTESTED As THAmmOrgAl.dat, but including additional P complexation. SWATRATESPHASEStunable default Nitrification, denitrification, organic matter decomposition following the SWAT2009 code (Neitsch et al., 2011). plantRATES5 default Plant nutrient uptake of individual elements. This file assumes plants prefer nitrate to ammonium. Pellets.dat default Pelletised fertiliser rate law (Ritger and Peppas, 1987; Sofyane et al., 2020). A copy of this file is replicated for each fertiliser in the current directory. See Section ??. *exchcoeffs Alternative parameterisations for the EXCHANGE block. TRACERS Defines an intert tracer and Ti species (these are now included in the databases THAmmOrg.dat). Fakhraei_Driscoll_2015_Table1_organalog Alternative organic acids from Fakhraei and Driscoll (2015). Files useful for defining ds.phasetable feedstock_phases_Lewis_THERMODDEM.csv Phases for six feedstocks analysed by Lewis et al. (2021). feedstock_phases_Lewis_SUPCRT.csv As feedstock_phases_Lewis_THERMODDEM.csv, but some thermodynamic data from SUPCRT (e.g., Goddéris et al., 2006, their Table 2). feedstock_phases_pure.csv Pure wollastonite and forsterite. secondary_phases.csv A selection of secondary phases (calcite, silica, gibbsite, nesquehonite, anatase, kaolinite, smectite). Other data files cropstoichiometry.csv default Defines the elemental stoichiometry of different crops. MMuptakeparams.csv default Michaelis-Menten uptake parameters for different nutrients (e.g., Roose et al., 2001). Palandri_Kharaka_2004_params.csv Weathering rate parameters from Palandri and Kharaka (2004) (not in format of ds.phasetable). periodictable.csv required Periodic table data, including valence, PHREEQC master species, and PHREEQC EXCHANGE species where relevant (read by getperiodictable.m). Chapter 3 Model parameters and settings 3.1 The optional defaultsettings.m script Parameter settings (e.g., Tables 3.1,3.2,3.4), directories for data and outputs, PHREEQC cells and other entities may be defined in an optional MATLAB script named “defaultsettings.m” in the current directory. All other input files for the preprocessor (Figure 1.1, Table 2.1) are plain-text .csv files which will be read from a specified directory (ds.csvdir) which may be defined in defaultsettings.m or as a preprocessor argument. The preprocessor executes the optional defaultsettings.m script prior to processing any command-line arguments or reading any input .csv files. Settings from defaultsettings.m can therefore be overridden by command-line arguments and data in the .csv files. Parameters and other settings can be entered in the defaultsettings.m file. Here are a few examples from the defaultsettings.m file used for Energy Farm example runs: ds.dbdir=[getenv("HOME") ’/soilmodel/phreeqcdbincludefiles/’]; % Databases ds.dbfile=’THAmmOrgAl.dat’; % Includes Erlandsson et al’s Al sorption stuff ds.outdir=[getenv("SCRATCH") ’/EF_CLM5/’]; % Write to this directory ds.MIXdir=’topdown’; % MIX cell waters from the top down during percolation ds.SAmethod=’traditional’; % Simplistic shrinking sphere outside PHREEQC ds.EXCHmethod=’comp’; % Set EXCHANGE block explicitly (do not equilibrate) ds.parm.CECopt=’CECcarboxylic’; % Assign some CEC to carboxylic organic sites ds.parm.setexaciditysp=’Al’; % 5 February 2025 Exchangeable acidity species 3.1.1 Surface area method Kinetics for minerals and amorphous solid phases are processed in the RATES block of the PHREEQC input file (or included file). Each phase has its own custom BASIC code in RATES for calculating the change in moles of the phase within PHREEQC if there is a corresponding entry in the KINETICS block of the PHREEQC input file. It is possible to include BASIC code for simplistic changes in the reactive surface area for each phase, but at present there is no coupling of PHREEQC with other software packages. For example, 19 20 Table 3.1: Parameters controlling feedstock applications and extending runs in the preprocessor (setupdatastruc.m). Note that all feedstocks can be excluded at runtime (getphreeqcrun.m) using the “nodust” keyword. parameter default units description fsname Name of feedstock (will replace all feedstocks in the timetable) repeattimetable Decimal year Vector with two years: Timetable entries between the first year and the end of the timetable will be repeated until the total run extends through the second year. Example: For a timetable with a three-year crop rotation and last entry in 2022, [2020.0 2071.0] will extend the timetable through 2071 by repeating all timetable entries from 2020.0 through 2022. nofeedstockyears Calendar year Array of years where feedstock will not be applied. The array can be in any order. Example: [2023:3:2070 2024:3:2070] will apply feedstock every third year starting in 2025. ceasefeedstockyear Calendar year There will be no feedstock applications after this year (reduces feedstock doses to zero and resets doseform to null strung for all feedstocks in timetable entries after this year). Table 3.2: Parameters controlling reactive surface area (RSA). parameter default units description RSAscale 1 Scaling factor for RSA of feedstock phases nativeRSAscale 1 Scaling factor for RSA of native soil minerals useBET 1 Use measured feedstock BET surface area instead of geometric P80 mm Particle diameter below which 80% of feedstock particles are found (for optional Rosin-Rammler PSD) RRspread Spread parameter for optional Rosin-Rammler PSD (applies to feedstocks) BETN2_m2g m g−1BET surface area for feedstocks (can also be given with each feedstock PSD BET_m2g m g−1BET surface area for feedstocks (alternative name) geo_m2g m g−1Geometric surface area for feedstocks (calculated from PSD or from BET/λif not given) lambda λ, BET/geometric surface area ratio for feedstocks Chapter 4 Energy Farm example 4.1 Input files and scripts Data files for running the Energy Farm example simulations are in the examples/EF_data/ subdirectory (Table 4.5). The remaining files for the five-year runs are in the examples/EF_CLM5/ subdirectory, while files for the long runs through 2070 are in the examples/EF_longruns/ subdirectory. Each of these directories has its own defaultsettings.m file where the databases and files to be used are specified along with settings such as tillage depth. Many of the settings are the same for short and long runs, but the file examples/EF_longruns/defaultsettings.m contains parameters specifically for the long runs which specify how to repeat data from the timetable and how often to save intermediate files. Each run has its own MATLAB script (Tables 4.1,4.2, and 4.3) with corresponding shell script for batch submission. For example, rBcompscript.m is the script for a baseline five-year run with the Blueridge metabasalt feedstock (“B”) where no processes are excluded, and runrBcompscript.sh is the corresponding shell script which would be submitted using the sbatch command. Two shell scripts allow submission of the calibration runs (submittuning) and the “main” set of runs (submitall). Note that these jobs must be submitted from the directories where they are found (examples/EF_CLM5/ or examples/EF_longruns). Table 4.1: MATLAB scripts and settings for Energy Farm calibration runs found in the subdirectory examples/EF_CLM5/. All use the Blueridge feedstock for the soy-maize-maize crop rotation. File RSAscale impscale Apatite Solver Note riBimp5s1p1script.m 1.1 5 implicit Runge-Kutta best riBimp5s1p3script.m 1.3 5 implicit Runge-Kutta rBimp1s1script.m 1.0 1 implicit Runge-Kutta rBimp1s2script.m 2.0 1 implicit Runge-Kutta rBimp5s1script.m 1.0 5 implicit Runge-Kutta rBrkaps1script.m 1.0 N/A kinetic Runge-Kutta crashed rBcvaps1script.m 1.0 N/A kinetic cvode The MATLAB scripts for the five-year runs must be run in the examples/EF_CLM5 directory. They were designed to allow easy changing of some of the key settings for the runs, employing a MATLAB function EFrunsettings.m (in the same subdirectory) to set inputs for the preprocessor. 27 28 Table 4.2: MATLAB scripts for the example five-year Energy Farm model runs, found in the subdirectory examples/EF_CLM5/. File Feedstock Option Description rBcompscript.m Blueridge baseline No processes excluded. rBnonitscript.m Blueridge No N-cycle No nitrification, denitrification or N-fixation. rBclosedscript.m Blueridge Closed system CO2limitation rBnososcript.m Blueridge No sorption Excludes sorption including cation exchange. rCcompscript.m control baseline No processes excluded. rCclosedscript.m control Closed system CO2limitation rCnososcript.m control No sorption Excludes sorption including cation exchange. rFcompscript.m Forsterite baseline No processes excluded. rFclosedscript.m Forsterite Closed system CO2limitation rFnososcript.m Forsterite No sorption Excludes sorption including cation exchange. rWcompscript.m Wollastonite baseline No processes excluded. rWclosedscript.m Wollastonite Closed system CO2limitation rWnososcript.m Wollastonite No sorption Excludes sorption including cation exchange. Table 4.3: MATLAB scripts for the example extended Energy Farm model runs through 2070, found in the subdirectory examples/EF_longruns/. All runs except the control run have annual feedstock applications from 2016–2020 inclusive; Nis the total number of feedstock applications. File Feedstock NDescription rlBrep4script.m Blueridge 4 No feedstock applications after 2020. rlBrep23script.m Blueridge 23 After 2020, applications every third year. rlBrep55script.m Blueridge 55 Annual applications throughout the run. rlCscript.m Control 0 Control run (not used for figures). 4.2 How to do an example Energy Farm run First, unpack the tarball containing the ARTEMIS code as described in Chapter 1. Make sure that the environment variables $HOME and $SCRATCH are defined and the MATLAB path definition in the startup.m file includes the main model directory (Chapter 1). The required MATLAB software listed in Table 1.1 must also be installed and included in the MATLAB path. On the command line, it is possible to simply call the MATLAB script (e.g., rBcompscript.m) on the MATLAB command line and both the preprocessor and the model should run, tying up MATLAB for at least half an hour: > rBcompscript > r=tweakruns(r,1); % Creates a few extra fields and changes drainage cell Here, the “>” symbol is the MATLAB prompt and should not be typed. The results of the preprocessor will be in structure “ds”, and the results of the main model will be in structure “r”. The tweakruns.m postprocessor was designed for the Energy Farm runs and is useful for making the figures and tables for the ARTEMIS GMD paper. It will calculate some overall fluxes, tweak the run labels, remove some irrelevant observational data, and reset the layer for comparison with the drainage water. Its arguments are the output structure from getphreeqcrun.m and a flag to add the base saturation to the run label (no longer used, so this setting makes no difference). Use of EFrunsettings.m makes the calling sequences for the preprocessor and the main model less obvious. To do the same run specifying the calling sequences: 29 > % setupdatastruc.m is the preprocessor > ds=setupdatastruc(’SAmethod’,’traditional’,’fsname=Blueridge’,... ’EXCHmethod’,’comp’,’parm=ignoreSws’,1,... ’parm=implicit’,{’Titanite_blueridge’,’Hydroxyapatite_blueridge’},... ’runname=rBcomp’,’runlab=Ap+Tit implicit’,’parm=RSAscale’,1.1,... ’parm=impscale’,5); > rBcomp=getphreeqcrun(ds,’nyears’,5); % main model (will take >30min to run) > rBcomp=tweakruns(rBcomp,1); % bespoke postprocessor for Energy Farm examples Here, the MATLAB function setupdatastruc.m is the preprocessor, which will print a number of things about the run setup. It returns a data structure “ds” which should be passed to the main model function getphreeqcrun.m; a copy of this data structure (possibly with runtime changes) will be included in the output from getphreeqcrun.m as a substructure. In the above example, rBcomp is the ARTEMIS output structure; rBcomp.ds contains the data structure from the preprocessor. In these runs (and indeed all runs included in the example), the soil moisture from the CLM5 model is excluded (ignoreSws) because it makes the soils very wet and ruins the nitrogen-cycling outputs. Also, the CLM5 moisture can greatly exceed field capacity and at times approach total saturation. ARTEMIS runs would likely be invalid given the limitations associated with gas diffusion in ARTEMIS. The main function getphreeqcrun.m will first print some information about the run setup (such as which processes are included and which are excluded), and it will then call PHREEQC for each day of the run, read the resulting selected-output file (.sel), add those results to the ARTEMIS output arrays, and print a line stating whether that call to PHREEQC was successful or not along with a timestamp and the number of PHREEQC warnings. At the end of the run, the number of days processed, the number of PHREEQC warnings and the elapsed time will be printed along with the name of the output .mat file (path for output files redacted): ... rBcomp iteration 1825 Day 1825 starting: 06-Sep-2025 19:20:38: PHREEQC OK 1 warnings (gettotalfields/getfsrelease) Titanite_blueridge is not a kinetic phase so will be excluded from feedstock release arrays (gettotalfields/getfsrelease) Hydroxyapatite_blueridge is not a kinetic phase so will be excluded from feedstock release arrays Runtime for rBcomp (1825 days, 1828 warnings): 37.2153 minutes savemat <X>/rBcomp_Blueridge_Jan2016_Dec2020_06-Sep-2025.mat at 06-Sep-2025 19:20:40 In this case, titanite and hydroxyapatite are lumped together as an implicit phase, which is included in the feedstock elemental release arrays. The 1825 lines for the PHREEQC iterations, along with the other printed statements, will soon fill up the MATLAB command window. If running in batch mode, the resulting slurm-xxx.out file will contain all those output lines which can be perused as desired. This is especially useful if the model crashes, as for the calibration run rBrkaps1 (path redacted): ... savemat <X>/rBrkaps1_Blueridge_Jan2016_Dec2020_06-Sep-2025.mat at 06-Sep-2025 19:22:02 <X>/rBrkaps1_Blueridge_Jan2016_Dec2020_06-Sep-2025.mat and 30 <X>/rBrkaps1_Blueridge_Jan2016_Dec2020_06-Sep-2025.dmp_iter1095 SAVED at iteration 1095 ... rBrkaps1 iteration 1269 Day 1269 starting: 06-Sep-2025 21:10:20: PHREEQC OK 6 warnings rBrkaps1 iteration 1270 Day 1270 starting: 06-Sep-2025 21:15:13: ERROR: Bad RK steps > 1000. Please decrease (time)step or increase -bad_step_max. 5 warnings Error processing rBrkaps1 at iteration 1270 runday=1270 Runtime for rBrkaps1 (ERROR at day 1270, 537 warnings): 153.299 minutes savemat <X>/rBrkaps1_Blueridge_Jan2016_Dec2020_06-Sep-2025.mat at 06-Sep-2025 21:18:53 The partial run can be loaded, examined and displayed with the other runs. One can see when the run crashed, when the last output file was saved, and the error thrown by PHREEQC. Troubleshooting PHREEQC errors can be tricky but ARTEMIS does include a few options for tweaking the PHREEQC numerical method. Here, the error suggests that PHREEQC reached the maximum number of attempts to integrate the system of kinetic equations during a reaction step. The PHREEQC default for the number of steps (PHREEQC setting -bad_step_max) is 500 steps. This parameter can be adjusted when calling getphreeqcrun using the ARTEMIS badstepmax keyword, and the number of evenly-divided steps for KINETICS (currently 3) can be adjusted with the ARTEMIS nkinsteps keyword. However, the Runge-Kutta method is not ideal for runs with kinetic apatite as the set of equations to be solved is too stiff (Parkhurst and Appelo, 2013, pages 110–111). Most parameters can be set as shown for RSAscale and impscale (Table 3.2). The run label is for the figure legends; output files from getphreeqcrun.m will have names starting with the given runname (Table 4.4). These will include the standard PHREEQC input, output, selected output and dump files (Parkhurst and Appelo, 2013) for the first and last days of the run, along with some intermediate days specified by the save interval (default: 365 days). To save space, these intermediate files have not been retained in the tarballs containing the other output files. The MATLAB structure returned by getphreeqcrun.m will be saved in the similarly-named file with a .mat extension, which can be loaded on the MATLAB command line. The .mat file will contain results covering the entire run. Table 4.4: Output files for the Energy Farm example run “rBcomp”. File Description rBcomp_Blueridge_Jan2016_Dec2020_03-Sep-2025.dmp_init PHREEQC dump file for the first day of the run. rBcomp_Blueridge_Jan2016_Dec2020_03-Sep-2025.dmp_iternPHREEQC dump file for day nof the run. rBcomp_Blueridge_Jan2016_Dec2020_03-Sep-2025.dmp PHREEQC dump file for the last day of the run. rBcomp_Blueridge_Jan2016_Dec2020_03-Sep-2025.in_init PHREEQC input file for the first day of the run. rBcomp_Blueridge_Jan2016_Dec2020_03-Sep-2025.in_iternPHREEQC input file for day nof the run. rBcomp_Blueridge_Jan2016_Dec2020_03-Sep-2025.in PHREEQC input file for the last day of the run. rBcomp_Blueridge_Jan2016_Dec2020_03-Sep-2025.out_outit PHREEQC output file for the first day of the run. rBcomp_Blueridge_Jan2016_Dec2020_03-Sep-2025.out_iternPHREEQC output file for day nof the run. rBcomp_Blueridge_Jan2016_Dec2020_03-Sep-2025.out PHREEQC output file for the last day of the run. rBcomp_Blueridge_Jan2016_Dec2020_03-Sep-2025.sel_init PHREEQC selected-output file for the first day of the run. rBcomp_Blueridge_Jan2016_Dec2020_03-Sep-2025.sel_iternPHREEQC selected-output file for day nof the run. rBcomp_Blueridge_Jan2016_Dec2020_03-Sep-2025.sel PHREEQC selected-output file for the last day of the run. rBcomp_Blueridge_Jan2016_Dec2020_28-Aug-2025.screenout PHREEQC screen output from the last day of the run. rBcomp_Blueridge_Jan2016_Dec2020_29-Aug-2025.mat MATLAB file containing ARTEMIS output for the entire run. 31 Other settings will be found in the MATLAB scripts. Originally, runs allowing PHREEQC to equilibrate the EXCHANGE block (EXCHmethod=’eq’) were done for comparison with runs specifying the EXCHANGE composition explicitly (EXCHmethod=’comp’). Runs excluding the UAN fertilisers and plant nutrient uptake were also undertaken but these processes proved to have minor effects on the key outputs of interest (Ca and Mg release from feedstock and solution chemistry). Table 4.5: Data files for the example Energy Farm model runs, found in the subdirectory examples/EF_data/. File Description EF_CLM5_2016-2070_timeseries_CO2atm.csv Atmospheric CO2 from CLM5 for long runs. EF_CLM5_2016-2070_timeseries_infil_ET.csv Hydrological forcings from CLM5 for long runs. EF_CLM5_2016-2070_timeseries_soilresp.csv Soil respiration from CLM5 for long runs. EF_CLM5_2016-2070_timeseries_soiltemp.csv Soil temperature from CLM5 for long runs. EF_CLM5_timeseries_CO2atm.csv Atmospheric CO2 from CLM5 for five-year runs. EF_CLM5_timeseries_hydro.csv Hydrological forcings from CLM5 for five-year runs. EF_CLM5_timeseries_infil_ET.csv Infiltration from CLM5 for five-year runs. EF_CLM5_timeseries_soilresp.csv Soil respiration from CLM5 for five-year runs. EF_CLM5_timeseries_soiltemp.csv Soil temperature from CLM5 for five-year runs. EF_maize_B0_timetable.csv Timetable of events at the Energy Farm. EF_obs_soil_drains_B0_SOTON.csv Drainage major ion chemistry for the large “B0” metabasalt-treated plot for the Energy Farm soy-maize-maize rotation. EF_obs_soil_drains_insitu_SOTON.csv Drainage pH for the large “B0” metabasalttreated plot for the Energy Farm soy-maizemaize rotation. EF_obs_soil_lysimeters_SOTON.csv Major ion chemistry measured in lysimeters for the large “B0” metabasalt-treated plot (25 cm depth) and the four small plots (50 cm depth) for the Energy Farm soy-maize-maize rotation. EF_obs_soil_pH_largeplots_fromIlsaKantola_perscomm_10Feb2021.csv Soil pH for the large “B0” metabasalt-treated plot for the Energy Farm soy-maize-maize rotation. EF_obs_total_cCaMgrelease_means_stdev_acrossallblocks_TiCat.csv Ca and Mg weathering release data from Beerling et al. (2024). initlayers_Flanagan_EF_0-183cm.csv Soil layer data for one of the most common soils in the Energy Farm soy-maize-maize rotation plots. This file contains metadata including the data sources. Lewis_feedstocks_PSD.csv Particle size distribution data for the six feedstocks analysed by Lewis et al. (2021), including the Blueridge metabasalt applied at the Energy Farm. rainwater_Champaign_NTN_2021_initialChemistry.csv Rainwater chemistry represent an average for the rainwater collected over the period starting 2020-11-17 15:15 and ending 2020-11-24 15:15 by the National Trends Network under the National Atmospheric Deposition Program for Champaign, IL, USA, downloaded on Friday 25 June 2021. If all given species are included, percent error in the charge balance of this solution is 0.0351%. Richland_loess_Batestown_till_phases_combined.csv Data for minerals associated with the parent materials at the Energy Farm. 32 4.3 Loading and displaying the Energy Farm results The first thing to do is to unpack the tarballs containing the model runs. These scripts assume that ARTEMIS_GMD_EF_CLM5_outputfiles_Sep-2025.tar.gz was unpacked in examples/EF_CLM5/ and ARTEMIS_GMD_EF_longruns_outputfiles_Sep-2025.tar.gz in examples/EF_longruns/, creating outputfiles subdirectories. There are two MATLAB scripts which I used to load the outputs for most of the five-year runs: EF_CLM5/loadtuningruns.m Loads the five-year-long runs done for calibration with the blueridge feedstock, creating MATLAB cell array of structures “rt”. EF_CLM5/loadmainruns.m Loads the five-year-long runs done with different feedstocks, creating MATLAB cell array of structures “r” There are five MATLAB scripts which were used to create the figures for the GMD paper, two of which also load ARTEMIS runs: EF_CLM5/Flanaganlayerdiagram.m Creates GMD location map and layer diagram (load the main runs with EF_CLM5/loadmainruns.m first as the layer data will be read from r1,1) EF_CLM5/tuningfigs.m GMD figures related to calibration with the blueridge feedstock and loaded with EF_CLM5/loadtuningruns.m EF_CLM5/feedstockfigs.m GMD figures related to runs done with different feedstocks and loaded with EF_CLM5/loadmainruns.m EF_CLM5/BRnonitfigs.m Loads runs and creates GMD figures related to the effect of the nitrogen cycle EF_CLM5/longtermfigs.m Loads runs covering the period 2016–2070 and creates GMD figures and tables related to lag times and effects of continued treatments with the blueridge feedstock Each of these scripts contains MATLAB switches which determine which figures will be made. These should be edited to choose which figures to look at. They also contain a variable “png” which determines whether the figures are saved or not; this is set to zero (do not save) if it does not exist. 4.4 Energy Farm calibration runs ARTEMIS was calibrated such that the Ca and Mg release from feedstock was within one standard error of release rates derived by Beerling et al. (2024) using the “TiCAT” method (Reershemius et al., 2023) whereby changes in cation content of soil samples are compared to an “immobile” element, in this case Ti. Little or no attempt was made to calibrate other ARTEMIS outputs, other than to change a soil moisture threshold for denitrification. As a first attempt, all feedstock phases with rate laws were modelled kinetically using rate laws from Palandri and Kharaka (2004), and one Ca-bearing phase with no kinetic or thermodynamic parameters (titanite) was ignored. PHREEQC has two numerical solvers (Parkhurst and Appelo, 33 2013). Using the default Runge-Kutta (RK) solver with six time subintervals and three kinetic steps, PHREEQC crashed with convergence errors before completing the five-year run. Another run with PHREEQC’s implicit stiff ordinary differential equation solver “cvode” did not crash, but took 258 minutes to run for 1825 days. Cumulative Mg release was acceptable in the kinetic cvode run, but Ca release was low compared to TiCAT (Fig. 4.1), suggesting a problem with the representation of one or more Ca-bearing feedstock phases. One such phase is titanite, which has no kinetic rate parameters. Another is apatite (Table 4.6), which contains phosphorus. In both kinetic runs, apatite weathered too quickly such that P was released and leached during the fallow seasons after feedstock application, leaving little or none for plant uptake during the following growing seasons (Section 4.4.1). ARTEMIS includes kinetic plant P uptake and P release from decomposing organic matter, and the PHREEQC database includes phosphate speciation and sorption to hydrous ferric oxides. However, P can also sorb to a variety of other surfaces (e.g., Jalali et al., 2022, and references in their Table S1) which could potentially reduce the apatite saturation index and increase its weathering rates, but would also retain more P in the system. P cycling in ARTEMIS could be improved with surfaces specifically developed and calibrated for the Energy Farm soils, but this is outside the scope of the present example. Table 4.6: Mineralogy of the “blueridge” metabasalt (Lewis et al., 2021) applied at the Energy Farm (2016–2019). All parameters for these phases are given in the input file “feedstock_phases_Lewis_THERMODDEM.csv”, included with the software and example files. Thermodynamic data are from the THERMODDEM database (BRGM, 2020; Blanc et al., 2012). Quartz kinetic parameters are from Bandstra and Brantley (2008) and Tester et al. (1994); other phases follow Palandri and Kharaka (2004). The “kinetics” column shows how each phase was treated in the final calibrated model run. Phase Formula Weight% kinetics Albite NaAlSi3O819.6 irreversible Ferroactinolite Ca2(Mg0.75Fe0.25)5Si8O22(OH)211.6 irreversible Epidote Ca2Al2(Al0.84Fe0.16)(SiO4)(Si2O7)O(OH) 25.6 irreversible Chlorite (Mg0.63Fe0.37)5Al(Si3Al)O10(OH)836.3 irreversible Quartz SiO25.2 reversible Titanite CaTiSiO51.7 implicit Muscovite KAl3Si3O10(OH)1.8F0.23.3345 irreversible Apatite Ca5(PO4)3(OH)11.9325 implicit As expected, a test run ignoring both apatite and titanite produced poor results, with cumulative Ca+Mg release over five tCO2ha−1too low (Fig. 4.1). Alternatively, apatite can be treated as an implicit phase along with titanite, so that they effectively dissolve in tandem with the kinetic feedstock phases. Using this implicit “phase” without any scaling, however, does not improve the RMSEs, and increasing the whole-feedstock initial reactive surface area using the RSAscale keyword (Table 3.2) improves model Ca release at the expense of Mg release. Instead, a combination of slightly increased feedstock RSA scaling by a factor of 1.1 and substantially increased “implicit” phase dissolution rate by a factor of five produced results deemed acceptable for the purposes of this example study (Fig. 4.1, blue curve). Runs where apatite was not a kinetic phase took approximately 40 minutes to run on the high-performance computer at the University of Sheffield. 34 Figure 4.1: Comparison between model calibration runs (solid lines) and TiCAT (blue dots, Beerling et al., 2024; Reershemius et al., 2023) Ca+Mg (a), Ca (b) and Mg (c) release rates from the “blueridge” feedstock at the Energy Farm (2016–2019). The thick blue solid line represents the best model run used in subsequent experiments in the current study. 4.4.1 Phosphorus dynamics during calibration runs The calibrating exercise served to underscore key limitations with P cycling in ARTEMIS, which can be observed in the topsoil model P pools (Fig. 4.2). Where treated as a kinetic phase, fast-weathering apatite released P during the fallow season. P was then leached before it could be taken up by plants, suggesting insufficient P sorption. Feedstock apatite was the main source of P in these runs, with only small amounts being released from organic matter. Unsurprisingly, if apatite was ignored, there was too little P in the system. With “implicit” apatite and titanite, plant P uptake improved, but P was released very slowly and apatite accumulated in the soil. This effect is partially mitigated if apatite (and titanite) dissolution was increased relative to dissolution of kinetically-weathering phases. The best results were obtained by increasing the implicit dissolution by a factor of five, together with a modest 1.1-fold increase in the initial reactive surface area of the whole feedstock (Fig. 4.2a). Plant P demand was however not quite satisfied during these five-year runs (Fig. 4.3), which also suggests insufficient P retention in the model soil. 35 Figure 4.2: Topsoil P pools during calibrating runs. Here, grey is the P remaining in feedstock, blue is P in solution, and green is P taken up by the plant. Apatite dissolves quickly enough to facilitate plant P uptake in the run with the best RMSE for Ca+Mg release from feedstock (a). Other runs with slower “implicit” apatite dissolution provide less P for plants and allow apatite to accumulate (b, c). Without apatite (c), or with apatite treated as a fast-weathering kinetic phase (d,e), there is little P available for plant uptake. 36 Figure 4.3: Total cumulative P uptake from all soil layers by plants during calibrating runs. Here, the thin blue line is total plant P demand, and the thick lines are P uptake curves for the calibrating runs. 43 Figure 5.3: The main function getphreeqcrun.m writes the first PHREEQC input file which follows the sequence of events shown. showstackedfields.m Creates a stack of plots, calling showrun, showCEC, showsurfaces, or showcelleltbalance for multiple cells or multiple runs. showAlk Shows the major ions contributing to the alkalinity balance. 44 Figure 5.4: The main function getphreeqcrun.m writes the PHREEQC input file which follows the sequence of events shown. Bibliography Bandstra, Joel Z and Susan L Brantley (2008). Data fitting techniques with applications to mineral dissolution kinetics. In: Kinetics of water-rock interaction. Springer, 211–257 (cit. on p. 33) Beerling, David J, Dimitar Z Epihov, Ilsa B Kantola, Michael D Masters, Tom Reershemius, Noah J Planavsky, Christopher T Reinhard, Jacob S Jordan, Sarah J Thorne, James Weber, et al. (2024). Enhanced weathering in the US Corn Belt delivers carbon removal with agronomic benefits. Proceedings of the National Academy of Sciences 121, e2319436121. doi:https://doi.org/10.1073/pnas.2319436121 (cit. on pp. 2,17,31,32,34) Beerling, David J, Euripides Kantzas, Mark R Lomas, Peter Wade, Rafael M Eufrasio, Phil Renforth, Binoy Sarkar, M Grace Andrews, Rachael H James, Christopher R Pearce, et al. (2020). Potential for large-scale CO2removal via enhanced rock weathering with croplands. Nature 583, 242–248. doi:https://doi.org/10.1038/s41586-0202448-9 (cit. on pp. 3,4,21,23) Beerling, David J, Euripides P Kantzas, Mark R Lomas, Lyla L Taylor, Shuang Zhang, Yoshiki Kanzaki, Rafael M Eufrasio, Phil Renforth, Jean-Francois Mecure, Hector Pollitt, et al. (2025). Transforming US agriculture for carbon removal with enhanced weathering. Nature, 1–10. doi:10.1038/s41586-024-08429-2 (cit. on p. 4) Ben-Noah, Ilan and Shmulik P Friedman (2018). Review and evaluation of root respiration and of natural and agricultural processes of soil aeration. Vadose Zone Journal 17, 1–47. doi:https://doi.org/10.2136/vzj2017.06.0119 (cit. on p. 22) Blanc, Ph, Arnault Lassin, Patrice Piantone, Mohamed Azaroual, Nicolas Jacquemet, A Fabbri, and Eric C Gaucher (2012). Thermoddem: A geochemical database focused on low temperature water/rock interactions and waste materials. Applied geochemistry 27, 2107–2116. doi:10.1016/j.apgeochem.2012.06.002 (cit. on pp. 17,33) BRGM (2020). Thermoddem PHREEQC database. Version 1.10, 15 December 2020 (cit. on p. 33) Brunauer, Stephen, Paul Hugh Emmett, and Edward Teller (1938). Adsorption of gases in multimolecular layers. Journal of the American chemical society 60, 309–319 (cit. on p. 13) Cerling, Thure E (1991). Carbon dioxide in the atmosphere: evidence from Cenozoic and Mesozoic paleosols. American Journal of Science;(United States) 291, 377–400. doi:https://doi.org/10.2475/ajs.291.4.377 (cit. on pp. 22, 25,38) Choudhury, Bhaskar J (2001). Modeling radiation-and carbon-use efficiencies of maize, sorghum, and rice. Agricultural and Forest Meteorology 106, 317–330. doi:https://doi.org/10.1016/S0168-1923(00)00217-3 (cit. on p. 25) Cosby, BJ, GM Hornberger, RB Clapp, and ToR Ginn (1984). A statistical exploration of the relationships of soil moisture characteristics to the physical properties of soils. Water Resources Research 20, 682–690. doi:https: //doi.org/10.1029/WR020i006p00682 (cit. on p. 11) Dufrene, E, R Ochs, and B Saugier (1990). Photosynthèse et production du palmier à huile en relation avec les facteurs climatiques. Oléagineux 45, 8–9 (cit. on p. 25) Dzombak, David A and Francois MM Morel (1990). Surface complexation modeling: hydrous ferric oxide. John Wiley & Sons (cit. on p. 21) Erlandsson, M, EH Oelkers, Kevin Bishop, H Sverdrup, S Belyazid, JLJ Ledesma, and SJ Köhler (2016). Spatial and temporal variations of base cation release from chemical weathering on a hillslope scale. Chemical Geology 441, 1–13. doi:https://doi.org/10.1016/j.chemgeo.2016.08.008 (cit. on pp. 17,18) 45 46 Fakhraei, Habibollah and Charles T Driscoll (2015). Proton and aluminum binding properties of organic acids in surface waters of the northeastern US. Environmental Science & Technology 49, 2939–2947. doi:https://doi. org/10.1021/es504024u (cit. on p. 18) Goddéris, Yves, Louis M François, Anne Probst, Jacques Schott, David Moncoulon, David Labat, and Daniel Viville (2006). Modelling weathering processes at the catchment scale: The WITCH numerical model. Geochimica et Cosmochimica Acta 70, 1128–1147. doi:10.1016/j.gca.2005.11.018 (cit. on p. 18) Hoffmann, Munir P, A Castaneda Vera, MT Van Wijk, Ken E Giller, Thomas Oberthuer, Christopher Donough, and Anthony M Whitbread (2014). Simulating potential growth and yield of oil palm (Elaeis guineensis) with PALMSIM: Model description, evaluation and application. Agricultural Systems 131, 1–10. doi:https://doi.org/10.1016/ j.agsy.2014.07.006 (cit. on p. 25) Jalali, Mohsen, Elham Amirabadi Farahani, and Mahdi Jalali (2022). Simulating phosphorus leaching from two agricultural soils as affected by different rates of phosphorus application based on the geochemical model PHREEQC. Environmental Monitoring and Assessment 194, 164. doi:10.1007/s10661-022-09828-6 (cit. on p. 33) Kantola, Ilsa B, Elena Blanc-Betes, Michael D Masters, Elliot Chang, Alison Marklein, Caitlin E Moore, Adam von Haden, Carl J Bernacchi, Adam Wolf, Dimitar Z Epihov, et al. (2023). Improved net carbon budgets in the US Midwest through direct measured impacts of enhanced weathering. Global Change Biology 29, 7012–7028. doi: https://doi.org/10.1111/gcb.16903 (cit. on p. 2) Kantola, Ilsa B, Michael D Masters, David J Beerling, Stephen P Long, and Evan H DeLucia (2017). Potential of global croplands and bioenergy crops for climate change mitigation through deployment for enhanced weathering. Biology letters 13, 20160714. doi:https://doi.org/10.1098/rsbl.2016.0714 (cit. on p. 2) Kantzas, E, M Val Martin, M Lomas, R Eufrasio, P Renforth, A Lewis, L Taylor, J-F Mecure, H Pollitt, P Vercoulen, N Vakilifard, P Holden, N Edwards, L Koh, N Pidgeon, S Banwart, and D Beerling (2022). Substantial carbon drawdown potential from enhanced rock weathering in the United Kingdom. Nature Geoscience 15, 382–389. doi: https://doi.org/10.1038/s41561-022-00925-2 (cit. on pp. 4,21,23) Kelland, Michael, Peter Wade, Amy Lewis, Lyla Taylor, Binoy Sarkar, Grace Andrews, Mark Lomas, Anne Cotton, Simon Kemp, Rachael James, Chris Pearce, Sue Hartley, Mark Hodson, Jonathan R Leake, Steve Banwart, and David Beerling (2020). Increased yield and CO2sequestration potential with the C4cereal Sorghum bicolor cultivated in basaltic rock dust-amended agricultural soil. Global Change Biology 26, 3658–3676. doi:10.1111/gcb.15089 (cit. on pp. 2,7) Kirschbaum, Miko UF and Rowena Mueller (2001). Net Ecosystem Exchange: Workshop Proceedings. Cooperative Research Centre for Greenhouse Accounting (cit. on p. 25) Lamade, Emmanuelle, Narcisse Djégui, and Philippe Leterme (1996). Estimation of carbon allocation to the roots from soil respiration measurements of oil palm. Plant and soil 181, 329–339. doi:https://doi.org/10.1007/ BF00012067 (cit. on p. 25) Lawrence, D, R Fisher, C Koven, K Oleson, S Swenson, and M Vertenstein (23 March 2020). Technical description of version 5.0 of the Community Land Model (CLM). National Center for Atmospheric Research. PO Box 3000, Boulder, CO 80307-300 (cit. on pp. 2,41) Lewis, Amy L, Binoy Sarkar, Peter Wade, Simon J Kemp, Mark E Hodson, Lyla L Taylor, Kok Loong Yeong, Kalu Davies, Paul N Nelson, Michael I Bird, et al. (2021). Effects of mineralogy, chemistry and physical properties of basalts on carbon capture potential and plant-nutrient element release via enhanced weathering. Applied Geochemistry 132, 105023. doi:10.1016/j.apgeochem.2021.105023 (cit. on pp. 18,31,33) Neitsch, Susan L, Jeffrey G Arnold, Jim R Kiniry, and Jimmy R Williams (2011). Soil and Water Assessment Tool theoretical documentation version 2009. Tech. rep. Texas Water Resources Institute (cit. on pp. 2,14,18,22) Noilhan, J and P Lacarrère (1995). GCM grid-scale evaporation from mesoscale modeling. Journal of Climate 8, 206– 223. doi:https://doi.org/10.1175/1520-0442(1995)008<0206:GGSEFM>2.0.CO;2 (cit. on p. 11) Palandri, James L and Yousif K Kharaka (2004). A compilation of rate parameters of water-mineral interaction kinetics for application to geochemical modeling. Tech. rep. US Geological Survey. doi:https://doi.org/10. 3133/ofr20041068 (cit. on pp. 17,18,32,33) Parkhurst, D.L. and C.A.J. Appelo (2013). Description of input and examples for PHREEQC version 3—A computer program for speciation, batch-reaction, one-dimensional transport, and inverse geochemical calculations. Techniques and Methods. U.S. Geological Survey. Chap. A43, 497 (cit. on pp. 1,2,24,30,32) 47 Raich, James W and William H Schlesinger (1992). The global carbon dioxide flux in soil respiration and its relationship to vegetation and climate. Tellus B 44, 81–99. doi:https://doi.org/10.1034/j.16000889.1992.t01-100001.x (cit. on p. 25) Reershemius, Tom, Mike E Kelland, Jacob S Jordan, Isabelle R Davis, Rocco D’Ascanio, Boriana Kalderon-Asael, Dan Asael, T Jesper Suhrhoff, Dimitar Z Epihov, David J Beerling, et al. (2023). Initial validation of a soil-based massbalance approach for empirical monitoring of enhanced rock weathering rates. Environmental Science & Technology 57, 19497–19507. doi:10.1021/acs.est.3c03609 (cit. on pp. 32,34) Renforth, Phil and Gideon Henderson (2017). Assessing ocean alkalinity for carbon sequestration. Reviews of Geophysics 55, 636–674. doi:https://doi.org/10.1002/2016RG000533 (cit. on p. 7) Ritger, Philip L and Nikolaos A Peppas (1987). A simple equation for description of solute release I. Fickian and non-fickian release from non-swellable devices in the form of slabs, spheres, cylinders or discs. Journal of Controlled Release 5, 23–36. doi:10.1016/0168-3659(87)90034-4 (cit. on p. 18) Roose, Tiina, AC Fowler, and PR Darrah (2001). A mathematical model of plant nutrient uptake. Journal of mathematical biology 42, 347–360. doi:10.1007/s002850000075 (cit. on p. 18) Sofyane, Asma, Mohammed Lahcini, Abdellatif El Meziane, Mehdi Khouloud, Abdelmalek Dahchour, Sylvain Caillol, and Mustapha Raihane (2020). Properties of coated controlled release diammonium phosphate fertilizers prepared with the use of bio-based amino oil. Journal of the American Oil Chemists’ Society 97, 751–763. doi:10.1002/ aocs.12360 (cit. on p. 18) Taylor, L, S. Banwart, P. Valdes, J. Leake, and D. Beerling (2012). Evaluating the effects of the co-evolution of terrestrial ecosystems, climate and CO2on continental weathering over geological time: a global-scale process-based approach. Philosophical Transactions of the Royal Society B 367, 565–582. doi:10.1098/rstb.2011.0251 (cit. on p. 2) Taylor, L, C Driscoll, P Groffman, G Rau, J Blum, and D Beerling (2021). Increased carbon capture by a silicate-treated forested watershed attenuated by acid deposition. Biogeosciences 8, 169–188. doi:10.5194/bg-18-169-2021 (cit. on p. 2) Taylor, L, M Lomas, Y Kanzaki, E Bolton, N Planavsky, C Reinhard, M Kelland, and D Beerling (2023). Rock-MIP: The Enhanced Rock Weathering Model Intercomparison Project. American Geophysical Union abstract GC51K0742, 15 December 2023, San Francisco (cit. on p. 2) Taylor, L, J Quirk, R Thorley, P Kharecha, J Hansen, A Ridgwell, M Lomas, S Banwart, and D Beerling (2016). Enhanced weathering strategies for stabilizing climate and averting ocean acidification. Nature Climate Change 6, 402–406. doi:10.1038/nclimate2882 (cit. on p. 2) Taylor, L, C Reinhard, M Lomas, Y Kanzaki, E Bolton, N Planavsky, and D Beerling (2024). Enhanced Rock Weathering Model Intercomparison Project Phase 1: Responses to Silicate Feedstocks and soil pCO2. American Geophysical Union abstract 1681700 talk B14B-03, 15 December 2024, Washington DC (cit. on p. 2) Tester, Jefferson W, W Gabriel Worley, Bruce A Robinson, Charles O Grigsby, and Jeffrey L Feerer (1994). Correlating quartz dissolution kinetics in pure water from 25 to 625 ◦C. Geochimica et cosmochimica acta 58, 2407–2420. doi: 10.1016/0016-7037(94)90020-5 (cit. on p. 33)