User Guide

Installation and usage

Installation

First, download and install Julia: https://julialang.org/downloads/. If you like you can also install julia-vscode.org to get a complete IDE. Both these steps are nicely summarized by https://techytok.com/julia-vscode/. ALternatively, you can also follow the official documentation of the VS Code Julia: https://www.julia-vscode.org/docs/dev/gettingstarted/#Installation-and-Configuration-1.

To install LWFBrook90.jl open a Julia REPL, enter the Pkg REPL by pressing ] and add the package:

(@v1.7) pkg> add LWFBrook90
(@v1.7) pkg> status

Dependencies of LWFBrook90.jl should automatically be installed. After that hit Ctrl-C to quit the Pkg REPL and return to the default Julia REPL. Now the installation can be tested with a single line using LWFBrook90; run_example().

Usage

Check out a step-by-step guide for a simulation in section Example Script 01

The steps in a typical simulation script are:

• load the package using LWFBrook90
• read input data model = loadSPAC()
  • optional arguments in loadSPAC() can be used to define model parameters (-> type ?loadSPAC to see documentation)
• pre-process model for simulation simulation = setup(model)
• compute simulation simulate!(simulation)
• post-process and plot simulation results, get_fluxes(simulation), get_states(simulation), get_forcing(simulation)

LWFBrook90.jl outputs quantities in daily resolution. Monthly or yearly quantities can be derived in the post processing. ~~For performance reasons, e.g. for Bayesian parameter estimation, computation of these additional quantities can also be deactivated during simulation, and they could be calculated in a post-processing step from the state vector.~~

Below is a code snippet providing a quick start. For more details refer to examples/example-script-01.jl.

# load problem defintion as SPAC object
input_path = "../../../examples/DAV2020-full/"; input_prefix = "DAV2020-full"; # make sure to get the path correct relative to `pwd()``
model = loadSPAC(input_path, input_prefix; simulate_isotopes = true);

# transform into a system of ODEs (space discretization etc...)
simulation     = setup(model)

# run simulation
simulate!(simulation)

# Get simulation output (states, fluxes, forcing) as DataFrames
get_states(simulation)
get_fluxes(simulation)
get_forcing(simulation)

using Plots, Measures; gr();
pl2 = plotisotopes(simulation, :d18O, (d18O = :auto, d2H = :auto), :showRWUcentroid)
pl2

pl1 = plotamounts(simulation, :above_and_belowground, :showRWUcentroid)
pl1

pl3 = plotforcingandstates(simulation)
pl3

Soil evaporation depth and isotope fractionation

Soil evaporation has two independent controls: the depth from which water is removed, and whether the evaporated vapor fractionates relative to the soil water.

Choosing the evaporation depth

The default SLVPDEPTH_m = 0.0 draws evaporation from the uppermost input soil layer. A value smaller than this layer's thickness also uses its full thickness. For example, with a 3 cm surface layer, both 0.0 and 0.01 use the top 3 cm; 0.10 uses the top 10 cm. The requested depth must fit within the soil profile.

For a deeper source zone, the discretization is refined at its lower boundary when needed. Evaporation is apportioned by layer thickness within this zone. Thus a 3 cm layer and a 7 cm layer in a 10 cm source zone provide 30% and 70% of the soil evaporation flux, respectively. Layers below the zone do not directly supply soil evaporation, although water and isotopes can still move between layers. Use a sufficiently fine surface grid to resolve the enrichment pattern of interest.

Enabling evaporation fractionation

With simulate_evaporation_fractionation = false, the evaporating flux has the same isotopic composition as the soil water in each contributing layer. Soil isotope values can still change through precipitation, transport, and mixing.

With simulate_evaporation_fractionation = true, the model computes the evaporated vapor's δ¹⁸O and δ²H using Craig–Gordon fractionation for each contributing layer. The calculation uses daytime air temperature, atmospheric humidity, and a kinetic fractionation exponent that varies with local soil moisture. Atmospheric vapor is assumed to be in isotopic equilibrium with the precipitation forcing; a separate measured vapor-isotope time series is not used. The option changes the isotopic composition of the evaporation flux; water loss is calculated by the hydrologic model. It applies to soil evaporation, rather than being a general switch for fractionation in interception or snow stores.

For a model loaded with simulate_isotopes = true, configure and run a simulation as follows. remakeSPAC returns a simulation ready for simulate!, so another setup call is unnecessary.

simulation = remakeSPAC(model;
    params = (SLVPDEPTH_m = 0.10,),
    solver_options = (simulate_evaporation_fractionation = true,),
    requested_tspan = model.tspan)
simulate!(simulation; save_everystep = false,
    saveat = range(model.tspan...; step = 1.0))

# Postprocessing depths are in millimeters, unlike SLVPDEPTH_m.
soil_isotopes = get_soil_([:d18O, :d2H], simulation;
    depths_to_read_out_mm = [15, 65, 150, 250])

Alternatively, enable fractionation while loading the model:

model = loadSPAC(input_path, input_prefix;
    simulate_isotopes = true,
    simulate_evaporation_fractionation = true)

Set SLVPDEPTH_m in param.csv, or change it afterward with remakeSPAC. Existing parameter files that omit it use the default 0.0 m.

Comparing the effects

Example Script 02 compares three simulations with identical forcing and initial conditions: 10 cm without fractionation, 10 cm with fractionation, and 30 cm with fractionation. The first pair isolates the fractionation option; the second pair isolates the evaporation depth. Heatmaps use a shared color scale for each isotope. Install Plots in the active Julia environment to run the plotting example.

During drying, a shallow source zone can enrich more rapidly near the surface, while a deeper source zone spreads the evaporative loss and enrichment over more soil. Rainfall, transport, and mixing can interrupt this pattern, so the meteorological example need not show monotonic enrichment. The default combination of 0.0 m and false preserves the previous soil-evaporation configuration.

Input data

Overview of input data

To run a simulation following input files are expected in a single folder

• soil_horizons.csv - containing the hydraulic parameters of the different soil horizons
• param.csv - containing scalar model parameters
• meteoveg.csv - containing daily values of meteorologic variables [and stand properties]
• meteoiso.csv - containing isotopic signatures of aggregate precipitation samples
• initial_conditions.csv - containing initial conditions of scalar state variables
• soil_discretization.csv - containing the initial conditions of the soil water status in the form of the soil matric potential (kPa) and the initial isotopic signatures (vector state variables); parameters of relative root density distributions defined for the soil domain; as well as the definition of the numerical discretization of the soil domain (nodes with upper and lower limits).
• meteo_storm_durations.csv - containing parameters of sub-daily storm/precipitation event patterns for each month

Of these only the CSV files for the parameters soil_horizons.csv and param.csv and the forcings meteoveg.csv, meteoiso.csv are absolutely needed. The remaing can be provided by the user as arguments to loadSPAC() in the Julia script:

• soil_discretizations.csv: not needed if Δz_thickness_m, root_distribution, and IC_soil provided.
• initial_conditions.csv: not needed if IC_scalar provided.
• meteo_storm_durations.csv: not needed if storm_durations_h provided.
• meteoveg.csv: columns with vegetation parameters are not needed if loadSPAC(canopy_evolution = ...) is provided.

The structure of the input CSV's is illustrated by the example input data sets isoBEA2010-18-* or DAV2020-bare-minimum or DAV2020-full located in the folder examples/, as well as below in this documentation. Please follow these examples closely when generating your own input files, including the exact column names and header lines containing the units.

As a special case, below canopy irrigation was implemented, avoiding interception losses. The example input data set DAV2020-irrigation illustrates the required inputs.

For convenience, input CSV files can be generated from a script that sets up a simulation with the R package LWFBrook90R (v0.4.3). Instead of running the simulation with run_LWFB90(), the same arguments can be used to generate the input files for LWFBrook90.jl using the R function provided in the file generate_LWFBrook90jl_Input.R. Note that the input file meteoiso.csv needs to be generated separately and the files containing the initial conditions (initial_conditions.csv and soil_discretization.csv) also need to be extended manually with the isotope values (see structure of these input files below).

To load input data and prepare a simulation follow the instructions in section Example Script 01 or alternatively use the sample script main_with_isotopes.jl. NOTE: these will be replaced with Jupyter-notebooks generated with Literate.jl.

In case you're unfamiliar with Julia, there are various ways to run a script such as main.jl: One possibility is to open the Julia REPL and run the script using include(“main.jl”). Alternatively, the editor VS Code in combination with the Julia extension (julia-vscode.org), provides a complete IDE for programming in Julia.

Structure of input data

soil_horizons.csv: (when using Mualem-VanGenuchten parametrization of the soil retention curve):

meteoveg.csv: contains time dependent parameters (meterologic variables and vegetation parameters):

If meteoveg.csv does not contain columns for DENSEF, HEIGHT, LAI, SAI, the parametrization of these should be provided to the function loadSPAC().

meteoiso.csv: contains time dependent isotopic signatures of precipiation. (Note that the dates contain the end dates of the collection interval of cumulative isotope samples. The values are backward interpolated with piecewise constants. The first row containng the NA is assumed to be the start date of the first collection interval.):

meteo_storm_durations.csv:

initial_conditions.csv:

soil_discretization.csv contains the initial conditions of the soil water status, root density distributions, and the definition of the numerical discretization of the soil domain (nodes with upper and lower limits):

param.csv contains scalar model parameters:

Calibration data (calibration not yet implemented)

Roadmap to include calibration data intends include of:

• Throughfall amounts (for parametrisation of interception)
• Soil moisture (volumetric water content, θ)
• Soil matric potential (ψ)

as well as

• Isotopic composition of soil water (δ_soil)
• Isotopic composition of xylem water (δ_xylem)

Below the example structure of data sets for psi.csv contains the soil matric potential (ψ) of different sensor series or as site average:

theta.csv contains the soil moisture (volumetric water content, θ):

delta_soil.csv contains the isotopic signature of soil water:

delta_xylem.csv contains the isotopic signature of xylem water: