EDCHM

library(EDCHM)
#> The legacy packages maptools, rgdal, and rgeos, underpinning the sp package,
#> which was just loaded, were retired in October 2023.
#> Please refer to R-spatial evolution reports for details, especially
#> https://r-spatial.org/r/2023/05/15/evolution4.html.
#> It may be desirable to make the sf package available;
#> package maintainers should consider adding sf to Suggests:.
#> Please note that 'maptools' will be retired during October 2023,
#> plan transition at your earliest convenience (see
#> https://r-spatial.org/r/2023/05/15/evolution4.html and earlier blogs
#> for guidance);some functionality will be moved to 'sp'.
#>  Checking rgeos availability: FALSE

Three Ways to Obtain One Model

EDCHM provides three methods for obtaining a hydrological model that simulates the hydrological process.

Use the Pre-Defined Model

EDCHM offers three pre-defined models: (Note: models with the _full suffix will return all variables in a list. Otherwise, the standard model returns only the discharge to speed up the calibration process.)

• EDCHM_GR4J() and EDCHM_GR4J_full(): original GR4J model, same as airGR::RunModel_GR4J().
• EDCHM_mini() and EDCHM_mini_full(): a minimal model structure with random process method.
• EDCHM_snow() and EDCHM_snow_full(): a snow-structure model with random process method.

The function arguments will show all the boundary (input) variable and parameters

Create a Model with a Pre-Defined Model Structure

The build_modell() function offers a standard structure that contains 15 processes and 34 sub-structures, which can generate more than ten millions models.

To build a model, you only need to:

1. Provide the method name for every process:

## select the process method
my_process_method <- c( 
atmosSnow = "atmosSnow_ThresholdT",
evatransPotential = "evatransPotential_TurcWendling",
evatransLand = "evatransActual_SupplyRatio",
evatransSoil = "evatransActual_SupplyRatio",
intercep = "intercep_Full",
snowMelt = "snowMelt_Kustas",
infilt = "infilt_AcceptRatio",
percola = "percola_SupplyRatio",
inteflow = "inteflow_Arno",
capirise = "capirise_HBV",
baseflow = "baseflow_SupplyRatio",
lateral = "lateral_GR4Jfix",
confluenLand = "confluenIUH_Kelly",
confluenSoil = "confluenIUH_Kelly",
confluenGround = "confluenIUH_Nash"
)

2. Create the model with build_modell(), this function will return the rang of all parameters in R-Space and a .cpp-file in a folder, which represents the model.

NOTE: Make sure that there is no model with the same name in the file path, when you give the file path.

## build
my_Model <- build_modell(my_process_method, "TestModel")

## check model code
my_Model$code_Model |> cat(sep = "\n")
#> #include <Rcpp.h>
#> // [[Rcpp::depends(EDCHM)]]
#> #include <EDCHM.h>
#> using namespace Rcpp;
#> using namespace EDCHM;
#> // [[Rcpp::export]]
#> 
#> 
#> 
#> NumericMatrix EDCHM_TestModel(
#> int n_time, 
#> int n_spat,
#> NumericMatrix atmos_netRadiat_MJ, 
#> NumericMatrix atmos_precipitation_mm, 
#> NumericMatrix atmos_solarRadiat_MJ, 
#> NumericMatrix atmos_temperature_Cel, 
#> NumericVector ground_capacity_mm, 
#> NumericVector ground_potentialLateral_mm, 
#> NumericVector ground_water_mm, 
#> NumericVector land_interceptCapacity_mm, 
#> NumericVector land_interceptWater_mm, 
#> NumericVector snow_ice_mm, 
#> NumericVector soil_capacity_mm, 
#> NumericVector soil_potentialCapirise_mm, 
#> NumericVector soil_potentialInteflow_mm, 
#> NumericVector soil_water_mm, 
#> NumericVector confluenGround_responseTime_TS, 
#> NumericVector confluenLand_responseTime_TS, 
#> NumericVector confluenSoil_responseTime_TS, 
#> NumericVector param_atmos_thr_Ts, 
#> NumericVector param_baseflow_sur_k, 
#> NumericVector param_confluenGround_nas_n, 
#> NumericVector param_confluenLand_kel_k, 
#> NumericVector param_confluenSoil_kel_k, 
#> NumericVector param_evatrans_sur_k, 
#> NumericVector param_evatrans_tur_k, 
#> NumericVector param_evatransLand_sur_k, 
#> NumericVector param_infilt_acr_k, 
#> NumericVector param_inteflow_arn_k, 
#> NumericVector param_inteflow_arn_thresh, 
#> NumericVector param_lateral_grf_gamma, 
#> NumericVector param_percola_sur_k, 
#> NumericVector param_snow_kus_fE, 
#> NumericVector param_snow_kus_fT
#> )
#> {
#> 
#> NumericVector atmos_potentialEvatrans_mm, atmos_snow_mm, confluenGround_iuh_1, confluenLand_iuh_1, confluenSoil_iuh_1, ground_capilarise_mm, ground_lateral_mm, land_evatrans_mm, land_intercept_mm, land_water_mm, snow_melt_mm, soil_evatrans_mm, soil_infilt_mm, soil_percolation_mm;
#> 
#> NumericMatrix land_runoff_mm(n_time, n_spat), ground_baseflow_mm(n_time, n_spat), confluen_streamflow_mm(n_time, n_spat), soil_interflow_mm(n_time, n_spat);
#> 
#> for (int i= 0; i < n_time; i++) {
#> 
#> atmos_snow_mm = atmosSnow_ThresholdT(atmos_precipitation_mm(i, _), atmos_temperature_Cel(i, _), param_atmos_thr_Ts);
#> atmos_precipitation_mm(i, _) = atmos_precipitation_mm(i, _) - atmos_snow_mm;
#> 
#> atmos_potentialEvatrans_mm = evatransPotential_TurcWendling(atmos_temperature_Cel(i, _), atmos_solarRadiat_MJ(i, _), param_evatrans_tur_k);
#> 
#> 
#> land_evatrans_mm = evatransActual_SupplyRatio(atmos_potentialEvatrans_mm, land_interceptWater_mm, land_interceptCapacity_mm, param_evatransLand_sur_k);
#> land_interceptWater_mm += - land_evatrans_mm;
#> 
#> soil_evatrans_mm = evatransActual_SupplyRatio(atmos_potentialEvatrans_mm, soil_water_mm, soil_capacity_mm, param_evatrans_sur_k);
#> soil_water_mm += - soil_evatrans_mm;
#> land_water_mm = atmos_precipitation_mm(i, _);
#> 
#> land_intercept_mm = intercep_Full(atmos_precipitation_mm(i, _), land_interceptWater_mm, land_interceptCapacity_mm);
#> land_interceptWater_mm += land_intercept_mm;
#> land_water_mm = atmos_precipitation_mm(i, _) - land_intercept_mm;
#> 
#> snow_melt_mm = snowMelt_Kustas(snow_ice_mm, atmos_temperature_Cel(i, _), atmos_netRadiat_MJ(i, _), param_snow_kus_fE, param_snow_kus_fT);
#> land_water_mm += snow_melt_mm;
#> snow_ice_mm += -snow_melt_mm;
#> snow_ice_mm += atmos_snow_mm;
#> 
#> soil_infilt_mm = infilt_AcceptRatio(land_water_mm, soil_water_mm, soil_capacity_mm, param_infilt_acr_k);
#> soil_water_mm += soil_infilt_mm;
#> land_runoff_mm(i, _) = land_water_mm - soil_infilt_mm;
#> 
#> soil_percolation_mm = percola_SupplyRatio(soil_water_mm, param_percola_sur_k);
#> ground_water_mm += soil_percolation_mm;
#> soil_water_mm += - soil_percolation_mm;
#> 
#> soil_interflow_mm(i, _) = inteflow_Arno(soil_water_mm, soil_capacity_mm, soil_potentialInteflow_mm, param_inteflow_arn_thresh, param_inteflow_arn_k);
#> soil_water_mm += - soil_interflow_mm(i, _);
#> 
#> ground_capilarise_mm = capirise_HBV(ground_water_mm, soil_water_mm, soil_capacity_mm, soil_potentialCapirise_mm);
#> ground_water_mm += - ground_capilarise_mm;
#> soil_water_mm += ground_capilarise_mm;
#> 
#> NumericVector baseflow_temp = ifelse(ground_water_mm < ground_capacity_mm, 0, ground_water_mm - ground_capacity_mm);
#> ground_water_mm = ifelse(ground_water_mm < ground_capacity_mm,ground_water_mm, ground_capacity_mm);
#> 
#> ground_baseflow_mm(i, _) = baseflow_SupplyRatio(ground_water_mm, param_baseflow_sur_k);
#> ground_water_mm += - ground_baseflow_mm(i, _);
#> ground_baseflow_mm(i, _) = ground_baseflow_mm(i, _) + baseflow_temp;
#> 
#> ground_lateral_mm = lateral_GR4Jfix(ground_water_mm, ground_capacity_mm, ground_potentialLateral_mm, param_lateral_grf_gamma);
#> ground_water_mm += ground_lateral_mm;
#> 
#> }
#> for (int j= 0; j < n_spat; j++) {
#> confluenLand_iuh_1 = confluenIUH_Kelly(confluenLand_responseTime_TS(j), param_confluenLand_kel_k(j));
#> confluenSoil_iuh_1 = confluenIUH_Kelly(confluenSoil_responseTime_TS(j), param_confluenSoil_kel_k(j));
#> confluenGround_iuh_1 = confluenIUH_Nash(confluenGround_responseTime_TS(j), param_confluenGround_nas_n(j));
#> 
#> confluen_streamflow_mm(_, j) = confluen_IUH3S(land_runoff_mm(_, j), soil_interflow_mm(_, j), ground_baseflow_mm(_, j), confluenLand_iuh_1, confluenSoil_iuh_1, confluenGround_iuh_1);
#> }
#> return confluen_streamflow_mm;
#> }

## check parameter range
my_Model$range_Parameter
#>                               min   max
#> param_atmos_thr_Ts         -1e+00 3.000
#> param_baseflow_sur_k        1e-02 1.000
#> param_confluenGround_nas_n  1e+00 8.000
#> param_confluenLand_kel_k    1e+00 4.000
#> param_confluenSoil_kel_k    1e+00 4.000
#> param_evatrans_sur_k        1e-01 1.000
#> param_evatrans_tur_k        6e-01 1.000
#> param_evatransLand_sur_k    1e-01 1.000
#> param_infilt_acr_k          1e-02 1.000
#> param_inteflow_arn_k        1e-01 1.000
#> param_inteflow_arn_thresh   1e-01 0.900
#> param_lateral_grf_gamma     1e-02 5.000
#> param_percola_sur_k         1e-02 1.000
#> param_snow_kus_fE           5e-04 0.003
#> param_snow_kus_fT           5e-02 1.000

The .cpp-file looks like the source code of EDCHM_mini, which will be shown in the next section on how to use it.

Manually build a model with EDCHM modules

Although build_modell() offers the possibility of generating over a million models, there may still be a requirement to create a model with a different structure. For example, it may involve adding or deleting some processes, or adding more layers.

As shown in the model source code of EDCHM_mini, it is straightforward to assemble a model with EDCHM modules. The steps are as follows:

1. Select the deciding processes and their methods, the Reference-Page shows all of the processes and process-methods.

2. Fill the arguments of the module (function) with specific variables and parameters. The original argument names are well-defined, so they can be used as variable names.

The function with arguments looks like:

infilt_UBC(
  land_water_mm,
  land_impermeableFrac_1,
  soil_water_mm,
  soil_capacity_mm,
  param_infilt_ubc_P0AGEN
)

The code line in the model source code looks like:

soil_infilt_mm = infilt_UBC(land_water_mm, land_impermeableFrac_1, soil_water_mm, soil_capacity_mm, param_infilt_ubc_P0AGEN);

In the most Situations the arguments names can be used directly.

3. Set the processing after the module, the relationship between the calculated process-variable and the standard variables.

soil_infilt_mm = infilt_UBC(land_water_mm, land_impermeableFrac_1, soil_water_mm, soil_capacity_mm, param_infilt_ubc_P0AGEN);
soil_water_mm += soil_infilt_mm;
land_runoff_mm(i, _) = land_water_mm - soil_infilt_mm;

The infiltration water will infiltrate from the land-layer to the soil-layer, while the remaining water will flow as runoff to the river or catchment outlet. So these two lines will be added in to the model.

4. Complete the boundary conditions (BC) variables (input, meteorological data) in matrix-form. The BC variables are always given as time series, so they must be in matrix form, for example, atmos_precipitation_mm and atmos_potentialEvatrans_mm.

soil_evatrans_mm = evatransActual_UBC(atmos_potentialEvatrans_mm(i, _), soil_water_mm, soil_capacity_mm, param_evatrans_ubc_gamma);
soil_water_mm += - soil_evatrans_mm;
land_water_mm = atmos_precipitation_mm(i, _);

Both variables is defined as matrices in input, so the time dimension (index i) must be specified.

5. Complete the outflow variables in matrix-form. The outflow variables refer to the variables that will flow to the catchment outlet through “routing”, for example, land_runoff_mm and ground_baseflow_mm.

land_runoff_mm(i, _) = land_water_mm - soil_infilt_mm;

ground_baseflow_mm(i, _) = ground_baseflow_mm(i, _) + baseflow_temp;

Both variables must be declared as matrices, so the time dimension (index i) must be specified.

6. Define all the parameters and input (BC) variables as arguments and declare all the other variables.

NumericMatrix EDCHM_mini(
// basci parameter
int n_time, 
int n_spat,
// BC data (Input)
NumericMatrix atmos_potentialEvatrans_mm, 
NumericMatrix atmos_precipitation_mm,
// parameters, that schould be estimated and innitial condition (IC) data
NumericVector ground_capacity_mm, 
NumericVector ground_water_mm, 
NumericVector land_impermeableFrac_1, 
NumericVector soil_capacity_mm, 
NumericVector soil_potentialPercola_mm, 
NumericVector soil_water_mm, 
// parameters, that need to calibrated
NumericVector confluenLand_responseTime_TS, 
NumericVector confluenGround_responseTime_TS, 
NumericVector param_baseflow_grf_gamma, 
NumericVector param_confluenLand_kel_k, 
NumericVector param_evatrans_ubc_gamma, 
NumericVector param_infilt_ubc_P0AGEN, 
NumericVector param_percola_arn_k, 
NumericVector param_percola_arn_thresh
)
{

// declare the vector-variable, that contain only on dimension time like confluenLand_iuh_1 and confluenGround_iuh_1
// or spat like land_water_mm, soil_evatrans_mm, soil_infilt_mm and soil_percolation_mm
NumericVector land_water_mm, soil_evatrans_mm, soil_infilt_mm, soil_percolation_mm, confluenLand_iuh_1, confluenGround_iuh_1;
// declare the matric-variable, that contain two dimensions time and spat
NumericMatrix land_runoff_mm(n_time, n_spat), ground_baseflow_mm(n_time, n_spat), confluen_streamflow_mm(n_time, n_spat);

return confluen_streamflow_mm;
}

With these six steps, we can build a new model, but it is important to ensure that the logic of relationships between variables is correct.

Run the model

Run the pre-defined model from EDCHM

The pre-defined model can be used like normal functions. You just need to provide the input data and the parameters in the correct format as specified in the function arguments. For example:

## reform the data into matrix
mat_EDCHM_TestData <- EDCHM_TestData |> as.matrix()

## Parameter initial
Param_mini <- c(land_impermeableFrac_1 = .5, 
                     soil_potentialPercola_mm = 5, 
                     confluenLand_responseTime_TS = 2, 
                     confluenGround_responseTime_TS = 4, 
                     Param_baseflow_grf_gamma = 3, 
                     Param_confluenLand_kel_k = 2, 
                     Param_evatrans_ubc_gamma = 1.5, 
                     Param_infilt_ubc_P0AGEN = 2, 
                     Param_percola_arn_k = .6, 
                     Param_percola_arn_thresh = .7)

## run the model
test_Q <- EDCHM_mini(
    180, # n_time, number of time step
    1, # n_spat, number of spatial unit
    mat_EDCHM_TestData[,"EC"] |> as.matrix(), # atmos_potentialEvatrans_mm, Input matrix
    mat_EDCHM_TestData[,"PA2"] |> as.matrix(), # atmos_precipitation_mm, Input matrix
    200., # ground_capacity_mm, (physical based) parameter
    100., # ground_water_mm, initial condition
    Param_mini["land_impermeableFrac_1"], 
    200., # soil_capacity_mm, (physical based) parameter
    Param_mini["soil_potentialPercola_mm"], 
    100., # soil_water_mm, initial condition
    Param_mini["confluenLand_responseTime_TS"], 
    Param_mini["confluenGround_responseTime_TS"], 
    Param_mini["Param_baseflow_grf_gamma"], 
    Param_mini["Param_confluenLand_kel_k"], 
    Param_mini["Param_evatrans_ubc_gamma"], 
    Param_mini["Param_infilt_ubc_P0AGEN"], 
    Param_mini["Param_percola_arn_k"], 
    Param_mini["Param_percola_arn_thresh"])

The input (boundary condition) data, including precipitation (P) and potential evapotranspiration (PET), as well as the simulated model result discharge (contribution in mm/TS, equivalent to precipitation) (Q), are shown below:

Run the model with cpp source code

If the model is still in CPP source code, you first need to compile it into a function using the Rcpp::sourceCpp() function. For example:

Rcpp::sourceCpp("path/to/your/model.cpp")

Once the model is compiled, you can use it as a normal function by providing the input data and parameters in the correct format.