Example Script 01

This example was generated 2024-04-05

Example data from Beatenberg is located in subfolder examples/. WSL is acknowledged for providing the input data (see section Acknowledgments).

To run your own simulations: follow the step-by-step instructions below.

Setup and run simulation

# How to work with a Julia script in a new environment
# create a new folder: "pgm"
# a) start `julia` -> `]` -> activate . -> add LWFBrook90
# b) in VSCode -> click at bottom and select "pgm" as Julia env
# c) run the code
# alternatively right-click on the folder in VSCode
#   -> Change to this Directory and Activate this environment

# Load packages:
using LWFBrook90

# `using LWFBrook90` brings the following exported functions into scope:
# - `loadSPAC()`                   # to define and run a simulation
# - `setup()`                      # to define and run a simulation
# - `simulate!()`                  # to define and run a simulation
# - `remakeSPAC()`                 # to define and run a simulation
# - `SPAC()`                       # to define and run a simulation
# - `DiscretizedSPAC()`            # to define and run a simulation
# - `get_states()`                 # to postprocess a simulation
# - `get_fluxes()`                 # to postprocess a simulation
# - `get_forcing()`                # to postprocess a simulation
# - `get_soil_()`                  # to postprocess a simulation
# - `get_water_partitioning()`     # to postprocess a simulation
# - `RelativeDaysFloat2DateTime()` # a helper function for time
# - `prepare_for_LWFBrook90R()`    # a helper function

Documentation can be accessed by typing: ?loadSPAC, ?setup, ?simulate!, ?plotisotopes, etc.

Define simulation model by reading in system definition and input data from input files. When printed, the generated SPAC model gives a summary.

# Read in input data
input_path = "../../../examples/DAV2020-full/"; input_prefix = "DAV2020-full"; # make sure to get the path correct relative to `pwd()``
# input_path = "examples/DAV2020-full/"; input_prefix = "DAV2020-full";
# input_path = "examples/isoBEAdense2010-18-reset-FALSE-input/"; input_prefix = "isoBEAdense2010-18-reset-FALSE";
model = loadSPAC(input_path, input_prefix; simulate_isotopes = false);
model = loadSPAC(input_path, input_prefix; simulate_isotopes = true);
model

SPAC model:
===== DATES:===============

===== METEO FORCING:===============
GLOBRAD (MJ/m2/day): avg:  13.90, range:   1.04 to  32.40
PREC       (mm/day): avg:   2.61, range:   0.00 to  35.40
TMAX           (°C): avg:   9.00, range:  -8.90 to  26.60
TMIN           (°C): avg:  -0.79, range: -18.80 to  12.90
VAPPRES       (kPa): avg:   0.64, range:   0.11 to   1.37
WIND          (m/s): avg:   2.59, range:   0.70 to   4.90
δ18O            (‰): avg: -13.26, range: -22.60 to  -6.87
δ2H             (‰): avg: -95.98, range:-171.53 to -44.68

===== CANOPY EVOLUTION:===============
model.pars.canopy_evolution was loaded form meteoveg.csv

===== INITIAL CONDITIONS:===============
Soil   IC: soil_discretization.csv
Scalar IC: 3×6 DataFrame
 Row │ u_GWAT_init_mm  u_INTS_init_mm  u_INTR_init_mm  u_SNOW_init_mm  u_CC_in ⋯
     │ Float64         Float64         Float64         Float64         Float64 ⋯
─────┼──────────────────────────────────────────────────────────────────────────
   1 │            1.0             0.0             0.0             0.0          ⋯
   2 │          -13.0           -13.0           -13.0           -13.0
   3 │          -95.0           -95.0           -95.0           -95.0
                                                               2 columns omitted

===== MODEL PARAMETRIZATION:===============
         VXYLEM_mm =>     20.0 |             RHOTP =>      2.0 |                TH =>     40.0 |
   DISPERSIVITY_mm =>     40.0 |                NN =>      2.5 |             MXKPL =>     15.6 |
           LAT_DEG =>     47.0 |          FRINTLAI =>      0.1 |            MXRTLN =>   3000.0 |
        ESLOPE_DEG =>      0.0 |          FSINTLAI =>      0.0 |          INITRLEN =>     12.0 |
        ASPECT_DEG =>      0.0 |          FRINTSAI =>      0.1 |          INITRDEP =>      0.2 |
               ALB =>      0.2 |          FSINTSAI =>      0.0 |          RGRORATE =>      0.0 |
             ALBSN =>      0.5 |            CINTRL =>      0.1 |           RGROPER =>     30.0 |
                C1 =>      0.2 |            CINTRS =>      0.1 |            FXYLEM =>      0.5 |
                C2 =>      0.5 |            CINTSL =>      0.6 |             PSICR =>     -1.0 |
                C3 =>      0.2 |            CINTSS =>      0.6 |             RTRAD =>      0.3 |
            WNDRAT =>      0.3 |            RSTEMP =>     -0.5 |            NOOUTF =>      1.0 |
             FETCH =>   5000.0 |            MELFAC =>      1.5 |          IDEPTH_m =>      0.4 |
               Z0W =>      0.0 |             CCFAC =>      0.3 |          QDEPTH_m =>      0.0 |
                ZW =>      2.0 |            LAIMLT =>      0.2 |              RSSA =>    795.3 |
            MAXLAI =>      3.0 |            SAIMLT =>      0.5 |              RSSB =>      1.0 |
  DENSEF_baseline_ =>      1.0 |            GRDMLT =>      0.3 |            INFEXP =>      0.5 |
     SAI_baseline_ =>      1.0 |            MAXLQF =>      0.1 |             BYPAR =>      0.0 |
  AGE_baseline_yrs =>    100.0 |             KSNVP =>      0.3 |             QFPAR =>      1.0 |
 HEIGHT_baseline_m =>     25.0 |            SNODEN =>      0.3 |              QFFC =>      0.0 |
            LWIDTH =>      0.1 |             GLMAX =>      0.0 |            IMPERV =>      0.0 |
               Z0G =>      0.0 |             GLMIN =>      0.0 |            DSLOPE =>      0.0 |
               Z0S =>      0.0 |                CR =>      0.5 |      LENGTH_SLOPE =>    200.0 |
               LPC =>      4.0 |                RM =>   1000.0 |             DRAIN =>      0.5 |
               CZS =>      0.1 |                R5 =>    287.0 |               GSC =>      0.0 |
               CZR =>      0.1 |              CVPD =>      2.0 |               GSP =>      0.0 |
                HS =>      1.0 |                TL =>      0.0 |                               |
                HR =>     10.0 |                T1 =>     10.0 |                               |
             ZMINH =>      2.0 |                T2 =>     30.0 |                              

===== SOIL DOMAIN:===============
Root distribution:       soil_discretization.csv
Soil layer properties:   3×3 DataFrame
 Row │ Upper_m  Lower_m  shp                                                                               
     │ Float64  Float64  MualemVa…                                                                         
─────┼─────────────────────────────────────────────────────────────────────────────────────────────────────
   1 │    0.0     -0.03  (θ from θr=0.000 to θs=0.379, Ks =  2854.9, STONEF=0.2, l=-3.3, n=1.2, α=   20.4)
   2 │   -0.03    -0.4   (θ from θr=0.000 to θs=0.379, Ks =  2854.9, STONEF=0.4, l=-3.3, n=1.2, α=   20.4)
   3 │   -0.4     -1.1   (θ from θr=0.000 to θs=0.379, Ks =  2854.9, STONEF=0.8, l=-3.3, n=1.2, α=   20.4)
Soil discretized into N=5 layers, Δz layers: (avg, min, max) = (0.220,0.030,0.400)m.
[-0.03, -0.1, -0.4, -0.8, -1.1]

===== SOLVER OPTIONS:===============
(compute_intermediate_quantities = true, simulate_isotopes = true, DTIMAX = 0.5, DSWMAX = 0.05, DPSIMAX = 0.0005)

If wanted, arguments can be supplied to modify the loaded SPAC (e.g. for parameter estimation), hence requiring less input CSVs. Warnings are output (see below) to signal to the user that existing input CSVs might be ignored.

# Read in providing arguments, requiring less CSVs
model_modified =
    loadSPAC(input_path, input_prefix;
        simulate_isotopes = true,
        Δz_thickness_m    = [fill(0.04, 5); # grid spacing (heterogenous), meter (N=21)
                             fill(0.05, 5); # write Δ in VSCode by typing \Delta and hit shift
                             fill(0.06, 5);
                             fill(0.07, 5);
                             fill(0.10, 1)],
        root_distribution = (beta = 0.98, z_rootMax_m=-0.5),
        IC_soil           = (PSIM_init_kPa = -7.0,
                            delta18O_init_permil = -9.0,
                            delta2H_init_permil = -11.0),
        canopy_evolution  = (DENSEF_rel = 100,
                             HEIGHT_rel = 100,
                             SAI_rel    = 100,
                             LAI_rel = (DOY_Bstart = 120,
                                     Bduration  = 21,
                                     DOY_Cstart = 270,
                                     Cduration  = 60,
                                     LAI_perc_BtoC = 100,
                                     LAI_perc_CtoB = 55)),
        storm_durations_h = [5.44, 5.44, 5.44, 5.44, 5.44, 5.44,
                             5.44, 5.44, 5.44, 5.44, 5.44, 5.44],
        IC_scalar         = (amount = (u_GWAT_init_mm = 1.,
                                       u_INTS_init_mm = 13.7,
                                       u_INTR_init_mm = 0.,
                                       u_SNOW_init_mm = 22.222,
                                       u_CC_init_MJ_per_m2 = 0.101010,
                                       u_SNOWLQ_init_mm =  0.),
                            d18O    = (u_GWAT_init_permil = -11.111,
                                       u_INTS_init_permil = -12.222,
                                       u_INTR_init_permil = -13.333,
                                       u_SNOW_init_permil = -14.444),
                            d2H     = (u_GWAT_init_permil = -95.111,
                                       u_INTS_init_permil = -95.222,
                                       u_INTR_init_permil = -95.333,
                                       u_SNOW_init_permil = -95.444)));

┌ Warning: Isotopic signature of precipitation is available for the period (2020-12-16T00:00:00 - 2021-11-04T00:00:00), i.e. (days -16.0 - 307.0)
│          Whereas the other meteorologic inputs are available for (2021-01-01T00:00:00 - 2021-12-31T00:00:00), i.e. (days 0.0 - 364.0).
│          Isotopic signatures will be extrapolated with constant values before and after the available period.
└ @ LWFBrook90 ~/work/LWFBrook90.jl/LWFBrook90.jl/src/func_read_inputData.jl:938
┌ Warning: Requested to overwrite storm durations. Values in ../../../examples/DAV2020-full/DAV2020-full_meteo_storm_durations.csv are ignored.
└ @ LWFBrook90 ~/work/LWFBrook90.jl/LWFBrook90.jl/src/func_read_inputData.jl:687
┌ Warning: Received canopy_evolution in loadSPAC(), overwriting values from `meteoveg.csv`.
└ @ LWFBrook90 ~/work/LWFBrook90.jl/LWFBrook90.jl/src/func_read_inputData.jl:142
┌ Warning: Requested to overwrite initial conditions. Values in ../../../examples/DAV2020-full/DAV2020-full_initial_conditions.csv are ignored.
└ @ LWFBrook90 ~/work/LWFBrook90.jl/LWFBrook90.jl/src/func_read_inputData.jl:172
┌ Warning: loadSPAC(...; Δz_thickness_m = ...) provided. Ignoring `soil_discretization.csv`, i.e. also overwriting root distribution and initial conditions.
└ @ LWFBrook90 ~/work/LWFBrook90.jl/LWFBrook90.jl/src/func_read_inputData.jl:202
┌ Warning: Spatial domain soil discretization is larger than the provided soil horizon information. Lowest soil horizon will be extended.
│ 
│ Lowest soil horizon from -0.4m to -1.1m will be extended down to -1.2000000000000006m.
└ @ LWFBrook90 ~/work/LWFBrook90.jl/LWFBrook90.jl/src/func_discretize_soil_domain.jl:84

Documentation can be accessed by typing: ?loadSPAC, ?setup, ?simulate!, ?plotisotopes, etc.

The objects model and model_modified are instances of a SPAC() (soil-plant-atmosphere continuum) and contain fully specified simulations. To run them they need to be transformed into a system of ODES (setup()) and solved (simulate()):

# Setup and run simulation
simulation     = setup(model)
simulation_mod = setup(model_modified)
# Inputs can be checked with:
# `plot_inputs(simulation)`. Note: NOT YET IMPLEMENTED

simulate!(simulation_mod)
simulate!(simulation_mod) # Run it a second time to showcase shorter runtime
# The computed solution is stored within the object.
# When using VSCode a progress bar is shown in the status bar.

┌ Info:   Start of simulation at 2024-04-05T13:46:40.622.
└         Saving intermediate results (`saveat=`) between: (0.0, 365.0) days
 runtime: 13.541669 seconds (20.97 M allocations: 1.394 GiB, 5.77% gc time, 94.15% compilation time)
[ Info:   Time steps for solving: 9088 (9088 accepted out of 9122 total)
[ Info:   End of simulation at 2024-04-05T13:46:54.297.
┌ Info:   Start of simulation at 2024-04-05T13:46:54.348.
└         Saving intermediate results (`saveat=`) between: (0.0, 365.0) days
 runtime: 0.755181 seconds (3.23 M allocations: 263.580 MiB, 9.08% gc time)
[ Info:   Time steps for solving: 9088 (9088 accepted out of 9122 total)
[ Info:   End of simulation at 2024-04-05T13:46:55.104.

Notice that the very first time a model is simulated, the runtime is longer because it includes compilation. This compilation is not needed for any subsequent run.

Postprocessing results

These simulations can then be post-processed with predefined functions as shown below. Alternatively the following variables contain the simuation results for user-defined post-processing:

simulation_mod.ODESolution;
simulation_mod.ODESolution_datetime;
typeof(simulation_mod.ODESolution);
propertynames(simulation_mod.ODESolution)
# HINT: To get these names, type `simulation_mod.` and wait for the autocomplete!

(:u, :u_analytic, :errors, :t, :k, :prob, :alg, :interp, :dense, :tslocation, :stats, :alg_choice, :retcode)

simulation_mod.ODESolution is a ODESolution object: documentation how to access under: https://docs.sciml.ai/DiffEqDocs/stable/basics/solution/ and more generally: https://docs.sciml.ai/Overview/stable/ simulationmod.ODESolution.t simulationmod.ODESolution.u simulationmod.ODESolution.retcode simulationmod.ODESolution.destats

Note that it is possible to use R Code from within Julia, e.g ggplot: https://stackoverflow.com/a/70073193/3915004

Postprocessing: Extract values

Below illustrative code how to extract states and fluxes:

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

The three functions: get_states(), get_fluxes(), and get_forcing() return DataFrames` with daily values of all the state variables, fluxes, and the input forcing.

get_states(simulation_mod)
get_fluxes(simulation_mod)
get_forcing(simulation_mod)

get_states(simulation_mod; days_to_read_out_d = 0:30:360) # specific days relative to simulation.parametrizedSPAC.reference_date can be specified

Below an illustrative code snippet how to export certain depths into a *.csv:

using CSV, DataFrames

# Single soil variables: How to get θ, or ψ, or δ18Os?
get_soil_(:θ, simulation_mod;
    depths_to_read_out_mm = nothing, days_to_read_out_d = nothing)
depth_to_read_out_mm = [10, 150, 500, 1000, 1150]
get_soil_(:θ, simulation_mod;
    depths_to_read_out_mm = depth_to_read_out_mm, days_to_read_out_d = nothing)
get_soil_(:θ, simulation_mod;
    depths_to_read_out_mm = depth_to_read_out_mm)
get_soil_(:ψ, simulation_mod;
    depths_to_read_out_mm = depth_to_read_out_mm)
get_soil_(:δ18O, simulation_mod;
    depths_to_read_out_mm = depth_to_read_out_mm)

# How to export θ as CSV?
# Only every day, provide days_to_read_out
days_to_read_out = range(simulation_mod.ODESolution.prob.tspan...)
dates_to_read_out = LWFBrook90.RelativeDaysFloat2DateTime.(
    days_to_read_out, simulation_mod.parametrizedSPAC.reference_date)
df_out_daily = get_soil_(:θ, simulation_mod;
    depths_to_read_out_mm = depth_to_read_out_mm, days_to_read_out_d = days_to_read_out)

insertcols!(df_out_daily, 1, :dates => dates_to_read_out);
show(df_out_daily)
##plot(df_out_daily[:,:dates], Matrix(df_out_daily[:,Not([:dates, :time])]))
##CSV.write(
#    joinpath(out_dir, fname * "_θ_depths_daily.csv"),
#    df_out_daily)

366×7 DataFrame
 Row │ dates                time     θ_m3m3_10mm  θ_m3m3_150mm  θ_m3m3_500mm   ⋯
     │ DateTime             Float64  Float64      Float64       Float64        ⋯
─────┼──────────────────────────────────────────────────────────────────────────
   1 │ 2021-01-01T00:00:00      0.0     0.200587      0.200587      0.200587   ⋯
   2 │ 2021-01-02T00:00:00      1.0     0.195189      0.199116      0.203679
   3 │ 2021-01-03T00:00:00      2.0     0.19362       0.197972      0.205324
   4 │ 2021-01-04T00:00:00      3.0     0.192807      0.197298      0.206038
   5 │ 2021-01-05T00:00:00      4.0     0.192393      0.19691       0.20635    ⋯
   6 │ 2021-01-06T00:00:00      5.0     0.192124      0.196648      0.206501
   7 │ 2021-01-07T00:00:00      6.0     0.192024      0.196533      0.206607
   8 │ 2021-01-08T00:00:00      7.0     0.191919      0.196427      0.206691
  ⋮  │          ⋮              ⋮          ⋮            ⋮             ⋮         ⋱
 360 │ 2021-12-26T00:00:00    359.0     0.192669      0.197482      0.210105   ⋯
 361 │ 2021-12-27T00:00:00    360.0     0.19261       0.1974        0.209935
 362 │ 2021-12-28T00:00:00    361.0     0.192623      0.197387      0.209798
 363 │ 2021-12-29T00:00:00    362.0     0.192725      0.19746       0.209726
 364 │ 2021-12-30T00:00:00    363.0     0.280244      0.284455      0.244632   ⋯
 365 │ 2021-12-31T00:00:00    364.0     0.25555       0.267545      0.302739
 366 │ 2022-01-01T00:00:00    365.0     0.245386      0.257208      0.295146
                                                  2 columns and 351 rows omitted

Postprocessing: Plotting (using provided functions)

Below an example script using the provided plot recipes that plot a) amounts, b) isotopes, or c) forcing and states as an additional internal check:

using Plots, Measures; gr();
pl1 = plotamounts(simulation_mod, :above_and_belowground, :showRWUcentroid)
pl1

pl2 = plotisotopes(simulation_mod, :d18O, (d18O = :auto, d2H = :auto), :showRWUcentroid)
pl2

pl3 = plotforcingandstates(simulation_mod)
pl3

# Use of additional arguments to `plotisotopes` is illustrated here.
# Note that typing `?plotisotopes` provides the documentation to the function.
# plotisotopes(simulation_mod, :d2H)
# plotisotopes(simulation_mod, :d18O)
# plotisotopes(simulation_mod;
#              xlim = (DateTime("2010-01-01"), DateTime("2010-03-31")))
# plotisotopes(simulation_mod, :d2H, (d18O = :auto, d2H = :auto), :showRWUcentroid;
#              xlim = (DateTime("2010-01-01"), DateTime("2010-08-30")))

Saving plots can be done with savefig

# using Dates: today;
# out_dir = joinpath("gitignored","out", string(today()))
# mkpath(out_dir)

# fname = input_prefix
# savefig(pl1, joinpath(out_dir, fname*"_plotRecipe_AMT.png"))
# savefig(pl2, joinpath(out_dir, fname*"_plotRecipe_ISO.png"))
# savefig(pl3, joinpath(out_dir, fname*"_plotRecipe_CHECK.png"))

Postprocessing: Plotting (using your own functions)

Below an example script using the manually written code to plot the simulation_mod.#

simstates = getstates(simulation) # this gives an error since it has not been simulated!

sim_states = get_states(simulation_mod)
sim_fluxes = get_fluxes(simulation_mod)

names(sim_states) # show column names
names(sim_fluxes) # show column names

53-element Vector{String}:
 "dates"
 "cum_d_prec"
 "cum_d_sfal"
 "cum_d_sthr"
 "cum_d_sint"
 "cum_d_rfal"
 "cum_d_rint"
 "cum_d_rthr"
 "cum_d_rsno"
 "cum_d_rnet"
 ⋮
 "TRANI_mmday_630mm"
 "TRANI_mmday_690mm"
 "TRANI_mmday_750mm"
 "TRANI_mmday_820mm"
 "TRANI_mmday_890mm"
 "TRANI_mmday_960mm"
 "TRANI_mmday_1030mm"
 "TRANI_mmday_1100mm"
 "TRANI_mmday_1200mm"

plot some of the scalar states

states_to_plot = ["INTS_mm", "INTR_mm", "SNOW_mm", "GWAT_mm", "SWAT_mm"]
sim_states_to_plot = sim_states[:, states_to_plot]
plot(sim_states[:,"dates"], Matrix(sim_states_to_plot),
    labels=permutedims(names(sim_states_to_plot)), legend=:topleft, ylabel = "-")

plot the fates of water partitioning

fluxes_to_plot = ["flow", "seep", "evap"]
sim_fluxes_to_plot = sim_fluxes[:, fluxes_to_plot]
plot(sim_fluxes[:,"dates"], Matrix(sim_fluxes_to_plot),
    labels=permutedims(names(sim_fluxes_to_plot)), legend=:topleft, ylabel = "mm/day")

plot internal fluxes

fluxes_to_plot = [
    "cum_d_prec", "cum_d_rnet", "cum_d_smlt",
    "cum_d_tran",
    "srfl", "slfl", "byfl", "dsfl", "gwfl", "vrfln"]
sim_fluxes_to_plot = sim_fluxes[:, fluxes_to_plot]
plot(sim_fluxes[:,"dates"], Matrix(sim_fluxes_to_plot),
    labels=permutedims(names(sim_fluxes_to_plot)), legend=:topleft, ylabel = "mm/day")

plot isotope signatures of fluxes

fluxes_to_plot = ["RWU_d18O", "RWU_d2H", "PREC_d18O", "PREC_d2H"]
sim_fluxes_to_plot = sim_fluxes[:, fluxes_to_plot]
plot(sim_fluxes[:,"dates"], Matrix(sim_fluxes_to_plot),
    labels=permutedims(names(sim_fluxes_to_plot)), legend=:topleft, ylabel = "mm/day")

Belowground quantities (θ,ψ,δ of soil water)

if you want some specific depths: use getsoil(): e.g depthtoreadoutmm = [150, 500, 800, 1500] dfδsoil = getsoil([:δ18O, :δ2H], simulationmod; depthstoreadoutmm = depthtoreadoutmm) dfθψ = getsoil([:θ, :ψ], simulationmod, depthstoreadoutmm = depthtoreadoutmm)

if we use depths that are already contained in sim_states we can use those:

@show propertynames(sim_states[:,r"(d18O)|(d2H)"]);
df_δsoil = sim_states[:, ["d18O_permil_160mm", "d18O_permil_510mm",
                          "d18O_permil_820mm", "d18O_permil_1200mm",
                          "d2H_permil_160mm", "d2H_permil_510mm",
                          "d2H_permil_820mm", "d2H_permil_1200mm"]]
df_θψ    = sim_states[:, ["θ_m3m3_160mm", "θ_m3m3_510mm",
                          "θ_m3m3_820mm", "θ_m3m3_1200mm",
                          "ψ_kPa_160mm", "ψ_kPa_510mm",
                          "ψ_kPa_820mm", "ψ_kPa_1200mm"]]
pl_θ = plot(sim_states.dates,
    Matrix(select(df_θψ, r"θ_")),
    labels = permutedims(names(select(df_θψ, r"θ_"))),
    xlabel = "Date",
    ylabel = "θ\n[-]",
    legend = :outerright)
pl_ψ = plot(sim_states.dates,
    # Matrix(select(df_θψ, r"ψ_")),
    log10.(-Matrix(select(df_θψ, r"ψ_"))), yflip = true,
    labels = permutedims(names(select(df_θψ, r"ψ_"))),
    xlabel = "Date",
    ylabel = "log10(ψ\n[kPa])",
    legend = :outerright);
pl_δ18O = plot(sim_states.dates,
    Matrix(select(df_δsoil, r"d18O_")),
    labels = permutedims(names(select(df_δsoil, r"d18O_"))),
    xlabel = "Date",
    ylabel = "δ¹⁸O soil\n[‰]",
    legend = :outerright);
pl_δ2H = plot(sim_states.dates,
    Matrix(select(df_δsoil, r"d2H_")),
    labels = permutedims(names(select(df_δsoil, r"d2H_"))),
    xlabel = "Date",
    ylabel = "δ²H soil\n[‰]",
    legend = :outerright);
# add precipitation to soil δ
PREC_color = :black
plot!(pl_δ18O, sim_fluxes.dates, sim_fluxes.PREC_d18O,
    labels = "PREC", color = PREC_color, linestyle = :dot);
plot!(pl_δ2H, sim_fluxes.dates, sim_fluxes.PREC_d2H,
    labels = "PREC", color = PREC_color, linestyle = :dot);


pl_PREC = plot(
    sim_fluxes.dates,
    sim_fluxes.cum_d_prec,
    t = :bar, color=PREC_color,
    legend = :outerright, labels = "PREC    ", # whitespace for hardcoded alignment of legend
    ylabel = "PREC\n[mm]");

plot(plot(pl_PREC, xlab = "", xticks = :none, topmargin = 5mm, bottommargin = 0mm),
    plot(pl_θ;     xlab = "", xticks = :none, topmargin = 0mm, bottommargin = 0mm),
    plot(pl_ψ;     xlab = "", xticks = :none, topmargin = 0mm, bottommargin = 0mm),
    plot(pl_δ18O;  xlab = "", xticks = :none, topmargin = 0mm, bottommargin = 0mm),
    plot(pl_δ2H;   xtick_direction=:out     , topmargin = 0mm, bottommargin = 5mm),
    link = :x,
    layout = grid(5, 1, heights=[0.1, 0.25 ,0.25, 0.2, 0.2]),
    size=(600,500), dpi = 300, margin = 5mm)

####################

This page was generated using Literate.jl.