Skip to content

Commit

Permalink
bug_fix_aet rootzone (#105)
Browse files Browse the repository at this point in the history 
* fixed bugs in the AET root zone related to flag option and memory.

* minor refactoring of the code, names consistency etc.
  • Loading branch information
@ajkhattak
ajkhattak authored Feb 16, 2024
1 parent 041c11c commit 44bb30c
Show file tree
Hide file tree
Showing 5 changed files with 62 additions and 57 deletions.
14 changes: 7 additions & 7 deletions configs/README.md
View file
Original file line number Diff line number Diff line change
Expand Up @@ -30,10 +30,10 @@ Example configuration files are provided in this directory. To build and run the
| *b_Xinanjiang_shape_parameter=1 | *double* | | | parameter_adjustable | direct runoff | when `surface_partitioning_scheme=Xinanjiang` |
| *x_Xinanjiang_shape_parameter=1 | *double* | | | parameter_adjustable | direct runoff | when `surface_partitioning_scheme=Xinanjiang` |
| urban_decimal_fraction | *double* | 0.0 - 1.0 | | parameter_adjustable | direct runoff | when `surface_partitioning_scheme=Xinanjiang` |
| aet_rootzone | *boolean* | True, true or 1 | | coupling parameter | `rootzone-based AET` | when `CFE coupled to SoilMoistureProfile` |
| max_root_zone_layer | *double* | | meters [m] | parameter_adjustable | AET | layer of the soil that is the maximum root zone depth. That is, the depth of the layer where the AET is drawn from |
| is_aet_rootzone | *boolean* | True, true or 1 | | coupling parameter | `rootzone-based AET` | when `CFE coupled to SoilMoistureProfile` |
| max_rootzone_layer | *double* | | meters [m] | parameter_adjustable | AET | layer of the soil that is the maximum root zone depth. That is, the depth of the layer where the AET is drawn from |
| soil_layer_depths | 1D array | | meters [m] | parameter_adjustable | AET | an array of depths from the surface. Example, soil_layer_depths=0.1,0.4,1.0,2.0
| sft_coupled | *boolean* | True, true or 1 | | coupling parameter | `ice-fraction based runoff` | when `CFE coupled to SoilFreezeThaw`|
| is_sft_coupled | *boolean* | True, true or 1 | | coupling parameter | `ice_fraction-based runoff` | when `CFE coupled to SoilFreezeThaw`|


## Direct runoff options in CFE
Expand All @@ -51,18 +51,18 @@ If the **Xinanjiang** scheme is choosen, four parameters need to be included in

## Rootzone-based Actual Evapotranspiration (AET)
The user has the option to turn ON and OFF rootzone-based AET, default option is OFF. To turn it ON, the following parameters need to be included in the configuration file.
1. `aet_rootzone=true`
1. `is_aet_rootzone=true`
2. `soil_layer_depths`
3. `max_root_zone_layer`
3. `max_rootzone_layer`

## CFE coupled to Soil freeze-thaw model (SFT)
The Soil Freeze-Thaw (SFT) model is a standalone model. For detailed information please refer to the [SFT repo](https://github.com/NOAA-OWP/SoilFreezeThaw). A few things to note when coupling CFE to SFT:
1. SFT model provides `ice fraction` to CFE runoff schemes (Schaake `ice_fraction_schaake` and Xinanjiang `ice_fraction_xinanjiang`)
2. To turn ON/OFF SFT set sft_coupled flag.
* `sft_coupled` : (type boolean) if `true`, SFT is turned ON. (options: True, true, 1).
* `is_sft_coupled` : (type boolean) if `true`, SFT is turned ON. (options: True, true, 1).
* If the runoff scheme is Xinanjiang, no additional parameters are needed in the CFE config files.
* If the runoff scheme is Schaake, the CFE config file will need an additional parameter, namely:
* `ice_content_threshold` : (type double, unit m). This represent the ice content above which soil is impermeable.


**Note:** By defualt `sft_coupled` and `aet_rootzone` are `OFF`, that means these changes do not effect the basic functionality of CFE.
**Note:** By defualt `is_sft_coupled` and `is_aet_rootzone` are set to `OFF`, that means these changes do not affect the basic functionality of CFE.
2 changes: 1 addition & 1 deletion configs/laramie_bmi_config_cfe_pass_aet_rz.txt
View file
Original file line number Diff line number Diff line change
Expand Up @@ -28,4 +28,4 @@ surface_partitioning_scheme=Schaake
#urban_decimal_fraction=0.0
aet_rootzone=true
soil_layer_depths=0.1,0.4,1.0,2.0
max_root_zone_layer=2
max_rootzone_layer=2
31 changes: 16 additions & 15 deletions include/conceptual_reservoir.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -15,27 +15,28 @@ struct conceptual_reservoir {
// iff is_exponential==TRUE, then it uses the exponential discharge function from the NWM V2.0 forumulation
// as the primary discharge with a zero threshold, and does not calculate a secondary discharge.
//--------------------------------------------------------------------------------------------------
int is_exponential; // set this true TRUE to use the exponential form of the discharge equation
double gw_storage; // Initial Storage - LKC: added since I need to keep track of it when changing parameters
double storage_max_m; // maximum storage in this reservoir
double storage_m; // state variable.
double storage_change_m; // storage change in the current step
double coeff_primary; // the primary outlet
int is_exponential; // set this true TRUE to use the exponential form of the discharge equation
double gw_storage; // Initial Storage - LKC: added since I need to keep track of it when changing parameters
double storage_max_m; // maximum storage in this reservoir
double storage_m; // state variable.
double storage_change_m; // storage change in the current step
double coeff_primary; // the primary outlet
double exponent_primary;
double storage_threshold_primary_m;
double storage_threshold_secondary_m;
double coeff_secondary;
double exponent_secondary;
double ice_fraction_schaake, ice_fraction_xinanjiang;
int is_sft_coupled; // boolean - true if SFT is ON otherwise OFF (default is OFF)
double ice_fraction_schaake;
double ice_fraction_xinanjiang;
int is_sft_coupled; // boolean - true if SFT is ON otherwise OFF (default is OFF)

//---Root zone adjusted AET development -rlm -ahmad -------------
double *smc_profile; //soil moisture content profile
int n_soil_layers; // number of soil layers
double *soil_layer_depths_m; // soil layer depths defined in the config file in units of [m]
int aet_root_zone; // boolean - true if aet_root_zone is ON otherwise OFF (default is OFF)
int max_root_zone_layer; // maximum root zone layer is used to identify the maximum layer to remove water from for AET
double *delta_soil_layer_depth_m; // used to calculate the total soil moisture in each layer
//---Root zone adjusted AET development -rlm -ajk -------------
double *smc_profile; //soil moisture content profile
int n_soil_layers; // number of soil layers
double *soil_layer_depths_m; // soil layer depths defined in the config file in units of [m]
int is_aet_rootzone; // boolean - true if aet_root_zone is ON otherwise OFF (default is OFF)
int max_rootzone_layer; // maximum root zone layer is used to identify the maximum layer to remove water from for AET
double *delta_soil_layer_depth_m; // used to calculate the total soil moisture in each layer
double soil_water_content_field_capacity; // water content [m/m] at field capacity. Used in AET routine

//---------------------------------------------------------------
Expand Down
52 changes: 28 additions & 24 deletions src/bmi_cfe.c
View file
Original file line number Diff line number Diff line change
Expand Up @@ -170,7 +170,7 @@ Variable var_info[] = {
// -------------------------------------------
{ 91, "soil_moisture_profile", "double", 1},
{ 92, "soil_layer_depths_m", "double", 1},
{ 93, "max_root_zone_layer", "int", 1},
{ 93, "max_rootzone_layer", "int", 1},
//--------------------------------------------
};

Expand Down Expand Up @@ -493,7 +493,7 @@ int read_init_config_cfe(const char* config_file, cfe_state_struct* model)
/* ------ Root zone adjusted AET development -rlm ------- */
int is_aet_rootzone_set = FALSE;
int is_soil_layer_depths_string_val_set = FALSE;
int is_max_root_zone_layer_set = FALSE;
int is_max_rootzone_layer_set = FALSE;
/*--------------------------------------------------------*/
// Default value
model->NWM_soil_params.refkdt = 3.0;
Expand Down Expand Up @@ -732,13 +732,14 @@ int read_init_config_cfe(const char* config_file, cfe_state_struct* model)

/*-------------------- Root zone AET development -rlm -----------------------*/
if (strcmp(param_key, "aet_rootzone") == 0) {
model->soil_reservoir.aet_root_zone = strtod(param_value, NULL);

if ( strcmp(param_value, "true")==0 || strcmp(param_value, "True")==0 || strcmp(param_value,"1")==0)
is_aet_rootzone_set = TRUE;

continue;
}

if (is_aet_rootzone_set == TRUE ) {
if (is_aet_rootzone_set == TRUE) {
if (strcmp(param_key, "soil_layer_depths") == 0) {
#if CFE_DEBUG >= 1
printf("Found configured soil depth values ('%s')\n", param_value);
Expand All @@ -747,9 +748,9 @@ int read_init_config_cfe(const char* config_file, cfe_state_struct* model)
is_soil_layer_depths_string_val_set = TRUE;
continue;
}
if (strcmp(param_key, "max_root_zone_layer") == 0) {
model->soil_reservoir.max_root_zone_layer = strtod(param_value, NULL);
is_max_root_zone_layer_set = TRUE;
if (strcmp(param_key, "max_rootzone_layer") == 0) {
model->soil_reservoir.max_rootzone_layer = strtod(param_value, NULL);
is_max_rootzone_layer_set = TRUE;
}
}

Expand Down Expand Up @@ -786,7 +787,7 @@ int read_init_config_cfe(const char* config_file, cfe_state_struct* model)

/* Ice fraction: if set to true and runoff scheme is Schaake, additional parameters are needed in the config file,
*//////////////////////////////////////////////////////////////////////////////
if (strcmp(param_key, "sft_coupled") == 0) {
if (strcmp(param_key, "is_sft_coupled") == 0) {
if ( strcmp(param_value, "true")==0 || strcmp(param_value, "True")==0 || strcmp(param_value,"1")==0)
is_sft_coupled_set = TRUE;

Expand Down Expand Up @@ -978,7 +979,7 @@ int read_init_config_cfe(const char* config_file, cfe_state_struct* model)

if(!is_ice_content_threshold_set) {
#if CFE_DEBUG >= 1
printf("sft_coupled and Schaake scheme are set to TRUE but param 'ice_fraction_threshold' not found in config file\n");
printf("is_sft_coupled and Schaake scheme are set to TRUE but param 'ice_fraction_threshold' not found in config file\n");
exit(-9);
#endif
}
Expand All @@ -993,9 +994,9 @@ int read_init_config_cfe(const char* config_file, cfe_state_struct* model)

/*------------------- Root zone AET development -rlm----------------------------- */
if (is_aet_rootzone_set == TRUE ) {
if (is_max_root_zone_layer_set == FALSE) {
if (is_max_rootzone_layer_set == FALSE) {
#if CFE_DEBUG >= 1
printf("Config param 'max_root_zone_layer' not found in config file\n");
printf("Config param 'max_rootzone_layer' not found in config file\n");
#endif
return BMI_FAILURE;
}
Expand Down Expand Up @@ -1041,10 +1042,11 @@ int read_init_config_cfe(const char* config_file, cfe_state_struct* model)

// Calculate the thickness (delta) of each soil layer
if (is_aet_rootzone_set == TRUE ) {
model->soil_reservoir.is_aet_rootzone = TRUE;
model->soil_reservoir.delta_soil_layer_depth_m = malloc(sizeof(double) * (model->soil_reservoir.n_soil_layers + 1));
double previous_depth = 0;
double current_depth = 0;
for (int i=1; i <= model->soil_reservoir.n_soil_layers; i++){
for (int i=1; i <= model->soil_reservoir.n_soil_layers; i++) {
current_depth = model->soil_reservoir.soil_layer_depths_m[i];
if (current_depth <= previous_depth)
printf("WARNING: soil depths may be out of order. One or more soil layer depths is less than or equal to the previous layer. Check CFE config file.\n");
Expand All @@ -1058,9 +1060,10 @@ int read_init_config_cfe(const char* config_file, cfe_state_struct* model)
model->soil_reservoir.soil_water_content_field_capacity = model->NWM_soil_params.smcmax * pow(base, exponent);
}
else {
model->soil_reservoir.aet_root_zone = FALSE;
model->soil_reservoir.n_soil_layers = 0;
model->soil_reservoir.is_aet_rootzone = FALSE;
model->soil_reservoir.n_soil_layers = 0;
}

/*--------------------END OF ROOT ZONE ADJUSTED AET DEVELOPMENT -rlm ------------------------------*/

#if CFE_DEBUG >= 1
Expand Down Expand Up @@ -1305,8 +1308,8 @@ static int Initialize (Bmi *self, const char *file)
cfe_bmi_data_ptr->soil_reservoir.ice_fraction_schaake = 0.0;
cfe_bmi_data_ptr->soil_reservoir.ice_fraction_xinanjiang = 0.0;

if (cfe_bmi_data_ptr->soil_reservoir.aet_root_zone == TRUE)
cfe_bmi_data_ptr->soil_reservoir.smc_profile = malloc(sizeof(double)*cfe_bmi_data_ptr->soil_reservoir.n_soil_layers);
if (cfe_bmi_data_ptr->soil_reservoir.is_aet_rootzone == TRUE)
cfe_bmi_data_ptr->soil_reservoir.smc_profile = malloc(sizeof(double)*cfe_bmi_data_ptr->soil_reservoir.n_soil_layers + 1);
else
cfe_bmi_data_ptr->soil_reservoir.smc_profile = malloc(sizeof(double)*1);

Expand Down Expand Up @@ -1906,8 +1909,8 @@ static int Get_value_ptr (Bmi *self, const char *name, void **dest)
return BMI_SUCCESS;
}

//--------------Root zone adjusted AET development -rlm -ahmad -------
if (strcmp (name, "soil_moisture_profile") == 0){
//--------------Root zone adjusted AET development -rlm -ajk -------
if (strcmp (name, "soil_moisture_profile") == 0) {
*dest = (void *) ((cfe_state_struct *)(self->data))->soil_reservoir.smc_profile;
return BMI_SUCCESS;
}
Expand Down Expand Up @@ -1980,12 +1983,13 @@ static int Set_value_at_indices (Bmi *self, const char *name, int * inds, int le
// For now, all variables are non-array scalar values, with only 1 item of type double

// Thus, there is only ever one value to return (len must be 1) and it must always be from index 0
//AJ: modifying it to work with soil moisture column for root zone depth based AET
if (strcmp(name, "soil_moisture_profile") == 0 || strcmp(name, "soil_layer_depths_m") == 0){ //Adding soil layer depths since they will be needed for root zone adjusted AET estimations -rlm
void *ptr = NULL;
// ptr = (double*) malloc (sizeof (double)*4);
status = Get_value_ptr(self, name, &ptr);
// ajk: modifying it to work with soil moisture column for rootzone depth based AET
if (strcmp(name, "soil_moisture_profile") == 0 || strcmp(name, "soil_layer_depths_m") == 0) { //Adding soil layer depths since they will be needed for root zone adjusted AET estimations -rlm

len = ((cfe_state_struct *)(self->data))->soil_reservoir.n_soil_layers + 1;
void *ptr = NULL; //(double*) malloc (sizeof (double)* len);
status = Get_value_ptr(self, name, &ptr);


if (status == BMI_FAILURE)
return BMI_FAILURE;
Expand Down Expand Up @@ -2381,7 +2385,7 @@ static int Get_state_var_ptrs (Bmi *self, void *ptr_list[])
// Root zone AET development -rlm
// ------------------------------------------------------------
ptr_list[90] = &(state->soil_reservoir.soil_layer_depths_m);
ptr_list[91] = &(state->soil_reservoir.max_root_zone_layer);
ptr_list[91] = &(state->soil_reservoir.max_rootzone_layer);
//-------------------------------------------------------------


Expand Down
20 changes: 10 additions & 10 deletions src/cfe.c
View file
Original file line number Diff line number Diff line change
Expand Up @@ -610,24 +610,24 @@ void et_from_soil(struct conceptual_reservoir *soil_res, struct evapotranspirati

et_struct->actual_et_from_soil_m_per_timestep = 0;

// if root_zone-based AET turned ON
if (soil_res->aet_root_zone) {
// Assuming the layer with the most moisture is the bottom layer of the root zone (max_root_zone_layer)
// if rootzone-based AET turned ON
if (soil_res->is_aet_rootzone) {
// Assuming the layer with the most moisture is the bottom layer of the root zone (max_rootzone_layer)
// Convert volumetric soil moisture from the max root zone layer to moisture content [m] (layer_storage_m)
double layer_storage_m = soil_res->smc_profile[soil_res->max_root_zone_layer] *
soil_res->delta_soil_layer_depth_m[soil_res->max_root_zone_layer];
double layer_storage_m = soil_res->smc_profile[soil_res->max_rootzone_layer] *
soil_res->delta_soil_layer_depth_m[soil_res->max_rootzone_layer];

// If the moisture content from the layer with the most moisture is less than the
// wilting point, actual et from this timestep is 0 and no water is removed.
if ( soil_res->smc_profile[soil_res->max_root_zone_layer] <= soil_parms->wltsmc ) {
if ( soil_res->smc_profile[soil_res->max_rootzone_layer] <= soil_parms->wltsmc ) {
et_struct->actual_et_from_soil_m_per_timestep = 0;
}
// calculate the amount of moisture removed by evapotranspiration for the bottom layer of the root zone
else if (soil_res->smc_profile[soil_res->max_root_zone_layer] >= soil_res->soil_water_content_field_capacity) {
else if (soil_res->smc_profile[soil_res->max_rootzone_layer] >= soil_res->soil_water_content_field_capacity) {
et_struct->actual_et_from_soil_m_per_timestep = min(et_struct->reduced_potential_et_m_per_timestep, layer_storage_m);
}
else {
Budyko_numerator = soil_res->smc_profile[soil_res->max_root_zone_layer] - soil_parms->wltsmc;
Budyko_numerator = soil_res->smc_profile[soil_res->max_rootzone_layer] - soil_parms->wltsmc;
Budyko_denominator = soil_res->soil_water_content_field_capacity - soil_parms->wltsmc;
Budyko = Budyko_numerator / Budyko_denominator;

Expand All @@ -636,8 +636,8 @@ void et_from_soil(struct conceptual_reservoir *soil_res, struct evapotranspirati

// Reduce remaining PET and remove moisture from soil profile equal to the calculated AET (actual_et_from_soil_m_per_timestep
et_struct->reduced_potential_et_m_per_timestep -= et_struct->actual_et_from_soil_m_per_timestep;
soil_res->smc_profile[soil_res->max_root_zone_layer] -= (et_struct->actual_et_from_soil_m_per_timestep /
soil_res->delta_soil_layer_depth_m[soil_res->max_root_zone_layer]);
soil_res->smc_profile[soil_res->max_rootzone_layer] -= (et_struct->actual_et_from_soil_m_per_timestep /
soil_res->delta_soil_layer_depth_m[soil_res->max_rootzone_layer]);
soil_res->storage_m -= et_struct->actual_et_from_soil_m_per_timestep;
}

Expand Down

0 comments on commit 44bb30c

Please sign in to comment.