Example Script 01

This example was generated 2025-02-15

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 2025-02-15T01:11:25.964.
└         Saving intermediate results (`saveat=`) between: (0.0, 365.0) days
 runtime: 13.772222 seconds (21.14 M allocations: 1.400 GiB, 5.57% gc time, 94.51% compilation time)
[ Info:   Time steps for solving: 9088 (9088 accepted out of 9122 total)
[ Info:   End of simulation at 2025-02-15T01:11:39.872.
┌ Info:   Start of simulation at 2025-02-15T01:11:39.924.
└         Saving intermediate results (`saveat=`) between: (0.0, 365.0) days
 runtime: 0.749511 seconds (3.23 M allocations: 263.581 MiB, 7.70% gc time)
[ Info:   Time steps for solving: 9088 (9088 accepted out of 9122 total)
[ Info:   End of simulation at 2025-02-15T01:11:40.673.

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)

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

13×102 DataFrame
 Row │ dates                GWAT_mm  INTS_mm   INTR_mm   SNOW_mm   CC_MJm2     ⋯
     │ DateTime             Float64  Float64   Float64   Float64   Float64     ⋯
─────┼──────────────────────────────────────────────────────────────────────────
   1 │ 2021-01-01T00:00:00      1.0  13.7      0.0        22.222   0.10101     ⋯
   2 │ 2021-01-31T00:00:00      1.0   1.09368  0.392421  158.295   0.0
   3 │ 2021-03-02T00:00:00      1.0   0.0      0.0       104.096   0.0
   4 │ 2021-04-01T00:00:00      1.0   0.0      0.0        71.5081  0.0
   5 │ 2021-05-01T00:00:00      1.0   0.0      0.0         0.0     0.0         ⋯
   6 │ 2021-05-31T00:00:00      1.0   0.0      0.0         0.0     0.0
   7 │ 2021-06-30T00:00:00      1.0   0.0      0.0         0.0     0.0
   8 │ 2021-07-30T00:00:00      1.0   0.0      0.0         0.0     0.0
   9 │ 2021-08-29T00:00:00      1.0   0.0      0.0         0.0     0.0         ⋯
  10 │ 2021-09-28T00:00:00      1.0   0.0      0.0         0.0     0.0
  11 │ 2021-10-28T00:00:00      1.0   0.0      0.0         0.0     0.0
  12 │ 2021-11-27T00:00:00      1.0   0.0      0.0         1.2592  0.00689034
  13 │ 2021-12-27T00:00:00      1.0   0.0      0.0        24.8801  0.0         ⋯
                                                              96 columns omitted

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();
pl2 = plotisotopes(simulation_mod, :d18O, (d18O = :auto, d2H = :auto), :showRWUcentroid)
pl2

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

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.#

# sim_states = get_states(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