Skip to content

Commit

Permalink
revised nash surface scheme: changes to infiltration based on soil de…
Browse files Browse the repository at this point in the history 
…ficit, tested satdk as well (but not appropriate for channel infiltration due to small fraction of wetted area), removed linear K_infiltration due to non-linearity issues, and c1 being 2-3 order of magnitude smaller than c0.
  • Loading branch information
@ajkhattak
ajkhattak committed May 16, 2024
1 parent a5871c8 commit a4f27c9
Show file tree
Hide file tree
Showing 4 changed files with 58 additions and 51 deletions.
9 changes: 4 additions & 5 deletions include/nash_cascade.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -11,18 +11,17 @@ struct nash_cascade_parameters {
double K_nash; // Fraction of storage per hour that moves from one reservoir to the next (time constant); [1/hour]
int nsubsteps; // the number of substeps that each dt is divided into
double *nash_storage; // storage array nash cascade reservoirs; [m]
double retention_depth; // parameter that represents retention depth process that is equivalent [m]
// to that used in WRF-Hydro
double retention_depth; /* parameter that represents retention depth process that is equivalent [m]
to that used in WRF-Hydro */
double runon_infiltration; // infiltration losses from surface runoff water to soil (or riparian groundwater) [m/hr]
int is_riparian_gw; // flag to turn on/off riparian groundwater (currently used in LASAM only)
double K_infiltration; // Fraction of storage per hour that moves from reservoirs to soil (time constant); [1/hour]
double Kinf_c0; // constant c0 in K_infiltration = c0 + c1 * S, where S is the reservoir storage
double Kinf_c1; // constant c1 in K_infiltration = c0 + c1 * S, where S is the reservoir storage
};

extern double nash_cascade(double flux_lat_m,int num_lateral_flow_nash_reservoirs,
double K_nash,double *nash_storage_arr);

double nash_cascade_surface_runoff(double runoff_m, struct nash_cascade_parameters *nash_params);
double nash_cascade_surface_runoff(double runoff_m, double soil_storage_deficit_m,
struct nash_cascade_parameters *nash_params);

#endif
27 changes: 6 additions & 21 deletions src/bmi_cfe.c
View file
Original file line number Diff line number Diff line change
Expand Up @@ -18,14 +18,14 @@
#define STATE_VAR_NAME_COUNT 94 // must match var_info array size


#define PARAM_VAR_NAME_COUNT 19
#define PARAM_VAR_NAME_COUNT 18
// NOTE: If you update the params, also update the unit test in ../test/main_unit_test_bmi.c
static const char *param_var_names[PARAM_VAR_NAME_COUNT] = {
"maxsmc", "satdk", "slope", "b", "Klf",
"Kn", "Cgw", "expon", "max_gw_storage",
"satpsi","wltsmc","alpha_fc","refkdt",
"a_Xinanjiang_inflection_point_parameter","b_Xinanjiang_shape_parameter","x_Xinanjiang_shape_parameter",
"Kinf_c0", "Kinf_c1",
"K_infiltration",
"N_nash"
};

Expand All @@ -34,7 +34,6 @@ static const char *param_var_types[PARAM_VAR_NAME_COUNT] = {
"double", "double", "double", "double",
"double", "double", "double", "double",
"double","double","double", "double",
"double",
"int"
};
//----------------------------------------------
Expand Down Expand Up @@ -809,10 +808,7 @@ int read_init_config_cfe(const char* config_file, cfe_state_struct* model)
continue;
}
if (strcmp(param_key, "Kinf_nash_surface") == 0) {
//model->nash_surface_params.K_infiltration = strtod(param_value, NULL);
// Kinf = C0 + C1 * S, so if Kinf is provided in the config file, set C0 and C1 as folows:
model->nash_surface_params.Kinf_c0 = strtod(param_value, NULL);
model->nash_surface_params.Kinf_c1 = 0.0;
model->nash_surface_params.K_infiltration = strtod(param_value, NULL);
is_K_infiltration_nash_surface_set = TRUE;
continue;
}
Expand Down Expand Up @@ -1110,10 +1106,7 @@ int read_init_config_cfe(const char* config_file, cfe_state_struct* model)
#if CFE_DEBUG >= 1
printf("Config param 'Kinf_nash_surface' not found in config file, default value is 0.05 [1/hr] \n");
#endif
// model->nash_surface_params.K_infiltration = 0.05; // used in the runon infiltration
// Kinf = C0 + C1 * S, so if Kinf is not provided the default values for C0 and C1 are set below,
model->nash_surface_params.Kinf_c0 = 0.05; // set C0 = value
model->nash_surface_params.Kinf_c1 = 0.0; // set C1 = 0.0
model->nash_surface_params.K_infiltration = 0.05; // used in the runon infiltration
}
if (is_retention_depth_nash_surface_set == FALSE) {
#if CFE_DEBUG >= 1
Expand Down Expand Up @@ -1919,21 +1912,13 @@ static int Get_value_ptr (Bmi *self, const char *name, void **dest)
return BMI_SUCCESS;
}

if (strcmp (name, "Kinf_c0") == 0) {
if (strcmp (name, "K_infiltraton") == 0) {
cfe_state_struct *cfe_ptr;
cfe_ptr = (cfe_state_struct *) self->data;
*dest = (void*)&cfe_ptr->nash_surface_params.Kinf_c0;
*dest = (void*)&cfe_ptr->nash_surface_params.K_infiltration;
return BMI_SUCCESS;
}

if (strcmp (name, "Kinf_c1") == 0) {
cfe_state_struct *cfe_ptr;
cfe_ptr = (cfe_state_struct *) self->data;
*dest = (void*)&cfe_ptr->nash_surface_params.Kinf_c1;
return BMI_SUCCESS;
}


/***********************************************************/
/*********** OUTPUT ***********************************/
/***********************************************************/
Expand Down
8 changes: 5 additions & 3 deletions src/cfe.c
View file
Original file line number Diff line number Diff line change
Expand Up @@ -244,10 +244,12 @@ extern void cfe(
direct_runoff_m = giuh_convolution_integral(infiltration_excess_m,num_giuh_ordinates,
giuh_ordinates_arr,runoff_queue_m_per_timestep_arr);
else if (surface_runoff_scheme == NASH_CASCADE)
direct_runoff_m = nash_cascade_surface_runoff(infiltration_excess_m, nash_surface_params);

direct_runoff_m = nash_cascade_surface_runoff(infiltration_excess_m, soil_reservoir_storage_deficit_m,
nash_surface_params);

soil_reservoir_struct->storage_m += nash_surface_params->runon_infiltration; // put the runon infiltrated water in the soil.

soil_reservoir_storage_deficit_m -= nash_surface_params->runon_infiltration; // update soil deficit

massbal_struct->vol_out_surface += direct_runoff_m;
massbal_struct->volout += direct_runoff_m;
massbal_struct->volout += flux_from_deep_gw_to_chan_m;
Expand Down
65 changes: 43 additions & 22 deletions src/nash_cascade.c
View file
Original file line number Diff line number Diff line change
Expand Up @@ -40,7 +40,8 @@ extern double nash_cascade(double flux_lat_m,int num_lateral_flow_nash_reservoir
//##############################################################
//############## NASH CASCADE SURFACE RUNOFF ################
//##############################################################
double nash_cascade_surface_runoff(double runoff_m, struct nash_cascade_parameters *nash_params)
double nash_cascade_surface_runoff(double runoff_m, double soil_storage_deficit_m,
struct nash_cascade_parameters *nash_params)
{
//##############################################################
// Solve for the flow through the Nash cascade to delay the
Expand All @@ -51,9 +52,7 @@ double nash_cascade_surface_runoff(double runoff_m, struct nash_cascade_paramete
int N_nash = nash_params->N_nash;
double K_nash = nash_params->K_nash;
double depth_r = nash_params->retention_depth;
//double K_infil = nash_params->K_infiltration;
double c0 = nash_params->Kinf_c0;
double c1 = nash_params->Kinf_c1;
double K_infil = nash_params->K_infiltration;

// local vars
double dt_h = 1.0; // model timestep [hour]
Expand All @@ -64,10 +63,11 @@ double nash_cascade_surface_runoff(double runoff_m, struct nash_cascade_paramete
double Q_r; // discharge from reservoir
double Q_out; // discharge at the outlet (the last reservoir) per subtimestep
double Q_infil; // discharge from reservoirs to soil
double K_infil; // K infiltration = c0 + c1 * S
double Q_to_channel_m = 0.0; // total outflow to channel per timestep
double Q_to_soil_m = 0.0; // runon infiltration (losses from surface runoff to soil)

double soil_deficit_m = soil_storage_deficit_m; // local variable to track the soil storage deficit
double infil_m;

nash_params->nash_storage[0] += runoff_m;

// Loop through number of sub-timesteps
Expand All @@ -76,9 +76,44 @@ double nash_cascade_surface_runoff(double runoff_m, struct nash_cascade_paramete
//Loop through reservoirs
for(int i = 0; i < N_nash; i++) {

// First: infiltration capacity should be satisfied before routing water through Nash reservoirs

// compute runon infiltration
Q_infil = K_infil * nash_params->nash_storage[i];
infil_m = Q_infil * subdt;

// case 1: soil_deficit greater than infil_m, all water can infiltrate, and soil_deficit decreases
// case 2: soil_deficit less than infil_m, portion of water can infiltrate, soil_deficit gets zero
// case 3: soil_deficit is zero, no infiltration

// determine the right amount of infiltrated water based on soil deficit
if (soil_deficit_m > 0.0 && soil_deficit_m > infil_m) {
soil_deficit_m -= infil_m; // update soil deficit for the next subtimstep iteration
}
else if (soil_deficit_m > 0.0) {
infil_m = soil_deficit_m;
soil_deficit_m = 0.0;
}
else {
infil_m = 0.0; // if got here, then soil is fully saturated, so no infiltration
}

// update nash storage (subtract infiltrated water)
if (nash_params->nash_storage[i] >= infil_m)
nash_params->nash_storage[i] = nash_params->nash_storage[i] - infil_m;
else {
infil_m = nash_params->nash_storage[i];
nash_params->nash_storage[i] = 0.0;
}


Q_to_soil_m += infil_m; // water volume that infiltrates per subtimestep from each reservoir

/*=========================================================================================*/
// Second: time to route water through Nash reservoirs
S = nash_params->nash_storage[i];

// first reservoir
// determine the amount of surface water available for routing from the first reservoir
if (i == 0 && S > depth_r)
S -= depth_r;
else if (i == 0)
Expand All @@ -92,21 +127,7 @@ double nash_cascade_surface_runoff(double runoff_m, struct nash_cascade_paramete
nash_params->nash_storage[i+1] += dS;
else
Q_out = Q_r;

//compute runon infiltration
K_infil = c0 + c1 * nash_params->nash_storage[i];
Q_infil = K_infil * nash_params->nash_storage[i];
dS_infil = Q_infil * subdt;

if (nash_params->nash_storage[i] >= dS_infil)
nash_params->nash_storage[i] = nash_params->nash_storage[i] - dS_infil;
else {
dS_infil = nash_params->nash_storage[i];
nash_params->nash_storage[i] = 0.0;
}

Q_to_soil_m += dS_infil; // water volume that infiltrates per subtimestep from each reservoir


}

Q_to_channel_m += Q_out * subdt; // Q_r at the end of N_nash loop is the discharge at the outlet
Expand Down

0 comments on commit a4f27c9

Please sign in to comment.