From 464ecdecf81d15e912a31449421b27e2f3b141cc Mon Sep 17 00:00:00 2001 From: Ahmad Jan Date: Thu, 13 Feb 2025 22:43:27 -0500 Subject: [PATCH] fix: AET should come from surface water retention depth before it is extracted from the soil. --- include/cfe.h | 5 +++++ src/bmi_cfe.c | 12 ++++++----- src/cfe.c | 56 ++++++++++++++++++++++++++++++++++++++++++--------- 3 files changed, 59 insertions(+), 14 deletions(-) diff --git a/include/cfe.h b/include/cfe.h index 66e608ff..045f7deb 100644 --- a/include/cfe.h +++ b/include/cfe.h @@ -91,6 +91,7 @@ struct evapotranspiration_structure { double potential_et_m_per_timestep; double reduced_potential_et_m_per_timestep; double actual_et_from_rain_m_per_timestep; + double actual_et_from_retention_depth_m_per_timestep; double actual_et_from_soil_m_per_timestep; double actual_et_m_per_timestep; }; @@ -122,6 +123,7 @@ struct massbal double volin ; double volout ; double volend ; + double vol_et_from_retention_depth; }; typedef struct massbal massbal; @@ -158,6 +160,9 @@ extern void Xinanjiang_partitioning_scheme(double water_input_depth_m, double fi extern void et_from_rainfall(double *timestep_rainfall_input_m, struct evapotranspiration_structure *et_struct); +extern void et_from_retention_depth(struct nash_cascade_parameters *nash_surface_params, + struct evapotranspiration_structure *et_struct); + extern void et_from_soil(struct conceptual_reservoir *soil_res, struct evapotranspiration_structure *et_struct, struct NWM_soil_parameters *soil_parms); diff --git a/src/bmi_cfe.c b/src/bmi_cfe.c index b7dbd8f4..0b1280d0 100644 --- a/src/bmi_cfe.c +++ b/src/bmi_cfe.c @@ -3227,6 +3227,7 @@ extern void initialize_volume_trackers(cfe_state_struct* cfe_ptr) { cfe_ptr->vol_struct.vol_soil_start = cfe_ptr->soil_reservoir.storage_m; cfe_ptr->vol_struct.vol_et_from_soil = 0; cfe_ptr->vol_struct.vol_et_from_rain = 0; + cfe_ptr->vol_struct.vol_et_from_retention_depth = 0; cfe_ptr->vol_struct.vol_et_to_atm = 0; } @@ -3325,11 +3326,12 @@ extern void mass_balance_check(cfe_state_struct* cfe_ptr){ giuh_residual = cfe_ptr->vol_struct.vol_runoff - cfe_ptr->vol_struct.vol_out_surface - cfe_ptr->vol_struct.vol_end_surface - cfe_ptr->vol_struct.vol_runon_infilt; printf("********************* SURFACE MASS BALANCE *********************\n"); - printf(" Volume into surface = %8.4lf m\n",cfe_ptr->vol_struct.vol_runoff); - printf(" Volume out surface = %8.4lf m\n",cfe_ptr->vol_struct.vol_out_surface); - printf(" Volume end surface = %8.4lf m\n",cfe_ptr->vol_struct.vol_end_surface); - printf(" Runon infiltration = %8.4lf m\n",cfe_ptr->vol_struct.vol_runon_infilt); - printf(" Surface residual = %6.4e m\n",giuh_residual); // should equal zero + printf(" Volume into surface = %8.4lf m\n",cfe_ptr->vol_struct.vol_runoff); + printf(" Volume out surface = %8.4lf m\n",cfe_ptr->vol_struct.vol_out_surface); + printf(" Volume end surface = %8.4lf m\n",cfe_ptr->vol_struct.vol_end_surface); + printf(" Runon infiltration = %8.4lf m\n",cfe_ptr->vol_struct.vol_runon_infilt); + printf(" Vol_et_from_retention_depth = %8.4lf m\n",cfe_ptr->vol_struct.vol_et_from_retention_depth); + printf(" Surface residual = %6.4e m\n", giuh_residual); // should equal zero if(!is_fabs_less_than_epsilon(giuh_residual,1.0e-12)) printf("WARNING: GIUH MASS BALANCE CHECK FAILED\n"); diff --git a/src/cfe.c b/src/cfe.c index 484538bc..0215b0fa 100644 --- a/src/cfe.c +++ b/src/cfe.c @@ -73,21 +73,40 @@ extern void cfe( evap_struct->reduced_potential_et_m_per_timestep = evap_struct->potential_et_m_per_s * time_step_size; evap_struct->actual_et_from_rain_m_per_timestep = 0; - if(timestep_rainfall_input_m > 0) {et_from_rainfall(×tep_rainfall_input_m,evap_struct);} + if (timestep_rainfall_input_m > 0) { + et_from_rainfall(×tep_rainfall_input_m,evap_struct); + } massbal_struct->vol_et_from_rain = massbal_struct->vol_et_from_rain + evap_struct->actual_et_from_rain_m_per_timestep; - massbal_struct->vol_et_to_atm = massbal_struct->vol_et_to_atm + evap_struct->actual_et_from_rain_m_per_timestep; - massbal_struct->volout=massbal_struct->volout+evap_struct->actual_et_from_rain_m_per_timestep; + massbal_struct->vol_et_to_atm = massbal_struct->vol_et_to_atm + evap_struct->actual_et_from_rain_m_per_timestep; + massbal_struct->volout = massbal_struct->volout + evap_struct->actual_et_from_rain_m_per_timestep; + + // AET from surface retention depth + evap_struct->actual_et_from_retention_depth_m_per_timestep = 0; + if (nash_surface_params->nash_storage[0] > 0.0 && evap_struct->reduced_potential_et_m_per_timestep > 0.0) { + et_from_retention_depth(nash_surface_params, evap_struct); + } - // LKC: Change this. Now evaporation happens before runoff calculation. This was creating issues since it modifies storage_m and not storage_deficit + massbal_struct->vol_et_from_retention_depth += evap_struct->actual_et_from_retention_depth_m_per_timestep; + massbal_struct->vol_et_to_atm += evap_struct->actual_et_from_retention_depth_m_per_timestep; + massbal_struct->volout += evap_struct->actual_et_from_retention_depth_m_per_timestep; + massbal_struct->vol_out_surface += evap_struct->actual_et_from_retention_depth_m_per_timestep; + + // LKC: Change this. Now evaporation happens before runoff calculation. This was creating issues since it modifies + // storage_m and not storage_deficit evap_struct->actual_et_from_soil_m_per_timestep = 0; - if(soil_reservoir_struct->storage_m > NWM_soil_params_struct.wilting_point_m) - {et_from_soil(soil_reservoir_struct, evap_struct, &NWM_soil_params_struct);} + if(soil_reservoir_struct->storage_m > NWM_soil_params_struct.wilting_point_m) { + et_from_soil(soil_reservoir_struct, evap_struct, &NWM_soil_params_struct); + } + massbal_struct->vol_et_from_soil = massbal_struct->vol_et_from_soil + evap_struct->actual_et_from_soil_m_per_timestep; - massbal_struct->vol_et_to_atm = massbal_struct->vol_et_to_atm + evap_struct->actual_et_from_soil_m_per_timestep; - massbal_struct->volout=massbal_struct->volout+evap_struct->actual_et_from_soil_m_per_timestep; + massbal_struct->vol_et_to_atm = massbal_struct->vol_et_to_atm + evap_struct->actual_et_from_soil_m_per_timestep; + massbal_struct->volout = massbal_struct->volout + evap_struct->actual_et_from_soil_m_per_timestep; + + evap_struct->actual_et_m_per_timestep = evap_struct->actual_et_from_rain_m_per_timestep + + evap_struct->actual_et_from_retention_depth_m_per_timestep + + evap_struct->actual_et_from_soil_m_per_timestep; - evap_struct->actual_et_m_per_timestep=evap_struct->actual_et_from_rain_m_per_timestep+evap_struct->actual_et_from_soil_m_per_timestep; // LKC: This needs to be calcualted here after et_from_soil since soil_reservoir_struct->storage_m changes soil_reservoir_storage_deficit_m=(NWM_soil_params_struct.smcmax*NWM_soil_params_struct.D-soil_reservoir_struct->storage_m); @@ -545,6 +564,25 @@ void et_from_rainfall(double *timestep_rainfall_input_m, struct evapotranspirati } } +//############################################################## +//########### ET FROM SURFACE RETENTION DEPTH ############## +//############################################################## +void et_from_retention_depth(struct nash_cascade_parameters *nash_surface_params, struct evapotranspiration_structure *et_struct) +{ + if (et_struct->reduced_potential_et_m_per_timestep >= nash_surface_params->nash_storage[0]) { + et_struct->actual_et_from_retention_depth_m_per_timestep = nash_surface_params->nash_storage[0]; + nash_surface_params->nash_storage[0] = 0.0; + } + else { + et_struct->actual_et_from_retention_depth_m_per_timestep = et_struct->reduced_potential_et_m_per_timestep; + //et_struct->reduced_potential_et_m_per_timestep = 0.0; + nash_surface_params->nash_storage[0] -= et_struct->actual_et_from_retention_depth_m_per_timestep; + } + + et_struct->reduced_potential_et_m_per_timestep -= et_struct->actual_et_from_retention_depth_m_per_timestep; + +} + //############################################################## //#################### ET FROM SOIL ######################## //##############################################################