Correct Snow precipitation

.gitignore

@@ -1 +1,7 @@
 .data/
+run_model_on_snellius/exe/includedSupportPackages.txt
+run_model_on_snellius/exe/mccExcludedFiles.log
+run_model_on_snellius/exe/readme.txt
+run_model_on_snellius/exe/requiredMCRProducts.txt
+run_model_on_snellius/exe/run_STEMMUS_SCOPE.sh
+run_model_on_snellius/exe/unresolvedSymbols.txt

CHANGELOG.md

@@ -1,12 +1,20 @@
 # Unreleased
 
-**Added:**
+**Changed:**
+- Change array shape of precipitation and runoff fluxes from time-series array to 1D array (for BMI purposes) discussed in [#280](https://github.com/EcoExtreML/STEMMUS_SCOPE/pull/280) and added in [#282](https://github.com/EcoExtreML/STEMMUS_SCOPE/pull/282) 
+
+**Fixed:**
+- Correct snow precipitation calcuations discussed in [#279](https://github.com/EcoExtreML/STEMMUS_SCOPE/issues/279) and fixed in [#282](https://github.com/EcoExtreML/STEMMUS_SCOPE/pull/282)
+
+<a name="1.6.2"></a>
+# [1.6.2](https://github.com/EcoExtreML/STEMMUS_SCOPE/releases/tag/1.6.2) - 17 Dec 2024
 
+The 1.6.2 release comes with minor changes as:
+
+**Added:**
 - Documentation using mkdocs and a github action workflow to publish the
   documentation in [#264](https://github.com/EcoExtreML/STEMMUS_SCOPE/pull/264)
 
-**Changed:**
-
 <a name="1.6.1"></a>
 # [1.6.1](https://github.com/EcoExtreML/STEMMUS_SCOPE/releases/tag/1.6.1) - 26 Sep 2024
 

src/+energy/calculateBoundaryConditions.m

@@ -1,5 +1,5 @@
 function [RHS, EnergyMatrices] = calculateBoundaryConditions(BoundaryCondition, EnergyMatrices, HBoundaryFlux, ForcingData, SoilVariables, ...
-                                                             Precip, Evap, Delt_t, r_a_SOIL, Rn_SOIL, RHS, L, KT, ModelSettings, GroundwaterSettings)
+                                                             Evap, Delt_t, r_a_SOIL, Rn_SOIL, RHS, L, KT, ModelSettings, GroundwaterSettings)
     %{
         Determine the boundary condition for solving the energy equation, see
         STEMMUS Technical Notes.
@@ -46,6 +46,6 @@
     else
         L_ts = L(n);
         SH = 0.1200 * Constants.c_a * (SoilVariables.T(n) - ForcingData.Ta_msr(KT)) / r_a_SOIL(KT);
-        RHS(n) = RHS(n) + 100 * Rn_SOIL(KT) / 1800 - Constants.RHOL * L_ts * Evap - SH + Constants.RHOL * Constants.c_L * (ForcingData.Ta_msr(KT) * Precip(KT) + BoundaryCondition.DSTOR0 * SoilVariables.T(n) / Delt_t);    % J cm-2 s-1
+        RHS(n) = RHS(n) + 100 * Rn_SOIL(KT) / 1800 - Constants.RHOL * L_ts * Evap - SH + Constants.RHOL * Constants.c_L * (ForcingData.Ta_msr(KT) * ForcingData.infiltration + BoundaryCondition.DSTOR0 * SoilVariables.T(n) / Delt_t);    % J cm-2 s-1
     end
 end

src/+energy/solveEnergyBalanceEquations.m

@@ -2,7 +2,7 @@
                                                                                         AirVariabes, VaporVariables, GasDispersivity, ThermalConductivityCapacity, ...
                                                                                         HBoundaryFlux, BoundaryCondition, ForcingData, DRHOVh, DRHOVT, KL_T, ...
                                                                                         Xah, XaT, Xaa, Srt, L_f, RHOV, RHODA, DRHODAz, L, Delt_t, P_g, P_gg, ...
-                                                                                        TOLD, Precip, EVAP, r_a_SOIL, Rn_SOIL, KT, CHK, ModelSettings, GroundwaterSettings)
+                                                                                        TOLD, EVAP, r_a_SOIL, Rn_SOIL, KT, CHK, ModelSettings, GroundwaterSettings)
     %{
         Solve the Energy balance equation with the Thomas algorithm to update
         the soil temperature 'SoilVariables.TT', the finite difference
@@ -19,7 +19,7 @@
     [RHS, EnergyMatrices, SAVE] = energy.assembleCoefficientMatrices(InitialValues, EnergyMatrices, SoilVariables, Delt_t, P_g, P_gg, ModelSettings, GroundwaterSettings);
 
     [RHS, EnergyMatrices] = energy.calculateBoundaryConditions(BoundaryCondition, EnergyMatrices, HBoundaryFlux, ForcingData, SoilVariables, ...
-                                                               Precip, EVAP, Delt_t, r_a_SOIL, Rn_SOIL, RHS, L, KT, ModelSettings, GroundwaterSettings);
+                                                               EVAP, Delt_t, r_a_SOIL, Rn_SOIL, RHS, L, KT, ModelSettings, GroundwaterSettings);
 
     [SoilVariables, CHK, RHS, EnergyMatrices] = energy.solveTridiagonalMatrixEquations(EnergyMatrices, SoilVariables, RHS, CHK, ModelSettings, GroundwaterSettings);
 

src/+groundwater/updateDunnianRunoff.m

@@ -1,13 +1,13 @@
-function [dunnianRunoff, update_Precip_msr] = updateDunnianRunoff(Precip_msr, groundWaterDepth)
+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
-    dunnianRunoff = zeros(size(Precip_msr));
+    runoffDunnian = zeros(size(Precip_msr));
     update_Precip_msr = Precip_msr;
 
     if groundWaterDepth <= 1.0
-        dunnianRunoff = Precip_msr;
+        runoffDunnian = Precip_msr;
         update_Precip_msr = zeros(size(Precip_msr));
     end
 

src/+init/defineInitialValues.m

@@ -171,7 +171,7 @@
 
     Nmsrmn = Dur_tot * 10; % Here, it is made as big as possible, in case a long simulation period containing many time step is defined.
     fields = {
-              'Precip', 'Ta', 'Ts', 'U', 'HR_a', 'Rns', 'Rnl', 'Rn', ...
+              'Ta', 'Ts', 'U', 'HR_a', 'Rns', 'Rnl', 'Rn', ...
               'h_SUR', 'SH', 'MO', 'Zeta_MO', 'TopPg'
              };
     structures{4} = helpers.createStructure(zeros(Nmsrmn, 1), fields);

src/+io/loadForcingData.m

@@ -20,7 +20,7 @@
     Precip_msr = Mdata{:, 6}; % (cm/sec)
 
     %%%%%%%%%% Adjust precipitation to get applied infiltration after removing: (1) saturation excess runoff   %%%%%%%%%%
-    %%%%%%%%%%                                                                  (2) infiltration excess runoff %%%%%%%%%%
+    %%%%%%%%%%           (2) infiltration excess runoff (updated in +soilmoisture/calculateBoundaryConditions) %%%%%%%%%%
     % Note: Code changes are related to the issue: https://github.com/EcoExtreML/STEMMUS_SCOPE/issues/232
     % Note: Adjusting the precipitation after the canopy interception is not implemented yet.
 
@@ -32,13 +32,9 @@
         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.R_Dunn = Precip_msr .* fsat; % Dunnian runoff (saturation excess runoff, in c/sec)
-        Precip_msr = Precip_msr .* (1 - fsat); % applied infiltration after removing Dunnian runoff
+        ForcingData.runoffDunnian = Precip_msr .* fsat; % Dunnian runoff (saturation excess runoff, in cm/sec)
     end
 
-    % (2) Infiltration excess runoff (Hortonian runoff)
-    ForcingData.R_Hort = zeros(size(Precip_msr)); % will be updated in +soilmoisture/calculateBoundaryConditions
-
     % replace negative values
     for jj = 1:Dur_tot
         if ForcingData.Ta_msr(jj) < -100
@@ -49,7 +45,4 @@
     % Outputs to be used by other functions
     ForcingData.Tmin = min(ForcingData.Ta_msr);
     ForcingData.Precip_msr = Precip_msr;
-    % Applied infiltration (= precipitation after removal of Dunnian runoff)
-    ForcingData.applied_inf = Precip_msr; % later will be updated in the ....
-    % +soilmoisture/calculateBoundaryConditions after removal of Hortonian runoff
 end

src/+soilmoisture/calculateBoundaryConditions.m

@@ -1,20 +1,14 @@
-function [AVAIL0, RHS, HeatMatrices, Precip, ForcingData] = calculateBoundaryConditions(BoundaryCondition, HeatMatrices, ForcingData, SoilVariables, InitialValues, ...
-                                                                                        TimeProperties, SoilProperties, RHS, hN, KT, Delt_t, Evap, ModelSettings, GroundwaterSettings)
+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.
     %}
 
-    % n is the index of n_th item
     n = ModelSettings.NN;
-
     C4 = HeatMatrices.C4;
     C4_a = HeatMatrices.C4_a;
 
-    Precip = InitialValues.Precip;
-    Precip_msr = ForcingData.Precip_msr;
-    Precipp = 0;
-
-    %  Apply the bottom boundary condition called for by BoundaryCondition.NBChB
+    %%%%%%  Apply the bottom boundary condition called for by BoundaryCondition.NBChB  %%%%%%
     if ~GroundwaterSettings.GroundwaterCoupling  % no Groundwater coupling
         if BoundaryCondition.NBChB == 1            %  Specify matric head at bottom to be ---BoundaryCondition.BChB;
             RHS(1) = BoundaryCondition.BChB;
@@ -46,7 +40,39 @@
         end
     end
 
-    %  Apply the surface boundary condition called for by BoundaryCondition.NBCh
+    %%%%%%  Apply the surface boundary condition called for by BoundaryCondition.NBCh  %%%%%%
+    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
+        Precip_snowAccum = 0; % initalize accumulated snow for first time step
+    else
+        Precip_snowAccum = ForcingData.Precip_snowAccum;
+    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;
+        if KIT == 1 % accumulate snow at one iteration only within the 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
+            Precip_liquid = ForcingData.Precip_liquid;
+        end
+        Precip_snow = 0;
+    end
+
+    % Infiltration = Precipitation - total runoff (Dunnian runoff + Hortonian runoff)
+    infiltration = Precip_liquid - runoffDunn; % removing Dunnian runoff (Hortonian runoff is removed below)
+
     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
@@ -65,41 +91,27 @@
             RHS(n) = RHS(n) - BoundaryCondition.BCh;   % a specified matric head (saturation or dryness) was applied;
         end
     else % (BoundaryCondition.NBCh == 3, Specified atmospheric forcing)
-
-        % Calculate applied infiltration and infiltration excess runoff (Hortonian runoff), modified by Mostafa
+        % Calculate infiltration excess runoff (Hortonian runoff) and update applied infiltration, 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)
-        % Check applied infiltration doesn't exceed infiltration capacity
+        % Calculate infiltration capacity
         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; % (cm/sec)
         infCap_min = min(Ks0, infCap);
 
-        % Infiltration excess runoff (Hortonian runoff). Note: Dunnian runoff is calculated in the +io/loadForcingData file
-        if Precip_msr(KT) > infCap_min
-            ForcingData.R_Hort(KT) = Precip_msr(KT) - infCap_min;
-        else
-            ForcingData.R_Hort(KT) = 0;
-        end
-
-        Precip(KT) = min(Precip_msr(KT), infCap_min);
-        ForcingData.applied_inf(KT) = Precip(KT); % applied infiltration after removing Hortonian runoff
-
-        if SoilVariables.Tss(KT) > 0
-            Precip(KT) = Precip(KT);
+        % Infiltration excess runoff (Hortonian runoff)
+        if infiltration > infCap_min
+            runoffHort = infiltration - infCap_min;
         else
-            Precip(KT) = Precip(KT);
-            Precipp = Precipp + Precip(KT);
-            Precip(KT) = 0;
+            runoffHort = 0;
         end
+        % Update applied infiltration after removing Hortonian runoff
+        infiltration = min(infiltration, infCap_min);
 
-        if SoilVariables.Tss(KT) > 0
-            AVAIL0 = Precip(KT) + Precipp + BoundaryCondition.DSTOR0 / Delt_t; % (cm/sec)
-            Precipp = 0;
-        else
-            AVAIL0 = Precip(KT) + BoundaryCondition.DSTOR0 / Delt_t;
-        end
+        % Add depression water to applied infiltration
+        AVAIL0 = infiltration + BoundaryCondition.DSTOR0 / Delt_t; % (cm/sec)
 
         if BoundaryCondition.NBChh == 1
             RHS(n) = hN;
@@ -111,6 +123,15 @@
             RHS(n) = RHS(n) + AVAIL0 - Evap;
         end
     end
+
+    % Outputs to be exported or used in other functions
     HeatMatrices.C4 = C4;
     HeatMatrices.C4_a = C4_a;
+    ForcingData.Precip = Precip;
+    ForcingData.Precip_liquid = Precip_liquid;
+    ForcingData.Precip_snow = Precip_snow;
+    ForcingData.Precip_snowAccum = Precip_snowAccum;
+    ForcingData.infiltration = infiltration;
+    ForcingData.runoffDunn = runoffDunn;
+    ForcingData.runoffHort = runoffHort;
 end

src/+soilmoisture/solveSoilMoistureBalance.m

@@ -1,5 +1,5 @@
-function [SoilVariables, HeatMatrices, HeatVariables, HBoundaryFlux, Rn_SOIL, Evap, Trap, r_a_SOIL, Srt, CHK, AVAIL0, Precip, RWUs, RWUg, ForcingData] = solveSoilMoistureBalance(SoilVariables, InitialValues, ForcingData, VaporVariables, GasDispersivity, TimeProperties, SoilProperties, ...
-                                                                                                                                                                                  BoundaryCondition, Delt_t, RHOV, DRHOVh, DRHOVT, D_Ta, hN, RWU, fluxes, KT, hOLD, Srt, P_gg, ModelSettings, GroundwaterSettings)
+function [SoilVariables, HeatMatrices, HeatVariables, HBoundaryFlux, Rn_SOIL, Evap, Trap, r_a_SOIL, Srt, CHK, AVAIL0, RWUs, RWUg, ForcingData] = solveSoilMoistureBalance(SoilVariables, InitialValues, ForcingData, VaporVariables, GasDispersivity, TimeProperties, SoilProperties, ...
+                                                                                                                                                                          BoundaryCondition, Delt_t, RHOV, DRHOVh, DRHOVT, D_Ta, hN, RWU, fluxes, KT, KIT, hOLD, Srt, P_gg, ModelSettings, GroundwaterSettings)
     %{
         Solve the soil moisture balance equation with the Thomas algorithm to
         update the soil matric potential 'hh', the finite difference
@@ -24,8 +24,8 @@
         r_a_SOIL = InitialValues.r_a_SOIL;
     end
 
-    [AVAIL0, RHS, HeatMatrices, Precip, ForcingData] = soilmoisture.calculateBoundaryConditions(BoundaryCondition, HeatMatrices, ForcingData, SoilVariables, InitialValues, ...
-                                                                                                TimeProperties, SoilProperties, RHS, hN, KT, Delt_t, Evap, ModelSettings, GroundwaterSettings);
+    [AVAIL0, RHS, HeatMatrices, ForcingData] = soilmoisture.calculateBoundaryConditions(BoundaryCondition, HeatMatrices, ForcingData, SoilVariables, InitialValues, ...
+                                                                                        TimeProperties, SoilProperties, RHS, hN, KT, KIT, Delt_t, Evap, ModelSettings, GroundwaterSettings);
 
     [CHK, hh, C4] = soilmoisture.solveTridiagonalMatrixEquations(HeatMatrices.C4, SoilVariables.hh, HeatMatrices.C4_a, RHS, ModelSettings, GroundwaterSettings);
 

src/STEMMUS_SCOPE.m

@@ -354,7 +354,7 @@
         GroundwaterSettings.gw_Dep = groundwater.calculateGroundWaterDepth(GroundwaterSettings.topLevel, GroundwaterSettings.headBotmLayer, ModelSettings.Tot_Depth);
 
         % Update Dunnian runoff and ForcingData.Precip_msr
-        [ForcingData.R_Dunn, ForcingData.Precip_msr] = groundwater.updateDunnianRunoff(ForcingData.Precip_msr, GroundwaterSettings.gw_Dep);
+        [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);
@@ -653,10 +653,10 @@
             GasDispersivity = conductivity.calculateGasDispersivity(InitialValues, SoilVariables, P_gg, k_g, ModelSettings);
 
             % Srt is both input and output
-            [SoilVariables, HeatMatrices, HeatVariables, HBoundaryFlux, Rn_SOIL, Evap, Trap, r_a_SOIL, Srt, CHK, AVAIL0, Precip, RWUs, RWUg, ForcingData] = ...
+            [SoilVariables, HeatMatrices, HeatVariables, HBoundaryFlux, Rn_SOIL, Evap, Trap, r_a_SOIL, Srt, CHK, AVAIL0, RWUs, RWUg, ForcingData] = ...
                                   soilmoisture.solveSoilMoistureBalance(SoilVariables, InitialValues, ForcingData, VaporVariables, GasDispersivity, ...
                                                                         TimeProperties, SoilProperties, BoundaryCondition, Delt_t, RHOV, DRHOVh, ...
-                                                                        DRHOVT, D_Ta, hN, RWU, fluxes, KT, hOLD, Srt, P_gg, ModelSettings, GroundwaterSettings);
+                                                                        DRHOVT, D_Ta, hN, RWU, fluxes, KT, KIT, hOLD, Srt, P_gg, ModelSettings, GroundwaterSettings);
 
             if BoundaryCondition.NBCh == 1
                 DSTOR = 0;
@@ -678,7 +678,7 @@
                 DSTOR = min(EXCESS, DSTMAX); % Depth of depression storage at end of current time step
                 % Next line is commented and Surface runoff is re-calculated using different approach (the following 3 lines)
                 % RS(KT) = (EXCESS - DSTOR) / Delt_t; % surface runoff, (unit conversion from cm/30mins to cm/sec)
-                RS(KT) = ForcingData.R_Hort(KT) + ForcingData.R_Dunn(KT); % total surface runoff (cm/sec)
+                ForcingData.runoff = ForcingData.runoffHort + ForcingData.runoffDunn; % total surface runoff (cm/sec)
             end
 
             if ModelSettings.Soilairefc == 1
@@ -701,7 +701,7 @@
                                                                                                       AirVariabes, VaporVariables, GasDispersivity, ThermalConductivityCapacity, ...
                                                                                                       HBoundaryFlux, BoundaryCondition, ForcingData, DRHOVh, DRHOVT, KL_T, ...
                                                                                                       Xah, XaT, Xaa, Srt, L_f, RHOV, RHODA, DRHODAz, L, Delt_t, P_g, P_gg, ...
-                                                                                                      TOLD, Precip, Evap, r_a_SOIL, Rn_SOIL, KT, CHK, ModelSettings, GroundwaterSettings);
+                                                                                                      TOLD, Evap, r_a_SOIL, Rn_SOIL, KT, CHK, ModelSettings, GroundwaterSettings);
             end
 
             if max(CHK) < 0.1
@@ -829,5 +829,5 @@
 % Calculate the total simulatiom time, added by mostafa
 end_time = clock;
 simtime_min = etime(end_time, start_time) / 60;
-simtime_hr = simtime_min / 24;
+simtime_hr = simtime_min / 60;
 disp(['Simulation time is : ' num2str(simtime_hr) ' hrs (' num2str(simtime_min) ' minutes)']);