Skip to content

Update cfe.c #139

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
ajkhattak merged 2 commits into from
Mar 3, 2025
Merged

Update cfe.c #139

Changes from all commits
Commits
Show all changes
2 commits
Select commit Hold shift + click to select a range
8d4578b
Update cfe.c
fred-ogden Mar 3, 2025
e68cf1e
fix: adjust white space for better review, fix struct parameter
hellkite500 Mar 3, 2025
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
156 changes: 94 additions & 62 deletions src/cfe.c
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -403,7 +403,7 @@ void Schaake_partitioning_scheme(double timestep_h, double field_capacity_m, dou

void Xinanjiang_partitioning_scheme(double water_input_depth_m, double field_capacity_m,
double max_soil_moisture_storage_m, double column_total_soil_water_m,
struct infiltration_excess_parameters_structure *parms,
struct infiltration_excess_parameters_structure *parms,
double *infiltration_excess_m, double *infiltration_depth_m,
double ice_fraction_xinanjiang)
{
Expand All @@ -425,6 +425,8 @@ void Xinanjiang_partitioning_scheme(double water_input_depth_m, double field_cap
//-------------------------------------------------------------------------
// Written by RLM May 2021
// Adapted by JMFrame September 2021 for new version of CFE
// Reviewed by FLO Feb. 2025, compared against Jayawardena & Zhou (2000) and NWM 3.0 code. Fixed bug and refactored.
// Compared against NWM 3.0 and found identical performance within approx. 1e-05 on calculated outputs
//-------------------------------------------------------------------------
// Inputs
// double water_input_depth_m amount of water input to soil surface this time step [m]
Expand All @@ -438,85 +440,115 @@ void Xinanjiang_partitioning_scheme(double water_input_depth_m, double field_cap
// double ice_fraction_xinanjiang fraction of top soil layer that is frozen [unitless decimal]
//
// Outputs
// double infiltration_excess_m amount of water partitioned to surface water this time step [m]
// double infiltration_excess_m amount of water partitioned to surface water this time step [m]
// double infiltration_depth_m amount of water partitioned as infiltration (soil water input) this time step [m]
//-------------------------------------------------------------------------

double tension_water_m, free_water_m, max_tension_water_m, max_free_water_m, pervious_runoff_m;

// local variables
double tension_water_m, free_water_m, max_tension_water_m, max_free_water_m, pervious_runoff_m,water_input_pervious_fraction_m;
double impervious_fraction, impervious_runoff_m;
double f_over_F; // notation from Jayawardena and Zhou (2000) see Fig 2.

//could move this if statement outside of both the Schaake and Xinanjiang subroutines edit FLO- moved to main().

// partition the total soil water in the column between free water and tension water
free_water_m = column_total_soil_water_m - field_capacity_m;
// first thing, check whether we can just rreturn without calculating anything to save compute
//---------------------- NWM variable name
*infiltration_excess_m = 0.0; // RUNSRF
*infiltration_depth_m = 0.0; // PDDUM

if(0.0 < free_water_m) { //edit FLO
tension_water_m = field_capacity_m;
} else {
free_water_m = 0.0;
tension_water_m = column_total_soil_water_m;
if(water_input_depth_m < 1.0e-08) { // zero or really close to zero so calculations not needed
return;
}

// estimate the maximum free water and tension water available in the soil column
max_free_water_m = max_soil_moisture_storage_m - field_capacity_m;
max_tension_water_m = field_capacity_m;

if (max_free_water_m <= 0 || max_tension_water_m <=0) { //edited by RLM; added logic block to handle parameter values of zero.
*infiltration_excess_m = 0.95 * water_input_depth_m;

// calculations required

// initialize // NWM variable name
tension_water_m = 0.0; // WM
max_tension_water_m = 0.0; // WM_MAX
free_water_m = 0.0; // SM
max_free_water_m = 0.0; // SM_MAX
impervious_runoff_m = 0.0; // IRUNOFF


// Partition the total soil water in the column between free water and tension water assuming that
// total pore space in the soil Vtot, given by (porosity * soil_thickness) also equals max_tension_water_m + max_free_water_m.
//
// Water input to soil fills up the tension water storage first. Once it is full, any additiona
// input of water goes into free water storage. This means that iff there is _any_ free water then
// the tension water storage is full.

if( column_total_soil_water_m - field_capacity_m > 0.0) { // soil moisture greater than field capacity
free_water_m = column_total_soil_water_m - field_capacity_m;
tension_water_m = field_capacity_m;
} else {
tension_water_m = column_total_soil_water_m;
}
max_tension_water_m = field_capacity_m;
max_free_water_m = max_soil_moisture_storage_m - field_capacity_m;

if(tension_water_m > max_tension_water_m) tension_water_m = max_tension_water_m; // as done in NWM - could cause massbal err.
if(free_water_m > max_free_water_m) free_water_m = max_free_water_m; // as done in NWM - could cause massbal err.

// Calculate the impervious runoff (see eq. 309 from Knoben et al).
// estimate the fraction of the modeled area that is impervious (impervious_fraction) based on
// urban classification (hard coded 95% [0.95] impervious) and frozen soils (passed to cfe
// from freeze-thaw model) using a weighted average.

impervious_fraction = (parms->urban_decimal_fraction * 0.95) + ((1.0 - parms->urban_decimal_fraction) * ice_fraction_xinanjiang);
impervious_runoff_m = impervious_fraction * water_input_depth_m;

water_input_pervious_fraction_m = water_input_depth_m - impervious_runoff_m; // from here on just deal with
// water that enters the soil
//edited by RLM; added logic block to handle what happens when porosity or field capacity = 0
//FLO changed from 0.95 to 1.0 because what happens to the other 5% ?
if (max_free_water_m <= 0.0 || max_tension_water_m <=0.0) {
*infiltration_excess_m = 1.0 * water_input_pervious_fraction_m + impervious_runoff_m;
*infiltration_excess_m += impervious_runoff_m; // must add impervious runoff back in
*infiltration_depth_m = 0.0; //added by FLO
return; // this is an unusual situation and should only happen when one or both parameters are zero
}

// if code gets to here, then field capacity and soil porosity were both nonzero, so there is some space in
// the soil ffor more water, AND impervious runoff has already been abstracted from precipitation as impervious_runoff_m.

// solve pervious surface runoff (m) based on Eq. 310
// Use notation from Knoben et al. (2019) p. 71

double a = parms->a_Xinanjiang_inflection_point_parameter; // contributing area curve inflection point
double b = parms->b_Xinanjiang_shape_parameter; // contributing area curve shape parameter
double Ex = parms->x_Xinanjiang_shape_parameter; // contributing area curve shape parameter ffor direct runoff
double W = tension_water_m;
double Wmax = max_tension_water_m;

if ( ( W / Wmax ) <= (0.5 - a)) {
// THIS LINE IN ORIGINAL CFE XINANJIANG CODE IS A BUG
// f_over_F = (pow((0.5 - a),(1.0 - b)) * pow((1.0 - (tension_water_m/max_tension_water_m)) , b));

// FLO correct comparing against Eqn. 2a in Jayawardena & Zhou (2000)
f_over_F = (pow((0.5 - a),(1.0 - b)) * pow( W / Wmax , b)); // Eqn. 2a in Jayawardena & Zhue

} else { // if ( (0.5-a) < W / Wmax )
f_over_F = 1.0 - pow((0.5 + a), (1.0 - b)) * pow((1.0 - W / Wmax) , b); // Eqn. 2b in Jayawardena & Zhue
}

// check that the free_water_m and tension_water_m do not exceed the maximum and if so, change to the max value
if(max_free_water_m < free_water_m) free_water_m = max_free_water_m;
if(max_tension_water_m < tension_water_m) tension_water_m = max_tension_water_m;

// estimate the fraction of the modeled area that is impervious (impervious_fraction) based on
// urban classification (hard coded 95% [0.95] impervious) and frozen soils (passed to cfe
// from freeze-thaw model) using a weighted average.
impervious_fraction = (parms->urban_decimal_fraction * 0.95) + ((1 - parms->urban_decimal_fraction) * ice_fraction_xinanjiang);

// Calculate the impervious runoff (see eq. 309 from Knoben et al).
impervious_runoff_m = impervious_fraction * water_input_depth_m;

// Calculate total estimated pervious runoff.
if ((tension_water_m/max_tension_water_m) <= (0.5 - parms->a_Xinanjiang_inflection_point_parameter)) {
pervious_runoff_m = (1 - impervious_fraction) * water_input_depth_m *
(pow((0.5 - parms->a_Xinanjiang_inflection_point_parameter),
(1.0 - parms->b_Xinanjiang_shape_parameter)) *
pow((1.0 - (tension_water_m/max_tension_water_m)),
parms->b_Xinanjiang_shape_parameter));
double R = water_input_pervious_fraction_m * f_over_F; // water moved from tension to free water storage

} else {
pervious_runoff_m = (1 - impervious_fraction) * water_input_depth_m * (1.0 -
pow((0.5 + parms->a_Xinanjiang_inflection_point_parameter),
(1.0 - parms->b_Xinanjiang_shape_parameter)) *
pow((1.0 - (tension_water_m/max_tension_water_m)),
(parms->b_Xinanjiang_shape_parameter)));
}
// Separate the surface water from the pervious runoff
// NOTE: If impervious runoff is added to this subroutine, impervious runoff should be added to
// the infiltration_excess_m.
*infiltration_excess_m = pervious_runoff_m * (1.0 - pow((1.0 - (free_water_m/max_free_water_m)),parms->x_Xinanjiang_shape_parameter));
}
double S = free_water_m;
double Smax = max_free_water_m;

// Separate the surface water from the pervious runoff
*infiltration_excess_m = pervious_runoff_m * (1.0 - pow((1.0 -
(free_water_m/max_free_water_m)),parms->x_Xinanjiang_shape_parameter)) + impervious_runoff_m;
*infiltration_excess_m = R * (1.0 - pow((1.0 - (S/Smax)), Ex)) + impervious_runoff_m;

// The surface runoff depth is bounded by a minimum of 0 and a maximum of the water input depth.
// Check that the estimated surface runoff is not less than 0.0 and if so, change the value to 0.0.
if(*infiltration_excess_m < 0.0) *infiltration_excess_m = 0.0;
// Check that the estimated surface runoff does not exceed the amount of water input to the soil surface. If it does,
// change the surface water runoff value to the water input depth.
if(*infiltration_excess_m > water_input_depth_m) *infiltration_excess_m = water_input_depth_m;
// Separate the infiltration from the total water input depth to the soil surface.
*infiltration_depth_m = water_input_depth_m - *infiltration_excess_m;

return;
*infiltration_depth_m = water_input_pervious_fraction_m - *infiltration_excess_m;

#ifdef DEBUG
if(fabs(water_input_depth_m - (*infiltration_depth_m) - (*infiltration_excess_m)) > 1.0e-06) printf("Massball err. warning: %f\n",
water_input_depth_m - (*infiltration_depth_m) - (*infiltration_excess_m));
#endif
return;
}



//##############################################################
//#################### ET FROM RAINFALL ####################
//##############################################################
Expand Down