diff --git a/CHANGELOG.md b/CHANGELOG.md
index 4430830f..38a972a5 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -16,9 +16,12 @@
 - The bug in the QVT calculations discussed in
   [#230](https://github.com/EcoExtreML/STEMMUS_SCOPE/issues/230) and fixed in
   [#234](https://github.com/EcoExtreML/STEMMUS_SCOPE/pull/234)
-- The bug in  activating the dry air calculations discussed in
+- The bug in activating the dry air calculations discussed in
   [#227](https://github.com/EcoExtreML/STEMMUS_SCOPE/issues/230) and fixed in
   [#234](https://github.com/EcoExtreML/STEMMUS_SCOPE/pull/227)
+- Defining the indices of the first four layers discussed in
+  [#237](https://github.com/EcoExtreML/STEMMUS_SCOPE/issues/237) and fixed in
+  [#238](https://github.com/EcoExtreML/STEMMUS_SCOPE/pull/238)
 
 
 <a name="1.5.0"></a>
diff --git a/run_model_on_snellius/exe/STEMMUS_SCOPE b/run_model_on_snellius/exe/STEMMUS_SCOPE
index b702486b..e95dfd9c 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/+soilmoisture/calculateBoundaryConditions.m b/src/+soilmoisture/calculateBoundaryConditions.m
index 4c841eb2..6d1fff71 100644
--- a/src/+soilmoisture/calculateBoundaryConditions.m
+++ b/src/+soilmoisture/calculateBoundaryConditions.m
@@ -72,7 +72,7 @@
         % Check applied infiltration doesn't exceed 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(51:54) * SoilVariables.Theta_UU(51:54, 1); % actual moisture 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);
 
diff --git a/src/calc_rsoil.m b/src/calc_rsoil.m
index 446c0a65..a37cc034 100644
--- a/src/calc_rsoil.m
+++ b/src/calc_rsoil.m
@@ -2,7 +2,7 @@
     % load Constants
     Constants = io.define_constants();
 
-    SMC = SoilVariables.Theta_LL(1:54, 2);
+    SMC = SoilVariables.Theta_LL(1:end - 1, 2);
     Se  = (SMC - VanGenuchten.Theta_r') ./ (VanGenuchten.Theta_s' - VanGenuchten.Theta_r');
     Ksoil = SoilVariables.Ks' .* Se.^Constants.l .* (1 - (1 - Se.^(1 ./ VanGenuchten.m')).^(VanGenuchten.m')).^2;
 
diff --git a/src/calc_sfactor.m b/src/calc_sfactor.m
index 866f13bf..4e58b213 100644
--- a/src/calc_sfactor.m
+++ b/src/calc_sfactor.m
@@ -1,5 +1,5 @@
 function [sfactor] = calc_sfactor(Rl, Theta_s, Theta_r, Theta_LL, bbx, Ta, Theta_f)
-    SMC = Theta_LL(1:54, 1); % soil surface moisture
+    SMC = Theta_LL(1:end - 1, 1); % soil surface moisture
     wfrac = 1 ./ (1 + exp((-100 .* Theta_s') .* (SMC - (Theta_f' + Theta_r') / 2))) .* bbx; % The soil water stress factor in each layer
     RL = Rl .* bbx;
     RLfrac = RL ./ (sum(sum(RL))); % root fraction in each layer