diff --git a/run_model_on_snellius/exe/STEMMUS_SCOPE b/run_model_on_snellius/exe/STEMMUS_SCOPE index f52bc8d1..93f93e8b 100755 Binary files a/run_model_on_snellius/exe/STEMMUS_SCOPE and b/run_model_on_snellius/exe/STEMMUS_SCOPE differ diff --git a/src/+groundwater/updateDunnianRunoff.m b/src/+groundwater/updateDunnianRunoff.m deleted file mode 100644 index 62557307..00000000 --- a/src/+groundwater/updateDunnianRunoff.m +++ /dev/null @@ -1,14 +0,0 @@ -function [runoffDunnian, update_Precip_msr] = updateDunnianRunoff(Precip_msr, groundWaterDepth) - % Dunnian runoff = Direct water input from precipitation + return flow - % (a) direct water input from precipitation when soil is fully saturated (depth to water table = 0) - % (b) Return flow (from groundwater exfiltration) calculated in MODFLOW and added to Dunnian runoff (through BMI) - % here approach (a) is implemented - runoffDunnian = zeros(size(Precip_msr)); - update_Precip_msr = Precip_msr; - - if groundWaterDepth <= 1.0 - runoffDunnian = Precip_msr; - update_Precip_msr = zeros(size(Precip_msr)); - end - -end diff --git a/src/+io/bin_to_csv.m b/src/+io/bin_to_csv.m index 5bc13136..70086824 100644 --- a/src/+io/bin_to_csv.m +++ b/src/+io/bin_to_csv.m @@ -12,6 +12,15 @@ function bin_to_csv(fnames, n_col, ns, options, SoilLayer, GroundwaterSettings, 'umol m-2 s-1', ' umol m-2 s-1', 'umol m-2 s-1', 'W m-2', 'umol m-2 s-1', 'W m-2', 'W m-2', 'mm s-1', 'mm s-1', 'mm s-1', 'Kg m-2 s-1', 'Kg m-2 s-1'}; write_output(flu_names, flu_units, fnames.flu_file, n_col.flu, ns); + %% Water balance fluxes + wbal_names = {'simulation_number', 'nu_iterations', 'year', 'DoY', 'Precip', 'Precip_liquid', 'Precip_snow', 'Precip_snowmelt', 'Precip_totalLiquid', ... + 'Eff_Precip', 'R_Hort', 'R_Dunn', 'Runoff', 'RWUs', 'RWUg', 'Trap', 'Evap', 'ET', 'botmFlux', 'botmFluxIn', 'botmFluxOut', 'deltaS', 'deltaS_in', ... + 'deltaS_out', 'correctedDeltaS', 'correctedDeltaS_in', 'correctedDeltaS_out', 'Total_inflow', 'Total_outflow', 'Residual', 'Error_org', 'Error_corrected'}; + wbal_units = {'', '', '', '', 'mm/30min', 'mm/30min', 'mm/30min', 'mm/30min', 'mm/30min', 'mm/30min', 'mm/30min', 'mm/30min', 'mm/30min', ... + 'mm/30min', 'mm/30min', 'mm/30min', 'mm/30min', 'mm/30min', 'mm/30min', 'mm/30min', 'mm/30min', 'mm/30min', 'mm/30min', ... + 'mm/30min', 'mm/30min', 'mm/30min', 'mm/30min', 'mm/30min', 'mm/30min', 'mm/30min', 'percentage', 'percentage'}; + write_output(wbal_names, wbal_units, fnames.wbal_file, n_col.wbal, ns); + %% surftemp surftemp_names = {'simulation_number', 'year', 'DoY', ... 'Ta', 'Ts(1)', 'Ts(2)', 'Tcave', 'Tsave'}; diff --git a/src/+io/constrainSoilVariables.m b/src/+io/constrainSoilVariables.m new file mode 100755 index 00000000..fae4008d --- /dev/null +++ b/src/+io/constrainSoilVariables.m @@ -0,0 +1,32 @@ +function [Theta] = constrainSoilVariables(Theta, VanGenuchten) + % Added by Prajwal and Mostafa + % Constrain soil moisture values within Van Genuchten bounds + % + % Args: + % Theta 2D array (e.g. SoilVariables.Theta_LL or SoilVariables.Theta_UU) + % VanGenuchten Structure with Theta_r (residual) and Theta_s (saturated) bounds + % + % Returns: + % SoilVariables: Updated structure with constrained mositure (Theta) + + % Input validation + if ~ismatrix(Theta) && ~ndims(Theta) == 2 + error('Theta must be a 2D matrix'); + end + + if ~isstruct(VanGenuchten) || ~isfield(VanGenuchten, 'Theta_r') || ~isfield(VanGenuchten, 'Theta_s') + error('VanGenuchten must be a structure with Theta_r and Theta_s fields'); + end + + % Apply constraints using vectorized operations + smConstrained = Theta; + smLowerBound = zeros(size(Theta)); + smLowerBound(1:end - 1, 1) = VanGenuchten.Theta_r' + 0.00000001; + smLowerBound(1:end - 1, 2) = VanGenuchten.Theta_r' + 0.00000001; + smUpperBound = zeros(size(Theta)); + smUpperBound(1:end - 1, 1) = VanGenuchten.Theta_s' - 0.00000001; + smUpperBound(1:end - 1, 2) = VanGenuchten.Theta_s' - 0.00000001; + smConstrained = max(min(smConstrained, smUpperBound), smLowerBound); + % Update soil mositure array with constrained boundaries + Theta = smConstrained; +end diff --git a/src/+io/create_output_files_binary.m b/src/+io/create_output_files_binary.m index 76695055..b94a0e4e 100644 --- a/src/+io/create_output_files_binary.m +++ b/src/+io/create_output_files_binary.m @@ -30,13 +30,14 @@ fnames.waterStressFactor_file = fullfile(Output_dir, 'waterStressFactor.bin'); fnames.waterPotential_file = fullfile(Output_dir, 'waterPotential.bin'); % leaf water potential fnames.Sim_hh_file = fullfile(Output_dir, 'Sim_hh.bin'); % soil matric potential - fnames.Sim_qlh_file = fullfile(Output_dir, 'qlh.bin'); % liquid flux due to matric potential gradient - fnames.Sim_qlt_file = fullfile(Output_dir, 'qlt.bin'); % liquid flux due to temprature gradient - fnames.Sim_qla_file = fullfile(Output_dir, 'qla.bin'); % liquid flux due to dry air pressure gradient - fnames.Sim_qvh_file = fullfile(Output_dir, 'qvh.bin'); % vapour flux due to matric potential gradient - fnames.Sim_qvt_file = fullfile(Output_dir, 'qvt.bin'); % vapour flux due to temprature gradient - fnames.Sim_qva_file = fullfile(Output_dir, 'qva.bin'); % vapour flux due to dry air pressure gradient - fnames.Sim_qtot_file = fullfile(Output_dir, 'qtot.bin'); % total flux (liquid + vapour) + fnames.Sim_qlh_file = fullfile(Output_dir, 'qLiquid_head.bin'); % liquid flux due to matric potential gradient + fnames.Sim_qlt_file = fullfile(Output_dir, 'qLiquid_temp.bin'); % liquid flux due to temprature gradient + fnames.Sim_qla_file = fullfile(Output_dir, 'qLiquid_dryair.bin'); % liquid flux due to dry air pressure gradient + fnames.Sim_qvh_file = fullfile(Output_dir, 'qVapour_head.bin'); % vapour flux due to matric potential gradient + fnames.Sim_qvt_file = fullfile(Output_dir, 'qVapour_temp.bin'); % vapour flux due to temprature gradient + fnames.Sim_qva_file = fullfile(Output_dir, 'qVapour_dryair.bin'); % vapour flux due to dry air pressure gradient + fnames.Sim_qtot_file = fullfile(Output_dir, 'qTotal.bin'); % total flux (liquid + vapour) + fnames.wbal_file = fullfile(Output_dir, 'waterBalance.bin'); % water balance fluxes and error if options.calc_ebal fnames.spectrum_obsdir_BlackBody_file = fullfile(Output_dir, 'spectrum_obsdir_BlackBody.bin'); % spectrum observation direction diff --git a/src/+io/loadForcingData.m b/src/+io/loadForcingData.m index c1dbdf87..a2115e55 100644 --- a/src/+io/loadForcingData.m +++ b/src/+io/loadForcingData.m @@ -19,17 +19,6 @@ ForcingData.G_msr = Mdata{:, 7} * 0.15; Precip_msr = Mdata{:, 6}; % (cm/sec) - % Calculate saturation excess runoff (Dunnian runoff) - if ~GroundwaterSettings.GroundwaterCoupling % Groundwater Coupling is not activated - % Concept is adopted from the CLM model (see issue 232, https://github.com/EcoExtreML/STEMMUS_SCOPE/issues/232) - % Check also the CLM documents (https://doi.org/10.5065/D6N877R0, https://doi.org/10.1029/2005JD006111) - wat_Dep = Tot_Depth / 100; % (m), this assumption water depth = total soil depth is not fully correct (to be improved) - fover = 0.5; % decay factor (fixed to 0.5 m-1) - fmax = SoilProperties.fmax; % potential maximum value of fsat - fsat = (fmax .* exp(-0.5 * fover * wat_Dep)); % fraction of saturated area (unitless) - ForcingData.runoffDunnian = Precip_msr .* fsat; % Dunnian runoff (saturation excess runoff, in cm/sec) - end % In case Groundwater Coupling is activated, Dunnian runoff is calculated in +groundwater/updateDunnianRunoff - % replace negative values for jj = 1:Dur_tot if ForcingData.Ta_msr(jj) < -100 diff --git a/src/+io/output_data_binary.m b/src/+io/output_data_binary.m index d5a918d4..2242085c 100644 --- a/src/+io/output_data_binary.m +++ b/src/+io/output_data_binary.m @@ -1,7 +1,7 @@ function n_col = output_data_binary(f, k, xyt, rad, canopy, ScopeParameters, vi, vmax, options, fluxes, meteo, iter, thermal, ... spectral, gap, profiles, Sim_Theta_U, Sim_Temp, Trap, Evap, WaterStress, WaterPotential, ... Sim_hh, Sim_qlh, Sim_qlt, Sim_qvh, Sim_qvt, Sim_qla, Sim_qva, Sim_qtot, ... - ForcingData, RS, RWUs, RWUg) + ForcingData, RS, RWUs, RWUg, wbal) %% OUTPUT DATA % author C. Van der Tol @@ -14,6 +14,17 @@ n_col.flu = length(flu_out); fwrite(f.flu_file, flu_out, 'double'); + %% Water balance fluxes + unitConv = 10 * 1800; % convert fluxes unit from cm/s to mm/30min + wbal_out = [k iter.counter xyt.year(k) xyt.t(k) ForcingData.Precip * unitConv ForcingData.Precip_liquid * unitConv ForcingData.Precip_snow * unitConv ... + ForcingData.Precip_snowmelt * unitConv ForcingData.Precip_totalLiquid * unitConv ForcingData.effectivePrecip * unitConv ForcingData.runoffHort * unitConv ... + ForcingData.runoffDunn * unitConv ForcingData.runoff * unitConv RWUs * unitConv RWUg * unitConv Trap * unitConv Evap * unitConv ... + (Trap + Evap) * unitConv wbal.botmFlux * unitConv wbal.botmFluxIn * unitConv wbal.botmFluxOut * unitConv wbal.deltaStorage * unitConv ... + wbal.deltaStorageIn * unitConv wbal.deltaStorageOut * unitConv wbal.correctedDeltaS * unitConv wbal.correctedDeltaSIn * unitConv ... + wbal.correctedDeltaSOut * unitConv wbal.totalInflow * unitConv wbal.totalOutflow * unitConv wbal.residual * unitConv wbal.errorInit wbal.error]; + n_col.wbal = length(wbal_out); + fwrite(f.wbal_file, wbal_out, 'double'); + %% surftemp surftemp_out = [k xyt.year(k) xyt.t(k) thermal.Ta thermal.Ts(1) thermal.Ts(2) thermal.Tcave thermal.Tsave]; n_col.surftemp = length(surftemp_out); diff --git a/src/+io/read_config.m b/src/+io/read_config.m index 60f33185..367042e4 100644 --- a/src/+io/read_config.m +++ b/src/+io/read_config.m @@ -1,4 +1,4 @@ -function [InputPath, OutputPath, InitialConditionPath, FullCSVfiles] = read_config(config_file) +function [InputPath, OutputPath, InitialConditionPath, FullCSVfiles, closeWaterBalance] = read_config(config_file) file_id = fopen(config_file); config = textscan(file_id, '%s %s', 'HeaderLines', 0, 'Delimiter', '='); @@ -24,3 +24,10 @@ else FullCSVfiles = str2num(config_paths{indx}); end + + indx = find(strcmp(config_vars, 'closeWaterBalance')); + if isempty(indx) + closeWaterBalance = 0; + else + closeWaterBalance = str2num(config_paths{indx}); + end diff --git a/src/+soilmoisture/calculateBoundaryConditions.m b/src/+soilmoisture/calculateBoundaryConditions.m index 8e8264b2..b6a0b690 100644 --- a/src/+soilmoisture/calculateBoundaryConditions.m +++ b/src/+soilmoisture/calculateBoundaryConditions.m @@ -1,14 +1,14 @@ function [AVAIL0, RHS, HeatMatrices, ForcingData] = calculateBoundaryConditions(BoundaryCondition, HeatMatrices, ForcingData, SoilVariables, InitialValues, ... TimeProperties, SoilProperties, RHS, hN, KT, KIT, Delt_t, Evap, ModelSettings, GroundwaterSettings) %{ - Determine the boundary condition for solving the soil moisture equation. + Determine boundary conditions for solving the soil moisture equation. %} n = ModelSettings.NN; C4 = HeatMatrices.C4; C4_a = HeatMatrices.C4_a; - %%%%%% Apply the bottom boundary condition called for by BoundaryCondition.NBChB %%%%%% + %%%%%% Apply the bottom boundary condition called for by BoundaryCondition.NBChB, modified by Mostafa %%%%%% if ~GroundwaterSettings.GroundwaterCoupling % no Groundwater coupling if BoundaryCondition.NBChB == 1 % Specify matric head at bottom to be ---BoundaryCondition.BChB; RHS(1) = BoundaryCondition.BChB; @@ -21,7 +21,7 @@ elseif BoundaryCondition.NBChB == 3 % BoundaryCondition.NBChB=3, Gravity drainage at bottom--specify flux= hydraulic conductivity; RHS(1) = RHS(1) - SoilVariables.KL_h(1, 1); end - else % Groundwater coupling is activated, added by Mostafa + else % Groundwater coupling is activated indxBotmLayer_R = GroundwaterSettings.indxBotmLayer_R; indxBotm = GroundwaterSettings.indxBotmLayer; % index of bottom boundary layer after neglecting the saturated layers (from bottom to top) soilThick = GroundwaterSettings.soilThick; % cumulative soil layers thickness @@ -40,45 +40,93 @@ end end - %%%%%% Apply the surface boundary condition called for by BoundaryCondition.NBCh %%%%%% + %%%%%% Prepare surface boundary conditions (precipitation and runoff) for BoundaryCondition.NBCh, modified by Mostafa %%%%%% Precip = ForcingData.Precip_msr(KT); % total precipitation (liquid + snow) - runoffDunn = ForcingData.runoffDunnian(KT); % Dunnian runoff (calculated in +io/loadForcingData file) - % Check if surface temperature is less than zero, then Precipitation is snow (modified by Mostafa) - if KT == 1 % see issue 279 (https://github.com/EcoExtreML/STEMMUS_SCOPE/issues/279) - Precip_snowAccum = 0; % initalize accumulated snow for first time step + % Check if surface temperature is <= 0 C; if so, precipitation is snow + % see issue 279 (https://github.com/EcoExtreML/STEMMUS_SCOPE/issues/279) + % Initialize new variables + Precip_snowmelt = 0; % snowmelt + Precip_liquid = 0; % only rainfall + Precip_totalLiquid = 0; % total liquid precipitation (rainfall + snowmelt) + if KT == 1 + Precip_snowAccum = 0; % accumulated snow at first time step else - Precip_snowAccum = ForcingData.Precip_snowAccum; + Precip_snowAccum = ForcingData.Precip_snowAccum; % accumulated snow from previous time step end if SoilVariables.Tss(KT) <= 0 % surface temperature is equal or less than zero Precip_snow = Precip; % snow precipitation - Precip_liquid = 0; % liquid precipitation (rainfall) - runoffDunn = 0; % update Dunnian runoff in case precpitation is snow - if KIT == 1 % accumulate snow at one iteration only within the time step + Precip_liquid = 0; + Precip_snowmelt = 0; + Precip_totalLiquid = 0; + if KIT == 1 % accumulate snow at one iteration only per time step Precip_snowAccum = Precip + Precip_snowAccum; - else - Precip_snowAccum = Precip_snowAccum; end - else % surface temperature is more than zero - if KIT == 1 % add accumulated snow of previous time steps to liquid precipitation at first time step when surface temperature > zero - Precip_liquid = Precip + Precip_snowAccum; - Precip_snowAccum = 0; - else + else % surface temperature > 0 C + if KIT == 1 % all accumulated snow of previous time steps becomes snowmelt + Precip_snowmelt = Precip_snowAccum; + Precip_liquid = Precip; + Precip_totalLiquid = Precip_liquid + Precip_snowmelt; + Precip_snowAccum = 0; % reset accumulated snow for following iterations + else % retrieve values from the previous iteration of the same time step Precip_liquid = ForcingData.Precip_liquid; + Precip_snowmelt = ForcingData.Precip_snowmelt; + Precip_totalLiquid = ForcingData.Precip_totalLiquid; end Precip_snow = 0; end %%% Calculate effective precipitation after removing canopy interception and total runoff %%% % effective precipitation = precipitation - canopy interception - (Dunnian runoff + Hortonian runoff) - % Currently, canopy interception is not implemented in the code yet - % (1) Remove saturation excess runoff (Dunnian runoff) - effectivePrecip = Precip_liquid - runoffDunn; % Hortonian runoff is removed below + % Note: currently canopy interception is not implemented + % 1.1. Calculate saturation excess runoff (Dunnian runoff) + if ~GroundwaterSettings.GroundwaterCoupling % Groundwater Coupling is not activated + % Concept is adopted from the CLM model (see issue 232, https://github.com/EcoExtreML/STEMMUS_SCOPE/issues/232) + % Check also the CLM documents (https://doi.org/10.5065/D6N877R0, https://doi.org/10.1029/2005JD006111) + wat_Dep = ModelSettings.Tot_Depth / 100; % (m), this assumption water depth = total soil depth is not fully correct (to be improved) + fover = 0.5; % decay factor (fixed to 0.5 m-1) + fmax = SoilProperties.fmax; % potential maximum value of fsat + fsat = (fmax .* exp(-0.5 * fover * wat_Dep)); % fraction of saturated area (unitless) + runoffDunn = Precip_totalLiquid * fsat; % Dunnian runoff (saturation excess runoff, in cm/sec) + else % Groundwater Coupling is activated + % Dunnian runoff = Direct water input from precipitation + return flow + % (a) Direct water input from precipitation when soil is fully saturated (water table depth = 0) + % (b) Return flow from groundwater exfiltration, calculated in MODFLOW and added through BMI + % Here approach (a) is implemented + runoffDunn = 0; % Dunnian runoff = 0 when soil is not fully saturated + if GroundwaterSettings.gw_Dep <= 1.0 + runoffDunn = Precip_totalLiquid; + end + end + % 1.2. Remove saturation excess runoff (Dunnian runoff) from precipitation + effectivePrecip = Precip_totalLiquid - runoffDunn; % Hortonian runoff is removed below + + % 2.1. Calculate infiltration excess runoff (Hortonian runoff) and update effective precipitation + Ks0 = SoilProperties.Ks0 / (3600 * 24); % saturated vertical hydraulic conductivity. unit conversion from cm/day to cm/sec + % Note: Ks0 is not adjusted by the fsat as in the CLM model (Check CLM document: https://doi.org/10.5065/D6N877R) + topThick = 5; % first 5 cm of the soil + satCap = SoilProperties.theta_s0 * topThick; % saturation capacity represented by saturated water content of the top 5 cm of the soil + actTheta = ModelSettings.DeltZ(end - 3:end) * SoilVariables.Theta_UU(end - 4:end - 1, 1); % actual moisture of the top 5 cm of the soil + infCap = (satCap - actTheta) / TimeProperties.DELT; % % infiltration capacity (cm/sec) + infCap_min = min(Ks0, infCap); % minimum infiltration capacity + + if effectivePrecip > infCap_min + runoffHort = effectivePrecip - infCap_min; % Hortonian runoff + else + runoffHort = 0; + end + % 2.2. Update effective precipitation after removing Hortonian runoff + effectivePrecip = min(effectivePrecip, infCap_min); + + % 2.3. Add depression water to effective precipitation + AVAIL0 = effectivePrecip + BoundaryCondition.DSTOR0 / Delt_t; % (cm/sec) + + %%%%%% Apply the bottom boundary condition called for by BoundaryCondition.NBChB %%%%%% if BoundaryCondition.NBCh == 1 % Specified matric head at surface---equal to hN; - % h_SUR: Observed matric potential at surface. This variable - % is not calculated anywhere! see issue 98, item 6 + % h_SUR: Observed matric potential at surface. + % Note: this variable is not calculated anywhere (see issue 98, item 6) RHS(n) = InitialValues.h_SUR(KT); C4(n, 1) = 1; RHS(n - 1) = RHS(n - 1) - C4(n - 1, 2) * RHS(n); @@ -91,29 +139,9 @@ RHS(n - 1) = RHS(n - 1) - C4(n - 1, 2) * RHS(n); C4(n - 1, 2) = 0; else - RHS(n) = RHS(n) - BoundaryCondition.BCh; % a specified matric head (saturation or dryness) was applied; + RHS(n) = RHS(n) - BoundaryCondition.BCh; % a specified matric head (saturation or dryness) is applied; end else % (BoundaryCondition.NBCh == 3, Specified atmospheric forcing) - % (2) Calculate infiltration excess runoff (Hortonian runoff) and update effective precpitation, modified by Mostafa - Ks0 = SoilProperties.Ks0 / (3600 * 24); % saturated vertical hydraulic conductivity. unit conversion from cm/day to cm/sec - % Note: Ks0 is not adjusted by the fsat as in the CLM model (Check CLM document: https://doi.org/10.5065/D6N877R) - topThick = 5; % first 5 cm of the soil - satCap = SoilProperties.theta_s0 * topThick; % saturation capacity represented by saturated water content of the top 5 cm of the soil - actTheta = ModelSettings.DeltZ(end - 3:end) * SoilVariables.Theta_UU(end - 4:end - 1, 1); % actual moisture of the top 5 cm of the soil - infCap = (satCap - actTheta) / TimeProperties.DELT; % infiltration capcaity (cm/sec) - infCap_min = min(Ks0, infCap); % minimum infiltration capcaity - - if effectivePrecip > infCap_min - runoffHort = effectivePrecip - infCap_min; % Hortonian runoff - else - runoffHort = 0; - end - % Update effective precipitation after removing Hortonian runoff - effectivePrecip = min(effectivePrecip, infCap_min); - - % Add depression water to effective precipitation - AVAIL0 = effectivePrecip + BoundaryCondition.DSTOR0 / Delt_t; % (cm/sec) - if BoundaryCondition.NBChh == 1 RHS(n) = hN; C4(n, 1) = 1; @@ -132,6 +160,8 @@ ForcingData.Precip_liquid = Precip_liquid; ForcingData.Precip_snow = Precip_snow; ForcingData.Precip_snowAccum = Precip_snowAccum; + ForcingData.Precip_snowmelt = Precip_snowmelt; + ForcingData.Precip_totalLiquid = Precip_totalLiquid; ForcingData.effectivePrecip = effectivePrecip; ForcingData.runoffDunn = runoffDunn; ForcingData.runoffHort = runoffHort; diff --git a/src/+soilmoisture/calculateEvapotranspiration.m b/src/+soilmoisture/calculateEvapotranspiration.m index c865ef66..a7a69fba 100644 --- a/src/+soilmoisture/calculateEvapotranspiration.m +++ b/src/+soilmoisture/calculateEvapotranspiration.m @@ -5,6 +5,9 @@ else % Groundwater Coupling is activated % index of bottom layer after neglecting saturated layers (from bottom to top) indxBotm = GroundwaterSettings.indxBotmLayer; + if indxBotm >= ModelSettings.NL - 2 % avoid model crash in case of fully saturated soil + indxBotm = ModelSettings.NL - 2; + end end Rn = (ForcingData.Rn_msr(KT)) * 8.64 / 24 / 100 * 1; @@ -39,7 +42,7 @@ RWU = RWU * 100; % unit conversion from m/sec to cm/sec if GroundwaterSettings.GroundwaterCoupling == 1 RWUs = sum(RWU(indxBotm:ModelSettings.NL)); - RWUg = sum(RWU(1:indxBotm - 1)); + RWUg = sum(RWU(1:indxBotm + 2)); % add +2 (2 layers to account for capillary frinage zone) else RWUs = sum(RWU(1:ModelSettings.NL)); RWUg = 0; diff --git a/src/STEMMUS_SCOPE.m b/src/STEMMUS_SCOPE.m index 4c2edbc9..3358c064 100644 --- a/src/STEMMUS_SCOPE.m +++ b/src/STEMMUS_SCOPE.m @@ -39,7 +39,7 @@ if strcmp(bmiMode, "initialize") || strcmp(runMode, "full") % Read the configPath file. Due to using MATLAB compiler, we cannot use run(CFG) disp (['Reading config from ', CFG]); - [InputPath, OutputPath, InitialConditionPath, FullCSVfiles] = io.read_config(CFG); + [InputPath, OutputPath, InitialConditionPath, FullCSVfiles, closeWaterBalance] = io.read_config(CFG); % Prepare forcing and soil data [SiteProperties, SoilProperties, TimeProperties] = io.prepareInputData(InputPath); @@ -242,7 +242,7 @@ atmfile = [path_input 'radiationdata/' char(F(4).FileName(1))]; atmo.M = helpers.aggreg(atmfile, spectral.SCOPEspec); - %% 13. create output files and + %% 13. create output files [Output_dir, fnames] = io.create_output_files_binary(parameter_file, SiteProperties.sitename, path_of_code, path_input, path_output, spectral, options); %% Initialize Temperature, Matric potential and soil air pressure. @@ -353,9 +353,6 @@ % Update GroundwaterSettings.headBotmLayer and GroundwaterSettings.tempBotm, from MODFLOW through BMI GroundwaterSettings.gw_Dep = groundwater.calculateGroundWaterDepth(GroundwaterSettings.topLevel, GroundwaterSettings.headBotmLayer, ModelSettings.Tot_Depth); - % Update Dunnian runoff and ForcingData.Precip_msr - [ForcingData.runoffDunnian, ForcingData.Precip_msr] = groundwater.updateDunnianRunoff(ForcingData.Precip_msr, GroundwaterSettings.gw_Dep); - % Calculate the index of the bottom layer level [GroundwaterSettings.indxBotmLayer, GroundwaterSettings.indxBotmLayer_R] = groundwater.calculateIndexBottomLayer(GroundwaterSettings.soilThick, GroundwaterSettings.gw_Dep, ModelSettings); @@ -621,6 +618,7 @@ SoilVariables.T = T; SoilVariables.h_frez = h_frez; SoilVariables.Tss(KT) = Tss; + Theta_UU_previous = SoilVariables.Theta_UU; for KIT = 1:ModelSettings.NIT [TT_CRIT, hh_frez] = HT_frez(SoilVariables.hh, ModelSettings.T0, Constants.g, L_f, SoilVariables.TT, ModelSettings.NN, SoilConstants.hd, ForcingData.Tmin); @@ -704,7 +702,22 @@ TOLD, Evap, r_a_SOIL, Rn_SOIL, KT, CHK, ModelSettings, GroundwaterSettings); end - if max(CHK) < 0.1 + % Recharge calculations, added by Mostafa + if GroundwaterSettings.GroundwaterCoupling == 1 % Groundwater coupling is enabled + gwfluxes = groundwater.calculateGroundwaterRecharge(EnergyVariables, SoilVariables, KT, gwfluxes, ModelSettings, GroundwaterSettings); + if GroundwaterSettings.gw_Dep <= 1 % soil is fully saturated + gwfluxes.recharge = 0; + end + else + gwfluxes.recharge = 0; + gwfluxes.indxRchrg = NaN; + end + + % Water balance calculations, added by Mostafa + [wbal, Theta_UU_corrected] = calculateWaterBalance(ForcingData, SoilVariables, TimeProperties, KT, Theta_UU_previous, VanGenuchten, Evap, Trap, gwfluxes, HBoundaryFlux, ModelSettings, GroundwaterSettings, closeWaterBalance); + SoilVariables.Theta_UU = Theta_UU_corrected; + + if max(CHK) < 0.1 && abs(wbal.error) < 1 break end hSAVE = SoilVariables.hh(ModelSettings.NN); @@ -716,11 +729,17 @@ KIT = 0; [TT_CRIT, hh_frez] = HT_frez(SoilVariables.hh, ModelSettings.T0, Constants.g, L_f, SoilVariables.TT, ModelSettings.NN, SoilConstants.hd, ForcingData.Tmin); - % updates inputs for UpdateSoilWaterContent + + % Update inputs for UpdateSoilWaterContent SoilVariables.TT_CRIT = TT_CRIT; SoilVariables.hh_frez = hh_frez; SoilVariables = UpdateSoilWaterContent(KIT, L_f, SoilVariables, VanGenuchten, ModelSettings); + % Update water balance (after UpdateSoilWaterContent) + [wbal, Theta_UU_corrected] = calculateWaterBalance(ForcingData, SoilVariables, TimeProperties, KT, Theta_UU_previous, VanGenuchten, Evap, Trap, gwfluxes, HBoundaryFlux, ModelSettings, GroundwaterSettings, closeWaterBalance); + SoilVariables.Theta_UU = Theta_UU_corrected; + Theta_UU_previous = SoilVariables.Theta_UU; % for next time step + if IRPT1 == 0 && IRPT2 == 0 if KT % In case last time step is not convergent and needs to be repeated. for i = 1:ModelSettings.NL @@ -775,17 +794,6 @@ end end - % Recharge calculations, added by Mostafa - if GroundwaterSettings.GroundwaterCoupling == 1 % Groundwater coupling is enabled - gwfluxes = groundwater.calculateGroundwaterRecharge(EnergyVariables, SoilVariables, KT, gwfluxes, ModelSettings, GroundwaterSettings); - if GroundwaterSettings.gw_Dep <= 1 % soil is fully saturated - gwfluxes.recharge = 0; - end - else - gwfluxes.recharge = 0; - gwfluxes.indxRchrg = NaN; - end - % set SoilVariables for the rest of the loop h = SoilVariables.h; hh = SoilVariables.hh; @@ -799,7 +807,7 @@ n_col = io.output_data_binary(file_ids, k, xyt, rad, canopy, ScopeParameters, vi, vmax, options, fluxes, ... meteo, iter, thermal, spectral, gap, profiles, Sim_Theta_U, Sim_Temp, Trap, ... Evap, WaterStressFactor, WaterPotential, Sim_hh, Sim_qlh, Sim_qlt, Sim_qvh, ... - Sim_qvt, Sim_qla, Sim_qva, Sim_qtot, ForcingData, RS, RWUs, RWUg); + Sim_qvt, Sim_qla, Sim_qva, Sim_qtot, ForcingData, RS, RWUs, RWUg, wbal); fclose("all"); end end diff --git a/src/calculateWaterBalance.m b/src/calculateWaterBalance.m new file mode 100755 index 00000000..1ebf93da --- /dev/null +++ b/src/calculateWaterBalance.m @@ -0,0 +1,176 @@ +function [wbal, Theta_UU_corrected] = calculateWaterBalance(ForcingData, SoilVariables, TimeProperties, KT, Theta_UU_previous, VanGenuchten, ... + Evap, Trap, gwfluxes, HBoundaryFlux, ModelSettings, GroundwaterSettings, correctWaterBalance) + %{ + Added by Mostafa + This function performs two tasks: 1) checks the original water balance of the unsaturated zone + 2) (optional) if water balance error is large -> correct water balance error + Outputs + wbal a structure that includes water balance fluxes + Theta_UU_corrected a 2D array of the corrected soil moisture profile + + Task 1: Check water balance + Water balance concept: Total inflow = Total outflow + residual + Total inflow = precipitation + Total outflow = evapotranspiration + runoff + bottom boundary flux + delta storage + Delta storage = (storage at time t+1 - storage at time t) / time step length + = (soil moisture(t+1) - soil moisture(t)) * soil thickness / time step length + Residual = total inflow - total outflow - delta storage + Error = (residual / total inflow) * 100% + + Task 2: Correct water balance + To close the water balance, one of the water fluxes (ET, runoff, bottom boundary flux, or delta storage) must be adjusted. + Since all fluxes depend on soil moisture, delta storage is chosen for correction. Why? + Because the adjustment is distributed across the entire soil moisture profile, so each soil layer + undergoes only a minor modification, minimizing potential disturbances. + This approach prevents large changes in other fluxes that could occur if they were corrected directly. + + %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% + ET One value Evapotranspiration [sum of Evaporation (Evap) + Transpiration (Trap)] + runoff One value Runoff [sum of Hortonian runoff + Dunnian runoff] + botmFluxIn One value Flux at bottom boundary of soil (input to unsaturated zone) + botmFluxOut One value Flux at bottom boundary of soil (output from unsaturated zone) + botmFlux One value Flux at bottom boundary of soil (difference between botmFluxIn and botmFluxOut) + Theta_UU_previous 2D array Soil moisture (liquid + ice) of the previous time step + Theta_UU_current 2D array Soil moisture (liquid + ice) of the current time step + Theta_UU_corrected 2D array Corrected soil moisture (liquid + ice) at the end of the current time step + deltaStorageIn One value Change in storage entering unsaturated zone (sum over all soil layers) + deltaStorageOut One value Change in storage leaving unsaturated zone (sum over all soil layers) + deltaStorage One value Net change in storage of the unsaturated zone (sum over all soil layers) + deltaStorageLay 1D array Change in storage of the unsaturated zone per soil layer + totalInflowInit One value Total inflow to unsaturated zone (before water balance correction) + totalOutflowInit One value Total outflow from unsaturated zone (before water balance correction) + residualInit One value Residual (total inflow - total outflow) + errorInit One value Error of water balance (before water balance correction) + thetaCorrection One value Corrected soil moisture value + thetaCorrectionLay 1D array Corrected soil moisture array + correctedDeltaSIn One value Corrected change in storage entering unsaturated zone (sum over all soil layers) + correctedDeltaSOut One value Corrected change in storage leaving unsaturated zone (sum over all soil layers) + correctedDeltaS One value Corrected net change in storage of the unsaturated zone + totalInflow One value Corrected total inflow to unsaturated zone (after water balance correction) + totalOutflow One value Corrected total outflow from unsaturated zone (after water balance correction) + error One value Error of water balance (after water balance correction) + %} + + %%%%%%%%%%%%%%%%%% Task 1: Check water balance error %%%%%%%%%%%%%%%%%% + % 1.1. Get water balance components (ET and runoff) + ET = Evap + Trap; + runoff = ForcingData.runoffHort + ForcingData.runoffDunn; + + % 1.2. Calculate delta storage + if KT == 1 % For the first time step only (for later steps, Theta_UU_previous is provided as input) + Theta_UU_previous = SoilVariables.Theta_UU; + end + Theta_UU_current = SoilVariables.Theta_UU; + Theta_UU_corrected = SoilVariables.Theta_UU; % will be corrected below + deltaStorageLay = (Theta_UU_current(1:end - 1, 1) - Theta_UU_previous(1:end - 1, 1))' .* ModelSettings.DeltZ; + deltaStorage = sum(deltaStorageLay) / TimeProperties.DELT; + deltaStorageIn = max(deltaStorage, 0); + deltaStorageOut = max(-deltaStorage, 0); + + % 1.3. Get flux at bottom boundary + if ~GroundwaterSettings.GroundwaterCoupling % no groundwater coupling + botmFlux = HBoundaryFlux.QMB; + else + botmFlux = gwfluxes.recharge; + end + botmFluxIn = max(botmFlux, 0); % input to soil zone (capillary rise flux in case of groundwater coupling) + botmFluxOut = max(-botmFlux, 0); % output from soil zone (recharge in case of groundwater coupling) + + % 1.4. Calculate initial total inflow, total outflow, residual, and error + totalInflowInit = ForcingData.Precip + botmFluxIn + deltaStorageIn; + totalOutflowInit = runoff + ET + botmFluxOut + deltaStorageOut; + residualInit = totalInflowInit - totalOutflowInit; + errorDenominator = mean([totalInflowInit + totalOutflowInit]); % avoid division by zero if total inflow = 0 + errorInit = residualInit / max(errorDenominator, eps) * 100; % use small value (eps) to avoid division by zero + errorInit = max(min(errorInit, 100), -100); % constrain extreme error values + + % Update for values in csv file (use only net values of delta storage and bottom flux) + totalInflowInit = ForcingData.Precip; + totalOutflowInit = runoff + ET - botmFlux - deltaStorage; + + % define outputs (used if correction below is skipped or error < error threshold) + correctedDeltaS = 0; + correctedDeltaSIn = 0; + correctedDeltaSOut = 0; + residual = residualInit; + error = errorInit; + totalInflow = totalInflowInit; + totalOutflow = totalOutflowInit; + Theta_UU_corrected = SoilVariables.Theta_UU; + + %%%%%%%%%%%%%%%%%% Task 2: Enforce closing water balance by correcting soil moisture profile %%%%%%%%%%%%%%%%%% + errorThreshold = 1; % unit is percentage + maxIterations = 30; + if correctWaterBalance && abs(errorInit) > errorThreshold + iteration = 0; + residual = 0; % starting value at each iteration + + % Loop to correct soil moisture until the water balance error is below the threshold + while abs(error) > errorThreshold && iteration < maxIterations + iteration = iteration + 1; + + % 2.1. Determine the total soil thickness + if GroundwaterSettings.GroundwaterCoupling + indxBotm = GroundwaterSettings.indxBotmLayer; + soil_depth = sum(ModelSettings.DeltZ(indxBotm + 1:ModelSettings.NL)); + else + indxBotm = 1; + soil_depth = sum(ModelSettings.DeltZ); + end + + % 2.2. Calculate correction value of soil moisture + thetaCorrection = residual / soil_depth * TimeProperties.DELT; + + % Apply max change in soil moisture (to avoid dramatic changes) + maxThetaChange = 0.02; + thetaCorrection = sign(thetaCorrection) * min(abs(thetaCorrection), maxThetaChange); + + % Add the correction value to the soil moisture profile + for i = indxBotm + 1:ModelSettings.NL + Theta_UU_corrected(i, 1) = Theta_UU_current(i, 1) + thetaCorrection; + Theta_UU_corrected(i, 2) = Theta_UU_current(i, 2) + thetaCorrection; + end + + % 2.3. Ensure corrected soil moisture values remain within residual and saturated water content limits + [Theta_UU_corrected] = io.constrainSoilVariables(Theta_UU_corrected, VanGenuchten); + + % Recalculate the correction value of soil moisture after checking residual and saturation bounds + thetaCorrectionLay = Theta_UU_corrected(1:end - 1, 1) - Theta_UU_previous(1:end - 1, 1); + + % 2.4. Recalculate corrected delta storage (based on the corrected soil moisture profile) + correctedDeltaS = sum(thetaCorrectionLay .* ModelSettings.DeltZ') / TimeProperties.DELT; + correctedDeltaSIn = max(-correctedDeltaS, 0); + correctedDeltaSOut = max(correctedDeltaS, 0); + + % 2.5. Calculate corrected total inflow, total outflow, residual and error + totalInflow = ForcingData.Precip + botmFluxIn + deltaStorageIn + correctedDeltaSIn; + totalOutflow = runoff + ET + botmFluxOut + deltaStorageOut + correctedDeltaSOut; + residual = totalInflow - totalOutflow; + errorDenominator = mean([totalInflow + totalOutflow]); % avoid division by zero if total inflow = 0 + error = residual / max(errorDenominator, eps) * 100; % use small value (eps) to avoid division by zero + error = max(min(error, 100), -100); % constrain extreme error values + end + + % Update values in csv file (use only net values of delta storage and bottom flux) + totalInflow = ForcingData.Precip; + totalOutflow = runoff + ET - botmFlux - deltaStorage + correctedDeltaS; + end + + % Outputs to be exported in csv or exposed to BMI + wbal.totalInflowInit = totalInflowInit; + wbal.totalOutflowInit = totalOutflowInit; + wbal.totalInflow = totalInflow; + wbal.totalOutflow = totalOutflow; + wbal.residual = residual; + wbal.errorInit = errorInit; + wbal.error = error; + wbal.deltaStorage = deltaStorage; + wbal.deltaStorageIn = deltaStorageIn; + wbal.deltaStorageOut = deltaStorageOut; + wbal.correctedDeltaS = correctedDeltaS; + wbal.correctedDeltaSIn = correctedDeltaSIn; + wbal.correctedDeltaSOut = correctedDeltaSOut; + wbal.botmFlux = botmFlux; + wbal.botmFluxIn = botmFluxIn; + wbal.botmFluxOut = botmFluxOut; +end