From c0a813a35ff28fd1947f2e6d5d4aaa024a3503fd Mon Sep 17 00:00:00 2001 From: Ryoko Araki <400mhs2@gmail.com> Date: Tue, 12 Apr 2022 14:45:50 -0700 Subject: [PATCH 1/6] debugged. '_m' were missing on primary&secondaryflux --- cfe.py | 7 ++++--- 1 file changed, 4 insertions(+), 3 deletions(-) diff --git a/cfe.py b/cfe.py index 33bef5d7..0e3feb47 100755 --- a/cfe.py +++ b/cfe.py @@ -2,6 +2,7 @@ import numpy as np import pandas as pd + class CFE(): def __init__(self): super(CFE, self).__init__() @@ -68,8 +69,8 @@ def run_cfe(self, cfe_state): self.conceptual_reservoir_flux_calc(cfe_state, cfe_state.soil_reservoir) # ________________________________________________ - cfe_state.flux_perc_m = cfe_state.primary_flux - cfe_state.flux_lat_m = cfe_state.secondary_flux + cfe_state.flux_perc_m = cfe_state.primary_flux_m + cfe_state.flux_lat_m = cfe_state.secondary_flux_m # ________________________________________________ cfe_state.gw_reservoir_storage_deficit_m = cfe_state.gw_reservoir['storage_max_m'] - cfe_state.gw_reservoir['storage_m'] @@ -98,7 +99,7 @@ def run_cfe(self, cfe_state): self.conceptual_reservoir_flux_calc(cfe_state, cfe_state.gw_reservoir) # ________________________________________________ - cfe_state.flux_from_deep_gw_to_chan_m = cfe_state.primary_flux + cfe_state.flux_from_deep_gw_to_chan_m = cfe_state.primary_flux_m cfe_state.vol_from_gw += cfe_state.flux_from_deep_gw_to_chan_m # ________________________________________________ From ef03a5eac0ea11f616bc08b3795c0ac6f27d17b2 Mon Sep 17 00:00:00 2001 From: Ryoko Araki <400mhs2@gmail.com> Date: Tue, 12 Apr 2022 14:50:20 -0700 Subject: [PATCH 2/6] debugged. '_m' was missing in primary&secondary fluxes --- cfe.py | 1 - 1 file changed, 1 deletion(-) diff --git a/cfe.py b/cfe.py index 0e3feb47..5444683a 100755 --- a/cfe.py +++ b/cfe.py @@ -2,7 +2,6 @@ import numpy as np import pandas as pd - class CFE(): def __init__(self): super(CFE, self).__init__() From 9fe18c9914d597673d971b43ee3c62b824c63232 Mon Sep 17 00:00:00 2001 From: Ryoko Araki <400mhs2@gmail.com> Date: Tue, 12 Apr 2022 16:00:48 -0700 Subject: [PATCH 3/6] an additioanl fix on secondary_flux --- cfe.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/cfe.py b/cfe.py index 5444683a..9ee4fc1c 100755 --- a/cfe.py +++ b/cfe.py @@ -102,7 +102,7 @@ def run_cfe(self, cfe_state): cfe_state.vol_from_gw += cfe_state.flux_from_deep_gw_to_chan_m # ________________________________________________ - if not self.is_fabs_less_than_epsilon(cfe_state.secondary_flux, 1.0e-09): + if not self.is_fabs_less_than_epsilon(cfe_state.secondary_flux_m, 1.0e-09): print("problem with nonzero flux point 1\n") # ________________________________________________ From 4579f27337223ca627659733c89446500272ab87 Mon Sep 17 00:00:00 2001 From: Ryoko Araki <400mhs2@gmail.com> Date: Tue, 12 Apr 2022 16:02:48 -0700 Subject: [PATCH 4/6] an additioal fix on secondary_flux --- cfe.py | 1 + 1 file changed, 1 insertion(+) diff --git a/cfe.py b/cfe.py index 9ee4fc1c..fa292577 100755 --- a/cfe.py +++ b/cfe.py @@ -2,6 +2,7 @@ import numpy as np import pandas as pd + class CFE(): def __init__(self): super(CFE, self).__init__() From d25cc5f1b716aa69d6ee7d4f40f26008052c0c35 Mon Sep 17 00:00:00 2001 From: Ryoko Araki <400mhs2@gmail.com> Date: Wed, 21 Sep 2022 10:28:46 -0700 Subject: [PATCH 5/6] Added ODE solver for the soil reservoir calculation --- bmi_cfe.py | 33 ++++++--- cfe.py | 214 +++++++++++++++++++++++++++++++++++++++++------------ 2 files changed, 191 insertions(+), 56 deletions(-) diff --git a/bmi_cfe.py b/bmi_cfe.py index 28a1de7e..e20ce21b 100755 --- a/bmi_cfe.py +++ b/bmi_cfe.py @@ -134,11 +134,16 @@ def initialize(self, current_time_step=0): self.primary_flux = 0 # temporary vars. self.secondary_flux = 0 # temporary vars. self.total_discharge = 0 + self.diff_infilt = 0 + self.diff_perc = 0 # ________________________________________________ # Evapotranspiration self.potential_et_m_per_timestep = 0 self.actual_et_m_per_timestep = 0 + self.reduced_potential_et_m_per_timestep = 0 + self.actual_et_from_rain_m_per_timestep = 0 + self.actual_et_from_soil_m_per_timestep = 0 # ________________________________________________________ # Set these values now that we have the information from the configuration file. @@ -283,9 +288,13 @@ def reset_volume_tracking(self): self.vol_soil_to_lat_flow = 0 self.vol_soil_to_gw = 0 self.vol_soil_end = 0 + self.vol_et_to_atm = 0 + self.vol_et_from_soil = 0 + self.vol_et_from_rain = 0 self.volin = 0 self.volout = 0 self.volend = 0 + self.vol_PET = 0 return #________________________________________________________ @@ -351,13 +360,14 @@ def finalize_mass_balance(self, verbose=True): self.vol_soil_end = self.soil_reservoir['storage_m'] self.global_residual = self.volstart + self.volin - self.volout - self.volend -self.vol_end_giuh - self.schaake_residual = self.volin - self.vol_sch_runoff - self.vol_sch_infilt + self.schaake_residual = self.volin - self.vol_et_from_rain - self.vol_sch_runoff - self.vol_sch_infilt self.giuh_residual = self.vol_sch_runoff - self.vol_out_giuh - self.vol_end_giuh - self.soil_residual = self.vol_soil_start + self.vol_sch_infilt - \ - self.vol_soil_to_lat_flow - self.vol_soil_end - self.vol_to_gw + self.soil_residual = self.vol_soil_start - self.vol_soil_end + self.vol_sch_infilt - \ + self.vol_soil_to_lat_flow - self.vol_soil_to_gw - self.vol_et_from_soil self.nash_residual = self.vol_in_nash - self.vol_out_nash - self.vol_in_nash_end self.gw_residual = self.vol_in_gw_start + self.vol_to_gw - self.vol_from_gw - self.vol_in_gw_end - if verbose: + + if verbose: print("\nGLOBAL MASS BALANCE") print(" initial volume: {:8.4f}".format(self.volstart)) print(" volume input: {:8.4f}".format(self.volin)) @@ -365,24 +375,30 @@ def finalize_mass_balance(self, verbose=True): print(" final volume: {:8.4f}".format(self.volend)) print(" residual: {:6.4e}".format(self.global_residual)) + print("\n AET & PET") + print(" volume PET: {:8.4f}".format(self.vol_PET)) + print(" volume AET: {:8.4f}".format(self.vol_et_to_atm)) print("\nSCHAAKE MASS BALANCE") + print(" volume input: {:8.4f}".format(self.volin)) + print("ET from rainfall: {:8.4f}".format(self.vol_et_from_rain)) print(" surface runoff: {:8.4f}".format(self.vol_sch_runoff)) print(" infiltration: {:8.4f}".format(self.vol_sch_infilt)) - print("schaake residual: {:6.4e}".format(self.schaake_residual)) + print("schaake residual: {:6.4e}".format(self.schaake_residual)) - print("\nGIUH MASS BALANCE"); + print("\nGIUH MASS BALANCE") print(" vol. into giuh: {:8.4f}".format(self.vol_sch_runoff)) print(" vol. out giuh: {:8.4f}".format(self.vol_out_giuh)) print(" vol. end giuh q: {:8.4f}".format(self.vol_end_giuh)) print(" giuh residual: {:6.4e}".format(self.giuh_residual)) print("\nSOIL WATER CONCEPTUAL RESERVOIR MASS BALANCE") - print(" init soil vol: {:8.4f}".format(self.vol_soil_start)) + print(" init soil vol: {:8.4f}".format(self.vol_soil_start)) print(" vol. into soil: {:8.4f}".format(self.vol_sch_infilt)) print("vol.soil2latflow: {:8.4f}".format(self.vol_soil_to_lat_flow)) print(" vol. soil to gw: {:8.4f}".format(self.vol_soil_to_gw)) - print(" final vol. soil: {:8.4f}".format(self.vol_soil_end)) + print(" vol. soil to ET: {:8.4f}".format(self.vol_et_from_soil)) + print(" final vol. soil: {:8.4f}".format(self.vol_soil_end)) print("vol. soil resid.: {:6.4e}".format(self.soil_residual)) @@ -399,7 +415,6 @@ def finalize_mass_balance(self, verbose=True): print(" vol from gw: {:8.4f}".format(self.vol_from_gw)) print("final gw.storage: {:8.4f}".format(self.vol_in_gw_end)) print(" gw. residual: {:6.4e}".format(self.gw_residual)) - return diff --git a/cfe.py b/cfe.py index fa292577..c03ba7ac 100755 --- a/cfe.py +++ b/cfe.py @@ -1,7 +1,36 @@ import time import numpy as np import pandas as pd +from scipy.integrate import odeint +import math +def conceptual_reservoir_flux_calc(t, S, storage_threshold_primary_m, storage_max_m, coeff_primary, coeff_secondary, PET, infilt, wltsmc): + # ODE of soil moisture reservoir mass balance + storage_above_threshold_m = S - storage_threshold_primary_m + storage_diff = storage_max_m - storage_threshold_primary_m + storage_ratio = np.minimum(storage_above_threshold_m / storage_diff, 1) + + perc_lat_switch = np.multiply(S - storage_threshold_primary_m > 0, 1) + ET_switch = np.multiply(S - wltsmc > 0, 1) + + storage_above_threshold_m_paw = S - wltsmc + storage_diff_paw = storage_threshold_primary_m - wltsmc + storage_ratio_paw = np.minimum(storage_above_threshold_m_paw/storage_diff_paw, 1) # Equation 11 (Ogden's document). + + dS = infilt -1 * perc_lat_switch * (coeff_primary + coeff_secondary) * storage_ratio - ET_switch * PET * storage_ratio_paw + return dS + +def jac(t, S, storage_threshold_primary_m, storage_max_m, coeff_primary, coeff_secondary, PET, infilt, wltsmc): + # Jacobian matrix for the conceptual_reservoir_flux_calc + storage_diff = storage_max_m - storage_threshold_primary_m + + perc_lat_switch = np.multiply(S - storage_threshold_primary_m > 0, 1) + ET_switch = np.multiply((S - wltsmc > 0) and (S - storage_threshold_primary_m < 0), 1) + + storage_diff_paw = storage_threshold_primary_m - wltsmc + + dfdS = -1 * perc_lat_switch * (coeff_primary + coeff_secondary) * 1/storage_diff - ET_switch * PET * 1/storage_diff_paw + return [dfdS] class CFE(): def __init__(self): @@ -16,11 +45,18 @@ def run_cfe(self, cfe_state): # ________________________________________________ cfe_state.potential_et_m_per_timestep = cfe_state.potential_et_m_per_s * cfe_state.time_step_size - + cfe_state.reduced_potential_et_m_per_timestep = cfe_state.potential_et_m_per_s * cfe_state.time_step_size + cfe_state.vol_PET += cfe_state.potential_et_m_per_timestep + # ________________________________________________ # SUBROUTINE # timestep_rainfall_input_m = f(timestep_rainfall_input_m, potential_et_m_per_timestep) + cfe_state.actual_et_from_rain_m_per_timestep = 0 self.et_from_rainfall(cfe_state) + + cfe_state.vol_et_from_rain += cfe_state.actual_et_from_rain_m_per_timestep + cfe_state.vol_et_to_atm += cfe_state.actual_et_from_rain_m_per_timestep + cfe_state.volout += cfe_state.actual_et_from_rain_m_per_timestep # ________________________________________________ cfe_state.soil_reservoir_storage_deficit_m = (cfe_state.soil_params['smcmax'] * \ @@ -53,34 +89,38 @@ def run_cfe(self, cfe_state): # ________________________________________________ if cfe_state.previous_flux_perc_m > cfe_state.soil_reservoir_storage_deficit_m: - diff = cfe_state.previous_flux_perc_m - cfe_state.soil_reservoir_storage_deficit_m + self.diff_infilt = cfe_state.previous_flux_perc_m - cfe_state.soil_reservoir_storage_deficit_m cfe_state.infiltration_depth_m = cfe_state.soil_reservoir_storage_deficit_m - cfe_state.vol_sch_runoff += diff - cfe_state.vol_sch_infilt -= diff - cfe_state.surface_runoff_depth_m += diff + cfe_state.vol_sch_runoff += self.diff_infilt + cfe_state.vol_sch_infilt -= self.diff_infilt + cfe_state.surface_runoff_depth_m += self.diff_infilt - # ________________________________________________ - cfe_state.vol_to_soil += cfe_state.infiltration_depth_m - cfe_state.soil_reservoir['storage_m'] += cfe_state.infiltration_depth_m - # ________________________________________________ # SUBROUTINE - # primary_flux, secondary_flux = f(reservoir) - self.conceptual_reservoir_flux_calc(cfe_state, cfe_state.soil_reservoir) + # Solve ODE for the soil reservoir + # Calculates primary_flux, secondary_flux, AET, and infiltration simultaneously + cfe_state.actual_et_from_soil_m_per_timestep = 0 + self.soil_reservoir_flux_calc(cfe_state, cfe_state.soil_reservoir) # ________________________________________________ cfe_state.flux_perc_m = cfe_state.primary_flux_m cfe_state.flux_lat_m = cfe_state.secondary_flux_m + cfe_state.vol_to_soil += cfe_state.infiltration_depth_m + cfe_state.vol_et_from_soil += cfe_state.actual_et_from_soil_m_per_timestep + cfe_state.vol_et_to_atm += cfe_state.actual_et_from_soil_m_per_timestep + cfe_state.volout += cfe_state.actual_et_from_soil_m_per_timestep # ________________________________________________ - cfe_state.gw_reservoir_storage_deficit_m = cfe_state.gw_reservoir['storage_max_m'] - cfe_state.gw_reservoir['storage_m'] - + cfe_state.gw_reservoir_storage_deficit_m = cfe_state.gw_reservoir['storage_max_m'] - cfe_state.gw_reservoir[ + 'storage_m'] + # ________________________________________________ if cfe_state.flux_perc_m > cfe_state.gw_reservoir_storage_deficit_m: - diff = cfe_state.flux_perc_m - cfe_state.gw_reservoir_storage_deficit_m + # (?) When the groundwater storage is full, the overflowing amount goes to direct runoff? + self.diff_perc = cfe_state.flux_perc_m - cfe_state.gw_reservoir_storage_deficit_m cfe_state.flux_perc_m = cfe_state.gw_reservoir_storage_deficit_m - cfe_state.vol_sch_runoff+=diff - cfe_state.vol_sch_infilt-=diff + cfe_state.vol_sch_runoff += self.diff_perc + cfe_state.vol_sch_infilt -= self.diff_perc # ________________________________________________ cfe_state.vol_to_gw += cfe_state.flux_perc_m @@ -90,7 +130,7 @@ def run_cfe(self, cfe_state): cfe_state.soil_reservoir['storage_m'] -= cfe_state.flux_perc_m cfe_state.soil_reservoir['storage_m'] -= cfe_state.flux_lat_m cfe_state.vol_soil_to_lat_flow += cfe_state.flux_lat_m #TODO add this to nash cascade as input - cfe_state.volout = cfe_state.volout + cfe_state.flux_lat_m; + cfe_state.volout += cfe_state.flux_lat_m # ________________________________________________ @@ -203,21 +243,102 @@ def et_from_rainfall(self,cfe_state): """ iff it is raining, take PET from rainfall first. Wet veg. is efficient evaporator. """ - - if cfe_state.timestep_rainfall_input_m >0.0: + if cfe_state.timestep_rainfall_input_m > 0.0: if cfe_state.timestep_rainfall_input_m > cfe_state.potential_et_m_per_timestep: - - cfe_state.actual_et_m_per_timestep = cfe_state.potential_et_m_per_timestep - cfe_state.timestep_rainfall_input_m -= cfe_state.actual_et_m_per_timestep - else: + cfe_state.actual_et_from_rain_m_per_timestep = cfe_state.potential_et_m_per_timestep + cfe_state.timestep_rainfall_input_m -= cfe_state.actual_et_from_rain_m_per_timestep + + else: + cfe_state.actual_et_from_rain_m_per_timestep = cfe_state.timestep_rainfall_input_m + cfe_state.timestep_rainfall_input_m = 0.0 + + cfe_state.reduced_potential_et_m_per_timestep = cfe_state.potential_et_m_per_timestep - cfe_state.actual_et_from_rain_m_per_timestep - cfe_state.potential_et_m_per_timestep -= cfe_state.timestep_rainfall_input_m - cfe_state.timestep_rainfall_input_m=0.0 return - + + def soil_reservoir_flux_calc(self, cfe_state, reservoir): + """ + This function solves ODE for a soil reservoir + :param cfe_state: + :param reservoir: + :return: reservoir['storage_m'] + """ + + # Initialization + y0 = [reservoir['storage_m']] + + t = np.array([0, 0.05, 0.15, 0.3, 0.6, 1.0]) + + # Solve ODE + sol = odeint( + conceptual_reservoir_flux_calc, + y0, + t, + args=( + reservoir['storage_threshold_primary_m'], + reservoir['storage_max_m'], + reservoir['coeff_primary'], + reservoir['coeff_secondary'], + cfe_state.reduced_potential_et_m_per_timestep, + cfe_state.infiltration_depth_m, + cfe_state.soil_params['wltsmc'] + ), + tfirst=True, + Dfun=jac + ) + + # Finalize the results + ts_concat = t + ys_concat = np.concatenate(sol, axis=0) + + # Calculate fluxes (fine-tuning the residuals by scaling) + t_proportion = np.diff(ts_concat) + ys_avg = np.convolve(ys_concat, np.ones(2), 'valid') / 2 + + lateral_flux = np.zeros(ys_avg.shape) + perc_lat_switch = ys_avg - reservoir['storage_threshold_primary_m'] > 0 + lateral_flux[perc_lat_switch] = reservoir['coeff_secondary'] * np.minimum( + (ys_avg[perc_lat_switch] - reservoir['storage_threshold_primary_m']) / ( + reservoir['storage_max_m'] - reservoir['storage_threshold_primary_m']), 1) + lateral_flux_frac = lateral_flux * t_proportion + + perc_flux = np.zeros(ys_avg.shape) + perc_flux[perc_lat_switch] = reservoir['coeff_primary'] * np.minimum( + (ys_avg[perc_lat_switch] - reservoir['storage_threshold_primary_m']) / ( + reservoir['storage_max_m'] - reservoir['storage_threshold_primary_m']), 1) + perc_flux_frac = perc_flux * t_proportion + + et_from_soil = np.zeros(ys_avg.shape) + ET_switch = ys_avg - cfe_state.soil_params['wltsmc'] > 0 + et_from_soil[ET_switch] = cfe_state.reduced_potential_et_m_per_timestep * np.minimum( + (ys_avg[ET_switch] - cfe_state.soil_params['wltsmc']) / (reservoir['storage_threshold_primary_m'] - cfe_state.soil_params['wltsmc']), 1) + et_from_soil_frac = et_from_soil * t_proportion + + infilt_to_soil = np.repeat(cfe_state.infiltration_depth_m, ys_avg.shape) + infilt_to_soil_frac = infilt_to_soil * t_proportion + + # Scale fluxes + sum_outflux = lateral_flux_frac + perc_flux_frac + et_from_soil_frac + if sum_outflux.any() == 0: + flux_scale = 0 + else: + flux_scale = np.zeros(infilt_to_soil_frac.shape) + flux_scale[sum_outflux != 0] = (np.diff(-ys_concat, axis=0)[sum_outflux != 0] + infilt_to_soil_frac[ + sum_outflux != 0]) / sum_outflux[sum_outflux != 0] + flux_scale[sum_outflux == 0] = 0 + scaled_lateral_flux = lateral_flux_frac * flux_scale + scaled_perc_flux = perc_flux_frac * flux_scale + scaled_et_flux = et_from_soil_frac * flux_scale + + # Return the results + cfe_state.primary_flux_m = math.fsum(scaled_perc_flux) + cfe_state.secondary_flux_m = math.fsum(scaled_lateral_flux) + cfe_state.actual_et_from_soil_m_per_timestep = math.fsum(scaled_et_flux) + reservoir['storage_m'] = ys_concat[-1] + # __________________________________________________________________________________________________________ ########## SINGLE OUTLET EXPONENTIAL RESERVOIR ############### ########## -or- ############### @@ -333,36 +454,35 @@ def Schaake_partitioning_scheme(self,cfe_state): # __________________________________________________________________________________________________________ - def et_from_soil(self,cfe_state): + def et_from_soil(self, cfe_state): """ take AET from soil moisture storage, using Budyko type curve to limit PET if wilting 0: - - print("this should not happen yet. Still debugging the other functions.") - + + if cfe_state.reduced_potential_et_m_per_timestep > 0: if cfe_state.soil_reservoir['storage_m'] >= cfe_state.soil_reservoir['storage_threshold_primary_m']: - - cfe_state.actual_et_m_per_timestep = np.min(cfe_state.potential_et_m_per_timestep, - cfe_state.soil_reservoir['storage_m']) - cfe_state.soil_reservoir['storage_m'] -= cfe_state.actual_et_m_per_timestep + cfe_state.actual_et_from_soil_m_per_timestep = np.minimum(cfe_state.reduced_potential_et_m_per_timestep, + cfe_state.soil_reservoir['storage_m']) - cfe_state.et_struct['potential_et_m_per_timestep'] = 0.0 - - elif (cfe_state.soil_reservoir['storage_m'] > cfe_state.soil_reservoir['wilting_point_m'] and + elif (cfe_state.soil_reservoir['storage_m'] > cfe_state.soil_params['wltsmc'] * cfe_state.soil_params[ + 'D'] and cfe_state.soil_reservoir['storage_m'] < cfe_state.soil_reservoir['storage_threshold_primary_m']): - - Budyko_numerator = cfe_state.soil_reservoir['storage_m'] - cfe_state.soil_reservoir['wilting_point_m'] + + Budyko_numerator = cfe_state.soil_reservoir['storage_m'] - cfe_state.soil_params['wltsmc'] * \ + cfe_state.soil_params['D'] Budyko_denominator = cfe_state.soil_reservoir['storage_threshold_primary_m'] - \ - cfe_state.soil_reservoir['wilting_point_m'] - Budyki = Budyko_numerator / Budyko_denominator - - cfe_state.actual_et_m_per_timestep = Budyko * cfe_state.potential_et_m_per_timestep - - cfe_state.soil_reservoir['storage_m'] -= cfe_state.actual_et_m_per_timestep + cfe_state.soil_params['wltsmc'] * cfe_state.soil_params['D'] + Budyko = Budyko_numerator / Budyko_denominator + + cfe_state.actual_et_from_soil_m_per_timestep = np.minimum( + Budyko * cfe_state.reduced_potential_et_m_per_timestep, + cfe_state.soil_reservoir['storage_m']) + + cfe_state.soil_reservoir['storage_m'] -= cfe_state.actual_et_from_soil_m_per_timestep + cfe_state.reduced_potential_et_m_per_timestep -= cfe_state.actual_et_from_soil_m_per_timestep + return From e9923cd008b5354096447ecfe9848e75ea4cffa8 Mon Sep 17 00:00:00 2001 From: Ryoko Araki <400mhs2@gmail.com> Date: Thu, 22 Sep 2022 21:42:39 -0700 Subject: [PATCH 6/6] Fixed the water balance --- cfe.py | 107 +++++++++++++++++++-------------------------------------- 1 file changed, 35 insertions(+), 72 deletions(-) diff --git a/cfe.py b/cfe.py index c03ba7ac..5cdf4b44 100755 --- a/cfe.py +++ b/cfe.py @@ -6,6 +6,28 @@ def conceptual_reservoir_flux_calc(t, S, storage_threshold_primary_m, storage_max_m, coeff_primary, coeff_secondary, PET, infilt, wltsmc): # ODE of soil moisture reservoir mass balance + """ + Case 1: S (Soil moisture storage ) > storage_threshold_primary_m + Interpretation: When the soil moisture is plenty, AET(=PET), percolation, and lateral flow are all active. Thus, + Equation: dS/dt = Infiltration - PET - (Klf+Kperc) * (S - storage_threshold_primary_m)/(storage_max_m - storage_threshold_primary_m) + Case 2: storage_threshold_primary_m > S (Soil moisture storage) > storage_threshold_primary_m - wltsmc + Interpretation: When the soil moisture is in the middle range, AET is active and proportional to the soil moisture storage ratio + Equation: dS/dt = Infiltration - PET * (S - wltsmc)/(storage_threshold_primary_m - wltsmc) + Case 3: wltsmc > S (Soil moisture storage) + Interpretation: When the soil moisture is depleted, no outflux is active + Equation: dS/dt = Infitlation + :param t: + :param S: + :param storage_threshold_primary_m: + :param storage_max_m: maximum soil moisture storage, i.e., porosity + :param coeff_primary: K_perc, percolation coefficient + :param coeff_secondary: K_lf, lateral flow coefficient + :param PET: potential evapotranspiration + :param infilt: infiltration + :param wltsmc: wilting point (in meter) + :return: dS + """ + storage_above_threshold_m = S - storage_threshold_primary_m storage_diff = storage_max_m - storage_threshold_primary_m storage_ratio = np.minimum(storage_above_threshold_m / storage_diff, 1) @@ -65,54 +87,32 @@ def run_cfe(self, cfe_state): # ________________________________________________ # SUBROUTINE - # Calculates the value for surface_runoff_depth_m + # Calculates the value for surface_runoff_depth_m (infiltration excess overland flow) self.Schaake_partitioning_scheme(cfe_state) - - # ________________________________________________ - # SUBROUTINE - self.et_from_soil(cfe_state) + cfe_state.vol_sch_runoff += cfe_state.surface_runoff_depth_m # ________________________________________________ if cfe_state.soil_reservoir_storage_deficit_m < cfe_state.infiltration_depth_m: - # put won't fit back into runoff - cfe_state.surface_runoff_depth_m += (cfe_state.infiltration_depth_m - soil_reservoir_storage_deficit_m) + # Calculates saturation excess overland flow + cfe_state.surface_runoff_depth_m += (cfe_state.infiltration_depth_m - cfe_state.soil_reservoir_storage_deficit_m) cfe_state.infiltration_depth_m = cfe_state.soil_reservoir_storage_deficit_m - cfe_state.soil_reservoir['storage_m'] = cfe_state.soil_reservoir['storage_max_m'] - # ________________________________________________ - cfe_state.vol_sch_runoff += cfe_state.surface_runoff_depth_m cfe_state.vol_sch_infilt += cfe_state.infiltration_depth_m # ________________________________________________ - if cfe_state.current_time_step == 0: - cfe_state.previous_flux_perc_m = cfe_state.flux_perc_m - - # ________________________________________________ - if cfe_state.previous_flux_perc_m > cfe_state.soil_reservoir_storage_deficit_m: - self.diff_infilt = cfe_state.previous_flux_perc_m - cfe_state.soil_reservoir_storage_deficit_m - cfe_state.infiltration_depth_m = cfe_state.soil_reservoir_storage_deficit_m - cfe_state.vol_sch_runoff += self.diff_infilt - cfe_state.vol_sch_infilt -= self.diff_infilt - cfe_state.surface_runoff_depth_m += self.diff_infilt - - # ________________________________________________ - # SUBROUTINE # Solve ODE for the soil reservoir # Calculates primary_flux, secondary_flux, AET, and infiltration simultaneously cfe_state.actual_et_from_soil_m_per_timestep = 0 self.soil_reservoir_flux_calc(cfe_state, cfe_state.soil_reservoir) - # ________________________________________________ cfe_state.flux_perc_m = cfe_state.primary_flux_m cfe_state.flux_lat_m = cfe_state.secondary_flux_m - cfe_state.vol_to_soil += cfe_state.infiltration_depth_m cfe_state.vol_et_from_soil += cfe_state.actual_et_from_soil_m_per_timestep cfe_state.vol_et_to_atm += cfe_state.actual_et_from_soil_m_per_timestep cfe_state.volout += cfe_state.actual_et_from_soil_m_per_timestep # ________________________________________________ - cfe_state.gw_reservoir_storage_deficit_m = cfe_state.gw_reservoir['storage_max_m'] - cfe_state.gw_reservoir[ - 'storage_m'] + cfe_state.gw_reservoir_storage_deficit_m = cfe_state.gw_reservoir['storage_max_m'] - cfe_state.gw_reservoir['storage_m'] # ________________________________________________ if cfe_state.flux_perc_m > cfe_state.gw_reservoir_storage_deficit_m: @@ -123,16 +123,12 @@ def run_cfe(self, cfe_state): cfe_state.vol_sch_infilt -= self.diff_perc # ________________________________________________ - cfe_state.vol_to_gw += cfe_state.flux_perc_m - cfe_state.vol_soil_to_gw += cfe_state.flux_perc_m - - cfe_state.gw_reservoir['storage_m'] += cfe_state.flux_perc_m - cfe_state.soil_reservoir['storage_m'] -= cfe_state.flux_perc_m - cfe_state.soil_reservoir['storage_m'] -= cfe_state.flux_lat_m - cfe_state.vol_soil_to_lat_flow += cfe_state.flux_lat_m #TODO add this to nash cascade as input - cfe_state.volout += cfe_state.flux_lat_m + cfe_state.gw_reservoir['storage_m'] += cfe_state.flux_perc_m + cfe_state.vol_to_gw += cfe_state.flux_perc_m + cfe_state.vol_soil_to_gw += cfe_state.flux_perc_m + cfe_state.vol_soil_to_lat_flow += cfe_state.flux_lat_m #TODO add this to nash cascade as input + cfe_state.volout += cfe_state.flux_lat_m - # ________________________________________________ # SUBROUTINE # primary_flux, secondary_flux = f(reservoir) @@ -284,7 +280,7 @@ def soil_reservoir_flux_calc(self, cfe_state, reservoir): reservoir['coeff_secondary'], cfe_state.reduced_potential_et_m_per_timestep, cfe_state.infiltration_depth_m, - cfe_state.soil_params['wltsmc'] + cfe_state.soil_params['wltsmc'] * cfe_state.soil_params['D'] ), tfirst=True, Dfun=jac @@ -312,9 +308,9 @@ def soil_reservoir_flux_calc(self, cfe_state, reservoir): perc_flux_frac = perc_flux * t_proportion et_from_soil = np.zeros(ys_avg.shape) - ET_switch = ys_avg - cfe_state.soil_params['wltsmc'] > 0 + ET_switch = ys_avg - cfe_state.soil_params['wltsmc'] * cfe_state.soil_params['D']> 0 et_from_soil[ET_switch] = cfe_state.reduced_potential_et_m_per_timestep * np.minimum( - (ys_avg[ET_switch] - cfe_state.soil_params['wltsmc']) / (reservoir['storage_threshold_primary_m'] - cfe_state.soil_params['wltsmc']), 1) + (ys_avg[ET_switch] - cfe_state.soil_params['wltsmc']* cfe_state.soil_params['D']) / (reservoir['storage_threshold_primary_m'] - cfe_state.soil_params['wltsmc']* cfe_state.soil_params['D']), 1) et_from_soil_frac = et_from_soil * t_proportion infilt_to_soil = np.repeat(cfe_state.infiltration_depth_m, ys_avg.shape) @@ -451,40 +447,7 @@ def Schaake_partitioning_scheme(self,cfe_state): cfe_state.infiltration_depth_m = 0.0 return - - - # __________________________________________________________________________________________________________ - def et_from_soil(self, cfe_state): - """ - take AET from soil moisture storage, - using Budyko type curve to limit PET if wilting 0: - if cfe_state.soil_reservoir['storage_m'] >= cfe_state.soil_reservoir['storage_threshold_primary_m']: - cfe_state.actual_et_from_soil_m_per_timestep = np.minimum(cfe_state.reduced_potential_et_m_per_timestep, - cfe_state.soil_reservoir['storage_m']) - - elif (cfe_state.soil_reservoir['storage_m'] > cfe_state.soil_params['wltsmc'] * cfe_state.soil_params[ - 'D'] and - cfe_state.soil_reservoir['storage_m'] < cfe_state.soil_reservoir['storage_threshold_primary_m']): - - Budyko_numerator = cfe_state.soil_reservoir['storage_m'] - cfe_state.soil_params['wltsmc'] * \ - cfe_state.soil_params['D'] - Budyko_denominator = cfe_state.soil_reservoir['storage_threshold_primary_m'] - \ - cfe_state.soil_params['wltsmc'] * cfe_state.soil_params['D'] - Budyko = Budyko_numerator / Budyko_denominator - - cfe_state.actual_et_from_soil_m_per_timestep = np.minimum( - Budyko * cfe_state.reduced_potential_et_m_per_timestep, - cfe_state.soil_reservoir['storage_m']) - - cfe_state.soil_reservoir['storage_m'] -= cfe_state.actual_et_from_soil_m_per_timestep - cfe_state.reduced_potential_et_m_per_timestep -= cfe_state.actual_et_from_soil_m_per_timestep - - return - # __________________________________________________________________________________________________________ def is_fabs_less_than_epsilon(self,a,epsilon):