EcoExtreML · MostafaGomaa93 · Sep 2, 2025 · Sep 2, 2025 · Sep 5, 2025 · Sep 5, 2025
diff --git a/run_model_on_snellius/exe/STEMMUS_SCOPE b/run_model_on_snellius/exe/STEMMUS_SCOPE
diff --git a/src/+groundwater/updateDunnianRunoff.m b/src/+groundwater/updateDunnianRunoff.m
diff --git 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
@@ -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
@@ -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
@@ -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
@@ -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
@@ -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
@@ -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;