Multi-run simulations

[yjaguilara]

Hi, I have been using Brook90.jl to run individual simulations and perform hydrologic analysis of soil water moisture and groundwater storage.
I would like to be able to run several simulations varying the values of some parameters in order to perform sensitivity analysis. Is there any method or function that allows me to run multiple simulations, only varying the values of certain parameters within certain established ranges?

Since I need to do more than 1000 simulations, it is tedious to have to prepare the csv for each of the simulations and run them individually; therefore any guide or help to be able to do multirunning would be super useful.

Thank you in advance for your response and help.

[fabern]

Hi Juliana,

Yes there is. You can do this with remakeSPAC() (documentation under: https://fabern.github.io/LWFBrook90.jl/stable/function-docs/#LWFBrook90.remakeSPAC-Tuple%7BDiscretizedSPAC%7D)

Sorry, for my brevity. I currently don‘t have access to a computer. So I can‘t provide you an example script. Let me know if this is enough to get you started…
The function call depends on which parameters you want to vary. Do you already know which parameters these are?

Best regards, Fabian

[yjaguilara]

Hello Fabian,
Thank you for your reply.
Unfortunately I have not yet managed to run multiple simulations to perform my sensitivity analysis.
I want to vary 41 parameters, and I already have the ranges for each of them.
8 of these 41, are parameters that I want to vary by soil layer (ths, thr, alpha, npar, ksat, tort, gravel).
Given the number of parameters I want to analyze, I would have to run at least 20,000 simulations.

Would it be possible to do this?
I would greatly appreciate any help you can give me.
Thank you in advance for your assistance.

[fabern]

Hi Juliana,

Yes, in general this is possible. It might be more or less difficult depending on the ranges you want to test, how many soil layers you have, and if you vary soil parameters freely.

Workflow: The function remakeSPAC is intended to modify parameters from an existing simulation to be used for e.g. parameter estimation. You can find some example in the code that @baaponter posted in this discussion. The overall workflow with parameter estimation would then be something along:

• define a base simulation
• generate your 20'000+ parameter sets
• loop over the parameter sets and
  • generate the modified simulation with remakeSPAC
  • run the simulation
  • evaluate a goodness of fit criteria e.g. against observations
  • store goodness of fit to some list (but probably not the full simulation results since for 20'000 runs this is large)
• select the best or a number of best simulations based on the goodness of fit.
• (potentially rerun these best simulations and storing their results)

Depending on the computer you have available it can be possible to parallelize the loop and/or you might also need to reduce how much data is saved/kept beyond the loop iteration due to memory (RAM) constraints.

Possible to do this? It depends:
The number of parameters is large and you need to think about how you sample this large, multidimensional space. Finding global optimums in high-dimensional parameter spaces is a complex, problem-specific task.

Anything that can reduce the complexity could be helpful. E.g. smaller ranges or reducing the number of parameters.

How many soil layers do you have?
Do you have any indication on the properties of the soil layers?
These parameters might be tricky: on one hand because they are many (8xNLayer) and because they can make the simulation much more computationally intensive:

When you generate them freely without any constraints, you might sample large differences between adjacent layers (that may be not reasonable representation of reality). These large gradients will most certainly make the numerical solution much more difficult. So you would have certain simulations that - due to representation of the soil parameters - might be too slow to finish your whole simulation period.
The question is then, if it is possible to exclude these kind of simulations beforehand. E.g. by defining how you sample the space of soil profiles that are explored with the 20'000 parameter sets. For my applications I could rely on data from dug soil profiles and simply fitted a multiplicative factor, but keeping their relative values constant.

[yjaguilara]

Hello Fabian,
thank you very much for your answer and help.
I have 4 layers of soil and I have already managed to run 20000 simulations varying most of the parameters I need, using the remakeSPAC() function. However I still have problems with a couple of parameters.

Context: My goal is to perform sensitivity analysis of my model using SAFEToolbox, but this tool is not available in Julia. For this reason what I do is to generate my X (the 20000 samples of each of my parameters) in R using the sampling methods required by SAFEToolbox and then import that X into Julia to run my model (I was previously using LWFBrook90.R for simplicity of the same language, however it is important for me to analyze my isotopic data so I really want to run with LWFBrook90.jl).

Currently I have not been able to vary the Rootden_, emergedur and winlaifrac parameters.

• Regarding Rootden_, I would like to be able to vary this parameter per layer using remakeSPAC(), but defining it directly but not calculating it from beta and z_rootMax_m. In other words, changing the values of Rootden_ for the 4 layers that I initially defined in the file soil_discretization.csv.

In R, I have managed to vary it by using the rootmethod = soilvar, that can be used if the root density depth distribution is specified in column 'rootden' in the soil-data.frame.
Is there any way to do something similar in my model in Julia, either by changing the method to calculate rootden_ or by defining it inside root_distribution, when using remakeSPAC()?

• Regarding emergedur (Leaf growth duration until maxlai is reached) and winlaifrac (Minimum LAI as a fraction of maxlai): Both are parameters defined in the LWFBrook90.R model, but I have not been able to find a way to vary them in my Julia model.
  Is there any way I can vary these 2 parameters, using remakeSPAC() and defining them inside params?

Ideally it would be something like this: See example code below

I greatly appreciate your help and your prompt response.
Thank you very much.

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

simulation = setup(model)

sim2 = remakeSPAC(simulation; 
    params = (
        emergedur = 40,
        winlaifrac = 0.325,
    ),
    root_distribution = (
        Rootden_ = [1.79,0.73,0.09,0.01],
    ),
    
)
simulate!(sim2, save_everystep = false, saveat = 1)

[baaponter]

Hello.
I hope it's ok to post here as I'm also working with multi-simulation.
Giving some context, I am using the remakeSPAC function to modify some parameters and most work fine, but I am having problems specifically with root_distribution and IC_soil:

A basic sample of the code:

model = loadSPAC(input_path, input_prefix; simulate_isotopes = true)
simulation = setup(model)
simulation_var = remakeSPAC(simulation; 
    params = (
        RSSA = 500,
        GLMAX = 0.0048,
    ),
    soil_horizons = (
        ths_ = [0.41,0.28,0.43,0.39],
        thr_ = [0.11,0.04,0.17,0.10]
    ),
    IC_soil = (
        PSIM_init_kPa = -6.1, 
    ),
    root_distribution = (
        beta = 0.91,
        z_rootMax_m = 0.0043, 
    ),
    
)
simulate!(simulation_var, save_everystep = false, saveat = 1, maxiters=1e10)

This generates the two errors
In IC_soil:

On root_distribution:

I have found that to modify root_distribution and IC_soil from remakeSPAC you must initially load from the model like this:

model = loadSPAC(input_path, input_prefix; simulate_isotopes = true,
    root_distribution = (beta = 0.95, z_rootMax_m = -1,),
    IC_soil = (PSIM_init_kPa = -6.3, delta18O_init_permil = -8.3, delta2H_init_permil = -58, ),
)

simulation = setup(model)

simulation_var = remakeSPAC(simulation; 
    params = (
        RSSA = 500,
        GLMAX = 0.0048,
    ),
    soil_horizons = (
        ths_ = [0.41,0.28,0.43,0.39],
        thr_ = [0.11,0.04,0.17,0.10]
    ),
    IC_soil = (
        PSIM_init_kPa = -6.1, 
    ),
    root_distribution = (
        beta = 0.91,
        z_rootMax_m = -0.0043, 
    ),
    
)
simulate!(simulation_var, save_everystep = false, saveat = 1, maxiters=1e10)

This no longer generates any error, but I see that when performing this action the soil_discretization.csv file is no longer taken into account since the Rootden parameter per layer is one of which I also want to modify (although I have not yet found how to do it) takes values very different from those in the file.

I would like to know if it is possible to adjust the root_distribution and IC_soil without having to define them in the model or if necessary what is the correct way to respect the other data in the file.
Also if there is the possibility to adjust Rootden by layer from remakeSPAC or from some other place in the code.
Thank you very much.

[fabern]

Hi baaponter,

Thanks for joining into the discussion and providing your sample code!
Yes, as you realized, the model needs to be defined differently to modify the root_distribution and initial conditions.
It needs to be given root distributions and initial conditions parametrically, and in that case ignores a potential soil_discretizations.csv input file.

You're issue is mainly related to what has (so far) been implemented as interfaces for the remakeSPAC()- function. Since my use cases involved estimating parameters of the root distribution, I only needed the remakeSPAC function to support the root distributions parametrized with the beta model. Other use cases can be added, I just need to ensure the code remains legible and maintainable and that there is a certain backwards compatibility.

So, if I understand correctly, you would like to have the possibility to use remakeSPAC to adjust the root_distribution and IC_soil by providing vectors of values for each soil horizon? So basically similar to replacing the file soil_discretization.csv but all from withing the script?

Let me know if this should be added to the remakeSPAC()-function. Since it can be ambiguous if one means soil horizons (as defined in soil_horizons.csv) or discretized layers (as defined in soil_discretization.csv), the most robust way to accommodate both would be to provide a piecewise continuous function of these values as a function of the depth z.

[baaponter]

Hello Fabian,
Thank you for your reply.

Indeed I want to do multi-simulations, varying the following parameters:

• Parameters contained in the soil_discretization.csv file:
  (With these I have problem at the moment)
  Rootden_ (Vary it by soil layer, assigning a vector with 4 specific values, considering a 4 layer model; but not calculating it from beta and maxrootdepth)
  uAux_PSIM_init_kPa (varying it in general, so that it is the same value for all 4 layers)

• Parameters contained in the soil_horizonts.csv file: thr, ths, alpha, npar, Ksat, tort, gravel
  (Currently remakeSPAC allows me to change these parameters without any problem, assigning vectors of 4 values, considering a 4 layers model).

The ideal would be to be able to modify the parameters (Rootden_ and uAux_PSIM_init_kPa) with the remakeSPAC function.
Something like that:

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

simulation = setup(model)

simulation_var = remakeSPAC(simulation; 
    params = (
        RSSA = 500,
        GLMAX = 0.0048,
    ),
    soil_horizons = (
        ths_ = [0.41,0.28,0.43,0.39],
        thr_ = [0.11,0.04,0.17,0.10]
    ),
    IC_soil = (
        PSIM_init_kPa = -6.1, 
    ),
    root_distribution = (
        Rootden_ = [3.19,1.29,0.35,0.032], ####### would it be possible ?
    ),
    
)
simulate!(simulation_var, save_everystep = false, saveat = 1, maxiters=1e10)

Or is there any other option to modify these 2 parameters, without having to modify the soil_discretization.csv file for each simulation?

Thanks in advance for your help

[fabern]

Hi baaponter,

Thanks for clarifying. No, at the moment there is no other option to modify this.
I will provide an update to include this use case. I'll try to have it ready by the end of this week...

[fabern]

Hi @yjaguilara and @baaponter,

Sorry that it takes a while. The input by the two of you requires some changes to the internal structure how and in what order LWFBrook90.jl processes the user input. This had wider ramifications than I initially thought. So it will take some more time. The progress of this major API rewrite is tracked by PR #92.

Thanks for having a little more patience.
Cheers, Fabian