diff --git a/docs/src/model_docs/params_lateral.md b/docs/src/model_docs/params_lateral.md
index 2b818c18d..05d89cea3 100644
--- a/docs/src/model_docs/params_lateral.md
+++ b/docs/src/model_docs/params_lateral.md
@@ -158,7 +158,7 @@ in bold represent model parameters that can be set through static input data (ne
 soil related parameters `f`, `soilthickness`, `z_exp`, `theta_s` and `theta_r` are derived from the
 vertical `SBM` concept (including unit conversion for `f`, `z_exp` and `soilthickness`), and
 can be listed in the TOML configuration file under `[input.vertical]`, to map the internal
-model parameter to the external netCDF variable. The internal slope model parameter `beta_l` is
+model parameter to the external netCDF variable. The internal slope model parameter `slope` is
 set through the TOML file as follows:
 
 ```toml
@@ -199,7 +199,7 @@ the `layered_exponential` profile `kv` is used and `z_exp` is required as part o
 | **`theta_s`** | saturated water content (porosity) | -  | 0.6 |
 | **`theta_r`** | residual water content  | -  | 0.01 |
 | `dt` | model time step | d  | - |
-| **`beta_l`** (`slope`) | slope | m m``^{-1}``  | - |
+| **`slope`** | slope | m m``^{-1}``  | - |
 | `dl` | drain length | m | - |
 | `dw` | drain width | m | - |
 | `zi` | pseudo-water table depth (top of the saturated zone) | m | - |
@@ -256,9 +256,9 @@ model parameter `mannings_n`, and is listed in the Table below between parenthes
 | `zb_max`    | maximum channel bed elevation | m | - |
 | `mannings_n_sq` | Manning's roughness squared at edge/link | (s m``^{-\frac{1}{3}}``)``^2`` | - |
 | `h`    | water depth | m | - |
-| `eta_max`    | maximum water elevation | m | - |
-| `eta_src`    | water elevation of source node of edge | m | - |
-| `eta_dst`    | water elevation of downstream node of edge | m | - |
+| `zs_max`    | maximum water elevation | m | - |
+| `zs_src`    | water elevation of source node of edge | m | - |
+| `zs_dst`    | water elevation of downstream node of edge | m | - |
 | `hf`    | water depth at edge/link | m | - |
 | `h_av`    | average water depth | m | - |
 | `dl`    | river length | m | - |
diff --git a/src/flextopo_model.jl b/src/flextopo_model.jl
index f0f18d07e..5575facba 100644
--- a/src/flextopo_model.jl
+++ b/src/flextopo_model.jl
@@ -517,9 +517,9 @@ function initialize_flextopo_model(config::Config)
         ldd = set_pit_ldd(pits_2d, ldd, inds)
     end
 
-    beta_l =
+    landslope =
         ncread(nc, config, "lateral.land.slope"; optional = false, sel = inds, type = Float)
-    clamp!(beta_l, 0.00001, Inf)
+    clamp!(landslope, 0.00001, Inf)
 
     dl = map(detdrainlength, ldd, xl, yl)
     dw = (xl .* yl) ./ dl
@@ -527,7 +527,7 @@ function initialize_flextopo_model(config::Config)
         nc,
         config,
         inds;
-        sl = beta_l,
+        sl = landslope,
         dl = dl,
         width = map(det_surfacewidth, dw, riverwidth, river),
         iterate = kinwave_it,
@@ -548,7 +548,7 @@ function initialize_flextopo_model(config::Config)
 
     # the indices of the river cells in the land(+river) cell vector
     index_river = filter(i -> !isequal(river[i], 0), 1:n)
-    frac_toriver = fraction_runoff_toriver(graph, ldd, index_river, beta_l, n)
+    frac_toriver = fraction_runoff_toriver(graph, ldd, index_river, landslope, n)
 
     rf = initialize_surfaceflow_river(
         nc,
diff --git a/src/flow.jl b/src/flow.jl
index 102745b15..a26cf310d 100644
--- a/src/flow.jl
+++ b/src/flow.jl
@@ -397,7 +397,7 @@ end
     theta_s::Vector{T} | "-"               # Saturated water content (porosity) [-]
     theta_r::Vector{T} | "-"               # Residual water content [-]
     dt::T | "d" | 0 | "none" | "none"      # model time step [d]
-    beta_l::Vector{T} | "m m-1"            # Slope [m m⁻¹]
+    slope::Vector{T} | "m m-1"            # Slope [m m⁻¹]
     dl::Vector{T} | "m"                    # Drain length [m]
     dw::Vector{T} | "m"                    # Flow width [m]
     zi::Vector{T} | "m"                    # Pseudo-water table depth [m] (top of the saturated zone)
@@ -446,7 +446,7 @@ function update(ssf::LateralSSF, network, frac_toriver, ksat_profile)
                         ssf.zi[v],
                         ssf.recharge[v],
                         ssf.kh_0[v],
-                        ssf.beta_l[v],
+                        ssf.slope[v],
                         ssf.theta_s[v] - ssf.theta_r[v],
                         ssf.f[v],
                         ssf.soilthickness[v],
@@ -465,7 +465,7 @@ function update(ssf::LateralSSF, network, frac_toriver, ksat_profile)
                         ssf.zi[v],
                         ssf.recharge[v],
                         ssf.kh[v],
-                        ssf.beta_l[v],
+                        ssf.slope[v],
                         ssf.theta_s[v] - ssf.theta_r[v],
                         ssf.soilthickness[v],
                         ssf.dt,
@@ -504,9 +504,9 @@ end
     mannings_n_sq::Vector{T} | "(s m-1/3)2" | _ | "edge"    # Manning's roughness squared at edge/link
     mannings_n::Vector{T} | "s m-1/3"                       # Manning's roughness at node
     h::Vector{T} | "m"                                      # water depth
-    eta_max::Vector{T} | "m" | _ | "edge"                   # maximum water elevation at edge
-    eta_src::Vector{T} | "m"                                # water elevation of source node of edge
-    eta_dst::Vector{T} | "m"                                # water elevation of downstream node of edge
+    zs_max::Vector{T} | "m" | _ | "edge"                    # maximum water elevation at edge
+    zs_src::Vector{T} | "m"                                 # water elevation of source node of edge
+    zs_dst::Vector{T} | "m"                                 # water elevation of downstream node of edge
     hf::Vector{T} | "m" | _ | "edge"                        # water depth at edge/link
     h_av::Vector{T} | "m"                                   # average water depth
     dl::Vector{T} | "m"                                     # river length
@@ -685,9 +685,9 @@ function initialize_shallowwater_river(
         mannings_n_sq = mannings_n_sq,
         mannings_n = n_river,
         h = h,
-        eta_max = zeros(_ne),
-        eta_src = zeros(_ne),
-        eta_dst = zeros(_ne),
+        zs_max = zeros(_ne),
+        zs_src = zeros(_ne),
+        zs_dst = zeros(_ne),
         hf = zeros(_ne),
         h_av = zeros(n),
         width = width,
@@ -736,11 +736,11 @@ function shallowwater_river_update(sw::ShallowWaterRiver, network, dt, doy, upda
         i = sw.active_e[j]
         i_src = nodes_at_link.src[i]
         i_dst = nodes_at_link.dst[i]
-        sw.eta_src[i] = sw.zb[i_src] + sw.h[i_src]
-        sw.eta_dst[i] = sw.zb[i_dst] + sw.h[i_dst]
+        sw.zs_src[i] = sw.zb[i_src] + sw.h[i_src]
+        sw.zs_dst[i] = sw.zb[i_dst] + sw.h[i_dst]
 
-        sw.eta_max[i] = max(sw.eta_src[i], sw.eta_dst[i])
-        sw.hf[i] = (sw.eta_max[i] - sw.zb_max[i])
+        sw.zs_max[i] = max(sw.zs_src[i], sw.zs_dst[i])
+        sw.hf[i] = (sw.zs_max[i] - sw.zb_max[i])
 
         sw.a[i] = sw.width_at_link[i] * sw.hf[i] # flow area (rectangular channel)
         sw.r[i] = sw.a[i] / (sw.width_at_link[i] + 2.0 * sw.hf[i]) # hydraulic radius (rectangular channel)
@@ -749,8 +749,8 @@ function shallowwater_river_update(sw::ShallowWaterRiver, network, dt, doy, upda
             sw.hf[i] > sw.h_thresh,
             local_inertial_flow(
                 sw.q0[i],
-                sw.eta_src[i],
-                sw.eta_dst[i],
+                sw.zs_src[i],
+                sw.zs_dst[i],
                 sw.hf[i],
                 sw.a[i],
                 sw.r[i],
@@ -770,7 +770,7 @@ function shallowwater_river_update(sw::ShallowWaterRiver, network, dt, doy, upda
         sw.q_av[i] += sw.q[i] * dt
     end
     if !isnothing(sw.floodplain)
-        @tturbo @. sw.floodplain.hf = max(sw.eta_max - sw.floodplain.zb_max, 0.0)
+        @tturbo @. sw.floodplain.hf = max(sw.zs_max - sw.floodplain.zb_max, 0.0)
 
         n = 0
         @inbounds for i in sw.active_e
@@ -830,8 +830,8 @@ function shallowwater_river_update(sw::ShallowWaterRiver, network, dt, doy, upda
                 sw.floodplain.a[i] > 1.0e-05,
                 local_inertial_flow(
                     sw.floodplain.q0[i],
-                    sw.eta_src[i],
-                    sw.eta_dst[i],
+                    sw.zs_src[i],
+                    sw.zs_dst[i],
                     sw.floodplain.hf[i],
                     sw.floodplain.a[i],
                     sw.floodplain.r[i],
@@ -1230,10 +1230,10 @@ function shallowwater_update(
         # for flow in x dir) and floodplain flow is not calculated.
         if xu <= sw.n && sw.ywidth[i] != T(0.0)
 
-            eta_x = sw.z[i] + sw.h[i]
-            eta_xu = sw.z[xu] + sw.h[xu]
-            eta_max = max(eta_x, eta_xu)
-            hf = (eta_max - sw.zx_max[i])
+            zs_x = sw.z[i] + sw.h[i]
+            zs_xu = sw.z[xu] + sw.h[xu]
+            zs_max = max(zs_x, zs_xu)
+            hf = (zs_max - sw.zx_max[i])
 
             if hf > sw.h_thresh
                 length = T(0.5) * (sw.xl[i] + sw.xl[xu]) # can be precalculated
@@ -1242,8 +1242,8 @@ function shallowwater_update(
                     sw.qx0[i],
                     sw.qx0[xd],
                     sw.qx0[xu],
-                    eta_x,
-                    eta_xu,
+                    zs_x,
+                    zs_xu,
                     hf,
                     sw.ywidth[i],
                     length,
@@ -1270,10 +1270,10 @@ function shallowwater_update(
         # for flow in y dir) and floodplain flow is not calculated.
         if yu <= sw.n && sw.xwidth[i] != T(0.0)
 
-            eta_y = sw.z[i] + sw.h[i]
-            eta_yu = sw.z[yu] + sw.h[yu]
-            eta_max = max(eta_y, eta_yu)
-            hf = (eta_max - sw.zy_max[i])
+            zs_y = sw.z[i] + sw.h[i]
+            zs_yu = sw.z[yu] + sw.h[yu]
+            zs_max = max(zs_y, zs_yu)
+            hf = (zs_max - sw.zy_max[i])
 
             if hf > sw.h_thresh
                 length = T(0.5) * (sw.yl[i] + sw.yl[yu]) # can be precalculated
@@ -1282,8 +1282,8 @@ function shallowwater_update(
                     sw.qy0[i],
                     sw.qy0[yd],
                     sw.qy0[yu],
-                    eta_y,
-                    eta_yu,
+                    zs_y,
+                    zs_yu,
                     hf,
                     sw.xwidth[i],
                     length,
diff --git a/src/groundwater/aquifer.jl b/src/groundwater/aquifer.jl
index 9909d4d8e..d650bf799 100644
--- a/src/groundwater/aquifer.jl
+++ b/src/groundwater/aquifer.jl
@@ -273,9 +273,9 @@ function conductance(
             connectivity.width[nzi],
         )
     elseif conductivity_profile == "uniform"
-        phi_i = aquifer.head[i]
-        phi_j = aquifer.head[j]
-        if phi_i >= phi_j
+        head_i = aquifer.head[i]
+        head_j = aquifer.head[j]
+        if head_i >= head_j
             saturation =
                 saturated_thickness(aquifer, i) / (aquifer.top[i] - aquifer.bottom[i])
         else
@@ -298,9 +298,9 @@ function flux!(Q, aquifer, connectivity, conductivity_profile)
         for nzi in connections(connectivity, i)
             # connection from i -> j
             j = connectivity.rowval[nzi]
-            delta_phi = aquifer.head[i] - aquifer.head[j]
+            delta_head = aquifer.head[i] - aquifer.head[j]
             cond = conductance(aquifer, i, j, nzi, conductivity_profile, connectivity)
-            Q[i] -= cond * delta_phi
+            Q[i] -= cond * delta_head
         end
     end
     return Q
diff --git a/src/groundwater/boundary_conditions.jl b/src/groundwater/boundary_conditions.jl
index fad10c526..06a04897e 100644
--- a/src/groundwater/boundary_conditions.jl
+++ b/src/groundwater/boundary_conditions.jl
@@ -25,16 +25,16 @@ end
 
 function flux!(Q, river::River, aquifer)
     for (i, index) in enumerate(river.index)
-        phi= aquifer.head[index]
+        head = aquifer.head[index]
         stage = river.stage[i]
-        if stage > phi
+        if stage > head
             cond = river.infiltration_conductance[i]
-            delta_phi = min(stage - river.bottom[i], stage - phi)
+            delta_head = min(stage - river.bottom[i], stage - head)
         else
             cond = river.exfiltration_conductance[i]
-            delta_phi = stage - phi
+            delta_head = stage - head
         end
-        river.flux[i] = check_flux(cond * delta_phi, aquifer, index)
+        river.flux[i] = check_flux(cond * delta_head, aquifer, index)
         Q[index] += river.flux[i]
     end
     return Q
@@ -53,8 +53,8 @@ end
 function flux!(Q, drainage::Drainage, aquifer)
     for (i, index) in enumerate(drainage.index)
         cond = drainage.conductance[i]
-        delta_phi = min(0, drainage.elevation[i] - aquifer.head[index])
-        drainage.flux[i] = check_flux(cond * delta_phi, aquifer, index)
+        delta_head = min(0, drainage.elevation[i] - aquifer.head[index])
+        drainage.flux[i] = check_flux(cond * delta_head, aquifer, index)
         Q[index] += drainage.flux[i]
     end
     return Q
@@ -72,8 +72,8 @@ end
 function flux!(Q, headboundary::HeadBoundary, aquifer)
     for (i, index) in enumerate(headboundary.index)
         cond = headboundary.conductance[i]
-        delta_phi = headboundary.head[i] - aquifer.head[index]
-        headboundary.flux[i] = check_flux(cond * delta_phi, aquifer, index)
+        delta_head = headboundary.head[i] - aquifer.head[index]
+        headboundary.flux[i] = check_flux(cond * delta_head, aquifer, index)
         Q[index] += headboundary.flux[i]
     end
     return Q
diff --git a/src/hbv_model.jl b/src/hbv_model.jl
index e42351b66..aab5cbbd7 100644
--- a/src/hbv_model.jl
+++ b/src/hbv_model.jl
@@ -249,9 +249,9 @@ function initialize_hbv_model(config::Config)
         ldd = set_pit_ldd(pits_2d, ldd, inds)
     end
 
-    beta_l =
+    landslope =
         ncread(nc, config, "lateral.land.slope"; optional = false, sel = inds, type = Float)
-    clamp!(beta_l, 0.00001, Inf)
+    clamp!(landslope, 0.00001, Inf)
 
     dl = map(detdrainlength, ldd, xl, yl)
     dw = (xl .* yl) ./ dl
@@ -259,7 +259,7 @@ function initialize_hbv_model(config::Config)
         nc,
         config,
         inds;
-        sl = beta_l,
+        sl = landslope,
         dl = dl,
         width = map(det_surfacewidth, dw, riverwidth, river),
         iterate = kinwave_it,
@@ -282,7 +282,7 @@ function initialize_hbv_model(config::Config)
 
     # the indices of the river cells in the land(+river) cell vector
     index_river = filter(i -> !isequal(river[i], 0), 1:n)
-    frac_toriver = fraction_runoff_toriver(graph, ldd, index_river, beta_l, n)
+    frac_toriver = fraction_runoff_toriver(graph, ldd, index_river, landslope, n)
 
     rf = initialize_surfaceflow_river(
         nc,
diff --git a/src/horizontal_process.jl b/src/horizontal_process.jl
index ea2f7fb9a..2d10333ce 100644
--- a/src/horizontal_process.jl
+++ b/src/horizontal_process.jl
@@ -69,28 +69,28 @@ function kin_wave!(Q, graph, toposort, Qold, q, alpha, beta, DCL, dt)
 end
 
 "Returns water table depth `zi` based on lateral subsurface flow `ssf` and hydraulic conductivity profile `ksat_profile`"
-function ssf_water_table_depth(ssf, kh_0, beta, f, d, dw, z_exp, ksat_profile)
+function ssf_water_table_depth(ssf, kh_0, slope, f, d, dw, z_exp, ksat_profile)
     if ksat_profile == "exponential"
-        zi = log((f * ssf) / (dw * kh_0 * beta) + exp(-f * d)) / -f
+        zi = log((f * ssf) / (dw * kh_0 * slope) + exp(-f * d)) / -f
     elseif ksat_profile == "exponential_constant"
-        ssf_constant = kh_0 * beta * exp(-f * z_exp) * (d - z_exp) * dw
+        ssf_constant = kh_0 * slope * exp(-f * z_exp) * (d - z_exp) * dw
         if ssf > ssf_constant
-            zi = log((f * (ssf - ssf_constant)) / (dw * kh_0 * beta) + exp(-f * z_exp)) / -f
+            zi = log((f * (ssf - ssf_constant)) / (dw * kh_0 * slope) + exp(-f * z_exp)) / -f
         else
-            zi = d - ssf / (dw * kh_0 * beta * exp(-f * z_exp))
+            zi = d - ssf / (dw * kh_0 * slope * exp(-f * z_exp))
         end
     end
     return zi
 end
 
 "Returns kinematic wave celecity `Cn` of lateral subsurface flow based on hydraulic conductivity profile `ksat_profile`"
-function ssf_celerity(zi, kh_0, beta, theta_e, f, z_exp, ksat_profile)
+function ssf_celerity(zi, kh_0, slope, theta_e, f, z_exp, ksat_profile)
     if ksat_profile == "exponential"
-        Cn = (kh_0 * exp(-f * zi) * beta) / theta_e
+        Cn = (kh_0 * exp(-f * zi) * slope) / theta_e
     elseif ksat_profile == "exponential_constant"
-        Cn_const = (kh_0 * exp(-f * z_exp) * beta) / theta_e
+        Cn_const = (kh_0 * exp(-f * z_exp) * slope) / theta_e
         if zi < z_exp
-            Cn = (kh_0 * exp(-f * zi) * beta) / theta_e + Cn_const
+            Cn = (kh_0 * exp(-f * zi) * slope) / theta_e + Cn_const
         else
             Cn = Cn_const
         end
@@ -99,7 +99,7 @@ function ssf_celerity(zi, kh_0, beta, theta_e, f, z_exp, ksat_profile)
 end
 
 """
-    kinematic_wave_ssf(ssfin, ssf_prev, zi_prev, r, kh_0, beta, theta_e, f, d, dt, dx, dw, ssfmax, z_exp, ksat_profile)
+    kinematic_wave_ssf(ssfin, ssf_prev, zi_prev, r, kh_0, slope, theta_e, f, d, dt, dx, dw, ssfmax, z_exp, ksat_profile)
 
 Kinematic wave for lateral subsurface flow for a single cell and timestep. An exponential
 decline of hydraulic conductivity at the soil surface `kh_0`, controllled by parameter `f`,
@@ -116,7 +116,7 @@ function kinematic_wave_ssf(
     zi_prev,
     r,
     kh_0,
-    beta,
+    slope,
     theta_e,
     f,
     d,
@@ -139,9 +139,9 @@ function kinematic_wave_ssf(
         count = 1
 
         # Estimate zi on the basis of the relation between subsurface flow and zi
-        zi = ssf_water_table_depth(ssf, kh_0, beta, f, d, dw, z_exp, ksat_profile)
+        zi = ssf_water_table_depth(ssf, kh_0, slope, f, d, dw, z_exp, ksat_profile)
         # Reciprocal of derivative delta Q/ delta z_i, constrained w.r.t. neff on the basis of the continuity equation)
-        Cn = ssf_celerity(zi, kh_0, beta, theta_e, f, z_exp, ksat_profile)
+        Cn = ssf_celerity(zi, kh_0, slope, theta_e, f, z_exp, ksat_profile)
         # Term of the continuity equation for Newton-Raphson iteration for iteration 1
         # because celerity Cn is depending on zi, the increase or decrease of zi is moved to the recharge term of the continuity equation
         # then (1./Cn)*ssf_prev can be replaced with (1./Cn)*ssf, and thus celerity and lateral flow rate ssf are then in line
@@ -161,9 +161,9 @@ function kinematic_wave_ssf(
         # Start while loop of Newton-Raphson iteration m until continuity equation approaches zero
         while true
             # Estimate zi on the basis of the relation between lateral flow rate and groundwater level
-            zi = ssf_water_table_depth(ssf, kh_0, beta, f, d, dw, z_exp, ksat_profile)
+            zi = ssf_water_table_depth(ssf, kh_0, slope, f, d, dw, z_exp, ksat_profile)
             # Reciprocal of derivative delta Q/ delta z_i, constrained w.r.t. neff on the basis of the continuity equation
-            Cn = ssf_celerity(zi, kh_0, beta, theta_e, f, z_exp, ksat_profile)
+            Cn = ssf_celerity(zi, kh_0, slope, theta_e, f, z_exp, ksat_profile)
             # Term of the continuity equation for given Newton-Raphson iteration m
             # because celerity Cn is depending on zi, the increase or decrease of zi is moved to the recharge term of the continuity equation
             # then (1./Cn)*ssf_prev can be replaced with (1./Cn)*ssf, and thus celerity and lateral flow rate ssf are then in line
@@ -199,7 +199,7 @@ function kinematic_wave_ssf(
 end
 
 """
-    kinematic_wave_ssf(ssfin, ssf_prev, zi_prev, r, kh, beta, theta_e, d, dt, dx, dw, ssfmax)
+    kinematic_wave_ssf(ssfin, ssf_prev, zi_prev, r, kh, slope, theta_e, d, dt, dx, dw, ssfmax)
 
 Kinematic wave for lateral subsurface flow for a single cell and timestep, based on
 (average) hydraulic conductivity `kh`.
@@ -207,7 +207,7 @@ Kinematic wave for lateral subsurface flow for a single cell and timestep, based
 Returns lateral subsurface flow `ssf`, water table depth `zi` and exfiltration rate
 `exfilt`.
 """
-function kinematic_wave_ssf(ssfin, ssf_prev, zi_prev, r, kh, beta, theta_e, d, dt, dx, dw, ssfmax)
+function kinematic_wave_ssf(ssfin, ssf_prev, zi_prev, r, kh, slope, theta_e, d, dt, dx, dw, ssfmax)
 
     epsilon = 1.0e-12
     max_iters = 3000
@@ -219,7 +219,7 @@ function kinematic_wave_ssf(ssfin, ssf_prev, zi_prev, r, kh, beta, theta_e, d, d
         ssf = (ssf_prev + ssfin) / 2.0
         count = 1
         # celerity (Cn)
-        Cn = (beta * kh) / theta_e
+        Cn = (slope * kh) / theta_e
         # constant term of the continuity equation for Newton-Raphson
         c = (dt / dx) * ssfin + (1.0 / Cn) * ssf_prev + r
         # continuity equation of which solution should be zero
@@ -355,15 +355,15 @@ function lateral_snow_transport!(snow, snowwater, slope, network)
 end
 
 """
-    local_inertial_flow(q0, eta0, eta1, hf, A, R, length, mannings_n, g, froude_limit, dt)
+    local_inertial_flow(q0, zs0, zs1, hf, A, R, length, mannings_n, g, froude_limit, dt)
 
 Local inertial approach for flow through area `A`. Returns the flow `q` between two adjacent
 river cells (nodes) for a single timestep.
 """
 function local_inertial_flow(
     q0,
-    eta0,
-    eta1,
+    zs0,
+    zs1,
     hf,
     A,
     R,
@@ -374,7 +374,7 @@ function local_inertial_flow(
     dt,
 )
 
-    slope = (eta1 - eta0) / length
+    slope = (zs1 - zs0) / length
     pow_R = cbrt(R * R * R * R)
     unit = one(hf)
     q = (
@@ -390,7 +390,7 @@ function local_inertial_flow(
 end
 
 """
-    local_inertial_flow(theta, q0, qd, qu, eta0, eta1, hf, width, length, mannings_n, g, froude_limit, dt)
+    local_inertial_flow(theta, q0, qd, qu, zs0, zs1, hf, width, length, mannings_n, g, froude_limit, dt)
 
 Local inertial approach for flow through a rectangular area. Returns the flow `q` between
 two adjacent cells (nodes) for a single timestep. Algorithm is based on de Almeida et al.
@@ -401,8 +401,8 @@ function local_inertial_flow(
     q0,
     qd,
     qu,
-    eta0,
-    eta1,
+    zs0,
+    zs1,
     hf,
     width,
     length,
@@ -412,7 +412,7 @@ function local_inertial_flow(
     dt,
 )
 
-    slope = (eta1 - eta0) / length
+    slope = (zs1 - zs0) / length
     unit = one(theta)
     half = oftype(theta, 0.5)
     pow_hf = cbrt(hf * hf * hf * hf * hf * hf * hf)
diff --git a/src/sbm_gwf_model.jl b/src/sbm_gwf_model.jl
index 01ed2777b..6a76b5bdf 100644
--- a/src/sbm_gwf_model.jl
+++ b/src/sbm_gwf_model.jl
@@ -99,9 +99,9 @@ function initialize_sbm_gwf_model(config::Config)
     end
 
     # overland flow (kinematic wave)
-    beta_l =
+    landslope =
         ncread(nc, config, "lateral.land.slope"; optional = false, sel = inds, type = Float)
-    clamp!(beta_l, 0.00001, Inf)
+    clamp!(landslope, 0.00001, Inf)
     ldd_2d = ncread(nc, config, "ldd"; optional = false, allow_missing = true)
 
     ldd = ldd_2d[inds]
@@ -116,14 +116,14 @@ function initialize_sbm_gwf_model(config::Config)
 
     # the indices of the river cells in the land(+river) cell vector
     index_river = filter(i -> !isequal(river[i], 0), 1:n)
-    frac_toriver = fraction_runoff_toriver(graph, ldd, index_river, beta_l, n)
+    frac_toriver = fraction_runoff_toriver(graph, ldd, index_river, landslope, n)
 
     if land_routing == "kinematic-wave"
         olf = initialize_surfaceflow_land(
             nc,
             config,
             inds;
-            sl = beta_l,
+            sl = landslope,
             dl,
             width = map(det_surfacewidth, dw, riverwidth, river),
             iterate = kinwave_it,
@@ -421,7 +421,7 @@ function initialize_sbm_gwf_model(config::Config)
             reverse_indices = rev_inds,
             xl = xl,
             yl = yl,
-            slope = beta_l,
+            slope = landslope,
             altitude = altitude,
         )
     elseif land_routing == "local-inertial"
@@ -432,7 +432,7 @@ function initialize_sbm_gwf_model(config::Config)
             reverse_indices = rev_inds,
             xl = xl,
             yl = yl,
-            slope = beta_l,
+            slope = landslope,
             altitude = altitude,
             index_river = index_river_nf,
             staggered_indices = indices,
diff --git a/src/sbm_model.jl b/src/sbm_model.jl
index a5cfaa206..bc238bab9 100644
--- a/src/sbm_model.jl
+++ b/src/sbm_model.jl
@@ -102,9 +102,9 @@ function initialize_sbm_model(config::Config)
         ldd = set_pit_ldd(pits_2d, ldd, inds)
     end
 
-    beta_l =
+    landslope =
         ncread(nc, config, "lateral.land.slope"; optional = false, sel = inds, type = Float)
-    clamp!(beta_l, 0.00001, Inf)
+    clamp!(landslope, 0.00001, Inf)
 
     dl = map(detdrainlength, ldd, xl, yl)
     dw = (xl .* yl) ./ dl
@@ -140,7 +140,7 @@ function initialize_sbm_model(config::Config)
             theta_s = sbm.theta_s,
             theta_r = sbm.theta_r,
             dt = dt / basetimestep,
-            beta_l = beta_l,
+            slope = landslope,
             dl = dl,
             dw = dw,
             exfiltwater = fill(mv, n),
@@ -180,14 +180,14 @@ function initialize_sbm_model(config::Config)
 
     # the indices of the river cells in the land(+river) cell vector
     index_river = filter(i -> !isequal(river[i], 0), 1:n)
-    frac_toriver = fraction_runoff_toriver(graph, ldd, index_river, beta_l, n)
+    frac_toriver = fraction_runoff_toriver(graph, ldd, index_river, landslope, n)
 
     if land_routing == "kinematic-wave"
         olf = initialize_surfaceflow_land(
             nc,
             config,
             inds;
-            sl = beta_l,
+            sl = landslope,
             dl,
             width = map(det_surfacewidth, dw, riverwidth, river),
             iterate = kinwave_it,
@@ -340,7 +340,7 @@ function initialize_sbm_model(config::Config)
         reverse_indices = rev_inds,
         xl,
         yl,
-        slope = beta_l,
+        slope = landslope,
     )
     if land_routing == "local-inertial"
         land = merge(land, (index_river = index_river_nf, staggered_indices = indices))
diff --git a/src/sediment.jl b/src/sediment.jl
index 2de64ecc2..dd708a32a 100644
--- a/src/sediment.jl
+++ b/src/sediment.jl
@@ -135,8 +135,7 @@ function initialize_landsed(nc, config, river, riverfrac, xl, yl, inds)
     cmax, _, canopygapfraction, sl, swood, kext = initialize_canopy(nc, config, inds)
 
     # Initialise parameters for the transport capacity part
-    beta_l = slope
-    clamp!(beta_l, 0.00001, Inf)
+    clamp!(slope, 0.00001, Inf)
     ldd_2d = ncread(nc, config, "ldd"; optional = false, allow_missing = true)
     ldd = ldd_2d[inds]
 
diff --git a/src/sediment_model.jl b/src/sediment_model.jl
index 5a91e0ed0..e2229a1bd 100644
--- a/src/sediment_model.jl
+++ b/src/sediment_model.jl
@@ -95,9 +95,9 @@ function initialize_sediment_model(config::Config)
     # River processes
 
     # the indices of the river cells in the land(+river) cell vector
-    beta_l =
+    landslope =
         ncread(nc, config, "lateral.land.slope"; optional = false, sel = inds, type = Float)
-    clamp!(beta_l, 0.00001, Inf)
+    clamp!(landslope, 0.00001, Inf)
 
     riverlength = riverlength_2d[inds_riv]
     riverwidth = riverwidth_2d[inds_riv]
@@ -108,7 +108,7 @@ function initialize_sediment_model(config::Config)
     graph_riv = flowgraph(ldd_riv, inds_riv, pcr_dir)
 
     index_river = filter(i -> !isequal(river[i], 0), 1:n)
-    frac_toriver = fraction_runoff_toriver(graph, ldd, index_river, beta_l, n)
+    frac_toriver = fraction_runoff_toriver(graph, ldd, index_river, landslope, n)
 
     rs = initialize_riversed(nc, config, riverwidth, riverlength, inds_riv)
 
diff --git a/src/utils.jl b/src/utils.jl
index 18c7095a9..57563b7b6 100644
--- a/src/utils.jl
+++ b/src/utils.jl
@@ -768,10 +768,10 @@ end
 function initialize_lateralssf_exp!(ssf::LateralSSF)
     for i in eachindex(ssf.ssf)
         ssf.ssfmax[i] =
-            ((ssf.kh_0[i] * ssf.beta_l[i]) / ssf.f[i]) *
+            ((ssf.kh_0[i] * ssf.slope[i]) / ssf.f[i]) *
             (1.0 - exp(-ssf.f[i] * ssf.soilthickness[i]))
         ssf.ssf[i] =
-            ((ssf.kh_0[i] * ssf.beta_l[i]) / ssf.f[i]) *
+            ((ssf.kh_0[i] * ssf.slope[i]) / ssf.f[i]) *
             (exp(-ssf.f[i] * ssf.zi[i]) - exp(-ssf.f[i] * ssf.soilthickness[i])) *
             ssf.dw[i]
     end
@@ -782,16 +782,16 @@ function initialize_lateralssf_exp_const!(ssf::LateralSSF)
     ssf_constant = @. ssf.khfrac *
        ssf.kh_0 *
        exp(-ssf.f * ssf.z_exp) *
-       ssf.beta_l *
+       ssf.slope *
        (ssf.soilthickness - ssf.z_exp)
     for i in eachindex(ssf.ssf)
         ssf.ssfmax[i] =
-            ((ssf.khfrac[i] * ssf.kh_0[i] * ssf.beta_l[i]) / ssf.f[i]) *
+            ((ssf.khfrac[i] * ssf.kh_0[i] * ssf.slope[i]) / ssf.f[i]) *
             (1.0 - exp(-ssf.f[i] * ssf.z_exp[i])) + ssf_constant[i]
         if ssf.zi[i] < ssf.z_exp[i]
             ssf.ssf[i] =
                 (
-                    ((ssf.kh_0[i] * ssf.beta_l[i]) / ssf.f[i]) *
+                    ((ssf.kh_0[i] * ssf.slope[i]) / ssf.f[i]) *
                     (exp(-ssf.f[i] * ssf.zi[i]) - exp(-ssf.f[i] * ssf.z_exp[i])) +
                     ssf_constant[i]
                 ) * ssf.dw[i]
@@ -799,7 +799,7 @@ function initialize_lateralssf_exp_const!(ssf::LateralSSF)
             ssf.ssf[i] =
                 ssf.kh_0[i] *
                 exp(-ssf.f[i] * ssf.zi[i]) *
-                ssf.beta_l[i] *
+                ssf.slope[i] *
                 (ssf.soilthickness[i] - ssf.zi[i]) *
                 ssf.dw[i]
         end
@@ -810,7 +810,7 @@ end
 function initialize_lateralssf_layered!(ssf::LateralSSF, sbm::SBM, ksat_profile)
     for i in eachindex(ssf.ssf)
         ssf.kh[i] = kh_layered_profile(sbm, ssf.khfrac[i], i, ksat_profile)
-        ssf.ssf[i] = ssf.kh[i] * (ssf.soilthickness[i] - ssf.zi[i]) * ssf.beta_l[i] * ssf.dw[i]
+        ssf.ssf[i] = ssf.kh[i] * (ssf.soilthickness[i] - ssf.zi[i]) * ssf.slope[i] * ssf.dw[i]
         kh_max = 0.0
         for j = 1:sbm.nlayers[i]
             if j <= sbm.nlayers_kv[i]
@@ -823,6 +823,6 @@ function initialize_lateralssf_layered!(ssf::LateralSSF, sbm::SBM, ksat_profile)
             end
         end
         kh_max = kh_max * ssf.khfrac[i] * 0.001 * 0.001
-        ssf.ssfmax[i] = kh_max * ssf.beta_l[i]
+        ssf.ssfmax[i] = kh_max * ssf.slope[i]
     end
 end
diff --git a/test/groundwater.jl b/test/groundwater.jl
index 074a85541..13d529ed1 100644
--- a/test/groundwater.jl
+++ b/test/groundwater.jl
@@ -434,10 +434,10 @@ end
             @test all(diff(gwf.aquifer.head) .> 0.0)
         end
 
-        phi_analytical = [
+        head_analytical = [
             transient_aquifer_1d(x, time, conductivity, specific_yield, aquifer_length, beta) for x in xc
         ]
-        difference = gwf.aquifer.head .- phi_analytical
+        difference = gwf.aquifer.head .- head_analytical
         # @test all(difference .< ?)  #TODO
     end
 
@@ -494,10 +494,10 @@ end
             @test all(diff(gwf.aquifer.head) .> 0.0)
         end
 
-        phi_analytical = [
+        head_analytical = [
             transient_aquifer_1d(x, time, conductivity, specific_yield, aquifer_length, beta) for x in xc
         ]
-        difference = gwf.aquifer.head .- phi_analytical
+        difference = gwf.aquifer.head .- head_analytical
         # @test all(difference .< ?)  #TODO
     end
 
@@ -555,18 +555,18 @@ end
         end
 
         # test for symmetry on x and y axes
-        phi = reshape(gwf.aquifer.head, shape)
-        @test phi[1:halfnrow, :] ≈ phi[end:-1:halfnrow+2, :]
-        @test phi[:, 1:halfnrow] ≈ phi[:, end:-1:halfnrow+2]
+        head = reshape(gwf.aquifer.head, shape)
+        @test head[1:halfnrow, :] ≈ head[end:-1:halfnrow+2, :]
+        @test head[:, 1:halfnrow] ≈ head[:, end:-1:halfnrow+2]
 
         # compare with analytical solution
         start = -0.5 * aquifer_length + 0.5 * cellsize
         stop = 0.5 * aquifer_length - 0.5 * cellsize
         X = collect(range(start, stop = stop, step = cellsize))
-        phi_analytical =
+        head_analytical =
             [drawdown_theis(x, time, discharge, transmissivity, storativity) for x in X] .+ 10.0
         # compare left-side, since it's symmetric anyway. Skip the well cell, and its first neighbor
-        difference = phi[1:halfnrow-1, halfnrow] - phi_analytical[1:halfnrow-1]
+        difference = head[1:halfnrow-1, halfnrow] - head_analytical[1:halfnrow-1]
         @test all(difference .< 0.02)
     end
 
diff --git a/test/horizontal_process.jl b/test/horizontal_process.jl
index 0ae29c089..3e9eb41ed 100644
--- a/test/horizontal_process.jl
+++ b/test/horizontal_process.jl
@@ -198,9 +198,9 @@ end
         mannings_n_sq = mannings_n_sq,
         mannings_n = n_river,
         h = h_init,
-        eta_max = zeros(_ne),
-        eta_src = zeros(_ne),
-        eta_dst = zeros(_ne),
+        zs_max = zeros(_ne),
+        zs_src = zeros(_ne),
+        zs_dst = zeros(_ne),
         hf = zeros(_ne),
         h_av = zeros(n),
         width = width,