Rm unused variables and equations

src/+energy/assembleCoefficientMatrices.m

@@ -59,11 +59,12 @@
             - (C1(n - 1, 2) / Delt_t + C4_a(n - 1)) * SV.hh(n - 1) - (C1(n, 1) / Delt_t + C4(n, 1)) * SV.hh(n) ...
             + (C3(n - 1, 2) / Delt_t) * P_g(n - 1) + (C3(n, 1) / Delt_t) * P_g(n) ...
             + (C1(n - 1, 2) / Delt_t) * SV.h(n - 1) + (C1(n, 1) / Delt_t) * SV.h(n);
+
     elseif ~ModelSettings.Soilairefc && ModelSettings.Thmrlefc
         RHS(indxBotm) = -C7(indxBotm) + (C2(indxBotm, 1) * SV.T(indxBotm) + C2(indxBotm, 2) * SV.T(indxBotm + 1)) / Delt_t ...
             - (C1(indxBotm, 1) / Delt_t + C4(indxBotm, 1)) * SV.hh(indxBotm) - (C1(indxBotm, 2) / Delt_t + C4(indxBotm, 2)) * SV.hh(indxBotm + 1) ...
             + (C1(indxBotm, 1) / Delt_t) * SV.h(indxBotm) + (C1(indxBotm, 2) / Delt_t) * SV.h(indxBotm + 1);
-        for i = 2:ModelSettings.NL
+        for i = indxBotm:ModelSettings.NL
             ARG4 = C1(i - 1, 2) / Delt_t;
             ARG5 = C1(i, 1) / Delt_t;
             ARG6 = C1(i, 2) / Delt_t;
@@ -76,6 +77,7 @@
         RHS(n) = -C7(n) + (C2(n - 1, 2) * SV.T(n - 1) + C2(n, 1) * SV.T(n)) / Delt_t ...
             - (C1(n - 1, 2) / Delt_t + C4(n - 1, 2)) * SV.hh(n - 1) - (C1(n, 1) / Delt_t + C4(n, 1)) * SV.hh(n) ...
             + (C1(n - 1, 2) / Delt_t) * SV.h(n - 1) + (C1(n, 1) / Delt_t) * SV.h(n);
+
     else
         RHS(indxBotm) = -C7(indxBotm) + (C2(indxBotm, 1) * SV.T(indxBotm) + C2(1, 2) * SV.T(indxBotm + 1)) / Delt_t;
         for i = indxBotm + 1:ModelSettings.NL

src/+energy/calculateBoundaryConditions.m

@@ -17,6 +17,10 @@
         % index of bottom layer after neglecting saturated layers (from bottom to top)
         indxBotm = GroundwaterSettings.indxBotmLayer;
         tempBotm = GroundwaterSettings.tempBotm; % groundwater temperature
+
+        if isnan(GroundwaterSettings.tempBotm) % no available groundwater temperature data
+            tempBotm = SoilVariables.TT(indxBotm);
+        end
     end
 
     % Apply the bottom boundary condition called for by BoundaryCondition.NBCTB

src/+energy/solveEnergyBalanceEquations.m

@@ -37,13 +37,6 @@
         end
     end
 
-    if GroundwaterSettings.GroundwaterCoupling
-        % Assign the groundwater temperature to the saturated layers
-        for i = indxBotm - 1:-1:1
-            SoilVariables.TT(i) = GroundwaterSettings.tempBotm; % groundwater temperature;
-        end
-    end
-
     for i = 1:ModelSettings.NN
         if SoilVariables.TT(i) <= -272
             SoilVariables.TT(i) = -272;

src/+energy/solveTridiagonalMatrixEquations.m

@@ -26,12 +26,13 @@
 
     if GroundwaterSettings.GroundwaterCoupling
         for i = indxBotm:-1:2
-            RHS(i - 1) = RHS(i); % assign temp bottom boundary for all saturated layers
+            RHS(i - 1) = RHS(i); % assign the groundwater temperature to the saturated layers
         end
+        % correct soil temp above groundwater table (avoid wrong temp behaviour in case RHS(indxBotm) > RHS(indxBotm + 1))
+        RHS(indxBotm + 1) = (RHS(indxBotm) + RHS(indxBotm + 2)) / 2;
     end
 
     for i = 1:ModelSettings.NN
-        % note: this for loop causes wrong values if i = indxBotm:ModelSettings.NN instead of i = 1:ModelSettings.NN
         CHK(i) = abs(RHS(i) - SoilVariables.TT(i));
         SoilVariables.TT(i) = RHS(i);
     end

src/+groundwater/calculateGroundWaterDepth.m

@@ -0,0 +1,13 @@
+function groundWaterDepth = calculateGroundWaterDepth(topLevel, headBotmLayer, Tot_Depth)
+    % water table depth: depth from top soil layer to groundwater level
+    groundWaterDepth = topLevel - headBotmLayer; % depth from top layer to groundwater level
+
+    % Check that the position of the water table is within the soil column
+    if groundWaterDepth <= 0
+        warning('The soil is fully saturated up to the land surface level!');
+        groundWaterDepth = 1.0; % to avoid model crashing, assign minimum groundWaterDepth value of 1 cm
+    elseif groundWaterDepth > Tot_Depth
+        warning('Groundwater table is below the end of the soil column!');
+    end
+
+end

src/+groundwater/calculateGroundwaterRecharge.m

@@ -31,13 +31,13 @@ STheta_L                soil moisture at the start of the current time step (top
         STheta_LL               soil moisture at the end of the current time step (top to bottom)
         soilThick               cumulative soil layers thickness (from top to bottom)
         indxAqLay               index of the MODFLOW aquifer that corresponds to each STEMMUS soil layer
-        aqLayers                elevation of top surface level and all bottom levels of aquifer layers, received from MODFLOW through BMI
+        aqlevels                elevation of top surface level and all bottom levels of aquifer layers, received from MODFLOW through BMI
     %}
 
     % Start Recharge calculations
     % (a) Define the upper and lower boundaries of the moving balancing domain
     % The moving balancing domain is located between depToGWT_strt and depToGWT_end
-    [depToGWT_end, indxGWLay_end] = groundwater.findPhreaticSurface(SoilVariables, KT, GroundwaterSettings);
+    [depToGWT_end, indxGWLay_end] = groundwater.findPhreaticSurface(SoilVariables.hh, KT, GroundwaterSettings.soilThick, GroundwaterSettings.indxBotmLayer_R);
     % indxRchrg and indxRchrgMax are the indices of the upper and lower levels of the moving boundary
     % Following the HYDRUS-MODFLOW paper and also STEMMUS-MODFLOW, indxRchrg and indxRchrgMax are defined as in the next two lines
     indxRchrgMax = max(indxGWLay_strt, indxGWLay_end) + 2; % the positive 2 is a user-specified value to define lower boundary of the moving boundary
@@ -64,18 +64,21 @@ soilThick               cumulative soil layers thickness (from top to bottom)
         recharge_init = Q_flip(indxRchrg); % mean([Q_flip(indxRchrg), Q_flip(indxRchrg_above)])
     end
 
+    %{
     % (d) Calculations of SY
     % Note: In the HYDRUS-MODFLOW paper, Sy (from MODFLOW) was used. In Lianyu STEMMUS_MODFLOW code, a combination of Sy and Ss was used
-    indxAqLay = GroundwaterSettings.indxAqLay; % index of MODFLOW aquifer layers for each STEMMUS soil layer
-    aqLayers = GroundwaterSettings.aqLayers; % elevation of top surface level and all bottom levels of aquifer layers
+    aqlevels = GroundwaterSettings.aqlevels; % elevation of top surface level and all bottom levels of aquifer layers
+    numAqL = GroundwaterSettings.numAqL; % number of MODFLOW aquifer layers
+    soilThick = GroundwaterSettings.soilThick; % cumulative soil layer thickness (from top to bottom)
+    indxAqLay = groundwater.calculateIndexAquifer(aqlevels, numAqL, soilThick); % index of MODFLOW aquifer layers for each STEMMUS soil layer
+
     K = indxAqLay(indxGWLay_end);
-    Thk = aqLayers(1) - aqLayers(K) - depToGWT_end;
+    Thk = aqlevels(1) - aqlevels(K) - depToGWT_end;
     SY = GroundwaterSettings.SY;
     SS = GroundwaterSettings.SS;
     S = (SY(K) - SS(K) * Thk) * (depToGWT_strt - depToGWT_end);
 
     % (e) Calculations of sy
-    soilThick = GroundwaterSettings.soilThick; % cumulative soil layer thickness (from top to bottom)
     ModelSettings = io.getModelSettings();
     NN = ModelSettings.NN; % Number of nodes
     NL = ModelSettings.NL; % Number of layers
@@ -93,8 +96,10 @@ soilThick               cumulative soil layers thickness (from top to bottom)
     end
 
     % (f) Aggregate c, d, and e to get recharge
-    % after couple of tests, it appears that the effect of S and sy is very minor, so they are removed but kept in the code for further investigation
     % recharge = recharge_init + S - sy;
+    %}
+
+    % after couple of tests, it appears that the effect of S and sy is very minor, so they are commented but kept in the code for further investigation
     gwfluxes.recharge = recharge_init; % Note: in STEMMUS +ve means up-flow direction and -ve means down (opposite of MODFLOW), so recharge sign needs to be converted in BMI
 
     if isnan(gwfluxes.recharge) || isinf(gwfluxes.recharge)

src/+groundwater/calculateIndexAquifer.m

@@ -0,0 +1,28 @@
+function indxAqLay = calculateIndexAquifer(aqlevels, numAqL, soilThick)
+    %{
+        Assign the index of the MODFLOW aquifer that corresponds to each STEMMUS soil layer.
+        aqlevels            elevation of top surface level and all bottom levels of aquifer layers
+        numAqL              number of MODFLOW aquifer layers, received from MODFLOW through BMI
+    %}
+
+    % Load model settings
+    ModelSettings = io.getModelSettings();
+
+    numAqN = numAqL + 1; % number of MODFLOW aquifer nodes
+    indxAqLay = zeros(ModelSettings.NN, 1);
+    indxAqLay(1) = 1;
+    for i = 2:ModelSettings.NN
+        for K = 2:numAqN
+            Z1 = aqlevels(K - 1);
+            Z0 = aqlevels(K);
+            ZZ = aqlevels(1) - soilThick(i);
+            if ZZ <= Z1 && ZZ > Z0
+                indxAqLay(i) = K - 1;
+                break
+            elseif ZZ == Z0 && K == numAqN
+                indxAqLay(i) = K - 1;
+                break
+            end
+        end
+    end
+end

src/+groundwater/calculateIndexBottomLayer.m

@@ -0,0 +1,31 @@
+function [indxBotmLayer, indxBotmLayer_R] = calculateIndexBottomLayer(soilThick, gw_Dep)
+    %{
+        Calculate the index of the bottom layer level using MODFLOW data
+
+        indxBotmLayer_R     index of the bottom layer that contains the current headBotmLayer (top to bottom)
+        indxBotmLayer       index of the bottom layer that contains the current headBotmLayer (bottom to top)
+
+    %}
+
+    indxBotmLayer_R = [];
+
+    % Load model settings
+    ModelSettings = io.getModelSettings();
+
+    for i = 1:ModelSettings.NL
+        midThick = (soilThick(i) + soilThick(i + 1)) / 2;
+        if gw_Dep >= soilThick(i) && gw_Dep < soilThick(i + 1)
+            if gw_Dep < midThick
+                indxBotmLayer_R = i;
+            elseif gw_Dep >= midThick
+                indxBotmLayer_R = i + 1;
+            end
+            break
+        elseif gw_Dep >= soilThick(i + 1)
+            continue
+        end
+    end
+
+    % Note: indxBotmLayer_R starts from top to bottom, opposite of STEMMUS (bottom to top)
+    indxBotmLayer = ModelSettings.NN - indxBotmLayer_R + 1; % index of bottom layer (from bottom to top)
+end

src/+groundwater/calculateSoilLayerThickness.m

@@ -0,0 +1,21 @@
+function soilThick = calculateSoilLayerThickness()
+    %{
+        Calculate soil layers thickness (cumulative layers thickness; e.g. 1, 2,
+        3, 5, 10, ......., 480, total soil depth)
+
+        soilThick           cumulative soil layers thickness (from top to bottom)
+
+    %}
+
+    % Load model settings
+    ModelSettings = io.getModelSettings();
+
+    soilThick = zeros(ModelSettings.NN, 1); % cumulative soil layers thickness
+    soilThick(1) = 0;
+
+    for i = 2:ModelSettings.NL
+        soilThick(i) = soilThick(i - 1) + ModelSettings.DeltZ_R(i - 1);
+    end
+    soilThick(ModelSettings.NN) = ModelSettings.Tot_Depth; % total soil depth
+
+end

src/+groundwater/findPhreaticSurface.m

@@ -1,4 +1,4 @@
-function [depToGWT, indxGWLay] = findPhreaticSurface(SoilVariables, KT, GroundwaterSettings)
+function [depToGWT, indxGWLay] = findPhreaticSurface(soilWaterPotential, KT, soilThick, indxBotmLayer_R)
     %{
         added by Mostafa, modified after Lianyu
         This function finds the soil layer that includes the phreatic surface (saturated zone) and is used in the groundwater recharge calculations
@@ -16,7 +16,6 @@ and depToGWT (e.g. difference = 1-10 cm  based on the thickness of the soil laye
                                 %%%%%%%%%% Variables definitions %%%%%%%%%%
         depToGWT            water table depth: depth from top soil layer to groundwater level, calculated from STEMMUS variables
         indxBotmLayer_R     index of the bottom layer that contains the current headBotmLayer (top to bottom)
-        indxBotmLayer       index of the bottom layer that contains the current headBotmLayer (bottom to top)
         indxGWLay           index of the soil layer that includes the phreatic surface, so the recharge is ....
                             extracted from that layer (index order is from top layer to bottom)
                             Note: indxGWLay must equal to indxBotmLayer_R
@@ -27,10 +26,9 @@ soilThick           cumulative soil layers thickness (from top to bottom)
     % Load model settings
     ModelSettings = io.getModelSettings();
     NN = ModelSettings.NN; % number of nodes
-    soilThick = GroundwaterSettings.soilThick;
 
     % Call the matric potential
-    Shh(1:1:NN) = SoilVariables.hh(NN:-1:1);
+    Shh(1:1:NN) = soilWaterPotential(NN:-1:1);
     depToGWT = 0; % starting value for initialization, updated below
     indxGWLay = NN - 2; % starting value for initialization, updated below
 
@@ -63,7 +61,7 @@ soilThick           cumulative soil layers thickness (from top to bottom)
             warning('The phreatic surface level is lower than the end of the soil column!');
         end
         % check 3
-        diff = abs(indxGWLay - GroundwaterSettings.indxBotmLayer_R);
+        diff = abs(indxGWLay - indxBotmLayer_R);
         if diff > 1
             warning('Index of the bottom layer boundary that includes groundwater head ~= index of the layer that has zero matric potential!');
         end

src/+groundwater/initializeGroundwaterSettings.m

@@ -0,0 +1,41 @@
+function GroundwaterSettings = initializeGroundwaterSettings()
+    %{
+        added by Mostafa,
+        The concept followed to couple STEMMUS to MODFLOW can be found at https://doi.org/10.5194/hess-23-637-2019
+        and also in (preprint) https://doi.org/10.5194/gmd-2022-221
+
+                                %%%%%%%%%% Important note %%%%%%%%%%
+        The index order for the soil layers in STEMMUS is from bottom to top (NN:1), where NN is the number of STEMMUS soil layers,
+        which is the opposite of MODFLOW (top to bottom). So, when converting information between STEMMUS and MODFLOW, indices need to be flipped
+
+                                %%%%%%%%%% Variables definitions %%%%%%%%%%
+
+        headBotmLayer       groundwater head (cm) at the bottom layer, received from MODFLOW through BMI
+        tempBotm            groundwater temperature (C) at the bottom layer, received from MODFLOW through BMI
+        topLevel            elevation of the top surface aquifer layer, received from MODFLOW through BMI
+        % numAqL            number of MODFLOW aquifer layers, received from MODFLOW through BMI
+        % aqBotmlevels      elevation of all bottom levels of aquifer layers, received from MODFLOW through BMI
+        % aqlevels          elevation of top surface level and all bottom levels of aquifer layers
+        % SS                specific storage of MODFLOW aquifers, default value = 0.05 (unitless)
+        % SY                specific yield of MODFLOW aquifers, default value = 1e-5 (1/m)
+    %}
+
+    % Activate/deactivate Groundwater coupling
+    GroundwaterSettings.GroundwaterCoupling = 0; % (value = 0 -> deactivate coupling, or = 1 -> activate coupling); default = 0, update value to = 1 -> through BMI
+
+    % Initialize the variables (head, temperature) at the bottom boundary (start of saturated zone)
+    GroundwaterSettings.headBotmLayer = 1950.0; % groundwater head (cm) at bottom layer
+    GroundwaterSettings.tempBotm = NaN; % groundwater temperature at bottom layer (C)
+
+    % Call MODFLOW layers information (number of aquifer layers and their elevations, etc)
+    % elevation of the top surface aquifer layer
+    GroundwaterSettings.topLevel = 2000.0;
+
+    % GroundwaterSettings.numAqL = 5; % number of MODFLOW aquifer layers
+    % GroundwaterSettings.aqBotmlevels = [1900.0  1800.0  1700.0  1600.0  1500.0]; % elevation of all bottom levels of aquifer layers
+    % GroundwaterSettings.aqlevels = [GroundwaterSettings.topLevel, GroundwaterSettings.aqBotmlevels]; % elevation of top surface level and all bottom levels of aquifer layers
+
+    % Define Specific yield (SY) and Specific storage (SS) with default values (otherwise received from MODFLOW through BMI)
+    % GroundwaterSettings.SY = [0.05  0.05  0.05  0.05  0.05]; % default SY = 0.05 (unitless)
+    % GroundwaterSettings.SS = [1e-7  1e-7  1e-7  1e-7  1e-7]; % default SS = 1e-5 1/m = 1e-7 1/cm
+end