Repository: Deltares/Wflow.jl (language: python)

Issue:
Unicode Character Usage and Consistency
### Feature type

Changing existing functionality

### Improvement Description

The current codebase utilizes Unicode characters inconsistently. The advantage of Unicode characters is that the implementation remains close to the model equations, making the code more readable and easier to understand (improving maintainability of the code). In addition, to make fully use of Unicode characters, these characters should also be included in the documentation (code notation matches notation in documentation).  Possible disadvantages of using Unicode characters are: more difficult to type, searching (example is required) and editor or text processor could have problems correctly displaying these characters (config TOML file). It's essential to evaluate the usage of Unicode characters, considering their advantages and disadvantages, and ensure consistent usage throughout the codebase.

### Implementation Description

Make a well considered choice (use consistent or not) and refactor the code (and change docs) based on the choice.

### Additional Context

_No response_


Discussion / hints:
We have decided to stop using unicode characters, and spell them out if necessary (for example `theta_s` instead of the Greek symbol). We aim for the variables to be as descriptive as possible, while ensuring the names won't get too long as that could limit readability of the code

Fix patch:
diff --git a/docs/src/changelog.md b/docs/src/changelog.md
index 3b96a025e..9ef507d46 100644
--- a/docs/src/changelog.md
+++ b/docs/src/changelog.md
@@ -57,6 +57,10 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
 - Always use fractions for the computation of potential evapotranspiration (interception and
   transpiration) and potential evaporation (open water and soil). Replaced variable
   `et_reftopot` by crop coefficient `kc` (use of `et_reftopot` has been deprecated).
+- For improved code readability it is now discouraged to use non-ASCII characters for the
+  names of variables, structs, functions and macros. Using the non-ASCII character for
+  built-in operators is still allowed. This change in naming convention is now in effect and
+  all invalid uses of non-ASCII characters have been replaced by ASCII equivalents.
 
 ### Added
 - Total water storage as an export variable for `SBM` concept. This is the total water stored
diff --git a/docs/src/index.md b/docs/src/index.md
index 3fa7c9632..3a66f391d 100644
--- a/docs/src/index.md
+++ b/docs/src/index.md
@@ -4,7 +4,7 @@ CurrentModule = Wflow
 
 # About wflow
 
-Wflow is Deltares’ solution for modelling hydrological processes, allowing users to account
+Wflow is Deltares' solution for modelling hydrological processes, allowing users to account
 for precipitation, interception, snow accumulation and melt, evapotranspiration, soil water,
 surface water and groundwater recharge in a fully distributed environment. Successfully
 applied worldwide for analyzing flood hazards, drought, climate change impacts and land use
diff --git a/docs/src/model_docs/lateral/kinwave.md b/docs/src/model_docs/lateral/kinwave.md
index 1fdde084a..3e69a83f6 100644
--- a/docs/src/model_docs/lateral/kinwave.md
+++ b/docs/src/model_docs/lateral/kinwave.md
@@ -16,7 +16,7 @@ where ``Q`` is the surface runoff in the kinematic wave [m``^3``/s], ``x`` is th
 the runoff pathway [m], ``A`` is the cross-section area of the runoff pathway [m``^{2}``],
 ``t`` is the integration timestep [s] and ``\alpha`` and ``\beta`` are coefficients.
 
-These equations are solved with a nonlinear scheme using Newton’s method and can also be
+These equations are solved with a nonlinear scheme using Newton's method and can also be
 iterated depending on the  model space and time resolution. By default, the iterations are
 performed until a stable solution is reached (``\epsilon < 10^{-12}``). For larger models,
 the number of iterations can also be fixed for to a specific sub-timestep (in seconds) for
@@ -176,4 +176,4 @@ zero is applied to the upstream flow `qin`.
 
 ## References
 + Chow, V., Maidment, D. and Mays, L., 1988, Applied Hydrology. McGraw-Hill Book Company,
-  New York.
\ No newline at end of file
+  New York.
diff --git a/docs/src/model_docs/lateral/sediment_flux.md b/docs/src/model_docs/lateral/sediment_flux.md
index ffde2841e..eedc5ac51 100644
--- a/docs/src/model_docs/lateral/sediment_flux.md
+++ b/docs/src/model_docs/lateral/sediment_flux.md
@@ -38,7 +38,7 @@ mobilize 5 classes of sediment:
 where ``CLA``, ``SIL`` and ``SAN`` are the primary clay, silt, sand fractions of the topsoil
 and ``PCL``, ``PSI``, ``PSA``, ``SAG`` and ``LAG`` are the clay, silt, sand, small and large
 aggregates fractions of the detached sediment respectively. The transport capacity of the
-flow using Yalin’s equation with particle differentiation, developed by Foster (1982), is:
+flow using Yalin's equation with particle differentiation, developed by Foster (1982), is:
 ```math
    TC_{i} = (P_{e})_{i}  (S_{g})_{i} \, \rho_{w} \,  g \, d_{i}  V_{*}
 ```
@@ -79,7 +79,7 @@ are much rarer than for soil loss and inland dynamics. The simpler models such a
 default sediment river model uses again the transport capacity of the ﬂow to determine if
 there is erosion or deposition (Neitsch et al., 2011).  A more physics-based approach
 (Partheniades, 1965) to determine river erosion is used by Liu et al. (2018) and in the new
-SWAT’s approach developed by Narasimhan et al. (2017). For wﬂow\_sediment, the new
+SWAT's approach developed by Narasimhan et al. (2017). For wﬂow\_sediment, the new
 physics-based model of SWAT was chosen for transport and erosion as it enables the use of
 parameter estimation for erosion of bed and bank of the channel and separates the suspended
 from the bed loads.
@@ -313,7 +313,7 @@ deposited sediment and the river bed/bank erosion.
 As sediments have a higher density than water, moving sediments in water can be deposited in
 the river bed. The deposition process depends on the mass of the sediment, but also on flow
 characteristics such as velocity. In wflow_sediment, as in SWAT, deposition is modelled with
-Einstein’s equation (Neitsch et al, 2011):
+Einstein's equation (Neitsch et al, 2011):
 ```math
    P_{dep}=\left(1-\dfrac{1}{e^{x}}\right)100
 ```
@@ -359,14 +359,14 @@ and shifted to the outlet cell. wflow\_sediment handles the lakes and reservoirs
 cell belongs to a lake/reservoir and is not the outlet then the model assumes that no
 erosion/deposition of sediments is happening and the sediments are only all transported to
 the lake/reservoir outlet. Once the sediments reach the outlet, then sediments are deposited
-in the lake/reservoir according to Camp’s model (1945) (Verstraeten et al, 2000):
+in the lake/reservoir according to Camp's model (1945) (Verstraeten et al, 2000):
 ```math
    TE = \dfrac{\omega_{s}}{u_{cr,res}} = \dfrac{A_{res}}{Q_{out,res}} \omega_{s}
 ```
 where ``TE`` is the trapping efficiency of the lake/reservoir (or the fraction of particles
 trapped), ``\omega_{s}`` is the particle velocity from Stokes (m s``^{-1}``), ``u_{cr,res}``
-is the reservoir’s critical settling velocity (m/s) which is equal to the reservoir’s
-outflow ``Q_{out,res}`` (m``^{3}`` s``^{-1}``) divided by the reservoir’s surface area
+is the reservoir's critical settling velocity (m/s) which is equal to the reservoir's
+outflow ``Q_{out,res}`` (m``^{3}`` s``^{-1}``) divided by the reservoir's surface area
 ``A_{res}`` (m``^{2}``).
 
 For reservoirs, coarse sediment particles from the bed load are also assumed to be trapped by the
@@ -439,4 +439,4 @@ transport in overland flow.
   the Total Environment, 538:855-875, 2015. 10.1016/j.scitotenv.2015.08.095
 + O. Vigiak, A. Malago, F. Bouraoui, M. Vanmaercke, F. Obreja, J. Poesen, H. Habersack, J.
   Feher, and S. Groselj. Modelling sediment fluxes in the Danube River Basin with SWAT.
-  Science of the Total Environment, 2017. 10.1016/j.scitotenv.2017.04.236
\ No newline at end of file
+  Science of the Total Environment, 2017. 10.1016/j.scitotenv.2017.04.236
diff --git a/docs/src/model_docs/params_lateral.md b/docs/src/model_docs/params_lateral.md
index d0df95f21..05d89cea3 100644
--- a/docs/src/model_docs/params_lateral.md
+++ b/docs/src/model_docs/params_lateral.md
@@ -13,7 +13,7 @@ internal model parameter `sl`, and is listed in the Table below between parenthe
 
 |  parameter  | description  	  | unit  | default |
 |:--------------- | ------------------| ----- | -------- |
-| `β`             |  constant in Manning's equation | - | - |
+| `beta`             |  constant in Manning's equation | - | - |
 | **`sl`** (`slope`) |  slope       | m m``^{-1}``| - |
 | **`n`**         | Manning's roughness | s m``^{-\frac{1}{3}}``| 0.036 |
 | **`dl`**        |  length      | m     |  -      |
@@ -27,12 +27,12 @@ internal model parameter `sl`, and is listed in the Table below between parenthe
 | `h`             | water level | m | - |
 | `h_av`          | average water level | m | - |
 | **`bankfull_depth`**   | bankfull river depth  | m | 1.0 |
-| `Δt`            | model time step | s | - |
+| `dt`            | model time step | s | - |
 | `its`           | number of fixed iterations | - | - |
 | **`width`**     |  width       | m          | - |
 | `alpha_pow`     | used in the power part of ``\alpha`` | - | - |
 | `alpha_term`    | term used in computation of ``\alpha`` | - | - |
-| `α`             | constant in momentum equation ``A = \alpha Q^{\beta}`` | s``^{\frac{3}{5}}`` m``^{\frac{1}{5}}`` | - |
+| `alpha`             | constant in momentum equation ``A = \alpha Q^{\beta}`` | s``^{\frac{3}{5}}`` m``^{\frac{1}{5}}`` | - |
 | `cel`           | celerity of kinematic wave | m s``^{-1}`` | - |
 | `reservoir_index`   |  map cell to 0 (no reservoir) or i (pick reservoir i in reservoir field) | - | - |
 | `lake_index`   |  map cell to 0 (no lake) or i (pick lake i in lake field) | - | - |
@@ -50,7 +50,7 @@ internal model parameter `sl`, and is listed in the Table below between parenthe
 
 |  parameter  | description  	  | unit  | default |
 |:--------------- | ------------------| ----- | -------- |
-| `β`             |  constant in Manning's equation | - | - |
+| `beta`             |  constant in Manning's equation | - | - |
 | **`sl`** (`slope`) |  slope       | m m``^{-1}``| - |
 | **`n`**         | Manning's roughness | s m``^{-\frac{1}{3}}``| 0.072 |
 | `dl`            |  length      | m     |  -      |
@@ -62,12 +62,12 @@ internal model parameter `sl`, and is listed in the Table below between parenthe
 | `volume`        | kinematic wave volume |m``^3``| - |
 | `h`             | water level | m | - |
 | `h_av`          | average water level | m | - |
-| `Δt`            | model time step | s | - |
+| `dt`            | model time step | s | - |
 | `its`           | number of fixed iterations | - | - |
 | `width`         |  width       | m          | - |
 | `alpha_pow`     | used in the power part of ``\alpha`` | - | - |
 | `alpha_term`    | term used in computation of ``\alpha`` | - | - |
-| `α`             | constant in momentum equation ``A = \alpha Q^{\beta}`` | s``^{\frac{3}{5}}`` m``^{\frac{1}{5}}`` | - |
+| `alpha`             | constant in momentum equation ``A = \alpha Q^{\beta}`` | s``^{\frac{3}{5}}`` m``^{\frac{1}{5}}`` | - |
 | `cel`           | celerity of kinematic wave | m s``^{-1}`` | - |
 | `to_river`      | part of overland flow that flows to the river | m``^3`` s``^{-1}`` | - |
 | `kinwave_it`    | boolean for kinematic wave iterations | - | false |
@@ -98,7 +98,7 @@ locs = "wflow_reservoirlocs"
 | **`targetfullfrac`** | target fraction full (of max storage)| -        | - |
 | **`targetminfrac`** | target minimum full fraction (of max storage) | -  | - |
 | `demandrelease`| minimum (environmental) flow released from reservoir  | m``^3`` s``^{-1}``| - |
-| `Δt`             | model time step     | s | - |
+| `dt`             | model time step     | s | - |
 | `volume`             | volume     | m``^3`` | - |
 | `inflow`             | total inflow into reservoir | m``^3`` | - |
 | `outflow`             | outflow into reservoir | m``^3`` s``^{-1}`` | - |
@@ -141,7 +141,7 @@ between parentheses.
 | **`lowerlake_ind`** (`linkedlakelocs`) | Index of lower lake (linked lakes) | - | 0 |
 | **`sh`** | data for storage curve | - | - |
 | **`hq`** | data rating curve | - | - |
-| `Δt`             | model time step     | s |  - |
+| `dt`             | model time step     | s |  - |
 | `inflow` | total inflow to the lake | m``^3``  | - |
 | `storage` | storage lake | m``^3``  | - |
 | `maxstorage`| maximum storage lake with rating curve type 1 | m``^3``  | - |
@@ -155,10 +155,10 @@ between parentheses.
 The Table below shows the parameters (fields) of struct `LateralSSF`, including a
 description of these parameters, the unit, and default value if applicable. The parameters
 in bold represent model parameters that can be set through static input data (netCDF). The
-soil related parameters `f`, `soilthickness`, `z_exp`, `θₛ` and `θᵣ` are derived from the
+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 `βₗ` is
+model parameter to the external netCDF variable. The internal slope model parameter `slope` is
 set through the TOML file as follows:
 
 ```toml
@@ -166,8 +166,8 @@ set through the TOML file as follows:
 slope = "Slope"
 ```
 
-The parameter `kh₀` is computed by multiplying the vertical hydraulic conductivity at the
-soil surface `kv₀` (including unit conversion) of the vertical `SBM` concept with the
+The parameter `kh_0` is computed by multiplying the vertical hydraulic conductivity at the
+soil surface `kv_0` (including unit conversion) of the vertical `SBM` concept with the
 internal parameter `khfrac` \[-\] (default value of 1.0). The internal model parameter
 `khfrac` is set through the TOML file as follows:
 
@@ -176,34 +176,34 @@ internal parameter `khfrac` \[-\] (default value of 1.0). The internal model par
 ksathorfrac = "KsatHorFrac"
 ```
 
-The `khfrac` parameter compensates for anisotropy, small scale `kv₀` measurements (soil
+The `khfrac` parameter compensates for anisotropy, small scale `kv_0` measurements (soil
 core) that do not represent larger scale hydraulic conductivity, and smaller flow length
 scales (hillslope) in reality, not represented by the model resolution.
 
 For the vertical [SBM](@ref params_sbm) concept different vertical hydraulic conductivity
 depth profiles are possible, and these also determine which `LateralSSF` parameters are used
 including the input requirements for the computation of lateral subsurface flow. For the
-`exponential` profile the model parameters `kh₀` and `f` are used. For the
-`exponential_constant` profile `kh₀` and `f` are used, and `z_exp` is required as part of
+`exponential` profile the model parameters `kh_0` and `f` are used. For the
+`exponential_constant` profile `kh_0` and `f` are used, and `z_exp` is required as part of
 `[input.vertical]`. For the `layered` profile, `SBM` model parameter `kv` is used, and for
 the `layered_exponential` profile `kv` is used and `z_exp` is required as part of
 `[input.vertical]`.
 
 | parameter | description  	        | unit | default |
 |:---------------| --------------- | ---------------------- | ----- |
-| `kh₀`          | horizontal hydraulic conductivity at soil surface  | m d``^{-1}`` |  3.0 |
-| **`f`** | a scaling parameter (controls exponential decline of `kh₀`) | m``^{-1}``  | 1.0 |
+| `kh_0`          | horizontal hydraulic conductivity at soil surface  | m d``^{-1}`` |  3.0 |
+| **`f`** | a scaling parameter (controls exponential decline of `kh_0`) | m``^{-1}``  | 1.0 |
 | `kh` | horizontal hydraulic conductivity | m d``^{-1}``  | - |
 | **`khfrac`** (`ksathorfrac`) | a muliplication factor applied to vertical hydraulic conductivity `kv` | -  | 100.0 |
 | **`soilthickness`** | soil thickness | m  | 2.0 |
-| **`θₛ`** | saturated water content (porosity) | -  | 0.6 |
-| **`θᵣ`** | residual water content  | -  | 0.01 |
-| `Δt` | model time step | d  | - |
-| **`βₗ`** (`slope`) | slope | m m``^{-1}``  | - |
+| **`theta_s`** | saturated water content (porosity) | -  | 0.6 |
+| **`theta_r`** | residual water content  | -  | 0.01 |
+| `dt` | model time step | d  | - |
+| **`slope`** | slope | m m``^{-1}``  | - |
 | `dl` | drain length | m | - |
 | `dw` | drain width | m | - |
 | `zi` | pseudo-water table depth (top of the saturated zone) | m | - |
-| **`z_exp`** | depth from soil surface for which exponential decline of `kh₀` is valid | m | - |
+| **`z_exp`** | depth from soil surface for which exponential decline of `kh_0` is valid | m | - |
 | `exfiltwater` | exfiltration (groundwater above surface level, saturated excess conditions) | m Δt⁻¹ | - |
 | `recharge` | net recharge to saturated store  | m``^2`` Δt⁻¹ | - |
 | `ssf` | subsurface flow  | m``^3`` d``{-1}``  | - |
@@ -247,18 +247,18 @@ model parameter `mannings_n`, and is listed in the Table below between parenthes
 | `active_n`    |  active nodes | - | - |
 | `active_e`    |  active edges | - | - |
 | `g`    |  acceleration due to gravity | m s``^{-2}`` | - |
-| `α`    |  stability coefficient (Bates et al., 2010) | - | 0.7 |
+| `alpha`    |  stability coefficient (Bates et al., 2010) | - | 0.7 |
 | `h_thresh`    |  depth threshold for calculating flow | m | 0.001 |
-| `Δt`    |  model time step | s | - |
+| `dt`    |  model time step | s | - |
 | `q`    |  river discharge (subgrid channel) | m``^3`` s``^{-1}`` | - |
 | `q_av`    |  average river channel (+ floodplain) discharge | m``^3`` s``^{-1}`` | - |
 | `q_channel_av` | average river channel discharge | m``^3`` s``^{-1}`` | - |
 | `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 | - |
-| `η_max`    | maximum water elevation | m | - |
-| `η_src`    | water elevation of source node of edge | m | - |
-| `η_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 | - |
@@ -356,10 +356,10 @@ internal model parameter `z`, and is listed in the Table below between parenthes
 | `xwidth`| effective flow width x direction (floodplain) | m | - |
 | `ywidth`| effective flow width y direction (floodplain) | m | - |
 | `g` | acceleration due to gravity | m s``^{-2}`` | - |
-| `θ` | weighting factor (de Almeida et al., 2012) | - | 0.8 |
-| `α`    |  stability coefficient (Bates et al., 2010) | - | 0.7 |
+| `theta` | weighting factor (de Almeida et al., 2012) | - | 0.8 |
+| `alpha`    |  stability coefficient (Bates et al., 2010) | - | 0.7 |
 | `h_thresh`    |  depth threshold for calculating flow | m | 0.001 |
-| `Δt`   |  model time step| s | - |
+| `dt`   |  model time step| s | - |
 | `qy0`  |  flow in y direction at previous time step| m``^3`` s``^{-1}`` | - |
 | `qx0`  |  flow in x direction at previous time step| m``^3`` s``^{-1}`` | - |
 | `qx`  |  flow in x direction | m``^3`` s``^{-1}`` | - |
@@ -412,11 +412,11 @@ altitude = "wflow_dem"
 ```
 
 The input parameter `conductivity` (listed under `[input.lateral.subsurface]`) is not equal
-to the internal model parameter `kh₀`, and is listed in the Table below between parentheses.
+to the internal model parameter `kh_0`, and is listed in the Table below between parentheses.
 
 |  parameter  | description  	        | unit  | default |
 |:--------------- | ------------------| ----- | -------|
-| **`kh₀`** (`conductivity`) | horizontal conductivity  | m d``^{-1}``s | - |
+| **`kh_0`** (`conductivity`) | horizontal conductivity  | m d``^{-1}``s | - |
 | **`specific_yield`**    | specific yield  | m m``^{-1}`` | - |
 | **`top`** (`altitude`)  | top groundwater layer  | m | - |
 | `bottom`     | bottom groundwater layer  | m | - |
@@ -592,7 +592,7 @@ slope = "RiverSlope"
 | **`cbagnold`**| Bagnold c coefficient  | - | - |
 | **`ebagnold`**| Bagnold exponent | - | - |
 | `n`             | number of cells     | - |  - |
-| `Δt`             | model time step     | s |  - |
+| `dt`             | model time step     | s |  - |
 | `ak`             | Kodatie coefficient `a`    | - |  - |
 | `bk`             | Kodatie coefficient `b`    | - |  - |
 | `ck`             | Kodatie coefficient `c`    | - |  - |
diff --git a/docs/src/model_docs/params_vertical.md b/docs/src/model_docs/params_vertical.md
index 840eee979..02477283c 100644
--- a/docs/src/model_docs/params_vertical.md
+++ b/docs/src/model_docs/params_vertical.md
@@ -25,8 +25,8 @@ through the TOML file. Below an example for the `exponential_constant` profile:
 ksat_profile = "exponential_constant"
 ```
 
-For the `exponential` profile the input parameters `kv₀` and `f` are used. For the
-`exponential_constant` profile `kv₀` and `f` are used, and `z_exp` is required as input. For
+For the `exponential` profile the input parameters `kv_0` and `f` are used. For the
+`exponential_constant` profile `kv_0` and `f` are used, and `z_exp` is required as input. For
 the `layered` profile, input parameter `kv` is used, and for the `layered_exponential`
 profile `kv` is used and `z_layered` is required as input.
 
@@ -44,12 +44,12 @@ profile `kv` is used and `z_layered` is required as input.
 | **`g_sifrac`**  | fraction of the snowpack on top of the glacier converted into ice  | Δt``^{-1}`` | 0.001 day``^{-1}`` |
 | **`glacierfrac`**  | fraction covered by a glacier | - | 0.0  |
 | **`glacierstore`**  | water within the glacier  | mm | 5500.0  |
-| **`θₛ`** (`theta_s`) | saturated water content (porosity) | - | 0.6  |
-| **`θᵣ`** (`theta_r`) | residual water content | - | 0.01  |
-| **`kv₀`** (`kv_0`) | Vertical hydraulic conductivity at soil surface | mm Δt``^{-1}`` | 3000.0 mm day``^{-1}``|
+| **`theta_s`** | saturated water content (porosity) | - | 0.6  |
+| **`theta_r`** | residual water content | - | 0.01  |
+| **`kv_0`** | Vertical hydraulic conductivity at soil surface | mm Δt``^{-1}`` | 3000.0 mm day``^{-1}``|
 | **`kv`** |  Vertical hydraulic conductivity per soil layer | mm Δt``^{-1}`` | 1000.0 mm day``^{-1}``|
-| **`f`** | scaling parameter (controls exponential decline of `kv₀`) | mm``^{-1}`` | 0.001  |
-| **`z_exp`** | Depth from soil surface for which exponential decline of `kv₀` is valid | mm | -  |
+| **`f`** | scaling parameter (controls exponential decline of `kv_0`) | mm``^{-1}`` | 0.001  |
+| **`z_exp`** | Depth from soil surface for which exponential decline of `kv_0` is valid | mm | -  |
 | **`z_layered`** | Depth from soil surface for which layered profile (of `layered_exponential`) is valid | mm | -  |
 | **`hb`** | air entry pressure of soil (Brooks-Corey) | cm | 10.0  |
 | **`soilthickness`** | soil thickness | mm | 2000.0  |
@@ -71,7 +71,7 @@ profile `kv` is used and `z_layered` is required as input.
 | **`cmax`** | maximum canopy storage | mm | 1.0 |
 | **`e_r`** (`eoverr`) | Gash interception model parameter | - | 0.1 |
 | **`canopygapfraction`** | canopy gap fraction | - | 0.1 | - |
-| `Δt`             | model time step     | s | - |
+| `dt`             | model time step     | s | - |
 | `maxlayers`      | maximum number of soil layers     | - | - |
 | `n`      |  number of grid cells    | - | - |
 | `nlayers`      |  number of soil layers    | - | - |
@@ -117,11 +117,11 @@ profile `kv` is used and `z_layered` is required as input.
 | `excesswaterpath` | excess water for compacted fraction | mm Δt``^{-1}`` | - |
 | `runoff` | total surface runoff from infiltration and saturation excess | mm Δt``^{-1}`` | - |
 | `net_runoff` | net surface runoff (`runoff` - `ae_openw_l`) | mm Δt``^{-1}`` | - |
-| `vwc` | volumetric water content per soil layer (including θᵣ and saturated zone)  | - | - |
-| `vwc_perc` | volumetric water content per soil layer (including θᵣ and saturated zone)  | % | - |
-| `rootstore` | root water storage in unsaturated and saturated zone (excluding θᵣ)  | mm| - |
-| `vwc_root` | volumetric water content in root zone (including θᵣ and saturated zone) | -| - |
-| `vwc_percroot` | volumetric water content in root zone (including θᵣ and saturated zone) | % | - |
+| `vwc` | volumetric water content per soil layer (including `theta_r` and saturated zone)  | - | - |
+| `vwc_perc` | volumetric water content per soil layer (including `theta_r` and saturated zone)  | % | - |
+| `rootstore` | root water storage in unsaturated and saturated zone (excluding `theta_r`)  | mm| - |
+| `vwc_root` | volumetric water content in root zone (including `theta_r` and saturated zone) | -| - |
+| `vwc_percroot` | volumetric water content in root zone (including `theta_r` and saturated zone) | % | - |
 | `ustoredepth` | total amount of available water in the unsaturated zone | mm | - |
 | `transfer` | downward flux from unsaturated to saturated zone | mm Δt``^{-1}`` | - |
 | `recharge` | net recharge to saturated zone | mm Δt``^{-1}`` | - |
@@ -377,4 +377,4 @@ specific_leaf = "Sl"
 | `TCsilt` | transport capacity of overland flow for particle class silt  | ton Δt``^{-1}`` | - |
 | `TCsand` | transport capacity of overland flow for particle class sand  | ton Δt``^{-1}`` | - |
 | `TCsagg` | transport capacity of overland flow for particle class small aggregates  | ton Δt``^{-1}`` | - |
-| `TClagg` | transport capacity of overland flow for particle class large aggregates  | ton Δt``^{-1}`` | - |
\ No newline at end of file
+| `TClagg` | transport capacity of overland flow for particle class large aggregates  | ton Δt``^{-1}`` | - |
diff --git a/docs/src/model_docs/shared_concepts.md b/docs/src/model_docs/shared_concepts.md
index c8edbde6e..ebda6144b 100644
--- a/docs/src/model_docs/shared_concepts.md
+++ b/docs/src/model_docs/shared_concepts.md
@@ -88,7 +88,7 @@ Both the Gash and Rutter models are available to estimate rainfall interception
 vegetation. The selection of an interception model depends on the simulation timestep.
 
 ### The analytical (Gash) model
-The analytical model of rainfall interception is based on Rutter’s numerical model. The
+The analytical model of rainfall interception is based on Rutter's numerical model. The
 simplifications that introduced allow the model to be applied on a daily basis, although a
 storm-based approach will yield better results in situations with more than one storm per
 day. The amount of water needed to completely saturate the canopy is defined as:
@@ -182,4 +182,4 @@ The extinction coefficient `kext` can be related to land cover.
 ## References
 + Seibert, J., Vis, M. J. P., Kohn, I., Weiler, M., and Stahl, K., 2018, Technical note:
   Representing glacier geometry changes in a semi-distributed hydrological model, Hydrol.
-  Earth Syst. Sci., 22, 2211–2224, https://doi.org/10.5194/hess-22-2211-2018.
\ No newline at end of file
+  Earth Syst. Sci., 22, 2211–2224, https://doi.org/10.5194/hess-22-2211-2018.
diff --git a/docs/src/model_docs/structures.md b/docs/src/model_docs/structures.md
index f325c3d40..78ca20a85 100644
--- a/docs/src/model_docs/structures.md
+++ b/docs/src/model_docs/structures.md
@@ -54,7 +54,7 @@ example with a part of the `SBM` struct:
 ```julia
 @get_units @with_kw struct SBM{T,N,M}
     # Model time step [s]
-    Δt::T | "s"
+    dt::T | "s"
     # Maximum number of soil layers
     maxlayers::Int | "-"
     # number of cells
@@ -66,11 +66,11 @@ example with a part of the `SBM` struct:
     # Fraction of river [-]
     riverfrac::Vector{T} | "-"
     # Saturated water content (porosity) [mm mm⁻¹]
-    θₛ::Vector{T} | "mm mm-1"
+    theta_s::Vector{T} | "mm mm-1"
     # Residual water content [mm mm⁻¹]
-    θᵣ::Vector{T} | "mm mm-1"
+    theta_r::Vector{T} | "mm mm-1"
     # Vertical hydraulic conductivity [mm Δt⁻¹] at soil surface
-    kv₀::Vector{T} | "mm Δt-1"
+    kv_0::Vector{T} | "mm dt-1"
     # Muliplication factor [-] applied to kv_z (vertical flow)
     kvfrac::Vector{SVector{N,T}} | "-"
 ```
@@ -83,9 +83,9 @@ refers to the maximum number of soil layers + 1. See also part of the following
 
 ```julia
 sbm = SBM{Float,maxlayers,maxlayers + 1}(
-    Δt = tosecond(Δt),
+    dt = tosecond(dt),
     maxlayers = maxlayers,
     n = n,
 ```
 
-For the other model concepts, we refer to the code to check these type parameters.
\ No newline at end of file
+For the other model concepts, we refer to the code to check these type parameters.
diff --git a/docs/src/model_docs/vertical/flextopo.md b/docs/src/model_docs/vertical/flextopo.md
index 2ea678769..114b54920 100644
--- a/docs/src/model_docs/vertical/flextopo.md
+++ b/docs/src/model_docs/vertical/flextopo.md
@@ -129,10 +129,10 @@ The shape of the beta function for various values of ``\beta`` [-] is shown belo
         ax = Axis(fig[1, 1], xlabel = "S/Smax [-]", ylabel = "Fraction of runoff [-]")      # hide
         x = 0:0.01:1                                                                        # hide
         betas = [0.3, 1.0, 3.0]                                                             # hide
-        for β in betas                                                                      # hide
-            lines!(ax, x, (1 .- (1 .- x).^β), label = @sprintf("β = %.1f", β))              # hide
+        for beta in betas                                                                   # hide
+            lines!(ax, x, (1 .- (1 .- x).^beta), label = @sprintf("beta = %.1f", beta))     # hide
         end                                                                                 # hide
-        Legend(fig[1, 2], ax, "β")                                                          # hide
+        Legend(fig[1, 2], ax, "beta")                                                       # hide
         fig                                                                                 # hide
     end                                                                                     # hide
 ```
diff --git a/docs/src/model_docs/vertical/sbm.md b/docs/src/model_docs/vertical/sbm.md
index 6a7744339..ca6146f76 100644
--- a/docs/src/model_docs/vertical/sbm.md
+++ b/docs/src/model_docs/vertical/sbm.md
@@ -37,7 +37,7 @@ Rutter model. The simulation timestep defines which interception model is used,
 Rutter model.
 
 ### The analytical (Gash) model (Gash, 1979)
-The analytical model of rainfall interception is based on Rutter’s numerical model. Simplifications 
+The analytical model of rainfall interception is based on Rutter's numerical model. Simplifications 
 allow the model to be applied on a daily basis, although a
 storm-based approach will yield better results in situations with more than one storm per
 day. The amount of water needed to completely saturate the canopy is defined as:
@@ -549,7 +549,7 @@ the vector that contains all active cells within the spatial model domain, and `
 the layer position):
 
 ```julia
-    usl[k] * (sbm.θₛ[i] - sbm.θᵣ[i]) - usld[k]
+    usl[k] * (sbm.theta_s[i] - sbm.theta_r[i]) - usld[k]
 ```
 
 where `usl` [mm] is the unsaturated layer thickness, `usld` is the `ustorelayerdepth` \[mm\]
@@ -565,7 +565,7 @@ spatial model domain, and `k` refers to the layer position):
     netcapflux = capflux
     for k = n_usl:-1:1
         toadd =
-            min(netcapflux, max(usl[k] * (sbm.θₛ[i] - sbm.θᵣ[i]) - usld[k], 0.0))
+            min(netcapflux, max(usl[k] * (sbm.theta_s[i] - sbm.theta_r[i]) - usld[k], 0.0))
         usld = setindex(usld, usld[k] + toadd, k)
         netcapflux = netcapflux - toadd
         actcapflux = actcapflux + toadd
diff --git a/docs/src/model_docs/vertical/sediment.md b/docs/src/model_docs/vertical/sediment.md
index 34d672b8e..c18981feb 100644
--- a/docs/src/model_docs/vertical/sediment.md
+++ b/docs/src/model_docs/vertical/sediment.md
@@ -97,12 +97,12 @@ is the soil cover-management factor from the USLE equation, ``K_{USLE}`` is the
 erodibility factor from the USLE equation, ``A_{i}`` is the area of the cell (m``^{2}``) and
 ``R_{i}`` is the rainfall intensity (here in mm min``^{-1}``). There are several methods
 available to estimate the ``C`` and ``K`` factors from the USLE. They can come from user
-input maps, for example maps resulting from Panagos & al.’s recent studies for Europe
+input maps, for example maps resulting from Panagos & al.'s recent studies for Europe
 (Panagos et al, 2015) (Ballabio et al, 2016). To get an estimate of the ``C`` factor
 globally, the other method is to estimate ``C`` values for the different land use type in
 from global land cover maps (e.g. GlobCover). An example is given for the global land cover
 map GlobCover, summed up in the table below, the values come from a literature study
-including Panagos et al.’s review (2015), Gericke & al. (2015), Mansoor & al. (2013), Chadli
+including Panagos et al.'s review (2015), Gericke & al. (2015), Mansoor & al. (2013), Chadli
 et al. (2016), de Vente et al. (2009), Borrelli et al. (2014), Yang et al. (2003) and Bosco
 et al. (2015).
 
@@ -220,4 +220,4 @@ This process is described in [Sediment Flux in overland flow](@ref).
   effect on soil productivity. Journal of Soil and Water Conservation, 38(5):381-383, 1983.
 + D. Yang, S. Kanae, T. Oki, T. Koike, and K. Musiake. Global potential soil erosion with
   reference to land use and climate changes. Hydrological Processes, 17(14):2913-2928, 2003.
-  10.1002/hyp.1441
\ No newline at end of file
+  10.1002/hyp.1441
diff --git a/docs/src/user_guide/model-setup.md b/docs/src/user_guide/model-setup.md
index 65b952463..2affe2e69 100644
--- a/docs/src/user_guide/model-setup.md
+++ b/docs/src/user_guide/model-setup.md
@@ -57,7 +57,7 @@ routing](@ref) and parameters that are part of this component are described in t
 subsurface flow](@ref) section of Model parameters. Input parameters for this component of
 the SBM + Kinematic wave model are derived from the SBM vertical concept and the land slope.
 One external parameter [`ksathorfrac`](@ref params_ssf) is used to calculate the horizontal
-hydraulic conductivity at the soil surface `kh₀`.
+hydraulic conductivity at the soil surface `kh_0`.
 
 There is also the option to use the local inertial model as part of the `sbm` model type:
 + for river flow, see also  [SBM + Local inertial river](@ref) model.
@@ -131,4 +131,4 @@ map stack.
 + Eilander, D., van Verseveld, W., Yamazaki, D., Weerts, A., Winsemius, H. C., and Ward, P.
   J.: A hydrography upscaling method for scale-invariant parametrization of distributed
   hydrological models, Hydrol. Earth Syst. Sci., 25, 5287–5313,
-  <https://doi.org/10.5194/hess-25-5287-2021>, 2021.
\ No newline at end of file
+  <https://doi.org/10.5194/hess-25-5287-2021>, 2021.
diff --git a/docs/src/user_guide/step2_settings_file.md b/docs/src/user_guide/step2_settings_file.md
index 021d5a50a..490679524 100644
--- a/docs/src/user_guide/step2_settings_file.md
+++ b/docs/src/user_guide/step2_settings_file.md
@@ -152,7 +152,7 @@ e_r = "EoverR"
 infiltcappath = "InfiltCapPath"
 infiltcapsoil = "InfiltCapSoil"
 kext = "Kext"
-"kv₀" = "KsatVer"
+"kv_0" = "KsatVer"
 leaf_area_index = "LAI"             # Cyclic variable
 m = "M"
 maxleakage = "MaxLeakage"
@@ -172,8 +172,8 @@ ttm = "TTM"
 w_soil = "wflow_soil"
 water_holding_capacity = "WHC"
 waterfrac = "WaterFrac"
-"θᵣ" = "thetaR"
-"θₛ" = "thetaS"
+"theta_r" = "thetaR"
+"theta_s" = "thetaS"
 
 [input.lateral.river]
 length = "wflow_riverlength"
@@ -430,4 +430,4 @@ Note that the mapping to the external netCDF variable listed under the section
 potential_evaporation = "PET" # forcing
 # temperature = "TEMP" # forcing
 precipitation = "P" # forcing
-```
\ No newline at end of file
+```
diff --git a/docs/src/user_guide/step4_running.md b/docs/src/user_guide/step4_running.md
index 5294ca805..4f162647b 100644
--- a/docs/src/user_guide/step4_running.md
+++ b/docs/src/user_guide/step4_running.md
@@ -61,7 +61,7 @@ bar, that gives an estimate of how much time is needed to finish the simulation:
 ┌ Info: Run information
 │   model_type = "sbm"
 │   starttime = CFTime.DateTimeStandard(2000-01-01T00:00:00)
-│   Δt = 86400 seconds
+│   dt = 86400 seconds
 │   endtime = CFTime.DateTimeStandard(2000-12-31T00:00:00)
 └   nthreads() = 4
 
diff --git a/src/Wflow.jl b/src/Wflow.jl
index c139fb371..15784641d 100644
--- a/src/Wflow.jl
+++ b/src/Wflow.jl
@@ -22,7 +22,7 @@ using Polyester
 using LoopVectorization
 using IfElse
 
-@metadata get_units "mm Δt-1" String
+@metadata get_units "mm dt-1" String
 # metadata for BMI grid
 @metadata exchange 1 Integer
 @metadata grid_type "unstructured" String
@@ -38,7 +38,7 @@ const version =
 mutable struct Clock{T}
     time::T
     iteration::Int
-    Δt::Second
+    dt::Second
 end
 
 function Clock(config)
@@ -46,8 +46,8 @@ function Clock(config)
     # been constructed before, the config is complete
     calendar = get(config, "calendar", "standard")::String
     starttime = cftime(config.starttime, calendar)
-    Δt = Second(config.timestepsecs)
-    Clock(starttime, 0, Δt)
+    dt = Second(config.timestepsecs)
+    Clock(starttime, 0, dt)
 end
 
 function Clock(config, reader)
@@ -59,13 +59,13 @@ function Clock(config, reader)
         timestepsecs = Dates.value(Second(nctimes[2] - nctimes[1]))
         config.timestepsecs = timestepsecs
     end
-    Δt = Second(timestepsecs)
+    dt = Second(timestepsecs)
 
     # if the config file does not have a start or endtime, follow the netCDF times
     # and add them to the config
     starttime = get(config, "starttime", nothing)
     if starttime === nothing
-        starttime = first(nctimes) - Δt
+        starttime = first(nctimes) - dt
         config.starttime = starttime
     end
     endtime = get(config, "endtime", nothing)
@@ -77,11 +77,11 @@ function Clock(config, reader)
     calendar = get(config, "calendar", "standard")::String
     fews_run = get(config, "fews_run", false)::Bool
     if fews_run
-        config.starttime = starttime + Δt
+        config.starttime = starttime + dt
     end
     starttime = cftime(config.starttime, calendar)
 
-    Clock(starttime, 0, Δt)
+    Clock(starttime, 0, dt)
 end
 
 include("io.jl")
@@ -220,15 +220,15 @@ function run(model::Model; close_files = true)
     calendar = get(config, "calendar", "standard")::String
     @warn string(
         "The definition of `starttime` has changed (equal to model state time).\n Please",
-        " update your settings TOML file by subtracting one model timestep Δt from the",
+        " update your settings TOML file by subtracting one model timestep dt from the",
         " `starttime`, if it was used with a Wflow version up to v0.6.3.",
     )
     starttime = clock.time
-    Δt = clock.Δt
+    dt = clock.dt
     endtime = cftime(config.endtime, calendar)
-    times = range(starttime + Δt, endtime, step = Δt)
+    times = range(starttime + dt, endtime, step = dt)
 
-    @info "Run information" model_type starttime Δt endtime nthreads()
+    @info "Run information" model_type starttime dt endtime nthreads()
     runstart_time = now()
     @progress for (i, time) in enumerate(times)
         @debug "Starting timestep." time i now()
diff --git a/src/bmi.jl b/src/bmi.jl
index def51be69..a2a7b9899 100644
--- a/src/bmi.jl
+++ b/src/bmi.jl
@@ -66,14 +66,14 @@ end
 function BMI.update_until(model::Model, time::Float64)
     @unpack clock, network, config = model
     t = BMI.get_current_time(model)
-    _div, _rem = divrem(time - t, model.clock.Δt.value)
+    _div, _rem = divrem(time - t, model.clock.dt.value)
     steps = Int(_div)
     if steps < 0
         error("The current model timestamp $t is larger than provided `time` $time")
     elseif abs(_rem) > eps()
         error_message = string(
             "Provided `time` $time minus the current model timestamp $t",
-            " is not an integer multiple of model time step $(model.clock.Δt.value)",
+            " is not an integer multiple of model time step $(model.clock.dt.value)",
         )
         error(error_message)
     end
@@ -86,7 +86,7 @@ end
 "Write state output to netCDF and close files."
 function BMI.finalize(model::Model)
     @unpack config, writer, clock = model
-    # it is possible that the state dataset has been closed by `save_state`  
+    # it is possible that the state dataset has been closed by `save_state`
     if !isnothing(writer.state_dataset) && isopen(writer.state_dataset)
         write_netcdf_timestep(model, writer.state_dataset, writer.state_parameters)
     end
@@ -282,7 +282,7 @@ end
     BMI.set_value(model::Model, name::String, src::Vector{T}) where T<:AbstractFloat
 
 Set a model variable `name` to the values in vector `src`, overwriting the current contents.
-The type and size of `src` must match the model’s internal array.
+The type and size of `src` must match the model's internal array.
 """
 function BMI.set_value(model::Model, name::String, src::Vector{T}) where {T<:AbstractFloat}
     BMI.get_value_ptr(model, name) .= src
diff --git a/src/flextopo.jl b/src/flextopo.jl
index ea2db00d4..d547a34b9 100644
--- a/src/flextopo.jl
+++ b/src/flextopo.jl
@@ -1,6 +1,6 @@
 @get_units @exchange @grid_type @grid_location @with_kw struct FLEXTOPO{T,N}
     # Model time step [s]
-    Δt::T | "s" | 0 | "none" | "none"
+    dt::T | "s" | 0 | "none" | "none"
     # Number of classes
     nclass::Int | "-" | 0 | "none" | "none"
     # Number of cells
@@ -25,7 +25,7 @@
     # Correction factor for precipitation [-]
     pcorr::Vector{T} | "-"
     # Degree-day factor [mm ᵒC⁻¹ Δt⁻¹]
-    cfmax::Vector{T} | "mm ᵒC-1 Δt-1"
+    cfmax::Vector{T} | "mm ᵒC-1 dt-1"
     # Threshold temperature for snowfall [ᵒC]
     tt::Vector{T} | "ᵒC"
     # Threshold temperature interval length [ᵒC]
@@ -44,9 +44,9 @@
     # Threshold temperature for snowfall above glacier [ᵒC]
     g_tt::Vector{T} | "ᵒC"
     # Degree-day factor [mm ᵒC⁻¹ Δt⁻¹] for glacier
-    g_cfmax::Vector{T} | "mm ᵒC-1 Δt-1"
+    g_cfmax::Vector{T} | "mm ᵒC-1 dt-1"
     # Fraction of the snowpack on top of the glacier converted into ice [Δt⁻¹]
-    g_sifrac::Vector{T} | "Δt-1"
+    g_sifrac::Vector{T} | "dt-1"
     # Water within the glacier [mm]
     glacierstore::Vector{T} | "mm"
     # Fraction covered by a glacier [-]
@@ -59,16 +59,16 @@
     ##HORTON
     # Maximum storage capacity in the hortonian ponding storage [mm]
     shmax::Vector{SVector{N,T}} | "mm"
-    #recession coefficient of the hortonian runoff storage [Δt-1]
-    khf::Vector{SVector{N,T}} | "Δt-1"
+    #recession coefficient of the hortonian runoff storage [dt-1]
+    khf::Vector{SVector{N,T}} | "dt-1"
     #maximum modelled accumulated frost resulting in shmin [ᵒC Δt]
     facc0::Vector{SVector{N,T}} | "ᵒC"
     #minimum modelled accumulated frost resulting in shmax [ᵒC Δt]
     facc1::Vector{SVector{N,T}} | "ᵒC"
     #exponent for the decline of infiltration capacity [-]
     fdec::Vector{SVector{N,T}} | "-"
-    #maximum infiltration capacity from horton ponding [mm Δt-1]
-    fmax::Vector{SVector{N,T}} | "mm Δt-1"
+    #maximum infiltration capacity from horton ponding [mm dt-1]
+    fmax::Vector{SVector{N,T}} | "mm dt-1"
     #minimum storage capacity in horton ponding (relative to shmax) [-]
     shmin::Vector{SVector{N,T}} | "-"
     #melt coefficient for melt of frozen topsoil [-]
@@ -81,19 +81,19 @@
     # Exponent in soil runoff generation equation [-]
     beta::Vector{SVector{N,T}} | "-"
     # maximum percolation rate [mm Δt⁻¹]
-    perc::Vector{SVector{N,T}} | "mm Δt-1"
+    perc::Vector{SVector{N,T}} | "mm dt-1"
     # maximum capillary rise rate [mm Δt⁻¹]
-    cap::Vector{SVector{N,T}} | "mm Δt-1"
+    cap::Vector{SVector{N,T}} | "mm dt-1"
     #FAST
     # Exponent for non linear recession [-]
     alfa::Vector{SVector{N,T}} | "-"
-    #recession coefficient of fast storage [Δt-1]
-    kf::Vector{SVector{N,T}} | "Δt-1"
+    #recession coefficient of fast storage [dt-1]
+    kf::Vector{SVector{N,T}} | "dt-1"
     # fraction of qrootzone to slowstorage (1-ds to faststorage)
     ds::Vector{SVector{N,T}} | "-"
     # SLOW
-    #recession coefficient of slow storage [Δt-1]
-    ks::Vector{T} | "Δt-1"
+    #recession coefficient of slow storage [dt-1]
+    ks::Vector{T} | "dt-1"
 
     ## STATES
     ##SNOW
@@ -130,96 +130,96 @@
     ## FLUXES
     #SNOW
     # Precipitation [mm Δt⁻¹]
-    precipitation::Vector{T} | "mm Δt-1"
+    precipitation::Vector{T} | "mm dt-1"
     # Temperature [ᵒC]
     temperature::Vector{T} | "ᵒC"
     # Potential evapotranspiration [mm Δt⁻¹]
-    potential_evaporation::Vector{T} | "mm Δt-1"
+    potential_evaporation::Vector{T} | "mm dt-1"
     # Potential evapotranspiration corrected [mm Δt⁻¹]
-    epotcorr::Vector{T} | "mm Δt-1"
+    epotcorr::Vector{T} | "mm dt-1"
     # Precipitation corrected [mm Δt⁻¹]
-    precipcorr::Vector{T} | "mm Δt-1"
+    precipcorr::Vector{T} | "mm dt-1"
     # Snow melt + precipitation as rainfall [mm]
-    rainfallplusmelt::Vector{T} | "mm Δt-1"
+    rainfallplusmelt::Vector{T} | "mm dt-1"
     # Snowfall [mm]
-    snowfall::Vector{T} | "mm Δt-1"
+    snowfall::Vector{T} | "mm dt-1"
     # Snowmelt [mm]
-    snowmelt::Vector{T} | "mm Δt-1"
+    snowmelt::Vector{T} | "mm dt-1"
     #INTERCEPTION
     # Potential soil evaporation [mm Δt⁻¹]
-    potsoilevap::Vector{SVector{N,T}} | "mm Δt-1"
+    potsoilevap::Vector{SVector{N,T}} | "mm dt-1"
     # Evaporation from interception storage [mm Δt⁻¹]
-    intevap::Vector{SVector{N,T}} | "mm Δt-1"
+    intevap::Vector{SVector{N,T}} | "mm dt-1"
     #effective precipitation [mm Δt⁻¹]
-    precipeffective::Vector{SVector{N,T}} | "mm Δt-1"
+    precipeffective::Vector{SVector{N,T}} | "mm dt-1"
     # Evaporation from interception sum classes [mm Δt⁻¹]
-    intevap_m::Vector{T} | "mm Δt-1"
+    intevap_m::Vector{T} | "mm dt-1"
     #HORTONPONDING
     # Evaporation from the hortonion ponding storage [-]
-    hortonevap::Vector{SVector{N,T}} | "mm Δt-1"
+    hortonevap::Vector{SVector{N,T}} | "mm dt-1"
     # Flux from the hortonian ponding storage to the hortonian runoff storage [mm Δt⁻¹]
-    qhortonpond::Vector{SVector{N,T}} | "mm Δt-1"
+    qhortonpond::Vector{SVector{N,T}} | "mm dt-1"
     # Flux from the hortonian ponding storage to the root zone storage [mm Δt⁻¹]
-    qhortonrootzone::Vector{SVector{N,T}} | "mm Δt-1"
+    qhortonrootzone::Vector{SVector{N,T}} | "mm dt-1"
     # modeled accumulated frost [ᵒC Δt]
     facc::Vector{SVector{N,T}} | "ᵒC Δt"
     # Evaporation from the hortonian sum classes [mm Δt⁻¹]
-    hortonevap_m::Vector{T} | "mm Δt-1"
+    hortonevap_m::Vector{T} | "mm dt-1"
     #HORTONRUNOFF
     # Flux from the hortonian runoff storage [mm Δt⁻¹]
-    qhortonrun::Vector{SVector{N,T}} | "mm Δt-1"
+    qhortonrun::Vector{SVector{N,T}} | "mm dt-1"
     #ROOTZONE
     # Evaporation from the root-zone storage [mm Δt⁻¹]
-    rootevap::Vector{SVector{N,T}} | "mm Δt-1"
+    rootevap::Vector{SVector{N,T}} | "mm dt-1"
     # Flux from the root-zone storage [mm Δt⁻¹]
-    qrootzone::Vector{SVector{N,T}} | "mm Δt-1"
+    qrootzone::Vector{SVector{N,T}} | "mm dt-1"
     # Pref. recharge to fast storage [mm Δt⁻¹]
-    qrootzonefast::Vector{SVector{N,T}} | "mm Δt-1"
+    qrootzonefast::Vector{SVector{N,T}} | "mm dt-1"
     # Pref. recharge to slow storage sum classes [mm Δt⁻¹]
-    qrootzoneslow_m::Vector{T} | "mm Δt-1"
+    qrootzoneslow_m::Vector{T} | "mm dt-1"
     # Capillary flux from the slow to the root-zone storage [mm Δt⁻¹]
-    qcapillary::Vector{SVector{N,T}} | "mm Δt-1"
+    qcapillary::Vector{SVector{N,T}} | "mm dt-1"
     # Capillary flux from the slow to the root-zone storage sum classes [mm Δt⁻¹]
-    qcapillary_m::Vector{T} | "mm Δt-1"
+    qcapillary_m::Vector{T} | "mm dt-1"
     # Percolation flux from the root-zone to the slow storage [mm Δt⁻¹]
-    qpercolation::Vector{SVector{N,T}} | "mm Δt-1"
+    qpercolation::Vector{SVector{N,T}} | "mm dt-1"
     # Percolation flux from the root-zone to the slow storage sum classes [mm Δt⁻¹]
-    qpercolation_m::Vector{T} | "mm Δt-1"
+    qpercolation_m::Vector{T} | "mm dt-1"
     # Evaporation from the root-zone storage, interception and hortonian [mm Δt⁻¹]
-    actevap::Vector{SVector{N,T}} | "mm Δt-1"
+    actevap::Vector{SVector{N,T}} | "mm dt-1"
     # Evaporation from the root-zone storage, interception and hortonian sum classes [mm Δt⁻¹]
-    actevap_m::Vector{T} | "mm Δt-1"
+    actevap_m::Vector{T} | "mm dt-1"
     # Evaporation from the root-zone storage sum classes [mm Δt⁻¹]
-    rootevap_m::Vector{T} | "mm Δt-1"
+    rootevap_m::Vector{T} | "mm dt-1"
     #FAST
     # runoff from fast reservoir [mm Δt⁻¹]
-    qfast::Vector{SVector{N,T}} | "mm Δt-1"
+    qfast::Vector{SVector{N,T}} | "mm dt-1"
     #SLOW
     # runoff from slow reservoir [mm Δt⁻¹]
-    qslow::Vector{T} | "mm Δt-1"
+    qslow::Vector{T} | "mm dt-1"
     #Total [mm Δt⁻¹]
-    runoff::Vector{T} | "mm Δt-1"
+    runoff::Vector{T} | "mm dt-1"
     # fast runoff sum classes [mm Δt⁻¹]
-    qfast_tot = Vector{T} | "mm Δt-1"
+    qfast_tot = Vector{T} | "mm dt-1"
 
 
     ## WATERBALANCES
     #water balance snow store
-    wb_snow::Vector{T} | "mm Δt-1"
+    wb_snow::Vector{T} | "mm dt-1"
     #water balance interception storage [mm Δt⁻¹]
-    wb_interception::Vector{SVector{N,T}} | "mm Δt-1"
+    wb_interception::Vector{SVector{N,T}} | "mm dt-1"
     #water balance hortonian ponding storage [mm Δt⁻¹]
-    wb_hortonponding::Vector{SVector{N,T}} | "mm Δt-1"
+    wb_hortonponding::Vector{SVector{N,T}} | "mm dt-1"
     #water balance hortonian runoff storage [mm Δt⁻¹]
-    wb_hortonrunoff::Vector{SVector{N,T}} | "mm Δt-1"
+    wb_hortonrunoff::Vector{SVector{N,T}} | "mm dt-1"
     #water balance root-zone storage [mm Δt⁻¹]
-    wb_rootzone::Vector{SVector{N,T}} | "mm Δt-1"
+    wb_rootzone::Vector{SVector{N,T}} | "mm dt-1"
     #water balance fast storage [mm Δt⁻¹]
-    wb_fast::Vector{SVector{N,T}} | "mm Δt-1"
+    wb_fast::Vector{SVector{N,T}} | "mm dt-1"
     #water balance slow storage [mm Δt⁻¹]
-    wb_slow::Vector{T} | "mm Δt-1"
+    wb_slow::Vector{T} | "mm dt-1"
     #total water balance [mm Δt⁻¹]
-    wb_tot::Vector{T} | "mm Δt-1"
+    wb_tot::Vector{T} | "mm dt-1"
 
 end
 
@@ -302,7 +302,7 @@ function common_glaciers(flextopo::FLEXTOPO, config)
                     flextopo.g_tt[i],
                     flextopo.g_cfmax[i],
                     flextopo.g_sifrac[i],
-                    Second(flextopo.Δt),
+                    Second(flextopo.dt),
                 )
             # Convert to mm per grid cell and add to snowmelt
             glaciermelt = glaciermelt * flextopo.glacierfrac[i]
diff --git a/src/flextopo_model.jl b/src/flextopo_model.jl
index c8b3975e7..5575facba 100644
--- a/src/flextopo_model.jl
+++ b/src/flextopo_model.jl
@@ -10,7 +10,7 @@ function initialize_flextopo_model(config::Config)
 
     reader = prepare_reader(config)
     clock = Clock(config, reader)
-    Δt = clock.Δt
+    dt = clock.dt
 
     do_reservoirs = get(config.model, "reservoirs", false)::Bool
     do_lakes = get(config.model, "lakes", false)::Bool
@@ -82,7 +82,7 @@ function initialize_flextopo_model(config::Config)
             defaults = 3.0,
             type = Float,
             fill = 0.0,
-        ) .* (Δt / basetimestep)
+        ) .* (dt / basetimestep)
 
     g_sifrac =
         ncread(
@@ -93,7 +93,7 @@ function initialize_flextopo_model(config::Config)
             defaults = 0.001,
             type = Float,
             fill = 0.0,
-        ) .* (Δt / basetimestep)
+        ) .* (dt / basetimestep)
     glacierfrac = ncread(
         nc,
         config,
@@ -122,7 +122,7 @@ function initialize_flextopo_model(config::Config)
             sel = inds,
             defaults = 3.75653,
             type = Float,
-        ) .* (Δt / basetimestep)
+        ) .* (dt / basetimestep)
 
     tt = ncread(nc, config, "vertical.tt"; sel = inds, defaults = -1.41934, type = Float)
 
@@ -141,7 +141,7 @@ function initialize_flextopo_model(config::Config)
     sfcf = ncread(nc, config, "vertical.sfcf"; sel = inds, defaults = 1.0, type = Float)
     ks =
         ncread(nc, config, "vertical.ks"; sel = inds, defaults = 0.006, type = Float) .*
-        (Δt / basetimestep)
+        (dt / basetimestep)
 
     #initialize parameters that differ per class
     hrufrac = ncread(
@@ -180,7 +180,7 @@ function initialize_flextopo_model(config::Config)
             defaults = 0.5,
             type = Float,
             dimname = :classes,
-        ) .* (Δt / basetimestep)
+        ) .* (dt / basetimestep)
     facc0 = ncread(
         nc,
         config,
@@ -217,7 +217,7 @@ function initialize_flextopo_model(config::Config)
             defaults = 2.0,
             type = Float,
             dimname = :classes,
-        ) .* (Δt / basetimestep)
+        ) .* (dt / basetimestep)
     shmin = ncread(
         nc,
         config,
@@ -272,7 +272,7 @@ function initialize_flextopo_model(config::Config)
             defaults = 0.30,
             type = Float,
             dimname = :classes,
-        ) .* (Δt / basetimestep)
+        ) .* (dt / basetimestep)
     cap =
         ncread(
             nc,
@@ -282,7 +282,7 @@ function initialize_flextopo_model(config::Config)
             defaults = 0.20,
             type = Float,
             dimname = :classes,
-        ) .* (Δt / basetimestep)
+        ) .* (dt / basetimestep)
     kf =
         ncread(
             nc,
@@ -292,7 +292,7 @@ function initialize_flextopo_model(config::Config)
             defaults = 0.1,
             type = Float,
             dimname = :classes,
-        ) .* (Δt / basetimestep)
+        ) .* (dt / basetimestep)
     alfa = ncread(
         nc,
         config,
@@ -354,7 +354,7 @@ function initialize_flextopo_model(config::Config)
 
 
     flextopo = FLEXTOPO{Float,nclass}(
-        Δt = Float(tosecond(Δt)),
+        dt = Float(tosecond(dt)),
         nclass = nclass,
         n = n,
         dic_function = dic_function,
@@ -493,7 +493,7 @@ function initialize_flextopo_model(config::Config)
     pits = zeros(Bool, modelsize_2d)
     if do_reservoirs
         reservoirs, resindex, reservoir, pits =
-            initialize_simple_reservoir(config, nc, inds_riv, nriv, pits, tosecond(Δt))
+            initialize_simple_reservoir(config, nc, inds_riv, nriv, pits, tosecond(dt))
     else
         reservoir = ()
         reservoirs = nothing
@@ -503,7 +503,7 @@ function initialize_flextopo_model(config::Config)
     # lakes
     if do_lakes
         lakes, lakeindex, lake, pits =
-            initialize_lake(config, nc, inds_riv, nriv, pits, tosecond(Δt))
+            initialize_lake(config, nc, inds_riv, nriv, pits, tosecond(dt))
     else
         lake = ()
         lakes = nothing
@@ -517,9 +517,9 @@ function initialize_flextopo_model(config::Config)
         ldd = set_pit_ldd(pits_2d, ldd, inds)
     end
 
-    βₗ =
+    landslope =
         ncread(nc, config, "lateral.land.slope"; optional = false, sel = inds, type = Float)
-    clamp!(βₗ, 0.00001, Inf)
+    clamp!(landslope, 0.00001, Inf)
 
     dl = map(detdrainlength, ldd, xl, yl)
     dw = (xl .* yl) ./ dl
@@ -527,12 +527,12 @@ function initialize_flextopo_model(config::Config)
         nc,
         config,
         inds;
-        sl = βₗ,
+        sl = landslope,
         dl = dl,
         width = map(det_surfacewidth, dw, riverwidth, river),
         iterate = kinwave_it,
         tstep = kw_land_tstep,
-        Δt = Δt,
+        dt = dt,
     )
 
     graph = flowgraph(ldd, inds, pcr_dir)
@@ -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, βₗ, n)
+    frac_toriver = fraction_runoff_toriver(graph, ldd, index_river, landslope, n)
 
     rf = initialize_surfaceflow_river(
         nc,
@@ -562,7 +562,7 @@ function initialize_flextopo_model(config::Config)
         lake = lakes,
         iterate = kinwave_it,
         tstep = kw_river_tstep,
-        Δt = Δt,
+        dt = dt,
     )
 
     # setup subdomains for the land and river kinematic wave domain, if nthreads = 1
diff --git a/src/flow.jl b/src/flow.jl
index d81873196..a26cf310d 100644
--- a/src/flow.jl
+++ b/src/flow.jl
@@ -3,7 +3,7 @@ abstract type SurfaceFlow end
 
 @get_units @exchange @grid_type @grid_location @with_kw struct SurfaceFlowRiver{T,R,L} <:
                                                                SurfaceFlow
-    β::T | "-" | 0 | "scalar"                    # constant in Manning's equation
+    beta::T | "-" | 0 | "scalar"                 # constant in Manning's equation
     sl::Vector{T} | "m m-1"                      # Slope [m m⁻¹]
     n::Vector{T} | "s m-1/3"                     # Manning's roughness [s m⁻⅓]
     dl::Vector{T} | "m"                          # Drain length [m]
@@ -18,12 +18,12 @@ abstract type SurfaceFlow end
     h::Vector{T} | "m"                           # Water level [m]
     h_av::Vector{T} | "m"                        # Average water level [m]
     bankfull_depth::Vector{T} | "m"              # Bankfull water level [m]
-    Δt::T | "s" | 0 | "none" | "none"            # Model time step [s]
+    dt::T | "s" | 0 | "none" | "none"            # Model time step [s]
     its::Int | "-" | 0 | "none" | "none"         # Number of fixed iterations
     width::Vector{T} | "m"                       # Flow width [m]
-    alpha_pow::T | "-" | 0 | "scalar"            # Used in the power part of α
-    alpha_term::Vector{T} | "-"                  # Term used in computation of α
-    α::Vector{T} | "s3/5 m1/5"                   # Constant in momentum equation A = αQᵝ, based on Manning's equation
+    alpha_pow::T | "-" | 0 | "scalar"            # Used in the power part of alpha
+    alpha_term::Vector{T} | "-"                  # Term used in computation of alpha
+    alpha::Vector{T} | "s3/5 m1/5"               # Constant in momentum equation A = alpha*Q^beta, based on Manning's equation
     cel::Vector{T} | "m s-1"                     # Celerity of the kinematic wave
     reservoir_index::Vector{Int} | "-"           # map cell to 0 (no reservoir) or i (pick reservoir i in reservoir field)
     lake_index::Vector{Int} | "-"                # map cell to 0 (no lake) or i (pick lake i in lake field)
@@ -40,7 +40,7 @@ end
 
 @get_units @exchange @grid_type @grid_location @with_kw struct SurfaceFlowLand{T} <:
                                                                SurfaceFlow
-    β::T | "-" | 0 | "scalar"                       # constant in Manning's equation
+    beta::T | "-" | 0 | "scalar"                       # constant in Manning's equation
     sl::Vector{T} | "m m-1"                         # Slope [m m⁻¹]
     n::Vector{T} | "s m-1/3"                        # Manning's roughness [s m⁻⅓]
     dl::Vector{T} | "m"                             # Drain length [m]
@@ -52,18 +52,18 @@ end
     volume::Vector{T} | "m3"                        # Kinematic wave volume [m³] (based on water level h)
     h::Vector{T} | "m"                              # Water level [m]
     h_av::Vector{T} | "m"                           # Average water level [m]
-    Δt::T | "s" | 0 | "none" | "none"               # Model time step [s]
+    dt::T | "s" | 0 | "none" | "none"               # Model time step [s]
     its::Int | "-" | 0 | "none" | "none"            # Number of fixed iterations
     width::Vector{T} | "m"                          # Flow width [m]
-    alpha_pow::T | "-" | 0 | "scalar"               # Used in the power part of α
-    alpha_term::Vector{T} | "-"                     # Term used in computation of α
-    α::Vector{T} | "s3/5 m1/5"                      # Constant in momentum equation A = αQᵝ, based on Manning's equation
+    alpha_pow::T | "-" | 0 | "scalar"               # Used in the power part of alpha
+    alpha_term::Vector{T} | "-"                     # Term used in computation of alpha
+    alpha::Vector{T} | "s3/5 m1/5"                  # Constant in momentum equation A = alpha * Q^beta, based on Manning's equation
     cel::Vector{T} | "m s-1"                        # Celerity of the kinematic wave
     to_river::Vector{T} | "m3 s-1"                  # Part of overland flow [m³ s⁻¹] that flows to the river
     kinwave_it::Bool | "-" | 0 | "none" | "none"    # Boolean for iterations kinematic wave
 end
 
-function initialize_surfaceflow_land(nc, config, inds; sl, dl, width, iterate, tstep, Δt)
+function initialize_surfaceflow_land(nc, config, inds; sl, dl, width, iterate, tstep, dt)
     @info "Kinematic wave approach is used for overland flow." iterate
     if tstep > 0
         @info "Using a fixed sub-timestep (seconds) $tstep for kinematic wave overland flow."
@@ -74,7 +74,7 @@ function initialize_surfaceflow_land(nc, config, inds; sl, dl, width, iterate, t
     n = length(inds)
 
     sf_land = SurfaceFlowLand(
-        β = Float(0.6),
+        beta = Float(0.6),
         sl = sl,
         n = n_land,
         dl = dl,
@@ -86,12 +86,12 @@ function initialize_surfaceflow_land(nc, config, inds; sl, dl, width, iterate, t
         volume = zeros(Float, n),
         h = zeros(Float, n),
         h_av = zeros(Float, n),
-        Δt = Float(tosecond(Δt)),
-        its = tstep > 0 ? Int(cld(tosecond(Δt), tstep)) : tstep,
+        dt = Float(tosecond(dt)),
+        its = tstep > 0 ? Int(cld(tosecond(dt), tstep)) : tstep,
         width = width,
         alpha_pow = Float((2.0 / 3.0) * 0.6),
         alpha_term = fill(mv, n),
-        α = fill(mv, n),
+        alpha = fill(mv, n),
         cel = zeros(Float, n),
         to_river = zeros(Float, n),
         kinwave_it = iterate,
@@ -112,7 +112,7 @@ function initialize_surfaceflow_river(
     lake,
     iterate,
     tstep,
-    Δt,
+    dt,
 )
     @info "Kinematic wave approach is used for river flow." iterate
     if tstep > 0
@@ -149,7 +149,7 @@ function initialize_surfaceflow_river(
     n = length(inds)
 
     sf_river = SurfaceFlowRiver(
-        β = Float(0.6),
+        beta = Float(0.6),
         sl = sl,
         n = n_river,
         dl = dl,
@@ -164,12 +164,12 @@ function initialize_surfaceflow_river(
         h = zeros(Float, n),
         h_av = zeros(Float, n),
         bankfull_depth = bankfull_depth,
-        Δt = Float(tosecond(Δt)),
-        its = tstep > 0 ? Int(cld(tosecond(Δt), tstep)) : tstep,
+        dt = Float(tosecond(dt)),
+        its = tstep > 0 ? Int(cld(tosecond(dt), tstep)) : tstep,
         width = width,
         alpha_pow = Float((2.0 / 3.0) * 0.6),
         alpha_term = fill(mv, n),
-        α = fill(mv, n),
+        alpha = fill(mv, n),
         cel = zeros(Float, n),
         reservoir_index = reservoir_index,
         lake_index = lake_index,
@@ -188,16 +188,16 @@ function update(sf::SurfaceFlowLand, network, frac_toriver)
 
     ns = length(subdomain_order)
 
-    @. sf.alpha_term = pow(sf.n / sqrt(sf.sl), sf.β)
+    @. sf.alpha_term = pow(sf.n / sqrt(sf.sl), sf.beta)
     # use fixed alpha value based flow width
-    @. sf.α = sf.alpha_term * pow(sf.width, sf.alpha_pow)
+    @. sf.alpha = sf.alpha_term * pow(sf.width, sf.alpha_pow)
     @. sf.qlat = sf.inwater / sf.dl
 
     sf.q_av .= 0.0
     sf.h_av .= 0.0
     sf.to_river .= 0.0
 
-    Δt, its = stable_timestep(sf)
+    dt, its = stable_timestep(sf)
     for _ = 1:its
         sf.qin .= 0.0
         for k = 1:ns
@@ -225,15 +225,15 @@ function update(sf::SurfaceFlowLand, network, frac_toriver)
                         sf.qin[v],
                         sf.q[v],
                         sf.qlat[v],
-                        sf.α[v],
-                        sf.β,
-                        Δt,
+                        sf.alpha[v],
+                        sf.beta,
+                        dt,
                         sf.dl[v],
                     )
 
                     # update h, only if surface width > 0.0
                     if sf.width[v] > 0.0
-                        crossarea = sf.α[v] * pow(sf.q[v], sf.β)
+                        crossarea = sf.alpha[v] * pow(sf.q[v], sf.beta)
                         sf.h[v] = crossarea / sf.width[v]
                     end
                     sf.q_av[v] += sf.q[v]
@@ -255,9 +255,9 @@ function update(sf::SurfaceFlowRiver, network, doy)
 
     ns = length(subdomain_order)
 
-    @. sf.alpha_term = pow(sf.n / sqrt(sf.sl), sf.β)
+    @. sf.alpha_term = pow(sf.n / sqrt(sf.sl), sf.beta)
     # use fixed alpha value based on 0.5 * bankfull_depth
-    @. sf.α = sf.alpha_term * pow(sf.width + sf.bankfull_depth, sf.alpha_pow)
+    @. sf.alpha = sf.alpha_term * pow(sf.width + sf.bankfull_depth, sf.alpha_pow)
 
     @. sf.qlat = sf.inwater / sf.dl
 
@@ -276,7 +276,7 @@ function update(sf::SurfaceFlowRiver, network, doy)
         sf.lake.actevap .= 0.0
     end
 
-    Δt, its = stable_timestep(sf)
+    dt, its = stable_timestep(sf)
     for _ = 1:its
         sf.qin .= 0.0
         for k = 1:ns
@@ -298,9 +298,9 @@ function update(sf::SurfaceFlowRiver, network, doy)
                         sf.qin[v],
                         sf.q[v],
                         sf.qlat[v] + inflow,
-                        sf.α[v],
-                        sf.β,
-                        Δt,
+                        sf.alpha[v],
+                        sf.beta,
+                        dt,
                         sf.dl[v],
                     )
 
@@ -308,7 +308,7 @@ function update(sf::SurfaceFlowRiver, network, doy)
                         # run reservoir model and copy reservoir outflow to inflow (qin) of
                         # downstream river cell
                         i = sf.reservoir_index[v]
-                        update(sf.reservoir, i, sf.q[v] + sf.inflow_wb[v], Δt)
+                        update(sf.reservoir, i, sf.q[v] + sf.inflow_wb[v], dt)
 
                         downstream_nodes = outneighbors(graph, v)
                         n_downstream = length(downstream_nodes)
@@ -329,7 +329,7 @@ function update(sf::SurfaceFlowRiver, network, doy)
                         # run lake model and copy lake outflow to inflow (qin) of downstream river
                         # cell
                         i = sf.lake_index[v]
-                        update(sf.lake, i, sf.q[v] + sf.inflow_wb[v], doy, Δt)
+                        update(sf.lake, i, sf.q[v] + sf.inflow_wb[v], doy, dt)
 
                         downstream_nodes = outneighbors(graph, v)
                         n_downstream = length(downstream_nodes)
@@ -348,7 +348,7 @@ function update(sf::SurfaceFlowRiver, network, doy)
                     end
 
                     # update h
-                    crossarea = sf.α[v] * pow(sf.q[v], sf.β)
+                    crossarea = sf.alpha[v] * pow(sf.q[v], sf.beta)
                     sf.h[v] = crossarea / sf.width[v]
                     sf.q_av[v] += sf.q[v]
                     sf.h_av[v] += sf.h[v]
@@ -372,8 +372,8 @@ function stable_timestep(sf::S) where {S<:SurfaceFlow}
             courant = zeros(n)
             for v = 1:n
                 if sf.q[v] > 0.0
-                    sf.cel[v] = 1.0 / (sf.α[v] * sf.β * pow(sf.q[v], (sf.β - 1.0)))
-                    courant[v] = (sf.Δt / sf.dl[v]) * sf.cel[v]
+                    sf.cel[v] = 1.0 / (sf.alpha[v] * sf.beta * pow(sf.q[v], (sf.beta - 1.0)))
+                    courant[v] = (sf.dt / sf.dl[v]) * sf.cel[v]
                 end
             end
             filter!(x -> x ≠ 0.0, courant)
@@ -384,26 +384,26 @@ function stable_timestep(sf::S) where {S<:SurfaceFlow}
     end
 
     # sub time step
-    Δt = sf.Δt / its
-    return Δt, its
+    dt = sf.dt / its
+    return dt, its
 end
 
 @get_units @exchange @grid_type @grid_location @with_kw struct LateralSSF{T}
-    kh₀::Vector{T} | "m d-1"               # Horizontal hydraulic conductivity at soil surface [m d⁻¹]
-    f::Vector{T} | "m-1"                   # A scaling parameter [m⁻¹] (controls exponential decline of kh₀)
+    kh_0::Vector{T} | "m d-1"               # Horizontal hydraulic conductivity at soil surface [m d⁻¹]
+    f::Vector{T} | "m-1"                   # A scaling parameter [m⁻¹] (controls exponential decline of kh_0)
     kh::Vector{T} | "m d-1"                # Horizontal hydraulic conductivity [m d⁻¹]
     khfrac::Vector{T} | "-"                # A muliplication factor applied to vertical hydraulic conductivity `kv` [-]
     soilthickness::Vector{T} | "m"         # Soil thickness [m]
-    θₛ::Vector{T} | "-"                     # Saturated water content (porosity) [-]
-    θᵣ::Vector{T} | "-"                    # Residual water content [-]
-    Δt::T | "d" | 0 | "none" | "none"      # model time step [d]
-    βₗ::Vector{T} | "m m-1"                 # Slope [m m⁻¹]
+    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]
+    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)
-    z_exp::Vector{T} | "m"                 # Depth [m] from soil surface for which exponential decline of kv₀ is valid
-    exfiltwater::Vector{T} | "m Δt-1"      # Exfiltration [m Δt⁻¹] (groundwater above surface level, saturated excess conditions)
-    recharge::Vector{T} | "m2 Δt-1"        # Net recharge to saturated store [m² Δt⁻¹]
+    z_exp::Vector{T} | "m"                 # Depth [m] from soil surface for which exponential decline of kv_0 is valid
+    exfiltwater::Vector{T} | "m dt-1"      # Exfiltration [m Δt⁻¹] (groundwater above surface level, saturated excess conditions)
+    recharge::Vector{T} | "m2 dt-1"        # Net recharge to saturated store [m² Δt⁻¹]
     ssf::Vector{T} | "m3 d-1"              # Subsurface flow [m³ d⁻¹]
     ssfin::Vector{T} | "m3 d-1"            # Inflow from upstream cells [m³ d⁻¹]
     ssfmax::Vector{T} | "m2 d-1"           # Maximum subsurface flow [m² d⁻¹]
@@ -445,12 +445,12 @@ function update(ssf::LateralSSF, network, frac_toriver, ksat_profile)
                         ssf.ssf[v],
                         ssf.zi[v],
                         ssf.recharge[v],
-                        ssf.kh₀[v],
-                        ssf.βₗ[v],
-                        ssf.θₛ[v] - ssf.θᵣ[v],
+                        ssf.kh_0[v],
+                        ssf.slope[v],
+                        ssf.theta_s[v] - ssf.theta_r[v],
                         ssf.f[v],
                         ssf.soilthickness[v],
-                        ssf.Δt,
+                        ssf.dt,
                         ssf.dl[v],
                         ssf.dw[v],
                         ssf.ssfmax[v],
@@ -465,10 +465,10 @@ function update(ssf::LateralSSF, network, frac_toriver, ksat_profile)
                         ssf.zi[v],
                         ssf.recharge[v],
                         ssf.kh[v],
-                        ssf.βₗ[v],
-                        ssf.θₛ[v] - ssf.θᵣ[v],
+                        ssf.slope[v],
+                        ssf.theta_s[v] - ssf.theta_r[v],
                         ssf.soilthickness[v],
-                        ssf.Δt,
+                        ssf.dt,
                         ssf.dl[v],
                         ssf.dw[v],
                         ssf.ssfmax[v],
@@ -480,8 +480,8 @@ function update(ssf::LateralSSF, network, frac_toriver, ksat_profile)
 end
 
 @get_units @exchange @grid_type @grid_location @with_kw struct GroundwaterExchange{T}
-    Δt::T | "d" | 0 | "none" | "none"   # model time step [d]
-    exfiltwater::Vector{T} | "m Δt-1"   # Exfiltration [m Δt⁻¹]  (groundwater above surface level, saturated excess conditions)
+    dt::T | "d" | 0 | "none" | "none"   # model time step [d]
+    exfiltwater::Vector{T} | "m dt-1"   # Exfiltration [m Δt⁻¹]  (groundwater above surface level, saturated excess conditions)
     zi::Vector{T} | "m"                 # Pseudo-water table depth [m] (top of the saturated zone)
     to_river::Vector{T} | "m3 d-1"      # Part of subsurface flow [m³ d⁻¹] that flows to the river
     ssf::Vector{T} | "m3 d-1"           # Subsurface flow [m³ d⁻¹]
@@ -493,9 +493,9 @@ end
     active_n::Vector{Int} | "-"                             # active nodes
     active_e::Vector{Int} | "-" | _ | "edge"                # active edges/links
     g::T | "m s-2" | 0 | "scalar"                           # acceleration due to gravity
-    α::T | "-" | 0 | "scalar"                               # stability coefficient (Bates et al., 2010)
+    alpha::T | "-" | 0 | "scalar"                           # stability coefficient (Bates et al., 2010)
     h_thresh::T | "m" | 0 | "scalar"                        # depth threshold for calculating flow
-    Δt::T | "s" | 0 | "none" | "none"                       # model time step [s]
+    dt::T | "s" | 0 | "none" | "none"                       # model time step [s]
     q::Vector{T} | "m3 s-1" | _ | "edge"                    # river discharge (subgrid channel)
     q0::Vector{T} | "m3 s-1" | _ | "edge"                   # river discharge (subgrid channel) at previous time step
     q_av::Vector{T} | "m3 s-1" | _ | "edge"                 # average river channel (+ floodplain) discharge [m³ s⁻¹]
@@ -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
-    η_max::Vector{T} | "m" | _ | "edge"                     # maximum water elevation at edge
-    η_src::Vector{T} | "m"                                  # water elevation of source node of edge
-    η_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
@@ -544,7 +544,7 @@ function initialize_shallowwater_river(
     reservoir,
     lake_index,
     lake,
-    Δt,
+    dt,
     floodplain,
 )
     # The local inertial approach makes use of a staggered grid (Bates et al. (2010)),
@@ -674,9 +674,9 @@ function initialize_shallowwater_river(
         active_n = active_index,
         active_e = active_index,
         g = 9.80665,
-        α = alpha,
+        alpha = alpha,
         h_thresh = h_thresh,
-        Δt = tosecond(Δt),
+        dt = tosecond(dt),
         q = zeros(_ne),
         q0 = zeros(_ne),
         q_av = q_av,
@@ -685,9 +685,9 @@ function initialize_shallowwater_river(
         mannings_n_sq = mannings_n_sq,
         mannings_n = n_river,
         h = h,
-        η_max = zeros(_ne),
-        η_src = zeros(_ne),
-        η_dst = zeros(_ne),
+        zs_max = zeros(_ne),
+        zs_src = zeros(_ne),
+        zs_dst = zeros(_ne),
         hf = zeros(_ne),
         h_av = zeros(n),
         width = width,
@@ -724,7 +724,7 @@ function get_inflow_waterbody(sw::ShallowWaterRiver, src_edge)
     return q_in
 end
 
-function shallowwater_river_update(sw::ShallowWaterRiver, network, Δt, doy, update_h)
+function shallowwater_river_update(sw::ShallowWaterRiver, network, dt, doy, update_h)
 
     @unpack nodes_at_link, links_at_node = network
 
@@ -736,11 +736,11 @@ function shallowwater_river_update(sw::ShallowWaterRiver, network, Δt, doy, upd
         i = sw.active_e[j]
         i_src = nodes_at_link.src[i]
         i_dst = nodes_at_link.dst[i]
-        sw.η_src[i] = sw.zb[i_src] + sw.h[i_src]
-        sw.η_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.η_max[i] = max(sw.η_src[i], sw.η_dst[i])
-        sw.hf[i] = (sw.η_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, Δt, doy, upd
             sw.hf[i] > sw.h_thresh,
             local_inertial_flow(
                 sw.q0[i],
-                sw.η_src[i],
-                sw.η_dst[i],
+                sw.zs_src[i],
+                sw.zs_dst[i],
                 sw.hf[i],
                 sw.a[i],
                 sw.r[i],
@@ -758,7 +758,7 @@ function shallowwater_river_update(sw::ShallowWaterRiver, network, Δt, doy, upd
                 sw.mannings_n_sq[i],
                 sw.g,
                 sw.froude_limit,
-                Δt,
+                dt,
             ),
             0.0,
         )
@@ -767,10 +767,10 @@ function shallowwater_river_update(sw::ShallowWaterRiver, network, Δt, doy, upd
         sw.q[i] = IfElse.ifelse(sw.h[i_src] <= 0.0, min(sw.q[i], 0.0), sw.q[i])
         sw.q[i] = IfElse.ifelse(sw.h[i_dst] <= 0.0, max(sw.q[i], 0.0), sw.q[i])
 
-        sw.q_av[i] += sw.q[i] * Δt
+        sw.q_av[i] += sw.q[i] * dt
     end
     if !isnothing(sw.floodplain)
-        @tturbo @. sw.floodplain.hf = max(sw.η_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, Δt, doy, upd
                 sw.floodplain.a[i] > 1.0e-05,
                 local_inertial_flow(
                     sw.floodplain.q0[i],
-                    sw.η_src[i],
-                    sw.η_dst[i],
+                    sw.zs_src[i],
+                    sw.zs_dst[i],
                     sw.floodplain.hf[i],
                     sw.floodplain.a[i],
                     sw.floodplain.r[i],
@@ -839,7 +839,7 @@ function shallowwater_river_update(sw::ShallowWaterRiver, network, Δt, doy, upd
                     sw.floodplain.mannings_n_sq[i],
                     sw.g,
                     sw.froude_limit,
-                    Δt,
+                    dt,
                 ),
                 0.0,
             )
@@ -858,7 +858,7 @@ function shallowwater_river_update(sw::ShallowWaterRiver, network, Δt, doy, upd
 
             sw.floodplain.q[i] =
                 IfElse.ifelse(sw.floodplain.q[i] * sw.q[i] < 0.0, 0.0, sw.floodplain.q[i])
-            sw.floodplain.q_av[i] += sw.floodplain.q[i] * Δt
+            sw.floodplain.q_av[i] += sw.floodplain.q[i] * dt
         end
     end
     # For reservoir and lake locations the local inertial solution is replaced by the
@@ -868,35 +868,35 @@ function shallowwater_river_update(sw::ShallowWaterRiver, network, Δt, doy, upd
         i = sw.reservoir_index[v]
 
         q_in = get_inflow_waterbody(sw, links_at_node.src[i])
-        update(sw.reservoir, v, q_in + sw.inflow_wb[i], Δt)
+        update(sw.reservoir, v, q_in + sw.inflow_wb[i], dt)
         sw.q[i] = sw.reservoir.outflow[v]
-        sw.q_av[i] += sw.q[i] * Δt
+        sw.q_av[i] += sw.q[i] * dt
     end
     for v in eachindex(sw.lake_index)
         i = sw.lake_index[v]
 
         q_in = get_inflow_waterbody(sw, links_at_node.src[i])
-        update(sw.lake, v, q_in + sw.inflow_wb[i], doy, Δt)
+        update(sw.lake, v, q_in + sw.inflow_wb[i], doy, dt)
         sw.q[i] = sw.lake.outflow[v]
-        sw.q_av[i] += sw.q[i] * Δt
+        sw.q_av[i] += sw.q[i] * dt
     end
     if update_h
         @batch per = thread minbatch = 2000 for i in sw.active_n
 
             q_src = sum_at(sw.q, links_at_node.src[i])
             q_dst = sum_at(sw.q, links_at_node.dst[i])
-            sw.volume[i] = sw.volume[i] + (q_src - q_dst + sw.inwater[i]) * Δt
+            sw.volume[i] = sw.volume[i] + (q_src - q_dst + sw.inwater[i]) * dt
 
             if sw.volume[i] < 0.0
                 sw.error[i] = sw.error[i] + abs(sw.volume[i])
                 sw.volume[i] = 0.0 # set volume to zero
             end
-            sw.volume[i] = max(sw.volume[i] + sw.inflow[i] * Δt, 0.0) # add external inflow
+            sw.volume[i] = max(sw.volume[i] + sw.inflow[i] * dt, 0.0) # add external inflow
 
             if !isnothing(sw.floodplain)
                 q_src = sum_at(sw.floodplain.q, links_at_node.src[i])
                 q_dst = sum_at(sw.floodplain.q, links_at_node.dst[i])
-                sw.floodplain.volume[i] = sw.floodplain.volume[i] + (q_src - q_dst) * Δt
+                sw.floodplain.volume[i] = sw.floodplain.volume[i] + (q_src - q_dst) * dt
                 # TODO check following approach:
                 # if floodplain volume negative, extract from river volume first
                 if sw.floodplain.volume[i] < 0.0
@@ -918,11 +918,11 @@ function shallowwater_river_update(sw::ShallowWaterRiver, network, Δt, doy, upd
                     sw.floodplain.h[i] = 0.0
                     sw.floodplain.volume[i] = 0.0
                 end
-                sw.floodplain.h_av[i] += sw.floodplain.h[i] * Δt
+                sw.floodplain.h_av[i] += sw.floodplain.h[i] * dt
             else
                 sw.h[i] = sw.volume[i] / (sw.dl[i] * sw.width[i])
             end
-            sw.h_av[i] += sw.h[i] * Δt
+            sw.h_av[i] += sw.h[i] * dt
         end
     end
 end
@@ -948,20 +948,20 @@ function update(sw::ShallowWaterRiver{T}, network, doy; update_h = true) where {
     sw.h_av .= 0.0
 
     t = T(0.0)
-    while t < sw.Δt
-        Δt = stable_timestep(sw)
-        if t + Δt > sw.Δt
-            Δt = sw.Δt - t
+    while t < sw.dt
+        dt = stable_timestep(sw)
+        if t + dt > sw.dt
+            dt = sw.dt - t
         end
-        shallowwater_river_update(sw, network, Δt, doy, update_h)
-        t = t + Δt
+        shallowwater_river_update(sw, network, dt, doy, update_h)
+        t = t + dt
     end
-    sw.q_av ./= sw.Δt
-    sw.h_av ./= sw.Δt
+    sw.q_av ./= sw.dt
+    sw.h_av ./= sw.dt
 
     if !isnothing(sw.floodplain)
-        sw.floodplain.q_av ./= sw.Δt
-        sw.floodplain.h_av ./= sw.Δt
+        sw.floodplain.q_av ./= sw.dt
+        sw.floodplain.h_av ./= sw.dt
         sw.q_channel_av .= sw.q_av
         sw.q_av .= sw.q_channel_av .+ sw.floodplain.q_av
     end
@@ -989,10 +989,10 @@ const dirs = (:yd, :xd, :xu, :yu)
     xwidth::Vector{T} | "m" | _ | "edge"                    # effective flow width x direction (floodplain)
     ywidth::Vector{T} | "m" | _ | "edge"                    # effective flow width y direction (floodplain)
     g::T | "m2 s-1" | 0 | "scalar"                          # acceleration due to gravity
-    θ::T | "-" | 0 | "scalar"                               # weighting factor (de Almeida et al., 2012)
-    α::T | "-" | 0 | "scalar"                               # stability coefficient (de Almeida et al., 2012)
+    theta::T | "-" | 0 | "scalar"                           # weighting factor (de Almeida et al., 2012)
+    alpha::T | "-" | 0 | "scalar"                           # stability coefficient (de Almeida et al., 2012)
     h_thresh::T | "m" | 0 | "scalar"                        # depth threshold for calculating flow
-    Δt::T | "s" | 0 | "none" | "none"                       # model time step [s]
+    dt::T | "s" | 0 | "none" | "none"                       # model time step [s]
     qy0::Vector{T} | "m3 s-1" | _ | "edge"                  # flow in y direction at previous time step
     qx0::Vector{T} | "m3 s-1" | _ | "edge"                  # flow in x direction at previous time step
     qx::Vector{T} | "m3 s-1" | _ | "edge"                   # flow in x direction
@@ -1025,7 +1025,7 @@ function initialize_shallowwater_land(
     inds_riv,
     river,
     waterbody,
-    Δt,
+    dt,
 )
     froude_limit = get(config.model, "froude_limit", true)::Bool # limit flow to subcritical according to Froude number
     alpha = get(config.model, "inertial_flow_alpha", 0.7)::Float64 # stability coefficient for model time step (0.2-0.7)
@@ -1106,10 +1106,10 @@ function initialize_shallowwater_land(
         xwidth = we_x,
         ywidth = we_y,
         g = 9.80665,
-        θ = theta,
-        α = alpha,
+        theta = theta,
+        alpha = alpha,
         h_thresh = h_thresh,
-        Δt = tosecond(Δt),
+        dt = tosecond(dt),
         qx0 = zeros(n + 1),
         qy0 = zeros(n + 1),
         qx = zeros(n + 1),
@@ -1137,30 +1137,30 @@ end
 
 Compute a stable timestep size for the local inertial approach, based on Bates et al. (2010).
 
-Δt = α * (Δx / sqrt(g max(h))
+dt = alpha * (Δx / sqrt(g max(h))
 """
 function stable_timestep(sw::ShallowWaterRiver{T})::T where {T}
-    Δtₘᵢₙ = T(Inf)
+    dt_min = T(Inf)
     @tturbo for i = 1:sw.n
-        Δt = sw.α * sw.dl[i] / sqrt(sw.g * sw.h[i])
-        Δtₘᵢₙ = Δt < Δtₘᵢₙ ? Δt : Δtₘᵢₙ
+        dt = sw.alpha * sw.dl[i] / sqrt(sw.g * sw.h[i])
+        dt_min = dt < dt_min ? dt : dt_min
     end
-    Δtₘᵢₙ = isinf(Δtₘᵢₙ) ? T(10.0) : Δtₘᵢₙ
-    return Δtₘᵢₙ
+    dt_min = isinf(dt_min) ? T(10.0) : dt_min
+    return dt_min
 end
 
 function stable_timestep(sw::ShallowWaterLand{T})::T where {T}
-    Δtₘᵢₙ = T(Inf)
+    dt_min = T(Inf)
     @tturbo for i = 1:sw.n
-        Δt = IfElse.ifelse(
+        dt = IfElse.ifelse(
             sw.rivercells[i] == 0,
-            sw.α * min(sw.xl[i], sw.yl[i]) / sqrt(sw.g * sw.h[i]),
+            sw.alpha * min(sw.xl[i], sw.yl[i]) / sqrt(sw.g * sw.h[i]),
             T(Inf),
         )
-        Δtₘᵢₙ = Δt < Δtₘᵢₙ ? Δt : Δtₘᵢₙ
+        dt_min = dt < dt_min ? dt : dt_min
     end
-    Δtₘᵢₙ = isinf(Δtₘᵢₙ) ? T(10.0) : Δtₘᵢₙ
-    return Δtₘᵢₙ
+    dt_min = isinf(dt_min) ? T(10.0) : dt_min
+    return dt_min
 end
 
 function update(
@@ -1188,27 +1188,27 @@ function update(
     sw.h_av .= 0.0
 
     t = T(0.0)
-    while t < swr.Δt
-        Δt_river = stable_timestep(swr)
-        Δt_land = stable_timestep(sw)
-        Δt = min(Δt_river, Δt_land)
-        if t + Δt > swr.Δt
-            Δt = swr.Δt - t
+    while t < swr.dt
+        dt_river = stable_timestep(swr)
+        dt_land = stable_timestep(sw)
+        dt = min(dt_river, dt_land)
+        if t + dt > swr.dt
+            dt = swr.dt - t
         end
-        shallowwater_river_update(swr, network.river, Δt, doy, update_h)
-        shallowwater_update(sw, swr, network, Δt)
-        t = t + Δt
+        shallowwater_river_update(swr, network.river, dt, doy, update_h)
+        shallowwater_update(sw, swr, network, dt)
+        t = t + dt
     end
-    swr.q_av ./= swr.Δt
-    swr.h_av ./= swr.Δt
-    sw.h_av ./= sw.Δt
+    swr.q_av ./= swr.dt
+    swr.h_av ./= swr.dt
+    sw.h_av ./= sw.dt
 end
 
 function shallowwater_update(
     sw::ShallowWaterLand{T},
     swr::ShallowWaterRiver{T},
     network,
-    Δt,
+    dt,
 ) where {T}
 
     indices = network.land.staggered_indices
@@ -1230,27 +1230,27 @@ function shallowwater_update(
         # for flow in x dir) and floodplain flow is not calculated.
         if xu <= sw.n && sw.ywidth[i] != T(0.0)
 
-            η_x = sw.z[i] + sw.h[i]
-            η_xu = sw.z[xu] + sw.h[xu]
-            η_max = max(η_x, η_xu)
-            hf = (η_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
                 sw.qx[i] = local_inertial_flow(
-                    sw.θ,
+                    sw.theta,
                     sw.qx0[i],
                     sw.qx0[xd],
                     sw.qx0[xu],
-                    η_x,
-                    η_xu,
+                    zs_x,
+                    zs_xu,
                     hf,
                     sw.ywidth[i],
                     length,
                     sw.mannings_n_sq[i],
                     sw.g,
                     sw.froude_limit,
-                    Δt,
+                    dt,
                 )
                 # limit qx in case water is not available
                 if sw.h[i] <= T(0.0)
@@ -1270,27 +1270,27 @@ function shallowwater_update(
         # for flow in y dir) and floodplain flow is not calculated.
         if yu <= sw.n && sw.xwidth[i] != T(0.0)
 
-            η_y = sw.z[i] + sw.h[i]
-            η_yu = sw.z[yu] + sw.h[yu]
-            η_max = max(η_y, η_yu)
-            hf = (η_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
                 sw.qy[i] = local_inertial_flow(
-                    sw.θ,
+                    sw.theta,
                     sw.qy0[i],
                     sw.qy0[yd],
                     sw.qy0[yu],
-                    η_y,
-                    η_yu,
+                    zs_y,
+                    zs_yu,
                     hf,
                     sw.xwidth[i],
                     length,
                     sw.mannings_n_sq[i],
                     sw.g,
                     sw.froude_limit,
-                    Δt,
+                    dt,
                 )
                 # limit qy in case water is not available
                 if sw.h[i] <= T(0.0)
@@ -1326,7 +1326,7 @@ function shallowwater_update(
                         sw.qx[i] + sw.qy[yd] - sw.qy[i] +
                         swr.inflow[inds_riv[i]] +
                         sw.runoff[i]
-                    ) * Δt
+                    ) * dt
                 if sw.volume[i] < T(0.0)
                     sw.error[i] = sw.error[i] + abs(sw.volume[i])
                     sw.volume[i] = T(0.0) # set volume to zero
@@ -1345,18 +1345,18 @@ function shallowwater_update(
                     sw.h[i] = T(0.0)
                     swr.volume[inds_riv[i]] = sw.volume[i]
                 end
-                swr.h_av[inds_riv[i]] += swr.h[inds_riv[i]] * Δt
+                swr.h_av[inds_riv[i]] += swr.h[inds_riv[i]] * dt
             end
         else
             sw.volume[i] +=
-                (sw.qx[xd] - sw.qx[i] + sw.qy[yd] - sw.qy[i] + sw.runoff[i]) * Δt
+                (sw.qx[xd] - sw.qx[i] + sw.qy[yd] - sw.qy[i] + sw.runoff[i]) * dt
             if sw.volume[i] < T(0.0)
                 sw.error[i] = sw.error[i] + abs(sw.volume[i])
                 sw.volume[i] = T(0.0) # set volume to zero
             end
             sw.h[i] = sw.volume[i] / (sw.xl[i] * sw.yl[i])
         end
-        sw.h_av[i] += sw.h[i] * Δt
+        sw.h_av[i] += sw.h[i] * dt
     end
 end
 
@@ -1433,8 +1433,8 @@ Compute floodplain flow area based on flow depth `h` and floodplain `depth`, `ar
 `width` of a floodplain profile.
 """
 function flow_area(width, area, depth, h)
-    Δh = h - depth  # depth at i1
-    area = area + (width * Δh) # area at i1, width at i2
+    dh = h - depth  # depth at i1
+    area = area + (width * dh) # area at i1, width at i2
     return area
 end
 
@@ -1445,8 +1445,8 @@ Compute floodplain wetted perimeter based on flow depth `h` and floodplain `dept
 wetted perimeter `p` of a floodplain profile.
 """
 function wetted_perimeter(p, depth, h)
-    Δh = h - depth # depth at i1
-    p = p + (2.0 * Δh) # p at i1
+    dh = h - depth # depth at i1
+    p = p + (2.0 * dh) # p at i1
     return p
 end
 
@@ -1454,8 +1454,8 @@ end
 function flood_depth(profile::FloodPlainProfile{T}, flood_volume, dl, i::Int)::T where {T}
     i1, i2 = interpolation_indices(flood_volume, @view profile.volume[:, i])
     ΔA = (flood_volume - profile.volume[i1, i]) / dl
-    Δh = ΔA / profile.width[i2, i]
-    flood_depth = profile.depth[i1] + Δh
+    dh = ΔA / profile.width[i2, i]
+    flood_depth = profile.depth[i1] + dh
     return flood_depth
 end
 
@@ -1625,7 +1625,7 @@ function set_river_inwater(model::Model{N,L,V,R,W,T}, ssf_toriver) where {N,L,V<
                 network.land.xl[inds] *
                 network.land.yl[inds] *
                 0.001
-            ) / vertical.Δt
+            ) / vertical.dt
         )
     )
 end
@@ -1658,7 +1658,7 @@ function set_land_inwater(model::Model{N,L,V,R,W,T}) where {N,L,V,R,W,T<:SbmGwfM
 
     lateral.land.inwater .=
         (vertical.net_runoff .* network.land.xl .* network.land.yl .* 0.001) ./
-        lateral.land.Δt .+ drainflux
+        lateral.land.dt .+ drainflux
 end
 
 """
@@ -1670,7 +1670,7 @@ function set_land_inwater(model::Model{N,L,V,R,W,T}) where {N,L,V,R,W,T<:SbmMode
     @unpack lateral, vertical, network = model
     lateral.land.inwater .=
         (vertical.net_runoff .* network.land.xl .* network.land.yl .* 0.001) ./
-        lateral.land.Δt
+        lateral.land.dt
 end
 
 """
@@ -1681,7 +1681,7 @@ Set `inwater` of the lateral land component, based on `runoff` of the `vertical`
 function set_land_inwater(model)
     @unpack lateral, vertical, network = model
     lateral.land.inwater .=
-        (vertical.runoff .* network.land.xl .* network.land.yl .* 0.001) ./ lateral.land.Δt
+        (vertical.runoff .* network.land.xl .* network.land.yl .* 0.001) ./ lateral.land.dt
 end
 
 # Computation of inflow from the lateral components `land` and `subsurface` to water bodies
@@ -1786,7 +1786,7 @@ function surface_routing(model; ssf_toriver = 0.0)
     # run river flow
     set_river_inwater(model, ssf_toriver)
     set_inflow_waterbody(model)
-    update(lateral.river, network.river, julian_day(clock.time - clock.Δt))
+    update(lateral.river, network.river, julian_day(clock.time - clock.dt))
 end
 
 """
@@ -1806,14 +1806,14 @@ function surface_routing(
     @unpack lateral, vertical, network, clock = model
 
     @. lateral.land.runoff = (
-        (vertical.net_runoff / 1000.0) * (network.land.xl * network.land.yl) / vertical.Δt +
+        (vertical.net_runoff / 1000.0) * (network.land.xl * network.land.yl) / vertical.dt +
         ssf_toriver +
         # net_runoff_river
         (
             (vertical.net_runoff_river * network.land.xl * network.land.yl * 0.001) /
-            vertical.Δt
+            vertical.dt
         )
     )
     set_inflow_waterbody(model)
-    update(lateral.land, lateral.river, network, julian_day(clock.time - clock.Δt))
+    update(lateral.land, lateral.river, network, julian_day(clock.time - clock.dt))
 end
diff --git a/src/groundwater/aquifer.jl b/src/groundwater/aquifer.jl
index 43aa05438..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"
-        ϕᵢ = aquifer.head[i]
-        ϕⱼ = aquifer.head[j]
-        if ϕᵢ >= ϕⱼ
+        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]
-            Δϕ = 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 * Δϕ
+            Q[i] -= cond * delta_head
         end
     end
     return Q
@@ -322,7 +322,7 @@ The following criterion can be found in Chu & Willis (1984)
 Δt * k * H / (Δx * Δy * S) <= 1/4
 """
 function stable_timestep(aquifer, conductivity_profile::String)
-    Δtₘᵢₙ = Inf
+    dt_min = Inf
     for i in eachindex(aquifer.head)
         if conductivity_profile == "exponential"
             zi = aquifer.top[i] - aquifer.head[i]
@@ -334,23 +334,23 @@ function stable_timestep(aquifer, conductivity_profile::String)
             value = aquifer.k[i] * saturated_thickness(aquifer, i)
         end
 
-        Δt = aquifer.area[i] * storativity(aquifer)[i] / value
-        Δtₘᵢₙ = Δt < Δtₘᵢₙ ? Δt : Δtₘᵢₙ
+        dt = aquifer.area[i] * storativity(aquifer)[i] / value
+        dt_min = dt < dt_min ? dt : dt_min
     end
-    return 0.25 * Δtₘᵢₙ
+    return 0.25 * dt_min
 end
 
 minimum_head(aquifer::ConfinedAquifer) = aquifer.head
 minimum_head(aquifer::UnconfinedAquifer) = max.(aquifer.head, aquifer.bottom)
 
 
-function update(gwf, Q, Δt, conductivity_profile)
+function update(gwf, Q, dt, conductivity_profile)
     Q .= 0.0  # TODO: Probably remove this when linking with other components
     flux!(Q, gwf.aquifer, gwf.connectivity, conductivity_profile)
     for boundary in gwf.boundaries
         flux!(Q, boundary, gwf.aquifer)
     end
-    gwf.aquifer.head .+= (Q ./ gwf.aquifer.area .* Δt ./ storativity(gwf.aquifer))
+    gwf.aquifer.head .+= (Q ./ gwf.aquifer.area .* dt ./ storativity(gwf.aquifer))
     # Set constant head (dirichlet) boundaries
     gwf.aquifer.head[gwf.constanthead.index] .= gwf.constanthead.head
     # Make sure no heads ends up below an unconfined aquifer bottom
diff --git a/src/groundwater/boundary_conditions.jl b/src/groundwater/boundary_conditions.jl
index f38e6d321..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)
-        ϕ = aquifer.head[index]
+        head = aquifer.head[index]
         stage = river.stage[i]
-        if stage > ϕ
+        if stage > head
             cond = river.infiltration_conductance[i]
-            Δϕ = min(stage - river.bottom[i], stage - ϕ)
+            delta_head = min(stage - river.bottom[i], stage - head)
         else
             cond = river.exfiltration_conductance[i]
-            Δϕ = stage - ϕ
+            delta_head = stage - head
         end
-        river.flux[i] = check_flux(cond * Δϕ, 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]
-        Δϕ = min(0, drainage.elevation[i] - aquifer.head[index])
-        drainage.flux[i] = check_flux(cond * Δϕ, 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]
-        Δϕ = headboundary.head[i] - aquifer.head[index]
-        headboundary.flux[i] = check_flux(cond * Δϕ, 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/groundwater/connectivity.jl b/src/groundwater/connectivity.jl
index 6b496066b..90a9bf6f2 100644
--- a/src/groundwater/connectivity.jl
+++ b/src/groundwater/connectivity.jl
@@ -39,19 +39,19 @@ connections(C::Connectivity, id::Int) = C.colptr[id]:(C.colptr[id+1]-1)
 
 
 """
-    connection_geometry(I, J, Δx, Δy)
+    connection_geometry(I, J, dx, dy)
 
 Compute geometrical properties of connections for structured input.
 """
-function connection_geometry(I, J, Δx, Δy)
+function connection_geometry(I, J, dx, dy)
     if I[1] != J[1]  # connection in y
-        length1 = 0.5 * Δy[I[1]]
-        length2 = 0.5 * Δy[J[1]]
-        width = Δx[I[2]]
+        length1 = 0.5 * dy[I[1]]
+        length2 = 0.5 * dy[J[1]]
+        width = dx[I[2]]
     elseif I[2] != J[2]  # connection in x
-        length1 = 0.5 * Δx[I[2]]
-        length2 = 0.5 * Δx[J[2]]
-        width = Δy[I[1]]
+        length1 = 0.5 * dx[I[2]]
+        length2 = 0.5 * dx[J[2]]
+        width = dy[I[1]]
     else
         # TODO: more specific exception? --> Martijn
         error("Inconsistent index")
@@ -69,7 +69,7 @@ const neighbors = (
 )
 
 # Constructor for the Connectivity structure for structured input
-function Connectivity(indices, reverse_indices, Δx::Vector{T}, Δy::Vector{T}) where {T}
+function Connectivity(indices, reverse_indices, dx::Vector{T}, dy::Vector{T}) where {T}
     # indices: These map from the 1D internal domain to the 2D external domain.
     # reverse_indices: from the 2D external domain to the 1D internal domain,
     # providing an Int which can be used as a linear index
@@ -92,7 +92,7 @@ function Connectivity(indices, reverse_indices, Δx::Vector{T}, Δy::Vector{T})
             J = I + neighbor
             if (1 <= J[1] <= nrow) && (1 <= J[2] <= ncol && reverse_indices[J] != 0) # Check if it's inbounds and neighbor is active
                 rowval[i] = reverse_indices[J]
-                length1[i], length2[i], width[i] = connection_geometry(I, J, Δx, Δy)
+                length1[i], length2[i], width[i] = connection_geometry(I, J, dx, dy)
                 i += 1
             end
         end
diff --git a/src/hbv.jl b/src/hbv.jl
index dd22bb309..038514622 100644
--- a/src/hbv.jl
+++ b/src/hbv.jl
@@ -1,15 +1,15 @@
 @get_units @exchange @grid_type @grid_location @with_kw struct HBV{T}
-    Δt::T | "s" | 0 | "none" | "none"   # Model time step [s]
+    dt::T | "s" | 0 | "none" | "none"   # Model time step [s]
     n::Int | "-" | 0 | "none" | "none"  # Number of cells
     fc::Vector{T} | "mm"                # Field capacity [mm]
     betaseepage::Vector{T} | "-"        # Exponent in soil runoff generation equation [-]
     lp::Vector{T} | "-"                 # Fraction of field capacity below which actual evaporation=potential evaporation [-]
     threshold::Vector{T} | "mm"         # Threshold soilwater storage above which AE=PE [mm]
-    k4::Vector{T} | "Δt-1"              # Recession constant baseflow [Δt⁻¹]
-    kquickflow::Vector{T} | "Δt-1"      # Recession constant upper reservoir [Δt⁻¹]
+    k4::Vector{T} | "dt-1"              # Recession constant baseflow [Δt⁻¹]
+    kquickflow::Vector{T} | "dt-1"      # Recession constant upper reservoir [Δt⁻¹]
     suz::Vector{T} | "mm"               # Level over which k0 is used [mm]
-    k0::Vector{T} | "Δt-1"              # Recession constant upper reservoir [Δt⁻¹]
-    khq::Vector{T} | "Δt-1"             # Recession rate at flow hq [Δt⁻¹]
+    k0::Vector{T} | "dt-1"              # Recession constant upper reservoir [Δt⁻¹]
+    khq::Vector{T} | "dt-1"             # Recession rate at flow hq [Δt⁻¹]
     hq::Vector{T}                       # High flow rate hq for which recession rate of upper reservoir is known [mm Δt⁻¹]
     alphanl::Vector{T}                  # Measure of non-linearity of upper reservoir
     perc::Vector{T}                     # Percolation from upper to lower zone [mm Δt⁻¹]
@@ -25,11 +25,11 @@
     tti::Vector{T} | "ᵒC"               # Critical temperature for snowmelt and refreezing [ᵒC]
     tt::Vector{T} | "ᵒC"                # Defines interval in which precipitation falls as rainfall and snowfall [ᵒC]
     ttm::Vector{T} | "ᵒC"               # Threshold temperature for snowmelt [ᵒC]
-    cfmax::Vector{T} | "mm ᵒC-1 Δt-1"   # Meltconstant in temperature-index [-]
+    cfmax::Vector{T} | "mm ᵒC-1 dt-1"   # Meltconstant in temperature-index [-]
     whc::Vector{T} | "-"                # Fraction of snow volume that can store water [-]
     g_tt::Vector{T} | "ᵒC"              # Threshold temperature for snowfall above glacier [ᵒC]
-    g_cfmax::Vector{T} | "mm ᵒC-1 Δt-1" # Degree-day factor [mm ᵒC⁻¹ Δt⁻¹] for glacier
-    g_sifrac::Vector{T} | "Δt-1"        # Fraction of the snowpack on top of the glacier converted into ice [Δt⁻¹]
+    g_cfmax::Vector{T} | "mm ᵒC-1 dt-1" # Degree-day factor [mm ᵒC⁻¹ Δt⁻¹] for glacier
+    g_sifrac::Vector{T} | "dt-1"        # Fraction of the snowpack on top of the glacier converted into ice [Δt⁻¹]
     glacierstore::Vector{T} | "mm"      # Water within the glacier [mm]
     glacierfrac::Vector{T} | "-"        # Fraction covered by a glacier [-]
     precipitation::Vector{T}            # Precipitation [mm Δt⁻¹]
@@ -136,7 +136,7 @@ function update_after_snow(hbv::HBV, config)
                 hbv.g_tt[i],
                 hbv.g_cfmax[i],
                 hbv.g_sifrac[i],
-                Second(hbv.Δt),
+                Second(hbv.dt),
             )
             # Convert to mm per grid cell and add to snowmelt
             glaciermelt = glaciermelt * hbv.glacierfrac[i]
diff --git a/src/hbv_model.jl b/src/hbv_model.jl
index 4577167a9..aab5cbbd7 100644
--- a/src/hbv_model.jl
+++ b/src/hbv_model.jl
@@ -10,7 +10,7 @@ function initialize_hbv_model(config::Config)
 
     reader = prepare_reader(config)
     clock = Clock(config, reader)
-    Δt = clock.Δt
+    dt = clock.dt
 
     do_reservoirs = get(config.model, "reservoirs", false)::Bool
     do_lakes = get(config.model, "lakes", false)::Bool
@@ -36,7 +36,7 @@ function initialize_hbv_model(config::Config)
             sel = inds,
             defaults = 3.75653,
             type = Float,
-        ) .* (Δt / basetimestep)
+        ) .* (dt / basetimestep)
     tt = ncread(nc, config, "vertical.tt"; sel = inds, defaults = -1.41934, type = Float)
     tti = ncread(nc, config, "vertical.tti"; sel = inds, defaults = 1.0, type = Float)
     ttm = ncread(nc, config, "vertical.ttm"; sel = inds, defaults = -1.41934, type = Float)
@@ -60,7 +60,7 @@ function initialize_hbv_model(config::Config)
             defaults = 3.0,
             type = Float,
             fill = 0.0,
-        ) .* (Δt / basetimestep)
+        ) .* (dt / basetimestep)
     g_sifrac =
         ncread(
             nc,
@@ -70,7 +70,7 @@ function initialize_hbv_model(config::Config)
             defaults = 0.001,
             type = Float,
             fill = 0.0,
-        ) .* (Δt / basetimestep)
+        ) .* (dt / basetimestep)
     glacierfrac = ncread(
         nc,
         config,
@@ -95,7 +95,7 @@ function initialize_hbv_model(config::Config)
     lp = ncread(nc, config, "vertical.lp"; sel = inds, defaults = 0.53, type = Float)
     k4 =
         ncread(nc, config, "vertical.k4"; sel = inds, defaults = 0.02307, type = Float) .*
-        (Δt / basetimestep)
+        (dt / basetimestep)
     kquickflow =
         ncread(
             nc,
@@ -104,29 +104,29 @@ function initialize_hbv_model(config::Config)
             sel = inds,
             defaults = 0.09880,
             type = Float,
-        ) .* (Δt / basetimestep)
+        ) .* (dt / basetimestep)
     suz = ncread(nc, config, "vertical.suz"; sel = inds, defaults = 100.0, type = Float)
     k0 =
         ncread(nc, config, "vertical.k0"; sel = inds, defaults = 0.30, type = Float) .*
-        (Δt / basetimestep)
+        (dt / basetimestep)
     khq =
         ncread(nc, config, "vertical.khq"; sel = inds, defaults = 0.09880, type = Float) .*
-        (Δt / basetimestep)
+        (dt / basetimestep)
     hq =
         ncread(nc, config, "vertical.hq"; sel = inds, defaults = 3.27, type = Float) .*
-        (Δt / basetimestep)
+        (dt / basetimestep)
     alphanl =
         ncread(nc, config, "vertical.alphanl"; sel = inds, defaults = 1.1, type = Float)
     perc =
         ncread(nc, config, "vertical.perc"; sel = inds, defaults = 0.4, type = Float) .*
-        (Δt / basetimestep)
+        (dt / basetimestep)
     cfr = ncread(nc, config, "vertical.cfr"; sel = inds, defaults = 0.05, type = Float)
     pcorr = ncread(nc, config, "vertical.pcorr"; sel = inds, defaults = 1.0, type = Float)
     rfcf = ncread(nc, config, "vertical.rfcf"; sel = inds, defaults = 1.0, type = Float)
     sfcf = ncread(nc, config, "vertical.sfcf"; sel = inds, defaults = 1.0, type = Float)
     cflux =
         ncread(nc, config, "vertical.cflux"; sel = inds, defaults = 2.0, type = Float) .*
-        (Δt / basetimestep)
+        (dt / basetimestep)
     icf = ncread(nc, config, "vertical.icf"; sel = inds, defaults = 2.0, type = Float)
     cevpf = ncread(nc, config, "vertical.cevpf"; sel = inds, defaults = 1.0, type = Float)
     epf = ncread(nc, config, "vertical.epf"; sel = inds, defaults = 1.0, type = Float)
@@ -145,7 +145,7 @@ function initialize_hbv_model(config::Config)
     threshold = fc .* lp
 
     hbv = HBV{Float}(
-        Δt = Float(tosecond(Δt)),
+        dt = Float(tosecond(dt)),
         n = n,
         fc = fc,
         betaseepage = betaseepage,
@@ -225,7 +225,7 @@ function initialize_hbv_model(config::Config)
     pits = zeros(Bool, modelsize_2d)
     if do_reservoirs
         reservoirs, resindex, reservoir, pits =
-            initialize_simple_reservoir(config, nc, inds_riv, nriv, pits, tosecond(Δt))
+            initialize_simple_reservoir(config, nc, inds_riv, nriv, pits, tosecond(dt))
     else
         reservoir = ()
         reservoirs = nothing
@@ -235,7 +235,7 @@ function initialize_hbv_model(config::Config)
     # lakes
     if do_lakes
         lakes, lakeindex, lake, pits =
-            initialize_lake(config, nc, inds_riv, nriv, pits, tosecond(Δt))
+            initialize_lake(config, nc, inds_riv, nriv, pits, tosecond(dt))
     else
         lake = ()
         lakes = nothing
@@ -249,9 +249,9 @@ function initialize_hbv_model(config::Config)
         ldd = set_pit_ldd(pits_2d, ldd, inds)
     end
 
-    βₗ =
+    landslope =
         ncread(nc, config, "lateral.land.slope"; optional = false, sel = inds, type = Float)
-    clamp!(βₗ, 0.00001, Inf)
+    clamp!(landslope, 0.00001, Inf)
 
     dl = map(detdrainlength, ldd, xl, yl)
     dw = (xl .* yl) ./ dl
@@ -259,12 +259,12 @@ function initialize_hbv_model(config::Config)
         nc,
         config,
         inds;
-        sl = βₗ,
+        sl = landslope,
         dl = dl,
         width = map(det_surfacewidth, dw, riverwidth, river),
         iterate = kinwave_it,
         tstep = kw_land_tstep,
-        Δt = Δt,
+        dt = dt,
     )
 
     graph = flowgraph(ldd, inds, pcr_dir)
@@ -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, βₗ, n)
+    frac_toriver = fraction_runoff_toriver(graph, ldd, index_river, landslope, n)
 
     rf = initialize_surfaceflow_river(
         nc,
@@ -296,7 +296,7 @@ function initialize_hbv_model(config::Config)
         lake = lakes,
         iterate = kinwave_it,
         tstep = kw_river_tstep,
-        Δt = Δt,
+        dt = dt,
     )
 
     # setup subdomains for the land and river kinematic wave domain, if nthreads = 1
diff --git a/src/horizontal_process.jl b/src/horizontal_process.jl
index 380edd39e..2d10333ce 100644
--- a/src/horizontal_process.jl
+++ b/src/horizontal_process.jl
@@ -24,19 +24,19 @@ function flowgraph(ldd::AbstractVector, inds::AbstractVector, pcr_dir::AbstractV
 end
 
 "Kinematic wave flow rate for a single cell and timestep"
-function kinematic_wave(Qin, Qold, q, α, β, Δt, Δx)
-    ϵ = 1.0e-12
+function kinematic_wave(Qin, Qold, q, alpha, beta, dt, dx)
+    epsilon = 1.0e-12
     max_iters = 3000
 
     if Qin + Qold + q ≈ 0.0
         return 0.0
     else
         # common terms
-        ab_pQ = α * β * pow(((Qold + Qin) / 2.0), (β - 1.0))
-        Δtx = Δt / Δx
-        C = Δtx * Qin + α * pow(Qold, β) + Δt * q
+        ab_pQ = alpha * beta * pow(((Qold + Qin) / 2.0), (beta - 1.0))
+        dtx = dt / dx
+        C = dtx * Qin + alpha * pow(Qold, beta) + dt * q
 
-        Qkx = (Δtx * Qin + Qold * ab_pQ + Δt * q) / (Δtx + ab_pQ)
+        Qkx = (dtx * Qin + Qold * ab_pQ + dt * q) / (dtx + ab_pQ)
         if isnan(Qkx)
             Qkx = 0.0
         end
@@ -44,11 +44,11 @@ function kinematic_wave(Qin, Qold, q, α, β, Δt, Δx)
         count = 1
 
         while true
-            fQkx = Δtx * Qkx + α * pow(Qkx, β) - C
-            dfQkx = Δtx + α * β * pow(Qkx, (β - 1.0))
+            fQkx = dtx * Qkx + alpha * pow(Qkx, beta) - C
+            dfQkx = dtx + alpha * beta * pow(Qkx, (beta - 1.0))
             Qkx = Qkx - fQkx / dfQkx
             Qkx = max(Qkx, 1.0e-30)
-            if (abs(fQkx) <= ϵ) || (count >= max_iters)
+            if (abs(fQkx) <= epsilon) || (count >= max_iters)
                 break
             end
             count += 1
@@ -59,38 +59,38 @@ function kinematic_wave(Qin, Qold, q, α, β, Δt, Δx)
 end
 
 "Kinematic wave flow rate over the whole network for a single timestep"
-function kin_wave!(Q, graph, toposort, Qold, q, α, β, DCL, Δt)
+function kin_wave!(Q, graph, toposort, Qold, q, alpha, beta, DCL, dt)
     for v in toposort
         upstream_nodes = inneighbors(graph, v)
         Qin = sum_at(Q, upstream_nodes)
-        Q[v] = kinematic_wave(Qin, Qold[v], q, α[v], β, Δt, DCL[v])
+        Q[v] = kinematic_wave(Qin, Qold[v], q, alpha[v], beta, dt, DCL[v])
     end
     return Q
 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₀, β, 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₀ * β) + exp(-f * d)) / -f
+        zi = log((f * ssf) / (dw * kh_0 * slope) + exp(-f * d)) / -f
     elseif ksat_profile == "exponential_constant"
-        ssf_constant = kh₀ * β * 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₀ * β) + 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₀ * β * 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₀, β, θₑ, 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₀ * exp(-f * zi) * β) / θₑ
+        Cn = (kh_0 * exp(-f * zi) * slope) / theta_e
     elseif ksat_profile == "exponential_constant"
-        Cn_const = (kh₀ * exp(-f * z_exp) * β) / θₑ
+        Cn_const = (kh_0 * exp(-f * z_exp) * slope) / theta_e
         if zi < z_exp
-            Cn = (kh₀ * exp(-f * zi) * β) / θₑ + Cn_const
+            Cn = (kh_0 * exp(-f * zi) * slope) / theta_e + Cn_const
         else
             Cn = Cn_const
         end
@@ -99,58 +99,58 @@ function ssf_celerity(zi, kh₀, β, θₑ, f, z_exp, ksat_profile)
 end
 
 """
-    kinematic_wave_ssf(ssfin, ssfₜ₋₁, ziₜ₋₁, r, kh₀, β, θₑ, f, d, Δt, Δx, 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₀`, controllled by parameter `f`,
+decline of hydraulic conductivity at the soil surface `kh_0`, controllled by parameter `f`,
 is assumed. The hydraulic conductivity profile `ksat_profile` is either `exponential` or
 `exponential_constant`, with `z_exp` the depth from the soil surface for which the
-exponential decline of `kh₀` is valid.
+exponential decline of `kh_0` is valid.
 
 Returns lateral subsurface flow `ssf`, water table depth `zi` and exfiltration rate
 `exfilt`.
 """
 function kinematic_wave_ssf(
     ssfin,
-    ssfₜ₋₁,
-    ziₜ₋₁,
+    ssf_prev,
+    zi_prev,
     r,
-    kh₀,
-    β,
-    θₑ,
+    kh_0,
+    slope,
+    theta_e,
     f,
     d,
-    Δt,
-    Δx,
+    dt,
+    dx,
     dw,
     ssfmax,
     z_exp,
     ksat_profile,
 )
 
-    ϵ = 1.0e-12
+    epsilon = 1.0e-12
     max_iters = 3000
 
-    if ssfin + ssfₜ₋₁ ≈ 0.0 && r <= 0.0
+    if ssfin + ssf_prev ≈ 0.0 && r <= 0.0
         return 0.0, d, 0.0
     else
         # initial estimate
-        ssf = (ssfₜ₋₁ + ssfin) / 2.0
+        ssf = (ssf_prev + ssfin) / 2.0
         count = 1
 
         # Estimate zi on the basis of the relation between subsurface flow and zi
-        zi = ssf_water_table_depth(ssf, kh₀, β, 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₀, β, θₑ, 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ₜ₋₁ can be replaced with (1./Cn)*ssf, and thus celerity and lateral flow rate ssf are then in line
-        c = (Δt / Δx) * ssfin + (1.0 / Cn) * ssf + (r - (ziₜ₋₁ - zi) * θₑ * dw)
+        # then (1./Cn)*ssf_prev can be replaced with (1./Cn)*ssf, and thus celerity and lateral flow rate ssf are then in line
+        c = (dt / dx) * ssfin + (1.0 / Cn) * ssf + (r - (zi_prev - zi) * theta_e * dw)
 
         # Continuity equation of which solution should be zero
-        fQ = (Δt / Δx) * ssf + (1.0 / Cn) * ssf - c
+        fQ = (dt / dx) * ssf + (1.0 / Cn) * ssf - c
         # Derivative of the continuity equation w.r.t. Q_out for iteration 1
-        dfQ = (Δt / Δx) + 1.0 / Cn
+        dfQ = (dt / dx) + 1.0 / Cn
         # Update lateral outflow estimate ssf (Q_out) for iteration 1
         ssf = ssf - (fQ / dfQ)
         if isnan(ssf)
@@ -161,25 +161,25 @@ 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₀, β, 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₀, β, θₑ, 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ₜ₋₁ can be replaced with (1./Cn)*ssf, and thus celerity and lateral flow rate ssf are then in line
-            c = (Δt / Δx) * ssfin + (1.0 / Cn) * ssf + (r - (ziₜ₋₁ - zi) * θₑ * dw)
+            # then (1./Cn)*ssf_prev can be replaced with (1./Cn)*ssf, and thus celerity and lateral flow rate ssf are then in line
+            c = (dt / dx) * ssfin + (1.0 / Cn) * ssf + (r - (zi_prev - zi) * theta_e * dw)
 
             # Continuity equation of which solution should be zero
-            fQ = (Δt / Δx) * ssf + (1.0 / Cn) * ssf - c
+            fQ = (dt / dx) * ssf + (1.0 / Cn) * ssf - c
             # Derivative of the continuity equation w.r.t. Q_out for iteration m+1
-            dfQ = (Δt / Δx) + 1.0 / Cn
+            dfQ = (dt / dx) + 1.0 / Cn
             # Update lateral outflow estimate ssf (Q_out) for iteration m+1
             ssf = ssf - (fQ / dfQ)
             if isnan(ssf)
                 ssf = 0.0
             end
             ssf = max(ssf, 1.0e-30)
-            if (abs(fQ) <= ϵ) || (count >= max_iters)
+            if (abs(fQ) <= epsilon) || (count >= max_iters)
                 break
             end
             count += 1
@@ -188,9 +188,9 @@ function kinematic_wave_ssf(
         # Constrain the lateral flow rate ssf
         ssf = min(ssf, (ssfmax * dw))
         # On the basis of the lateral flow rate, estimate the amount of groundwater level above surface (saturation excess conditions), then rest = negative
-        rest = ziₜ₋₁ - (ssfin * Δt + r * Δx - ssf * Δt) / (dw * Δx) / θₑ
+        rest = zi_prev - (ssfin * dt + r * dx - ssf * dt) / (dw * dx) / theta_e
         # In case the groundwater level lies above surface (saturation excess conditions, rest = negative), calculate the exfiltration rate and set groundwater back to zero.
-        exfilt = min(rest, 0.0) * -θₑ
+        exfilt = min(rest, 0.0) * -theta_e
         zi = clamp(zi, 0.0, d)
 
         return ssf, zi, exfilt
@@ -199,7 +199,7 @@ function kinematic_wave_ssf(
 end
 
 """
-    kinematic_wave_ssf(ssfin, ssfₜ₋₁, ziₜ₋₁, r, kh, β, θₑ, d, Δt, Δx, 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,25 +207,25 @@ 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ₜ₋₁, ziₜ₋₁, r, kh, β, θₑ, d, Δt, Δx, dw, ssfmax)
+function kinematic_wave_ssf(ssfin, ssf_prev, zi_prev, r, kh, slope, theta_e, d, dt, dx, dw, ssfmax)
 
-    ϵ = 1.0e-12
+    epsilon = 1.0e-12
     max_iters = 3000
 
-    if ssfin + ssfₜ₋₁ ≈ 0.0 && r <= 0.0
+    if ssfin + ssf_prev ≈ 0.0 && r <= 0.0
         return 0.0, d, 0.0
     else
         # initial estimate
-        ssf = (ssfₜ₋₁ + ssfin) / 2.0
+        ssf = (ssf_prev + ssfin) / 2.0
         count = 1
         # celerity (Cn)
-        Cn = (β * kh) / θₑ
+        Cn = (slope * kh) / theta_e
         # constant term of the continuity equation for Newton-Raphson
-        c = (Δt / Δx) * ssfin + (1.0 / Cn) * ssfₜ₋₁ + r
+        c = (dt / dx) * ssfin + (1.0 / Cn) * ssf_prev + r
         # continuity equation of which solution should be zero
-        fQ = (Δt / Δx) * ssf + (1.0 / Cn) * ssf - c
+        fQ = (dt / dx) * ssf + (1.0 / Cn) * ssf - c
         # Derivative of the continuity equation
-        dfQ = (Δt / Δx) + 1.0 / Cn
+        dfQ = (dt / dx) + 1.0 / Cn
         # Update lateral subsurface flow estimate ssf
         ssf = ssf - (fQ / dfQ)
         if isnan(ssf)
@@ -235,14 +235,14 @@ function kinematic_wave_ssf(ssfin, ssfₜ₋₁, ziₜ₋₁, r, kh, β, θₑ,
 
         # Start while loop of Newton-Raphson iteration
         while true
-            fQ = (Δt / Δx) * ssf + (1.0 / Cn) * ssf - c
-            dfQ = (Δt / Δx) + 1.0 / Cn
+            fQ = (dt / dx) * ssf + (1.0 / Cn) * ssf - c
+            dfQ = (dt / dx) + 1.0 / Cn
             ssf = ssf - (fQ / dfQ)
             if isnan(ssf)
                 ssf = 0.0
             end
             ssf = max(ssf, 1.0e-30)
-            if (abs(fQ) <= ϵ) || (count >= max_iters)
+            if (abs(fQ) <= epsilon) || (count >= max_iters)
                 break
             end
             count += 1
@@ -251,12 +251,12 @@ function kinematic_wave_ssf(ssfin, ssfₜ₋₁, ziₜ₋₁, r, kh, β, θₑ,
         # Constrain the lateral subsurface flow rate ssf
         ssf = min(ssf, (ssfmax * dw))
         # On the basis of the lateral flow rate, estimate the amount of groundwater level above surface (saturation excess conditions), then rest = negative
-        zi = ziₜ₋₁ - (ssfin * Δt + r * Δx - ssf * Δt) / (dw * Δx) / θₑ
+        zi = zi_prev - (ssfin * dt + r * dx - ssf * dt) / (dw * dx) / theta_e
         if zi > d
-            ssf = max(ssf - (dw * Δx) * θₑ * (zi - d), 1.0e-30)
+            ssf = max(ssf - (dw * dx) * theta_e * (zi - d), 1.0e-30)
         end
         # Exfiltration rate and set zi to zero.
-        exfilt = min(zi, 0.0) * -θₑ
+        exfilt = min(zi, 0.0) * -theta_e
         zi = clamp(zi, 0.0, d)
 
         return ssf, zi, exfilt
@@ -355,15 +355,15 @@ function lateral_snow_transport!(snow, snowwater, slope, network)
 end
 
 """
-    local_inertial_flow(q0, η0, η1, hf, A, R, length, mannings_n, g, froude_limit, Δt)
+    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,
-    η0,
-    η1,
+    zs0,
+    zs1,
     hf,
     A,
     R,
@@ -371,14 +371,14 @@ function local_inertial_flow(
     mannings_n_sq,
     g,
     froude_limit,
-    Δt,
+    dt,
 )
 
-    slope = (η1 - η0) / length
+    slope = (zs1 - zs0) / length
     pow_R = cbrt(R * R * R * R)
     unit = one(hf)
     q = (
-        (q0 - g * A * Δt * slope) / (unit + g * Δt * mannings_n_sq * abs(q0) / (pow_R * A))
+        (q0 - g * A * dt * slope) / (unit + g * dt * mannings_n_sq * abs(q0) / (pow_R * A))
     )
 
     # if froude number > 1.0, limit flow
@@ -390,36 +390,36 @@ function local_inertial_flow(
 end
 
 """
-    local_inertial_flow(θ, q0, qd, qu, η0, η1, hf, width, length, mannings_n, g, froude_limit, Δt)
+    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.
 (2012).
 """
 function local_inertial_flow(
-    θ,
+    theta,
     q0,
     qd,
     qu,
-    η0,
-    η1,
+    zs0,
+    zs1,
     hf,
     width,
     length,
     mannings_n_sq,
     g,
     froude_limit,
-    Δt,
+    dt,
 )
 
-    slope = (η1 - η0) / length
-    unit = one(θ)
-    half = oftype(θ, 0.5)
+    slope = (zs1 - zs0) / length
+    unit = one(theta)
+    half = oftype(theta, 0.5)
     pow_hf = cbrt(hf * hf * hf * hf * hf * hf * hf)
 
     q = (
-        ((θ * q0 + half * (unit - θ) * (qu + qd)) - g * hf * width * Δt * slope) /
-        (unit + g * Δt * mannings_n_sq * abs(q0) / (pow_hf * width))
+        ((theta * q0 + half * (unit - theta) * (qu + qd)) - g * hf * width * dt * slope) /
+        (unit + g * dt * mannings_n_sq * abs(q0) / (pow_hf * width))
     )
     # if froude number > 1.0, limit flow
     if froude_limit
diff --git a/src/io.jl b/src/io.jl
index 77d6cdf45..16d89baed 100644
--- a/src/io.jl
+++ b/src/io.jl
@@ -335,7 +335,7 @@ function update_cyclic!(model)
     @unpack cyclic_dataset, cyclic_times, cyclic_parameters = reader
 
     # pick up the data that is valid for the past model time step
-    month_day = monthday(clock.time - clock.Δt)
+    month_day = monthday(clock.time - clock.dt)
 
     is_first_timestep = clock.iteration == 1
 
@@ -1345,19 +1345,19 @@ function reset_clock!(clock::Clock, config)
     # we want this method to be mutating
     clock.time = new_clock.time
     clock.iteration = new_clock.iteration
-    clock.Δt = new_clock.Δt
+    clock.dt = new_clock.dt
     return clock
 end
 
 function advance!(clock)
     clock.iteration += 1
-    clock.time += clock.Δt
+    clock.time += clock.dt
     return clock
 end
 
 function rewind!(clock)
     clock.iteration -= 1
-    clock.time -= clock.Δt
+    clock.time -= clock.dt
     return clock
 end
 
diff --git a/src/reservoir_lake.jl b/src/reservoir_lake.jl
index 1daf39d69..0b6acf768 100644
--- a/src/reservoir_lake.jl
+++ b/src/reservoir_lake.jl
@@ -1,5 +1,5 @@
 @get_units @exchange @grid_type @grid_location @with_kw struct SimpleReservoir{T}
-    Δt::T | "s" | 0 | "none" | "none"                   # Model time step [s]
+    dt::T | "s" | 0 | "none" | "none"                   # Model time step [s]
     maxvolume::Vector{T} | "m3"                         # maximum storage (above which water is spilled) [m³]
     area::Vector{T} | "m2"                              # reservoir area [m²]
     maxrelease::Vector{T} | "m3 s-1"                    # maximum amount that can be released if below spillway [m³ s⁻¹]
@@ -23,7 +23,7 @@
 end
 
 
-function initialize_simple_reservoir(config, nc, inds_riv, nriv, pits, Δt)
+function initialize_simple_reservoir(config, nc, inds_riv, nriv, pits, dt)
     # read only reservoir data if reservoirs true
     # allow reservoirs only in river cells
     # note that these locations are only the reservoir outlet pixels
@@ -133,7 +133,7 @@ function initialize_simple_reservoir(config, nc, inds_riv, nriv, pits, Δt)
     n = length(resarea)
     @info "Read `$n` reservoir locations."
     reservoirs = SimpleReservoir{Float}(
-        Δt = Δt,
+        dt = dt,
         demand = resdemand,
         maxrelease = resmaxrelease,
         maxvolume = resmaxvolume,
@@ -170,9 +170,9 @@ element rather than all at once.
 function update(res::SimpleReservoir, i, inflow, timestepsecs)
 
     # limit lake evaporation based on total available volume [m³]
-    precipitation = 0.001 * res.precipitation[i] * (timestepsecs / res.Δt) * res.area[i]
+    precipitation = 0.001 * res.precipitation[i] * (timestepsecs / res.dt) * res.area[i]
     available_volume = res.volume[i] + inflow * timestepsecs + precipitation
-    evap = 0.001 * res.evaporation[i] * (timestepsecs / res.Δt) * res.area[i]
+    evap = 0.001 * res.evaporation[i] * (timestepsecs / res.dt) * res.area[i]
     actevap = min(available_volume, evap) # [m³/timestepsecs]
 
     vol = res.volume[i] + (inflow * timestepsecs) + precipitation - actevap
@@ -205,7 +205,7 @@ function update(res::SimpleReservoir, i, inflow, timestepsecs)
 end
 
 @get_units @exchange @grid_type @grid_location @with_kw struct Lake{T}
-    Δt::T | "s" | 0 | "none" | "none"           # Model time step [s]
+    dt::T | "s" | 0 | "none" | "none"           # Model time step [s]
     lowerlake_ind::Vector{Int} | "-"            # Index of lower lake (linked lakes)
     area::Vector{T} | "m2"                      # lake area [m²]
     maxstorage::Vector{Union{T,Missing}} | "m3" # lake maximum storage from rating curve 1 [m³]
@@ -231,7 +231,7 @@ end
     end
 end
 
-function initialize_lake(config, nc, inds_riv, nriv, pits, Δt)
+function initialize_lake(config, nc, inds_riv, nriv, pits, dt)
     # read only lake data if lakes true
     # allow lakes only in river cells
     # note that these locations are only the lake outlet pixels
@@ -399,7 +399,7 @@ function initialize_lake(config, nc, inds_riv, nriv, pits, Δt)
     end
     n = length(lakearea)
     lakes = Lake{Float}(
-        Δt = Δt,
+        dt = dt,
         lowerlake_ind = lowerlake_ind,
         area = lakearea,
         maxstorage = maximum_storage(lake_storfunc, lake_outflowfunc, lakearea, sh, hq),
@@ -487,9 +487,9 @@ function update(lake::Lake, i, inflow, doy, timestepsecs)
     has_lowerlake = lo != 0
 
     # limit lake evaporation based on total available volume [m³]
-    precipitation = 0.001 * lake.precipitation[i] * (timestepsecs / lake.Δt) * lake.area[i]
+    precipitation = 0.001 * lake.precipitation[i] * (timestepsecs / lake.dt) * lake.area[i]
     available_volume = lake.storage[i] + inflow * timestepsecs + precipitation
-    evap = 0.001 * lake.evaporation[i] * (timestepsecs / lake.Δt) * lake.area[i]
+    evap = 0.001 * lake.evaporation[i] * (timestepsecs / lake.dt) * lake.area[i]
     actevap = min(available_volume, evap) # [m³/timestepsecs]
 
     ### Modified Puls Approach (Burek et al., 2013, LISFLOOD) ###
@@ -530,18 +530,18 @@ function update(lake::Lake, i, inflow, doy, timestepsecs)
         else
             if diff_wl >= 0.0
                 if lake.waterlevel[i] > lake.threshold[i]
-                    Δh = lake.waterlevel[i] - lake.threshold[i]
-                    outflow = lake.b[i] * pow(Δh, lake.e[i])
-                    maxflow = (Δh * lake.area[i]) / timestepsecs
+                    dh = lake.waterlevel[i] - lake.threshold[i]
+                    outflow = lake.b[i] * pow(dh, lake.e[i])
+                    maxflow = (dh * lake.area[i]) / timestepsecs
                     outflow = min(outflow, maxflow)
                 else
                     outflow = Float(0)
                 end
             else
                 if lake.waterlevel[lo] > lake.threshold[i]
-                    Δh = lake.waterlevel[lo] - lake.threshold[i]
-                    outflow = -1.0 * lake.b[i] * pow(Δh, lake.e[i])
-                    maxflow = (Δh * lake.area[lo]) / timestepsecs
+                    dh = lake.waterlevel[lo] - lake.threshold[i]
+                    outflow = -1.0 * lake.b[i] * pow(dh, lake.e[i])
+                    maxflow = (dh * lake.area[lo]) / timestepsecs
                     outflow = max(outflow, -maxflow)
                 else
                     outflow = Float(0)
diff --git a/src/sbm.jl b/src/sbm.jl
index e0282cda3..352ab5c08 100644
--- a/src/sbm.jl
+++ b/src/sbm.jl
@@ -1,6 +1,6 @@
 @get_units @exchange @grid_type @grid_location @with_kw struct SBM{T,N,M}
     # Model time step [s]
-    Δt::T | "s" | 0 | "none" | "none"
+    dt::T | "s" | 0 | "none" | "none"
     # Maximum number of soil layers
     maxlayers::Int | "-" | 0 | "none" | "none"
     # number of cells
@@ -14,11 +14,11 @@
     # Fraction of river [-]
     riverfrac::Vector{T} | "-"
     # Saturated water content (porosity) [-]
-    θₛ::Vector{T} | "-"
+    theta_s::Vector{T} | "-"
     # Residual water content [-]
-    θᵣ::Vector{T} | "-"
+    theta_r::Vector{T} | "-"
     # Vertical hydraulic conductivity [mm Δt⁻¹] at soil surface
-    kv₀::Vector{T}
+    kv_0::Vector{T}
     # Vertical hydraulic conductivity [mm Δt⁻¹] per soil layer
     kv::Vector{SVector{N,T}} | "-"
     # Muliplication factor [-] applied to kv_z (vertical flow)
@@ -57,9 +57,9 @@
     stemflow::Vector{T}
     # Throughfall [mm Δt⁻¹]
     throughfall::Vector{T}
-    # A scaling parameter [mm⁻¹] (controls exponential decline of kv₀)
+    # A scaling parameter [mm⁻¹] (controls exponential decline of kv_0)
     f::Vector{T} | "mm-1"
-    # Depth [mm] from soil surface for which exponential decline of kv₀ is valid
+    # Depth [mm] from soil surface for which exponential decline of kv_0 is valid
     z_exp::Vector{T} | "mm"
     # Depth [mm] from soil surface for which layered profile is valid
     z_layered::Vector{T} | "mm"
@@ -140,15 +140,15 @@
     runoff::Vector{T}
     # Net surface runoff (surface runoff - actual open water evaporation) [mm Δt⁻¹]
     net_runoff::Vector{T}
-    # Volumetric water content [-] per soil layer (including θᵣ and saturated zone)
+    # Volumetric water content [-] per soil layer (including theta_r and saturated zone)
     vwc::Vector{SVector{N,T}} | "-"
-    # Volumetric water content [%] per soil layer (including θᵣ and saturated zone)
+    # Volumetric water content [%] per soil layer (including theta_r and saturated zone)
     vwc_perc::Vector{SVector{N,T}} | "%"
-    # Root water storage [mm] in unsaturated and saturated zone (excluding θᵣ)
+    # Root water storage [mm] in unsaturated and saturated zone (excluding theta_r)
     rootstore::Vector{T} | "mm"
-    # Volumetric water content [-] in root zone (including θᵣ and saturated zone)
+    # Volumetric water content [-] in root zone (including theta_r and saturated zone)
     vwc_root::Vector{T} | "-"
-    # Volumetric water content [%] in root zone (including θᵣ and saturated zone)
+    # Volumetric water content [%] in root zone (including theta_r and saturated zone)
     vwc_percroot::Vector{T} | "%"
     # Amount of available water in the unsaturated zone [mm]
     ustoredepth::Vector{T} | "mm"
@@ -160,7 +160,7 @@
     actleakage::Vector{T}
     ### Snow parameters ###
     # Degree-day factor [mm ᵒC⁻¹ Δt⁻¹]
-    cfmax::Vector{T} | "mm ᵒC-1 Δt-1"
+    cfmax::Vector{T} | "mm ᵒC-1 dt-1"
     # Threshold temperature for snowfall [ᵒC]
     tt::Vector{T} | "ᵒC"
     # Threshold temperature interval length [ᵒC]
@@ -182,9 +182,9 @@
     # Threshold temperature for snowfall above glacier [ᵒC]
     g_tt::Vector{T} | "ᵒC"
     # Degree-day factor [mm ᵒC⁻¹ Δt⁻¹] for glacier
-    g_cfmax::Vector{T} | "mm ᵒC-1 Δt-1"
+    g_cfmax::Vector{T} | "mm ᵒC-1 dt-1"
     # Fraction of the snowpack on top of the glacier converted into ice [Δt⁻¹]
-    g_sifrac::Vector{T} | "Δt-1"
+    g_sifrac::Vector{T} | "dt-1"
     # Water within the glacier [mm]
     glacierstore::Vector{T} | "mm"
     # Fraction covered by a glacier [-]
@@ -262,7 +262,7 @@ end
 
 function initialize_sbm(nc, config, riverfrac, inds)
 
-    Δt = Second(config.timestepsecs)
+    dt = Second(config.timestepsecs)
     config_thicknesslayers = get(config.model, "thicknesslayers", Float[])
     ksat_profile = get(config.input.vertical, "ksat_profile", "exponential")::String
     if length(config_thicknesslayers) > 0
@@ -283,7 +283,7 @@ function initialize_sbm(nc, config, riverfrac, inds)
             sel = inds,
             defaults = 3.75653,
             type = Float,
-        ) .* (Δt / basetimestep)
+        ) .* (dt / basetimestep)
     tt = ncread(nc, config, "vertical.tt"; sel = inds, defaults = 0.0, type = Float)
     tti = ncread(nc, config, "vertical.tti"; sel = inds, defaults = 1.0, type = Float)
     ttm = ncread(nc, config, "vertical.ttm"; sel = inds, defaults = 0.0, type = Float)
@@ -296,7 +296,7 @@ function initialize_sbm(nc, config, riverfrac, inds)
             sel = inds,
             defaults = 0.1125,
             type = Float,
-        ) .* (Δt / basetimestep)
+        ) .* (dt / basetimestep)
     cf_soil =
         ncread(nc, config, "vertical.cf_soil"; sel = inds, defaults = 0.038, type = Float)
     # glacier parameters
@@ -318,7 +318,7 @@ function initialize_sbm(nc, config, riverfrac, inds)
             defaults = 3.0,
             type = Float,
             fill = 0.0,
-        ) .* (Δt / basetimestep)
+        ) .* (dt / basetimestep)
     g_sifrac =
         ncread(
             nc,
@@ -328,7 +328,7 @@ function initialize_sbm(nc, config, riverfrac, inds)
             defaults = 0.001,
             type = Float,
             fill = 0.0,
-        ) .* (Δt / basetimestep)
+        ) .* (dt / basetimestep)
     glacierfrac = ncread(
         nc,
         config,
@@ -348,34 +348,31 @@ function initialize_sbm(nc, config, riverfrac, inds)
         fill = 0.0,
     )
     # soil parameters
-    θₛ = ncread(
+    theta_s = ncread(
         nc,
         config,
         "vertical.theta_s";
-        alias = "vertical.θₛ",
         sel = inds,
         defaults = 0.6,
         type = Float,
     )
-    θᵣ = ncread(
+    theta_r = ncread(
         nc,
         config,
         "vertical.theta_r";
-        alias = "vertical.θᵣ",
         sel = inds,
         defaults = 0.01,
         type = Float,
     )
-    kv₀ =
+    kv_0 =
         ncread(
             nc,
             config,
             "vertical.kv_0";
-            alias = "vertical.kv₀",
             sel = inds,
             defaults = 3000.0,
             type = Float,
-        ) .* (Δt / basetimestep)
+        ) .* (dt / basetimestep)
     f = ncread(nc, config, "vertical.f"; sel = inds, defaults = 0.001, type = Float)
     hb = ncread(nc, config, "vertical.hb"; sel = inds, defaults = 10.0, type = Float)
     soilthickness = ncread(
@@ -394,7 +391,7 @@ function initialize_sbm(nc, config, riverfrac, inds)
             sel = inds,
             defaults = 10.0,
             type = Float,
-        ) .* (Δt / basetimestep)
+        ) .* (dt / basetimestep)
     infiltcapsoil =
         ncread(
             nc,
@@ -403,7 +400,7 @@ function initialize_sbm(nc, config, riverfrac, inds)
             sel = inds,
             defaults = 100.0,
             type = Float,
-        ) .* (Δt / basetimestep)
+        ) .* (dt / basetimestep)
     maxleakage =
         ncread(
             nc,
@@ -412,7 +409,7 @@ function initialize_sbm(nc, config, riverfrac, inds)
             sel = inds,
             defaults = 0.0,
             type = Float,
-        ) .* (Δt / basetimestep)
+        ) .* (dt / basetimestep)
 
     c = ncread(
         nc,
@@ -487,10 +484,10 @@ function initialize_sbm(nc, config, riverfrac, inds)
 
     cmax, e_r, canopygapfraction, sl, swood, kext = initialize_canopy(nc, config, inds)
 
-    θₑ = θₛ .- θᵣ
-    soilwatercapacity = soilthickness .* θₑ
+    theta_e = theta_s .- theta_r
+    soilwatercapacity = soilthickness .* theta_e
     satwaterdepth = 0.85 .* soilwatercapacity # cold state value for satwaterdepth
-    zi = max.(0.0, soilthickness .- satwaterdepth ./ θₑ) # cold state value for zi
+    zi = max.(0.0, soilthickness .- satwaterdepth ./ theta_e) # cold state value for zi
 
     # these are filled in the loop below
     # TODO see if we can replace this approach
@@ -545,7 +542,7 @@ function initialize_sbm(nc, config, riverfrac, inds)
                 defaults = 1000.0,
                 type = Float,
                 dimname = :layer,
-            ) .* (Δt / basetimestep)
+            ) .* (dt / basetimestep)
         if size(kv, 1) != maxlayers
             parname = param(config.input.vertical, "kv")
             size1 = size(kv, 1)
@@ -579,16 +576,16 @@ function initialize_sbm(nc, config, riverfrac, inds)
     end
 
     sbm = SBM{Float,maxlayers,maxlayers + 1}(
-        Δt = tosecond(Δt),
+        dt = tosecond(dt),
         maxlayers = maxlayers,
         n = n,
         nlayers = nlayers,
         n_unsatlayers = n_unsatlayers,
         nlayers_kv = nlayers_kv,
         riverfrac = riverfrac,
-        θₛ = θₛ,
-        θᵣ = θᵣ,
-        kv₀ = kv₀,
+        theta_s = theta_s,
+        theta_r = theta_r,
+        kv_0 = kv_0,
         kv = svectorscopy(kv, Val{maxlayers}()),
         kvfrac = svectorscopy(kvfrac, Val{maxlayers}()),
         hb = hb,
@@ -714,7 +711,7 @@ function update_until_snow(sbm::SBM, config)
 
         canopy_potevap =
             sbm.kc[i] * sbm.potential_evaporation[i] * (1.0 - canopygapfraction)
-        if Second(sbm.Δt) >= Hour(23)
+        if Second(sbm.dt) >= Hour(23)
             throughfall, interception, stemflow, canopystorage = rainfall_interception_gash(
                 cmax,
                 e_r,
@@ -794,7 +791,7 @@ function update_until_recharge(sbm::SBM, config)
                     sbm.g_tt[i],
                     sbm.g_cfmax[i],
                     sbm.g_sifrac[i],
-                    Second(sbm.Δt),
+                    Second(sbm.dt),
                 )
                 # Convert to mm per grid cell and add to snowmelt
                 glaciermelt = glaciermelt * sbm.glacierfrac[i]
@@ -869,13 +866,13 @@ function update_until_recharge(sbm::SBM, config)
                     sbm.satwaterdepth[i],
                     kv_z,
                     usl[1],
-                    sbm.θₛ[i],
-                    sbm.θᵣ[i],
+                    sbm.theta_s[i],
+                    sbm.theta_r[i],
                 )
                 usld = setindex(usld, ustorelayerdepth, 1)
             else
                 for m = 1:n_usl
-                    l_sat = usl[m] * (sbm.θₛ[i] - sbm.θᵣ[i])
+                    l_sat = usl[m] * (sbm.theta_s[i] - sbm.theta_r[i])
                     kv_z = hydraulic_conductivity_at_depth(sbm, z[m], i, m, ksat_profile)
                     ustorelayerdepth =
                         m == 1 ? sbm.ustorelayerdepth[i][m] + infiltsoilpath :
@@ -901,12 +898,12 @@ function update_until_recharge(sbm::SBM, config)
                     # Check if groundwater level lies below the surface
                     soilevapunsat =
                         potsoilevap *
-                        min(1.0, usld[1] / (sbm.zi[i] * (sbm.θₛ[i] - sbm.θᵣ[i])))
+                        min(1.0, usld[1] / (sbm.zi[i] * (sbm.theta_s[i] - sbm.theta_r[i])))
                 else
                     # In case first layer contains no saturated storage
                     soilevapunsat =
                         potsoilevap *
-                        min(1.0, usld[1] / (usl[1] * ((sbm.θₛ[i] - sbm.θᵣ[i]))))
+                        min(1.0, usld[1] / (usl[1] * ((sbm.theta_s[i] - sbm.theta_r[i]))))
                 end
             end
             # Ensure that the unsaturated evaporation rate does not exceed the
@@ -927,7 +924,7 @@ function update_until_recharge(sbm::SBM, config)
                     min(1.0, (sbm.act_thickl[i][1] - sbm.zi[i]) / sbm.act_thickl[i][1])
                 soilevapsat = min(
                     soilevapsat,
-                    (sbm.act_thickl[i][1] - sbm.zi[i]) * (sbm.θₛ[i] - sbm.θᵣ[i]),
+                    (sbm.act_thickl[i][1] - sbm.zi[i]) * (sbm.theta_s[i] - sbm.theta_r[i]),
                 )
             else
                 soilevapsat = 0.0
@@ -953,8 +950,8 @@ function update_until_recharge(sbm::SBM, config)
                 actevapustore,
                 sbm.c[i][k],
                 usl[k],
-                sbm.θₛ[i],
-                sbm.θᵣ[i],
+                sbm.theta_s[i],
+                sbm.theta_r[i],
                 sbm.hb[i],
                 ust,
             )
@@ -964,7 +961,7 @@ function update_until_recharge(sbm::SBM, config)
         # check soil moisture balance per layer
         du = 0.0
         for k = n_usl:-1:1
-            du = max(0.0, usld[k] - usl[k] * (sbm.θₛ[i] - sbm.θᵣ[i]))
+            du = max(0.0, usld[k] - usl[k] * (sbm.theta_s[i] - sbm.theta_r[i]))
             usld = setindex(usld, usld[k] - du, k)
             if k > 1
                 usld = setindex(usld, usld[k-1] + du, k - 1)
@@ -1006,7 +1003,7 @@ function update_until_recharge(sbm::SBM, config)
             netcapflux = capflux
             for k = n_usl:-1:1
                 toadd =
-                    min(netcapflux, max(usl[k] * (sbm.θₛ[i] - sbm.θᵣ[i]) - usld[k], 0.0))
+                    min(netcapflux, max(usl[k] * (sbm.theta_s[i] - sbm.theta_r[i]) - usld[k], 0.0))
                 usld = setindex(usld, usld[k] + toadd, k)
                 netcapflux = netcapflux - toadd
                 actcapflux = actcapflux + toadd
@@ -1074,7 +1071,7 @@ function update_after_subsurfaceflow(sbm::SBM, zi, exfiltsatwater)
         exfiltustore = 0.0
         for k = sbm.n_unsatlayers[i]:-1:1
             if k <= n_usl
-                exfiltustore = max(0, usld[k] - usl[k] * (sbm.θₛ[i] - sbm.θᵣ[i]))
+                exfiltustore = max(0, usld[k] - usl[k] * (sbm.theta_s[i] - sbm.theta_r[i]))
             else
                 exfiltustore = usld[k]
             end
@@ -1100,13 +1097,13 @@ function update_after_subsurfaceflow(sbm::SBM, zi, exfiltsatwater)
             if k <= n_usl
                 vwc = setindex(
                     vwc,
-                    (usld[k] + (sbm.act_thickl[i][k] - usl[k]) * (sbm.θₛ[i] - sbm.θᵣ[i])) / sbm.act_thickl[i][k] + sbm.θᵣ[i],
+                    (usld[k] + (sbm.act_thickl[i][k] - usl[k]) * (sbm.theta_s[i] - sbm.theta_r[i])) / sbm.act_thickl[i][k] + sbm.theta_r[i],
                     k,
                 )
             else
-                vwc = setindex(vwc, sbm.θₛ[i], k)
+                vwc = setindex(vwc, sbm.theta_s[i], k)
             end
-            vwc_perc = setindex(vwc_perc, (vwc[k] / sbm.θₛ[i]) * 100.0, k)
+            vwc_perc = setindex(vwc_perc, (vwc[k] / sbm.theta_s[i]) * 100.0, k)
         end
 
         rootstore_unsat = 0
@@ -1117,12 +1114,12 @@ function update_after_subsurfaceflow(sbm::SBM, zi, exfiltsatwater)
                 usld[k]
         end
 
-        rootstore_sat = max(0.0, sbm.rootingdepth[i] - zi[i]) * (sbm.θₛ[i] - sbm.θᵣ[i])
+        rootstore_sat = max(0.0, sbm.rootingdepth[i] - zi[i]) * (sbm.theta_s[i] - sbm.theta_r[i])
         rootstore = rootstore_sat + rootstore_unsat
-        vwc_root = rootstore / sbm.rootingdepth[i] + sbm.θᵣ[i]
-        vwc_percroot = (vwc_root / sbm.θₛ[i]) * 100.0
+        vwc_root = rootstore / sbm.rootingdepth[i] + sbm.theta_r[i]
+        vwc_percroot = (vwc_root / sbm.theta_s[i]) * 100.0
 
-        satwaterdepth = (sbm.soilthickness[i] - zi[i]) * (sbm.θₛ[i] - sbm.θᵣ[i])
+        satwaterdepth = (sbm.soilthickness[i] - zi[i]) * (sbm.theta_s[i] - sbm.theta_r[i])
 
         # update the outputs and states
         sbm.n_unsatlayers[i] = n_usl
diff --git a/src/sbm_gwf_model.jl b/src/sbm_gwf_model.jl
index 429708436..6a76b5bdf 100644
--- a/src/sbm_gwf_model.jl
+++ b/src/sbm_gwf_model.jl
@@ -21,7 +21,7 @@ function initialize_sbm_gwf_model(config::Config)
 
     reader = prepare_reader(config)
     clock = Clock(config, reader)
-    Δt = clock.Δt
+    dt = clock.dt
 
     do_reservoirs = get(config.model, "reservoirs", false)::Bool
     do_lakes = get(config.model, "lakes", false)::Bool
@@ -81,7 +81,7 @@ function initialize_sbm_gwf_model(config::Config)
     pits = zeros(Bool, modelsize_2d)
     if do_reservoirs
         reservoirs, resindex, reservoir, pits =
-            initialize_simple_reservoir(config, nc, inds_riv, nriv, pits, tosecond(Δt))
+            initialize_simple_reservoir(config, nc, inds_riv, nriv, pits, tosecond(dt))
     else
         reservoir = ()
         reservoirs = nothing
@@ -91,7 +91,7 @@ function initialize_sbm_gwf_model(config::Config)
     # lakes
     if do_lakes
         lakes, lakeindex, lake, pits =
-            initialize_lake(config, nc, inds_riv, nriv, pits, tosecond(Δt))
+            initialize_lake(config, nc, inds_riv, nriv, pits, tosecond(dt))
     else
         lake = ()
         lakes = nothing
@@ -99,9 +99,9 @@ function initialize_sbm_gwf_model(config::Config)
     end
 
     # overland flow (kinematic wave)
-    βₗ =
+    landslope =
         ncread(nc, config, "lateral.land.slope"; optional = false, sel = inds, type = Float)
-    clamp!(βₗ, 0.00001, Inf)
+    clamp!(landslope, 0.00001, Inf)
     ldd_2d = ncread(nc, config, "ldd"; optional = false, allow_missing = true)
 
     ldd = ldd_2d[inds]
@@ -116,19 +116,19 @@ 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, βₗ, 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 = βₗ,
+            sl = landslope,
             dl,
             width = map(det_surfacewidth, dw, riverwidth, river),
             iterate = kinwave_it,
             tstep = kw_land_tstep,
-            Δt,
+            dt,
         )
     elseif land_routing == "local-inertial"
         index_river_nf = rev_inds_riv[inds] # not filtered (with zeros)
@@ -146,7 +146,7 @@ function initialize_sbm_gwf_model(config::Config)
             inds_riv,
             river,
             waterbody = !=(0).(resindex + lakeindex),
-            Δt,
+            dt,
         )
     end
 
@@ -169,7 +169,7 @@ function initialize_sbm_gwf_model(config::Config)
             lake = lakes,
             iterate = kinwave_it,
             tstep = kw_river_tstep,
-            Δt = Δt,
+            dt = dt,
         )
     elseif river_routing == "local-inertial"
         rf, nodes_at_link = initialize_shallowwater_river(
@@ -184,7 +184,7 @@ function initialize_sbm_gwf_model(config::Config)
             reservoir = reservoirs,
             lake_index = lakeindex,
             lake = lakes,
-            Δt = Δt,
+            dt = dt,
             floodplain = floodplain_1d,
         )
     else
@@ -247,7 +247,7 @@ function initialize_sbm_gwf_model(config::Config)
 
     # reset zi and satwaterdepth with groundwater head from unconfined aquifer
     sbm.zi .= (altitude .- initial_head) .* 1000.0
-    sbm.satwaterdepth .= (sbm.soilthickness .- sbm.zi) .* (sbm.θₛ .- sbm.θᵣ)
+    sbm.satwaterdepth .= (sbm.soilthickness .- sbm.zi) .* (sbm.theta_s .- sbm.theta_r)
 
     # river boundary of unconfined aquifer
     infiltration_conductance = ncread(
@@ -421,7 +421,7 @@ function initialize_sbm_gwf_model(config::Config)
             reverse_indices = rev_inds,
             xl = xl,
             yl = yl,
-            slope = βₗ,
+            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 = βₗ,
+            slope = landslope,
             altitude = altitude,
             index_river = index_river_nf,
             staggered_indices = indices,
@@ -512,20 +512,20 @@ function update(model::Model{N,L,V,R,W,T}) where {N,L,V,R,W,T<:SbmGwfModel}
     # determine stable time step for groundwater flow
     conductivity_profile =
         get(config.input.lateral.subsurface, "conductivity_profile", "uniform")
-    Δt_gw = stable_timestep(lateral.subsurface.flow.aquifer, conductivity_profile) # time step in day (Float64)
-    Δt_sbm = (vertical.Δt / tosecond(basetimestep)) # vertical.Δt is in seconds (Float64)
-    if Δt_gw < Δt_sbm
+    dt_gw = stable_timestep(lateral.subsurface.flow.aquifer, conductivity_profile) # time step in day (Float64)
+    dt_sbm = (vertical.dt / tosecond(basetimestep)) # vertical.dt is in seconds (Float64)
+    if dt_gw < dt_sbm
         @warn(
-            "stable time step Δt $Δt_gw for groundwater flow is smaller than sbm Δt $Δt_sbm"
+            "stable time step dt $dt_gw for groundwater flow is smaller than sbm dt $dt_sbm"
         )
     end
 
     Q = zeros(vertical.n)
     # exchange of recharge between vertical sbm concept and groundwater flow domain
     # recharge rate groundwater is required in units [m d⁻¹]
-    lateral.subsurface.recharge.rate .= vertical.recharge ./ 1000.0 .* (1.0 / Δt_sbm)
+    lateral.subsurface.recharge.rate .= vertical.recharge ./ 1000.0 .* (1.0 / dt_sbm)
     # update groundwater domain
-    update(lateral.subsurface.flow, Q, Δt_sbm, conductivity_profile)
+    update(lateral.subsurface.flow, Q, dt_sbm, conductivity_profile)
 
     # determine excess water depth [m] (exfiltwater) in groundwater domain (head > surface)
     # and reset head
@@ -545,7 +545,7 @@ function update(model::Model{N,L,V,R,W,T}) where {N,L,V,R,W,T<:SbmGwfModel}
     )
 
     ssf_toriver = zeros(vertical.n)
-    ssf_toriver[inds_riv] = -lateral.subsurface.river.flux ./ lateral.river.Δt
+    ssf_toriver[inds_riv] = -lateral.subsurface.river.flux ./ lateral.river.dt
     surface_routing(model, ssf_toriver = ssf_toriver)
 
     return model
diff --git a/src/sbm_model.jl b/src/sbm_model.jl
index 626c979a5..bc238bab9 100644
--- a/src/sbm_model.jl
+++ b/src/sbm_model.jl
@@ -14,7 +14,7 @@ function initialize_sbm_model(config::Config)
 
     reader = prepare_reader(config)
     clock = Clock(config, reader)
-    Δt = clock.Δt
+    dt = clock.dt
 
     do_reservoirs = get(config.model, "reservoirs", false)::Bool
     do_lakes = get(config.model, "lakes", false)::Bool
@@ -78,7 +78,7 @@ function initialize_sbm_model(config::Config)
     pits = zeros(Bool, modelsize_2d)
     if do_reservoirs
         reservoirs, resindex, reservoir, pits =
-            initialize_simple_reservoir(config, nc, inds_riv, nriv, pits, tosecond(Δt))
+            initialize_simple_reservoir(config, nc, inds_riv, nriv, pits, tosecond(dt))
     else
         reservoir = ()
         reservoirs = nothing
@@ -88,7 +88,7 @@ function initialize_sbm_model(config::Config)
     # lakes
     if do_lakes
         lakes, lakeindex, lake, pits =
-            initialize_lake(config, nc, inds_riv, nriv, pits, tosecond(Δt))
+            initialize_lake(config, nc, inds_riv, nriv, pits, tosecond(dt))
     else
         lake = ()
         lakes = nothing
@@ -102,9 +102,9 @@ function initialize_sbm_model(config::Config)
         ldd = set_pit_ldd(pits_2d, ldd, inds)
     end
 
-    βₗ =
+    landslope =
         ncread(nc, config, "lateral.land.slope"; optional = false, sel = inds, type = Float)
-    clamp!(βₗ, 0.00001, Inf)
+    clamp!(landslope, 0.00001, Inf)
 
     dl = map(detdrainlength, ldd, xl, yl)
     dw = (xl .* yl) ./ dl
@@ -122,25 +122,25 @@ function initialize_sbm_model(config::Config)
             type = Float,
         )
 
-        # unit for lateral subsurface flow component is [m³ d⁻¹], sbm.kv₀ [mm Δt⁻¹]
-        kh₀ = khfrac .* sbm.kv₀ .* 0.001 .* (basetimestep / Δt)
+        # unit for lateral subsurface flow component is [m³ d⁻¹], sbm.kv_0 [mm Δt⁻¹]
+        kh_0 = khfrac .* sbm.kv_0 .* 0.001 .* (basetimestep / dt)
         f = sbm.f .* 1000.0
         zi = sbm.zi .* 0.001
         soilthickness = sbm.soilthickness .* 0.001
         z_exp = sbm.z_exp .* 0.001
 
         ssf = LateralSSF{Float}(
-            kh₀ = kh₀,
+            kh_0 = kh_0,
             f = f,
             kh = fill(mv, n),
             khfrac = khfrac,
             zi = zi,
             z_exp = z_exp,
             soilthickness = soilthickness,
-            θₛ = sbm.θₛ,
-            θᵣ = sbm.θᵣ,
-            Δt = Δt / basetimestep,
-            βₗ = βₗ,
+            theta_s = sbm.theta_s,
+            theta_r = sbm.theta_r,
+            dt = dt / basetimestep,
+            slope = landslope,
             dl = dl,
             dw = dw,
             exfiltwater = fill(mv, n),
@@ -163,7 +163,7 @@ function initialize_sbm_model(config::Config)
         # when the SBM model is coupled (BMI) to a groundwater model, the following
         # variables are expected to be exchanged from the groundwater model.
         ssf = GroundwaterExchange{Float}(
-            Δt = Δt / basetimestep,
+            dt = dt / basetimestep,
             exfiltwater = fill(mv, n),
             zi = fill(mv, n),
             to_river = fill(mv, n),
@@ -180,19 +180,19 @@ 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, βₗ, 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 = βₗ,
+            sl = landslope,
             dl,
             width = map(det_surfacewidth, dw, riverwidth, river),
             iterate = kinwave_it,
             tstep = kw_land_tstep,
-            Δt,
+            dt,
         )
     elseif land_routing == "local-inertial"
         index_river_nf = rev_inds_riv[inds] # not filtered (with zeros)
@@ -210,7 +210,7 @@ function initialize_sbm_model(config::Config)
             inds_riv,
             river,
             waterbody = !=(0).(resindex + lakeindex),
-            Δt,
+            dt,
         )
     end
 
@@ -231,7 +231,7 @@ function initialize_sbm_model(config::Config)
             lake = lakes,
             iterate = kinwave_it,
             tstep = kw_river_tstep,
-            Δt = Δt,
+            dt = dt,
         )
     elseif river_routing == "local-inertial"
         rf, nodes_at_link = initialize_shallowwater_river(
@@ -246,7 +246,7 @@ function initialize_sbm_model(config::Config)
             reservoir = reservoirs,
             lake_index = lakeindex,
             lake = lakes,
-            Δt = Δt,
+            dt = dt,
             floodplain = floodplain_1d,
         )
     else
@@ -340,7 +340,7 @@ function initialize_sbm_model(config::Config)
         reverse_indices = rev_inds,
         xl,
         yl,
-        slope = βₗ,
+        slope = landslope,
     )
     if land_routing == "local-inertial"
         land = merge(land, (index_river = index_river_nf, staggered_indices = indices))
@@ -518,7 +518,7 @@ function set_states(
             max.(
                 0.0,
                 vertical.soilthickness .-
-                vertical.satwaterdepth ./ (vertical.θₛ .- vertical.θᵣ),
+                vertical.satwaterdepth ./ (vertical.theta_s .- vertical.theta_r),
             )
         vertical.zi .= zi
         if land_routing == "kinematic-wave"
diff --git a/src/sediment.jl b/src/sediment.jl
index f77ec9039..dd708a32a 100644
--- a/src/sediment.jl
+++ b/src/sediment.jl
@@ -14,9 +14,9 @@
     # Depth of overland flow [m]
     h_land::Vector{T} | "m"
     # Canopy interception [mm Δt⁻¹]
-    interception::Vector{T} | "mm Δt⁻¹"
+    interception::Vector{T} | "mm dt-1"
     # Precipitation [mm Δt⁻¹]
-    precipitation::Vector{T} | "mm Δt⁻¹"
+    precipitation::Vector{T} | "mm dt-1"
     # Overland flow [m3/s]
     q_land::Vector{T} | "m3 s-1"
     # Canopy height [m]
@@ -38,17 +38,17 @@
     # USLE soil erodibility factor [-]
     usleK::Vector{T} | "-"
     # Sediment eroded by rainfall [ton Δt⁻¹]
-    sedspl::Vector{T} | "t Δt⁻¹"
+    sedspl::Vector{T} | "t dt-1"
     # Sediment eroded by overland flow [ton Δt⁻¹]
-    sedov::Vector{T} | "t Δt⁻¹"
+    sedov::Vector{T} | "t dt-1"
     # Total eroded soil [ton Δt⁻¹]
-    soilloss::Vector{T} | "t Δt⁻¹"
+    soilloss::Vector{T} | "t dt-1"
     # Eroded soil per particle class [ton Δt⁻¹]
-    erosclay::Vector{T} | "t Δt⁻¹"  # clay
-    erossilt::Vector{T} | "t Δt⁻¹"  # silt
-    erossand::Vector{T} | "t Δt⁻¹"  # sand
-    erossagg::Vector{T} | "t Δt⁻¹"  # small aggregates
-    eroslagg::Vector{T} | "t Δt⁻¹"  # large aggregates
+    erosclay::Vector{T} | "t dt-1"  # clay
+    erossilt::Vector{T} | "t dt-1"  # silt
+    erossand::Vector{T} | "t dt-1"  # sand
+    erossagg::Vector{T} | "t dt-1"  # small aggregates
+    eroslagg::Vector{T} | "t dt-1"  # large aggregates
     ## Interception related to leaf_area_index climatology ###
     # Specific leaf storage [mm]
     sl::Vector{T} | "mm"
@@ -85,13 +85,13 @@
     # Filter with river cells
     rivcell::Vector{T} | "-"
     # Total transport capacity of overland flow [ton Δt⁻¹]
-    TCsed::Vector{T} | "t Δt⁻¹"
+    TCsed::Vector{T} | "t dt-1"
     # Transport capacity of overland flow per particle class [ton Δt⁻¹]
-    TCclay::Vector{T} | "t Δt⁻¹"
-    TCsilt::Vector{T} | "t Δt⁻¹"
-    TCsand::Vector{T} | "t Δt⁻¹"
-    TCsagg::Vector{T} | "t Δt⁻¹"
-    TClagg::Vector{T} | "t Δt⁻¹"
+    TCclay::Vector{T} | "t dt-1"
+    TCsilt::Vector{T} | "t dt-1"
+    TCsand::Vector{T} | "t dt-1"
+    TCsagg::Vector{T} | "t dt-1"
+    TClagg::Vector{T} | "t dt-1"
 
     function LandSediment{T}(args...) where {T}
         equal_size_vectors(args)
@@ -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
-    βₗ = slope
-    clamp!(βₗ, 0.00001, Inf)
+    clamp!(slope, 0.00001, Inf)
     ldd_2d = ncread(nc, config, "ldd"; optional = false, allow_missing = true)
     ldd = ldd_2d[inds]
 
@@ -305,8 +304,8 @@ function update_until_ols(eros::LandSediment, config)
     # Options from config
     do_lai = haskey(config.input.vertical, "leaf_area_index")
     rainerosmethod = get(config.model, "rainerosmethod", "answers")::String
-    Δt = Second(config.timestepsecs)
-    ts = Float(Δt.value)
+    dt = Second(config.timestepsecs)
+    ts = Float(dt.value)
 
     for i = 1:eros.n
 
@@ -401,8 +400,8 @@ function update_until_oltransport(ols::LandSediment, config::Config)
 
     do_river = get(config.model, "runrivermodel", false)::Bool
     tcmethod = get(config.model, "landtransportmethod", "yalinpart")::String
-    Δt = Second(config.timestepsecs)
-    ts = Float(Δt.value)
+    dt = Second(config.timestepsecs)
+    ts = Float(dt.value)
 
     for i = 1:ols.n
         sinslope = sin(atan(ols.slope[i]))
@@ -598,35 +597,35 @@ end
     # Filter with river cells
     rivcell::Vector{T} | "-"
     # Total eroded soil [ton Δt⁻¹]
-    soilloss::Vector{T} | "t Δt⁻¹"
+    soilloss::Vector{T} | "t dt-1"
     # Eroded soil per particle class [ton Δt⁻¹]
-    erosclay::Vector{T} | "t Δt⁻¹"
-    erossilt::Vector{T} | "t Δt⁻¹"
-    erossand::Vector{T} | "t Δt⁻¹"
-    erossagg::Vector{T} | "t Δt⁻¹"
-    eroslagg::Vector{T} | "t Δt⁻¹"
+    erosclay::Vector{T} | "t dt-1"
+    erossilt::Vector{T} | "t dt-1"
+    erossand::Vector{T} | "t dt-1"
+    erossagg::Vector{T} | "t dt-1"
+    eroslagg::Vector{T} | "t dt-1"
     # Total transport capacity of overland flow [ton Δt⁻¹]
-    TCsed::Vector{T} | "t Δt⁻¹"
+    TCsed::Vector{T} | "t dt-1"
     # Transport capacity of overland flow per particle class [ton Δt⁻¹]
-    TCclay::Vector{T} | "t Δt⁻¹"
-    TCsilt::Vector{T} | "t Δt⁻¹"
-    TCsand::Vector{T} | "t Δt⁻¹"
-    TCsagg::Vector{T} | "t Δt⁻¹"
-    TClagg::Vector{T} | "t Δt⁻¹"
-    # Sediment flux in overland flow [ton Δt⁻¹]
-    olsed::Vector{T} | "t Δt⁻¹"
-    olclay::Vector{T} | "t Δt⁻¹"
-    olsilt::Vector{T} | "t Δt⁻¹"
-    olsand::Vector{T} | "t Δt⁻¹"
-    olsagg::Vector{T} | "t Δt⁻¹"
-    ollagg::Vector{T} | "t Δt⁻¹"
-    # Sediment reaching the river with overland flow [ton Δt⁻¹]
-    inlandsed::Vector{T} | "t Δt⁻¹"
-    inlandclay::Vector{T} | "t Δt⁻¹"
-    inlandsilt::Vector{T} | "t Δt⁻¹"
-    inlandsand::Vector{T} | "t Δt⁻¹"
-    inlandsagg::Vector{T} | "t Δt⁻¹"
-    inlandlagg::Vector{T} | "t Δt⁻¹"
+    TCclay::Vector{T} | "t dt-1"
+    TCsilt::Vector{T} | "t dt-1"
+    TCsand::Vector{T} | "t dt-1"
+    TCsagg::Vector{T} | "t dt-1"
+    TClagg::Vector{T} | "t dt-1"
+    # Sediment flux in overland flow [ton dt-1]
+    olsed::Vector{T} | "t dt-1"
+    olclay::Vector{T} | "t dt-1"
+    olsilt::Vector{T} | "t dt-1"
+    olsand::Vector{T} | "t dt-1"
+    olsagg::Vector{T} | "t dt-1"
+    ollagg::Vector{T} | "t dt-1"
+    # Sediment reaching the river with overland flow [ton dt-1]
+    inlandsed::Vector{T} | "t dt-1"
+    inlandclay::Vector{T} | "t dt-1"
+    inlandsilt::Vector{T} | "t dt-1"
+    inlandsand::Vector{T} | "t dt-1"
+    inlandsagg::Vector{T} | "t dt-1"
+    inlandlagg::Vector{T} | "t dt-1"
 
 
     function OverlandFlowSediment{T}(args...) where {T}
@@ -683,7 +682,7 @@ end
     # number of cells
     n::Int | "-" | 0 | "none" | "none"
     # Timestep [s]
-    Δt::T | "s" | 0 | "none" | "none"
+    dt::T | "s" | 0 | "none" | "none"
     # River geometry (slope [-], length [m], width [m])
     sl::Vector{T} | "m"
     dl::Vector{T} | "m"
@@ -725,12 +724,12 @@ end
     # River discharge [m3/s]
     q_riv::Vector{T} | "m3 s-1"
     # Sediment input from land erosion [ton Δt⁻¹]
-    inlandclay::Vector{T} | "t Δt⁻¹"
-    inlandsilt::Vector{T} | "t Δt⁻¹"
-    inlandsand::Vector{T} | "t Δt⁻¹"
-    inlandsagg::Vector{T} | "t Δt⁻¹"
-    inlandlagg::Vector{T} | "t Δt⁻¹"
-    inlandsed::Vector{T} | "t Δt⁻¹"
+    inlandclay::Vector{T} | "t dt-1"
+    inlandsilt::Vector{T} | "t dt-1"
+    inlandsand::Vector{T} | "t dt-1"
+    inlandsagg::Vector{T} | "t dt-1"
+    inlandlagg::Vector{T} | "t dt-1"
+    inlandsed::Vector{T} | "t dt-1"
     # Sediment / particle left in the cell [ton]
     sedload::Vector{T} | "t"
     clayload::Vector{T} | "t"
@@ -740,21 +739,21 @@ end
     laggload::Vector{T} | "t"
     gravload::Vector{T} | "t"
     # Sediment / particle stored on the river bed after deposition [ton Δt⁻¹]
-    sedstore::Vector{T} | "t Δt⁻¹"
-    claystore::Vector{T} | "t Δt⁻¹"
-    siltstore::Vector{T} | "t Δt⁻¹"
-    sandstore::Vector{T} | "t Δt⁻¹"
-    saggstore::Vector{T} | "t Δt⁻¹"
-    laggstore::Vector{T} | "t Δt⁻¹"
-    gravstore::Vector{T} | "t Δt⁻¹"
+    sedstore::Vector{T} | "t dt-1"
+    claystore::Vector{T} | "t dt-1"
+    siltstore::Vector{T} | "t dt-1"
+    sandstore::Vector{T} | "t dt-1"
+    saggstore::Vector{T} | "t dt-1"
+    laggstore::Vector{T} | "t dt-1"
+    gravstore::Vector{T} | "t dt-1"
     # Sediment / particle flux [ton Δt⁻¹]
-    outsed::Vector{T} | "t Δt⁻¹"
-    outclay::Vector{T} | "t Δt⁻¹"
-    outsilt::Vector{T} | "t Δt⁻¹"
-    outsand::Vector{T} | "t Δt⁻¹"
-    outsagg::Vector{T} | "t Δt⁻¹"
-    outlagg::Vector{T} | "t Δt⁻¹"
-    outgrav::Vector{T} | "t Δt⁻¹"
+    outsed::Vector{T} | "t dt-1"
+    outclay::Vector{T} | "t dt-1"
+    outsilt::Vector{T} | "t dt-1"
+    outsand::Vector{T} | "t dt-1"
+    outsagg::Vector{T} | "t dt-1"
+    outlagg::Vector{T} | "t dt-1"
+    outgrav::Vector{T} | "t dt-1"
     # Total sediment concentrations (SSconc + Bedconc) [g/m3]
     Sedconc::Vector{T} | "g m-3"
     # Suspended load concentration [g/m3]
@@ -762,15 +761,15 @@ end
     # Bed load concentration [g/m3]
     Bedconc::Vector{T} | "g m-3"
     # River transport capacity
-    maxsed::Vector{T} | "t Δt⁻¹"
+    maxsed::Vector{T} | "t dt-1"
     # Eroded sediment (total, bank and bed)
-    erodsed::Vector{T} | "t Δt⁻¹"
-    erodsedbank::Vector{T} | "t Δt⁻¹"
-    erodsedbed::Vector{T} | "t Δt⁻¹"
+    erodsed::Vector{T} | "t dt-1"
+    erodsedbank::Vector{T} | "t dt-1"
+    erodsedbed::Vector{T} | "t dt-1"
     # Deposited sediment
-    depsed::Vector{T} | "t Δt⁻¹"
+    depsed::Vector{T} | "t dt-1"
     # Sediment in
-    insed::Vector{T} | "t Δt⁻¹"
+    insed::Vector{T} | "t dt-1"
     # Reservoir and lakes
     wbcover::Vector{T} | "-"
     wblocs::Vector{T} | "-"
@@ -789,7 +788,7 @@ function initialize_riversed(nc, config, riverwidth, riverlength, inds_riv)
     nriv = length(inds_riv)
     # River flow transport capacity method: ["bagnold", "engelund", "yang", "kodatie", "molinas"]
     tcmethodriv = get(config.model, "rivtransportmethod", "bagnold")::String
-    Δt = Second(config.timestepsecs)
+    dt = Second(config.timestepsecs)
     # Reservoir / lakes
     do_reservoirs = get(config.model, "doreservoir", false)::Bool
     do_lakes = get(config.model, "dolake", false)::Bool
@@ -1053,7 +1052,7 @@ function initialize_riversed(nc, config, riverwidth, riverlength, inds_riv)
 
     rs = RiverSediment(
         n = nriv,
-        Δt = Float(Δt.value),
+        dt = Float(dt.value),
         # Parameters
         sl = riverslope,
         dl = riverlength,
@@ -1211,7 +1210,7 @@ function update(rs::RiverSediment, network, config)
             maxsed = rs.ak[v] * vmean^rs.bk[v] * rs.h_riv[v]^rs.ck[v] * rs.sl[v]^rs.dk[v]
             # Transport capacity [tons/m3]
             maxsed =
-                ifelse(rs.q_riv[v] > 0.0, maxsed * rs.width[v] / (rs.q_riv[v] * rs.Δt), 0.0)
+                ifelse(rs.q_riv[v] > 0.0, maxsed * rs.width[v] / (rs.q_riv[v] * rs.dt), 0.0)
         elseif tcmethod == "yang"
             ws = 411 * rs.d50[v]^2 / 3600
             vshear = (9.81 * hydrad * rs.sl[v])^0.5
@@ -1272,7 +1271,7 @@ function update(rs::RiverSediment, network, config)
         # 1285 g/L: boundary between streamflow and debris flow (Costa, 1988)
         maxsed = min(maxsed, 1.285)
         # Transport capacity [ton]
-        maxsed = maxsed * (rs.h_riv[v] * rs.width[v] * rs.dl[v] + rs.q_riv[v] * rs.Δt)
+        maxsed = maxsed * (rs.h_riv[v] * rs.width[v] * rs.dl[v] + rs.q_riv[v] * rs.dt)
         rs.maxsed[v] = maxsed
 
         ### River erosion ###
@@ -1318,7 +1317,7 @@ function update(rs::RiverSediment, network, config)
             #(assuming only one bank is eroding)
             Tex = max(TEffbank - rs.TCrbank[v], 0.0)
             # 1.4 is bank default bulk density
-            ERbank = max(0.0, rs.kdbank[v] * Tex * rs.dl[v] * rs.h_riv[v] * 1.4 * rs.Δt)
+            ERbank = max(0.0, rs.kdbank[v] * Tex * rs.dl[v] * rs.h_riv[v] * 1.4 * rs.dt)
             # 1.5 is bed default bulk density
             ERbed = max(
                 0.0,
@@ -1327,7 +1326,7 @@ function update(rs::RiverSediment, network, config)
                 rs.dl[v] *
                 rs.width[v] *
                 1.5 *
-                rs.Δt,
+                rs.dt,
             )
             # Relative potential erosion rates of the bed and the bank [-]
             RTEbank = ifelse(ERbank + ERbed > 0.0, ERbank / (ERbank + ERbed), 0.0)
@@ -1508,7 +1507,7 @@ function update(rs::RiverSediment, network, config)
         # Sediment transported out of the cell during the timestep [ton]
         # 0 in case all sediment are deposited in the cell
         # Reduce the fraction so that there is still some sediment staying in the river cell
-        fwaterout = min(rs.q_riv[v] * rs.Δt / (rs.h_riv[v] * rs.width[v] * rs.dl[v]), 1.0)
+        fwaterout = min(rs.q_riv[v] * rs.dt / (rs.h_riv[v] * rs.width[v] * rs.dl[v]), 1.0)
         rs.outsed[v] = fwaterout * (insed + erodsed - depsed)
         rs.outclay[v] = fwaterout * (inclay + erodclay - depclay)
         rs.outsilt[v] = fwaterout * (insilt + erodsilt - depsilt)
@@ -1528,7 +1527,7 @@ function update(rs::RiverSediment, network, config)
 
         ### Concentrations and suspended sediments ###
         # Conversion from load [ton] to concentration for rivers [mg/L]
-        toconc = ifelse(rs.q_riv[v] > 0.0, 1e6 / (rs.q_riv[v] * rs.Δt), 0.0)
+        toconc = ifelse(rs.q_riv[v] > 0.0, 1e6 / (rs.q_riv[v] * rs.dt), 0.0)
         rs.Sedconc[v] = rs.outsed[v] * toconc
 
         # Differentiation of bed and suspended load using Rouse number for suspension
diff --git a/src/sediment_model.jl b/src/sediment_model.jl
index 4dcbd5a1c..e2229a1bd 100644
--- a/src/sediment_model.jl
+++ b/src/sediment_model.jl
@@ -14,7 +14,7 @@ function initialize_sediment_model(config::Config)
 
     reader = prepare_reader(config)
     clock = Clock(config, reader)
-    Δt = clock.Δt
+    dt = clock.dt
 
     do_river = get(config.model, "runrivermodel", false)::Bool
 
@@ -95,9 +95,9 @@ function initialize_sediment_model(config::Config)
     # River processes
 
     # the indices of the river cells in the land(+river) cell vector
-    βₗ =
+    landslope =
         ncread(nc, config, "lateral.land.slope"; optional = false, sel = inds, type = Float)
-    clamp!(βₗ, 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, βₗ, 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 53167f887..57563b7b6 100644
--- a/src/utils.jl
+++ b/src/utils.jl
@@ -665,12 +665,12 @@ concept `SBM` (at index `i`) based on hydraulic conductivity profile `ksat_profi
 """
 function hydraulic_conductivity_at_depth(sbm::SBM, z, i, n, ksat_profile)
     if ksat_profile == "exponential"
-        kv_z = sbm.kvfrac[i][n] * sbm.kv₀[i] * exp(-sbm.f[i] * z)
+        kv_z = sbm.kvfrac[i][n] * sbm.kv_0[i] * exp(-sbm.f[i] * z)
     elseif ksat_profile == "exponential_constant"
         if z < sbm.z_exp[i]
-            kv_z = sbm.kvfrac[i][n] * sbm.kv₀[i] * exp(-sbm.f[i] * z)
+            kv_z = sbm.kvfrac[i][n] * sbm.kv_0[i] * exp(-sbm.f[i] * z)
         else
-            kv_z = sbm.kvfrac[i][n] * sbm.kv₀[i] * exp(-sbm.f[i] * sbm.z_exp[i])
+            kv_z = sbm.kvfrac[i][n] * sbm.kv_0[i] * exp(-sbm.f[i] * sbm.z_exp[i])
         end
     elseif ksat_profile == "layered"
         kv_z = sbm.kvfrac[i][n] * sbm.kv[i][n]
@@ -697,7 +697,7 @@ table depth `z` [mm] and hydraulic conductivity profile `ksat_profile`.
 function kh_layered_profile(sbm::SBM, khfrac, i, ksat_profile)
 
     m = sbm.nlayers[i]
-    t_factor = (tosecond(basetimestep) / sbm.Δt)
+    t_factor = (tosecond(basetimestep) / sbm.dt)
     if (sbm.soilthickness[i] - sbm.zi[i]) > 0.0
         transmissivity = 0.0
         sumlayers = @view sbm.sumlayers[i][2:end]
@@ -768,10 +768,10 @@ end
 function initialize_lateralssf_exp!(ssf::LateralSSF)
     for i in eachindex(ssf.ssf)
         ssf.ssfmax[i] =
-            ((ssf.kh₀[i] * ssf.βₗ[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₀[i] * ssf.βₗ[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
@@ -780,26 +780,26 @@ end
 "Initialize lateral subsurface variables `ssf` and `ssfmax` with `ksat_profile` `exponential_constant`"
 function initialize_lateralssf_exp_const!(ssf::LateralSSF)
     ssf_constant = @. ssf.khfrac *
-       ssf.kh₀ *
+       ssf.kh_0 *
        exp(-ssf.f * ssf.z_exp) *
-       ssf.βₗ *
+       ssf.slope *
        (ssf.soilthickness - ssf.z_exp)
     for i in eachindex(ssf.ssf)
         ssf.ssfmax[i] =
-            ((ssf.khfrac[i] * ssf.kh₀[i] * ssf.βₗ[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₀[i] * ssf.βₗ[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]
         else
             ssf.ssf[i] =
-                ssf.kh₀[i] *
+                ssf.kh_0[i] *
                 exp(-ssf.f[i] * ssf.zi[i]) *
-                ssf.βₗ[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.βₗ[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.βₗ[i]
+        ssf.ssfmax[i] = kh_max * ssf.slope[i]
     end
 end
diff --git a/src/vertical_process.jl b/src/vertical_process.jl
index 971cc6221..9801b1e2e 100644
--- a/src/vertical_process.jl
+++ b/src/vertical_process.jl
@@ -120,7 +120,7 @@ function rainfall_interception_modrut(
 end
 
 """
-    acttransp_unsat_sbm(rootingdepth, ustorelayerdepth, sumlayer, restpotevap, sum_actevapustore, c, usl, θₛ, θᵣ, hb, ust::Bool = false)
+    acttransp_unsat_sbm(rootingdepth, ustorelayerdepth, sumlayer, restpotevap, sum_actevapustore, c, usl, theta_s, theta_r, hb, ust::Bool = false)
 
 Compute actual transpiration for unsaturated zone.
 If `ust` is `true`, the whole unsaturated store is available for transpiration.
@@ -133,8 +133,8 @@ If `ust` is `true`, the whole unsaturated store is available for transpiration.
 - `sum_actevapustore` (cumulative actual transpiration (more than one unsaturated layers))
 - `c` (Brooks-Corey coefficient)
 - `usl` (thickness of unsaturated zone)
-- `θₛ`
-- `θᵣ`
+- `theta_s`
+- `theta_r`
 - `hb` (air entry pressure)
 - `ust`
 
@@ -151,8 +151,8 @@ function acttransp_unsat_sbm(
     sum_actevapustore,
     c,
     usl,
-    θₛ,
-    θᵣ,
+    theta_s,
+    theta_r,
     hb,
     ust::Bool = false,
 )
@@ -186,7 +186,7 @@ function acttransp_unsat_sbm(
         vwc = 0.0
     end
     vwc = max(vwc, 0.0000001)
-    head = hb / (pow(((vwc) / (θₛ - θᵣ)), (1.0 / par_lambda)))  # Note that in the original formula, thetaR is extracted from vwc, but thetaR is not part of the numerical vwc calculation
+    head = hb / (pow(((vwc) / (theta_s - theta_r)), (1.0 / par_lambda)))  # Note that in the original formula, thetaR is extracted from vwc, but thetaR is not part of the numerical vwc calculation
     head = max(head, hb)
 
     # Transform h to a reduction coefficient value according to Feddes et al. (1978).
@@ -293,7 +293,7 @@ function unsatzone_flow_layer(usd, kv_z, l_sat, c)
 end
 
 """
-    unsatzone_flow_sbm(ustorelayerdepth, soilwatercapacity, satwaterdepth, kv_z, usl, θₛ, θᵣ)
+    unsatzone_flow_sbm(ustorelayerdepth, soilwatercapacity, satwaterdepth, kv_z, usl, theta_s, theta_r)
 
 The transfer of water from the unsaturated store `ustorelayerdepth` to the saturated store `satwaterdepth`
 is controlled by the vertical saturated hydraulic conductivity `kv_z` at the water table and the ratio between
@@ -307,15 +307,15 @@ function unsatzone_flow_sbm(
     satwaterdepth,
     kv_z,
     usl,
-    θₛ,
-    θᵣ,
+    theta_s,
+    theta_r,
 )
 
     sd = soilwatercapacity - satwaterdepth
     if sd <= 0.00001
         ast = 0.0
     else
-        st = kv_z * min(ustorelayerdepth, usl * (θₛ - θᵣ)) / sd
+        st = kv_z * min(ustorelayerdepth, usl * (theta_s - theta_r)) / sd
         ast = min(st, ustorelayerdepth)
         ustorelayerdepth = ustorelayerdepth - ast
     end
@@ -402,7 +402,7 @@ function snowpack_hbv(
 end
 
 """
-    glacier_hbv(glacierfrac, glacierstore, snow, temperature, tt, cfmax, g_sifrac, Δt)
+    glacier_hbv(glacierfrac, glacierstore, snow, temperature, tt, cfmax, g_sifrac, dt)
 
 HBV-light type of glacier modelling.
 First, a fraction of the snowpack is converted into ice using the HBV-light
@@ -418,7 +418,7 @@ occurs if the snow cover < 10 mm.
 - `tt` temperature threshold for ice melting [°C]
 - `cfmax` ice degree-day factor in [mm/(°C/day)]
 - `g_sifrac` fraction of the snow turned into ice [-]
-- `Δt` model timestep [s]
+- `dt` model timestep [s]
 
 # Output
 - `snow`
@@ -427,14 +427,14 @@ occurs if the snow cover < 10 mm.
 - `glaciermelt`
 
 """
-function glacier_hbv(glacierfrac, glacierstore, snow, temperature, tt, cfmax, g_sifrac, Δt)
+function glacier_hbv(glacierfrac, glacierstore, snow, temperature, tt, cfmax, g_sifrac, dt)
 
     # Fraction of the snow transformed into ice (HBV-light model)
     snow2glacier = g_sifrac * snow
     snow2glacier = glacierfrac > 0.0 ? snow2glacier : 0.0
 
     # Max conversion to 8mm/day
-    snow2glacier = min(snow2glacier, 8.0 * (Δt / basetimestep))
+    snow2glacier = min(snow2glacier, 8.0 * (dt / basetimestep))
 
     snow = snow - (snow2glacier * glacierfrac)
     glacierstore = glacierstore + snow2glacier


Test patch:
diff --git a/server/test/client.jl b/server/test/client.jl
index 9c06e66fd..d5eaa70a9 100644
--- a/server/test/client.jl
+++ b/server/test/client.jl
@@ -36,7 +36,7 @@ end
     @test request((fn = "get_output_item_count",)) == Dict("output_item_count" => 193)
     to_check = [
         "vertical.nlayers",
-        "vertical.θᵣ",
+        "vertical.theta_r",
         "lateral.river.q",
         "lateral.river.reservoir.outflow",
     ]
@@ -52,7 +52,7 @@ vwc_1_size = 0
     @test request((fn = "get_var_itemsize", name = "lateral.subsurface.ssf")) ==
           Dict("var_itemsize" => sizeof(Wflow.Float))
     @test request((fn = "get_var_type", name = "vertical.n"))["status"] == "ERROR"
-    @test request((fn = "get_var_units", name = "vertical.θₛ")) == Dict("var_units" => "-")
+    @test request((fn = "get_var_units", name = "vertical.theta_s")) == Dict("var_units" => "-")
     @test request((fn = "get_var_location", name = "lateral.river.q")) ==
           Dict("var_location" => "node")
     zi_nbytes = request((fn = "get_var_nbytes", name = "vertical.zi"))["var_nbytes"]
@@ -67,7 +67,7 @@ vwc_1_size = 0
     @test request((fn = "get_var_grid", name = "lateral.river.h")) == Dict("var_grid" => 3)
     msg = (fn = "get_value", name = "vertical.zi", dest = fill(0.0, zi_size))
     @test mean(request(msg)["value"]) ≈ 277.8362161218515
-    msg = (fn = "get_value_ptr", name = "vertical.θₛ")
+    msg = (fn = "get_value_ptr", name = "vertical.theta_s")
     @test mean(request(msg)["value_ptr"]) ≈ 0.4409211971535584
     msg = (
         fn = "get_value_at_indices",
diff --git a/test/bmi.jl b/test/bmi.jl
index dfd1a5577..db8c95462 100644
--- a/test/bmi.jl
+++ b/test/bmi.jl
@@ -24,7 +24,7 @@ tomlpath = joinpath(@__DIR__, "sbm_config.toml")
             @test BMI.get_output_item_count(model) == 193
             to_check = [
                 "vertical.nlayers",
-                "vertical.θᵣ",
+                "vertical.theta_r",
                 "lateral.river.q",
                 "lateral.river.reservoir.outflow",
             ]
@@ -35,12 +35,12 @@ tomlpath = joinpath(@__DIR__, "sbm_config.toml")
         end
 
         @testset "variable information functions" begin
-            @test BMI.get_var_grid(model, "vertical.θₛ") == 6
+            @test BMI.get_var_grid(model, "vertical.theta_s") == 6
             @test BMI.get_var_grid(model, "lateral.river.h") == 3
             @test BMI.get_var_grid(model, "lateral.river.reservoir.inflow") == 0
             @test_throws ErrorException BMI.get_var_grid(model, "lateral.river.lake.volume")
             @test BMI.get_var_type(model, "lateral.river.reservoir.inflow") == "$Float"
-            @test BMI.get_var_units(model, "vertical.θₛ") == "-"
+            @test BMI.get_var_units(model, "vertical.theta_s") == "-"
             @test BMI.get_var_itemsize(model, "lateral.subsurface.ssf") == sizeof(Float)
             @test BMI.get_var_nbytes(model, "lateral.river.q") ==
                   length(model.lateral.river.q) * sizeof(Float)
diff --git a/test/groundwater.jl b/test/groundwater.jl
index 528d472e0..13d529ed1 100644
--- a/test/groundwater.jl
+++ b/test/groundwater.jl
@@ -4,14 +4,14 @@ function initial_head(x)
 end
 
 """
-    transient_aquifer_1d(x, time, conductivity, specific_yield, aquifer_length, Β)
+    transient_aquifer_1d(x, time, conductivity, specific_yield, aquifer_length, beta)
 
 Non-steady flow in an unconfined rectangular aquifer, with Dirichlet h(0, t) = 0
 on the left edge, and a Neumann Boundary Condition (dh/dx = 0) on the right.
 """
-function transient_aquifer_1d(x, time, conductivity, specific_yield, aquifer_length, Β)
+function transient_aquifer_1d(x, time, conductivity, specific_yield, aquifer_length, beta)
     initial_head(x) / 1.0 +
-    (Β * conductivity * initial_head(aquifer_length) * time) /
+    (beta * conductivity * initial_head(aquifer_length) * time) /
     (specific_yield * aquifer_length * aquifer_length)
 end
 
@@ -23,7 +23,7 @@ Non-steady flow in a confined aquifer, using the well function of Theis.
 """
 function drawdown_theis(distance, time, discharge, transmissivity, storativity)
     u = (storativity * distance^2) / (4 * transmissivity * time)
-    return discharge / (4 * π * transmissivity) * expint(u)
+    return discharge / (4 * pi * transmissivity) * expint(u)
 end
 
 
@@ -31,10 +31,10 @@ function homogenous_aquifer(nrow, ncol)
     shape = (nrow, ncol)
     # Domain, geometry
     domain = ones(Bool, shape)
-    Δx = fill(10.0, ncol)
-    Δy = fill(10.0, nrow)
+    dx = fill(10.0, ncol)
+    dy = fill(10.0, nrow)
     indices, reverse_indices = Wflow.active_indices(domain, false)
-    connectivity = Wflow.Connectivity(indices, reverse_indices, Δx, Δy)
+    connectivity = Wflow.Connectivity(indices, reverse_indices, dx, dy)
     ncell = connectivity.ncell
 
     conf_aqf = Wflow.ConfinedAquifer(
@@ -65,8 +65,8 @@ end
     ncol = 2
     nrow = 3
     shape = (ncol, nrow)
-    Δx = [10.0, 20.0]
-    Δy = [5.0, 15.0, 25.0]
+    dx = [10.0, 20.0]
+    dy = [5.0, 15.0, 25.0]
     collect_connections(con, cell_id) =
         [con.rowval[nzi] for nzi in Wflow.connections(con, cell_id)]
 
@@ -74,15 +74,15 @@ end
         @testset "connection_geometry: y" begin
             I = CartesianIndex(1, 1)
             J = CartesianIndex(2, 1)
-            @test Wflow.connection_geometry(I, J, Δx, Δy) == (2.5, 7.5, 10.0)
-            @test_throws Exception Wflow.connection_geometry(I, I, Δx, Δy)
+            @test Wflow.connection_geometry(I, J, dx, dy) == (2.5, 7.5, 10.0)
+            @test_throws Exception Wflow.connection_geometry(I, I, dx, dy)
         end
 
         @testset "connection_geometry: x" begin
             I = CartesianIndex(1, 1)
             J = CartesianIndex(1, 2)
-            @test Wflow.connection_geometry(I, J, Δx, Δy) == (5.0, 10.0, 5.0)
-            @test_throws Exception Wflow.connection_geometry(I, I, Δx, Δy)
+            @test Wflow.connection_geometry(I, J, dx, dy) == (5.0, 10.0, 5.0)
+            @test_throws Exception Wflow.connection_geometry(I, I, dx, dy)
         end
 
         @testset "Connectivity 1D(x)" begin
@@ -91,7 +91,7 @@ end
             # +---+---+
             domain = ones(Bool, (1, 2))
             indices, reverse_indices = Wflow.active_indices(domain, false)
-            conn = Wflow.Connectivity(indices, reverse_indices, Δx, [5.0])
+            conn = Wflow.Connectivity(indices, reverse_indices, dx, [5.0])
             @test conn.ncell == 2
             @test conn.nconnection == 2
             @test conn.length1 == [5.0, 10.0]
@@ -113,7 +113,7 @@ end
             # +---+
             domain = ones(Bool, (3, 1))
             indices, reverse_indices = Wflow.active_indices(domain, false)
-            conn = Wflow.Connectivity(indices, reverse_indices, [10.0], Δy)
+            conn = Wflow.Connectivity(indices, reverse_indices, [10.0], dy)
             @test conn.ncell == 3
             @test conn.nconnection == 4
             @test conn.length1 == [2.5, 7.5, 7.5, 12.5]
@@ -136,7 +136,7 @@ end
             # +---+---+
             domain = ones(Bool, (nrow, ncol))
             indices, reverse_indices = Wflow.active_indices(domain, false)
-            conn = Wflow.Connectivity(indices, reverse_indices, Δx, Δy)
+            conn = Wflow.Connectivity(indices, reverse_indices, dx, dy)
             @test conn.ncell == 6
             @test conn.nconnection == 14
             @test conn.colptr == [1, 3, 6, 8, 10, 13, 15]
@@ -160,7 +160,7 @@ end
             domain = ones(Bool, (nrow, ncol))
             domain[2, 1] = false
             indices, reverse_indices = Wflow.active_indices(domain, false)
-            conn = Wflow.Connectivity(indices, reverse_indices, Δx, Δy)
+            conn = Wflow.Connectivity(indices, reverse_indices, dx, dy)
             @test conn.ncell == 5
             @test conn.nconnection == 8
             @test conn.colptr == [1, 2, 3, 5, 7, 9]
@@ -351,9 +351,9 @@ end
         gwf.aquifer.head[gwf.constanthead.index] .= gwf.constanthead.head
 
         Q = zeros(3)
-        Δt = 0.25 # days
+        dt = 0.25 # days
         for _ = 1:50
-            Wflow.update(gwf, Q, Δt, conductivity_profile)
+            Wflow.update(gwf, Q, dt, conductivity_profile)
         end
 
         @test gwf.aquifer.head ≈ [2.0, 3.0, 4.0]
@@ -373,9 +373,9 @@ end
         gwf.aquifer.head[gwf.constanthead.index] .= gwf.constanthead.head
 
         Q = zeros(3)
-        Δt = 0.25 # days
+        dt = 0.25 # days
         for _ = 1:50
-            Wflow.update(gwf, Q, Δt, conductivity_profile)
+            Wflow.update(gwf, Q, dt, conductivity_profile)
         end
 
         @test gwf.aquifer.head ≈ [2.0, 3.0, 4.0]
@@ -390,17 +390,17 @@ end
         bottom = 0.0
         specific_yield = 0.15
         cellsize = 500.0
-        Β = 1.12
+        beta = 1.12
         aquifer_length = cellsize * ncol
         gwf_f = 3.0
         conductivity_profile = "uniform"
 
         # Domain, geometry
         domain = ones(Bool, shape)
-        Δx = fill(cellsize, ncol)
-        Δy = fill(cellsize, nrow)
+        dx = fill(cellsize, ncol)
+        dy = fill(cellsize, nrow)
         indices, reverse_indices = Wflow.active_indices(domain, false)
-        connectivity = Wflow.Connectivity(indices, reverse_indices, Δx, Δy)
+        connectivity = Wflow.Connectivity(indices, reverse_indices, dx, dy)
         ncell = connectivity.ncell
         xc = collect(range(0.0, stop = aquifer_length - cellsize, step = cellsize))
         aquifer = Wflow.UnconfinedAquifer(
@@ -422,22 +422,22 @@ end
             Wflow.AquiferBoundaryCondition[],
         )
 
-        Δt = Wflow.stable_timestep(gwf.aquifer, conductivity_profile)
+        dt = Wflow.stable_timestep(gwf.aquifer, conductivity_profile)
         Q = zeros(ncell)
         time = 20.0
-        nstep = Int(ceil(time / Δt))
-        time = nstep * Δt
+        nstep = Int(ceil(time / dt))
+        time = nstep * dt
 
         for i = 1:nstep
-            Wflow.update(gwf, Q, Δt, conductivity_profile)
+            Wflow.update(gwf, Q, dt, conductivity_profile)
             # Gradient dh/dx is positive, all flow to the left
             @test all(diff(gwf.aquifer.head) .> 0.0)
         end
 
-        ϕ_analytical = [
-            transient_aquifer_1d(x, time, conductivity, specific_yield, aquifer_length, Β) for x in xc
+        head_analytical = [
+            transient_aquifer_1d(x, time, conductivity, specific_yield, aquifer_length, beta) for x in xc
         ]
-        difference = gwf.aquifer.head .- ϕ_analytical
+        difference = gwf.aquifer.head .- head_analytical
         # @test all(difference .< ?)  #TODO
     end
 
@@ -450,17 +450,17 @@ end
         bottom = 0.0
         specific_yield = 0.15
         cellsize = 500.0
-        Β = 1.12
+        beta = 1.12
         aquifer_length = cellsize * ncol
         gwf_f = 3.0
         conductivity_profile = "exponential"
 
         # Domain, geometry
         domain = ones(Bool, shape)
-        Δx = fill(cellsize, ncol)
-        Δy = fill(cellsize, nrow)
+        dx = fill(cellsize, ncol)
+        dy = fill(cellsize, nrow)
         indices, reverse_indices = Wflow.active_indices(domain, false)
-        connectivity = Wflow.Connectivity(indices, reverse_indices, Δx, Δy)
+        connectivity = Wflow.Connectivity(indices, reverse_indices, dx, dy)
         ncell = connectivity.ncell
         xc = collect(range(0.0, stop = aquifer_length - cellsize, step = cellsize))
         aquifer = Wflow.UnconfinedAquifer(
@@ -482,22 +482,22 @@ end
             Wflow.AquiferBoundaryCondition[],
         )
 
-        Δt = Wflow.stable_timestep(gwf.aquifer, conductivity_profile)
+        dt = Wflow.stable_timestep(gwf.aquifer, conductivity_profile)
         Q = zeros(ncell)
         time = 20.0
-        nstep = Int(ceil(time / Δt))
-        time = nstep * Δt
+        nstep = Int(ceil(time / dt))
+        time = nstep * dt
 
         for i = 1:nstep
-            Wflow.update(gwf, Q, Δt, conductivity_profile)
+            Wflow.update(gwf, Q, dt, conductivity_profile)
             # Gradient dh/dx is positive, all flow to the left
             @test all(diff(gwf.aquifer.head) .> 0.0)
         end
 
-        ϕ_analytical = [
-            transient_aquifer_1d(x, time, conductivity, specific_yield, aquifer_length, Β) for x in xc
+        head_analytical = [
+            transient_aquifer_1d(x, time, conductivity, specific_yield, aquifer_length, beta) for x in xc
         ]
-        difference = gwf.aquifer.head .- ϕ_analytical
+        difference = gwf.aquifer.head .- head_analytical
         # @test all(difference .< ?)  #TODO
     end
 
@@ -521,10 +521,10 @@ end
 
         # Domain, geometry
         domain = ones(Bool, shape)
-        Δx = fill(cellsize, ncol)
-        Δy = fill(cellsize, nrow)
+        dx = fill(cellsize, ncol)
+        dy = fill(cellsize, nrow)
         indices, reverse_indices = Wflow.active_indices(domain, false)
-        connectivity = Wflow.Connectivity(indices, reverse_indices, Δx, Δy)
+        connectivity = Wflow.Connectivity(indices, reverse_indices, dx, dy)
         ncell = connectivity.ncell
         aquifer = Wflow.ConfinedAquifer(
             fill(startinghead, ncell),
@@ -544,29 +544,29 @@ end
         well = Wflow.Well([discharge], [0.0], [reverse_indices[wellrow, wellrow]])
         gwf = Wflow.GroundwaterFlow(aquifer, connectivity, constanthead, [well])
 
-        Δt = Wflow.stable_timestep(gwf.aquifer, conductivity_profile)
+        dt = Wflow.stable_timestep(gwf.aquifer, conductivity_profile)
         Q = zeros(ncell)
         time = 20.0
-        nstep = Int(ceil(time / Δt))
-        time = nstep * Δt
+        nstep = Int(ceil(time / dt))
+        time = nstep * dt
 
         for i = 1:nstep
-            Wflow.update(gwf, Q, Δt, conductivity_profile)
+            Wflow.update(gwf, Q, dt, conductivity_profile)
         end
 
         # test for symmetry on x and y axes
-        ϕ = reshape(gwf.aquifer.head, shape)
-        @test ϕ[1:halfnrow, :] ≈ ϕ[end:-1:halfnrow+2, :]
-        @test ϕ[:, 1:halfnrow] ≈ ϕ[:, 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))
-        ϕ_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 = ϕ[1:halfnrow-1, halfnrow] - ϕ_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 851429dd7..3e9eb41ed 100644
--- a/test/horizontal_process.jl
+++ b/test/horizontal_process.jl
@@ -1,4 +1,4 @@
-const Δt_sec = 86400.0
+const dt_sec = 86400.0
 const ldd_mv = 255
 
 # read the staticmaps into memory
@@ -29,14 +29,14 @@ sink = toposort[end]
 
 # calculate parameters of kinematic wave
 const q = 0.000001
-const β = 0.6
-const AlpPow = (2.0 / 3.0) * β
-AlpTermR = (N ./ sqrt.(slope)) .^ β
+const beta = 0.6
+const AlpPow = (2.0 / 3.0) * beta
+AlpTermR = (N ./ sqrt.(slope)) .^ beta
 P = Bw + (2.0 * waterlevel)
-α = AlpTermR .* P .^ AlpPow
+alpha = AlpTermR .* P .^ AlpPow
 
 Q = zeros(n)
-Q = Wflow.kin_wave!(Q, graph, toposort, Qold, q, α, β, DCL, Δt_sec)
+Q = Wflow.kin_wave!(Q, graph, toposort, Qold, q, alpha, beta, DCL, dt_sec)
 
 @testset "flow rate" begin
     @test sum(Q) ≈ 2.957806043289641e6
@@ -175,7 +175,7 @@ end
     )
 
     alpha = 0.7
-    Δt = 1.0
+    dt = 1.0
     h_thresh = 1.0e-03
     froude_limit = true
     h_init = zeros(n - 1)
@@ -187,9 +187,9 @@ end
         active_n = collect(1:n-1),
         active_e = collect(1:_ne),
         g = 9.80665,
-        α = alpha,
+        alpha = alpha,
         h_thresh = h_thresh,
-        Δt = Δt,
+        dt = dt,
         q0 = zeros(_ne),
         q = zeros(_ne),
         q_av = zeros(_ne),
@@ -198,9 +198,9 @@ end
         mannings_n_sq = mannings_n_sq,
         mannings_n = n_river,
         h = h_init,
-        η_max = zeros(_ne),
-        η_src = zeros(_ne),
-        η_dst = zeros(_ne),
+        zs_max = zeros(_ne),
+        zs_src = zeros(_ne),
+        zs_dst = zeros(_ne),
         hf = zeros(_ne),
         h_av = zeros(n),
         width = width,
@@ -227,14 +227,14 @@ end
     )
 
     # run until steady state is reached
-    ϵ = 1.0e-12
+    epsilon = 1.0e-12
     while true
         sw_river.inwater[1] = 20.0
         h0 = mean(sw_river.h)
-        Δt = Wflow.stable_timestep(sw_river)
-        Wflow.shallowwater_river_update(sw_river, network, Δt, 0.0, true)
+        dt = Wflow.stable_timestep(sw_river)
+        Wflow.shallowwater_river_update(sw_river, network, dt, 0.0, true)
         d = abs(h0 - mean(sw_river.h))
-        if d <= ϵ
+        if d <= epsilon
             break
         end
     end
diff --git a/test/io.jl b/test/io.jl
index f6fbfc329..a033230a0 100644
--- a/test/io.jl
+++ b/test/io.jl
@@ -24,14 +24,14 @@ config = Wflow.Config(tomlpath)
     @test config.output isa Wflow.Config
     @test collect(keys(config.output)) == ["lateral", "vertical", "path"]
 
-    # theta_s can also be provided under the alias θₛ
-    @test Wflow.get_alias(config.input.vertical, "theta_s", "θₛ", nothing) == "thetaS"
+    # theta_s can also be provided under the alias theta_s
+    @test Wflow.get_alias(config.input.vertical, "theta_s", "theta_s", nothing) == "thetaS"
     val = pop!(config.input.vertical, "theta_s")
-    config.input.vertical["θₛ"] = val
-    @test Wflow.get_alias(config.input.vertical, "theta_s", "θₛ", nothing) == "thetaS"
+    config.input.vertical["theta_s"] = val
+    @test Wflow.get_alias(config.input.vertical, "theta_s", "theta_s", nothing) == "thetaS"
 
     # modifiers can also be applied
-    kvconf = Wflow.get_alias(config.input.vertical, "kv_0", "kv₀", nothing)
+    kvconf = Wflow.get_alias(config.input.vertical, "kv_0", "kv_0", nothing)
     @test kvconf isa Wflow.Config
     ncname, modifier = Wflow.ncvar_name_modifier(kvconf, config = config)
     @test ncname === "KsatVer"
@@ -73,7 +73,7 @@ end
 
     @test clock.time == DateTimeProlepticGregorian(2000, 1, 1)
     @test clock.iteration == 0
-    @test clock.Δt == Second(Day(1))
+    @test clock.dt == Second(Day(1))
     # test that the missing keys have been added to the config
     @test config.starttime == DateTime(2000, 1, 1)
     @test config.endtime == DateTime(2001, 1, 1)
@@ -88,7 +88,7 @@ end
     clock = Wflow.Clock(config, reader)
     @test clock.time == DateTimeStandard(2003, 4, 5)
     @test clock.iteration == 0
-    @test clock.Δt == Second(Hour(1))
+    @test clock.dt == Second(Hour(1))
 
     close(ds)
     config = Wflow.Config(tomlpath)  # restore the config
@@ -97,58 +97,58 @@ end
 @testset "Clock{DateTimeStandard}" begin
     # 29 days in this February due to leap year
     starttime = DateTimeStandard(2000, 2, 28)
-    Δt = Day(1)
-    clock = Wflow.Clock(starttime, 0, Second(Δt))
+    dt = Day(1)
+    clock = Wflow.Clock(starttime, 0, Second(dt))
 
     Wflow.advance!(clock)
     Wflow.advance!(clock)
     @test clock.time == DateTimeStandard(2000, 3, 1)
     @test clock.iteration == 2
-    @test clock.Δt == Δt
+    @test clock.dt == dt
 
     Wflow.rewind!(clock)
     @test clock.time == DateTimeStandard(2000, 2, 29)
     @test clock.iteration == 1
-    @test clock.Δt == Δt
+    @test clock.dt == dt
 
     config = Wflow.Config(
-        Dict("starttime" => starttime, "timestepsecs" => Dates.value(Second(Δt))),
+        Dict("starttime" => starttime, "timestepsecs" => Dates.value(Second(dt))),
     )
     Wflow.reset_clock!(clock, config)
     @test clock.time == starttime
     @test clock.iteration == 0
-    @test clock.Δt == Δt
+    @test clock.dt == dt
 end
 
 @testset "Clock{DateTime360Day}" begin
     # 30 days in each month
     starttime = DateTime360Day(2000, 2, 29)
-    Δt = Day(1)
-    clock = Wflow.Clock(starttime, 0, Second(Δt))
+    dt = Day(1)
+    clock = Wflow.Clock(starttime, 0, Second(dt))
 
     Wflow.advance!(clock)
     Wflow.advance!(clock)
     @test clock.time == DateTime360Day(2000, 3, 1)
     @test clock.iteration == 2
-    @test clock.Δt == Δt
+    @test clock.dt == dt
 
     Wflow.rewind!(clock)
     @test clock.time == DateTime360Day(2000, 2, 30)
     @test clock.iteration == 1
-    @test clock.Δt == Δt
+    @test clock.dt == dt
 
     config = Wflow.Config(
         Dict(
             "starttime" => "2020-02-29",
             "calendar" => "360_day",
-            "timestepsecs" => Dates.value(Second(Δt)),
+            "timestepsecs" => Dates.value(Second(dt)),
         ),
     )
     Wflow.reset_clock!(clock, config)
     @test clock.time isa DateTime360Day
     @test string(clock.time) == "2020-02-29T00:00:00"
     @test clock.iteration == 0
-    @test clock.Δt == Δt
+    @test clock.dt == dt
 end
 
 @testset "CFTime" begin
@@ -436,8 +436,8 @@ end
     @test clock.time == DateTimeNoLeap(2000, 1, 1)
 
     starttime = DateTimeNoLeap(2000, 2, 28)
-    Δt = Day(1)
-    clock = Wflow.Clock(starttime, 0, Second(Δt))
+    dt = Day(1)
+    clock = Wflow.Clock(starttime, 0, Second(dt))
     Wflow.advance!(clock)
     @test clock.time == DateTimeNoLeap(2000, 3, 1)
 end
diff --git a/test/reservoir_lake.jl b/test/reservoir_lake.jl
index 9570ccb63..bc2446f92 100644
--- a/test/reservoir_lake.jl
+++ b/test/reservoir_lake.jl
@@ -1,6 +1,6 @@
 @testset "reservoir simple" begin
     res = Wflow.SimpleReservoir{Float64}(
-        Δt = 86400.0,
+        dt = 86400.0,
         demand = [52.523],
         maxrelease = [420.184],
         maxvolume = [25_000_000.0],
@@ -31,7 +31,7 @@ end
 
 @testset "lake" begin
     lake = Wflow.Lake{Float64}(
-        Δt = 86400.0,
+        dt = 86400.0,
         lowerlake_ind = [0],
         area = [180510409.0],
         maxstorage = Wflow.maximum_storage([1], [3], [180510409.0], [missing], [missing]),
@@ -72,7 +72,7 @@ sh = [
     @test typeof(values(sh[1])) == Tuple{Vector{Float},Vector{Float}}
 
     lake = Wflow.Lake{Float}(
-        Δt = 86400.0,
+        dt = 86400.0,
         lowerlake_ind = [2, 0],
         area = [472461536.0, 60851088.0],
         maxstorage = Wflow.maximum_storage(
@@ -124,7 +124,7 @@ end
 
 @testset "overflowing lake with sh and hq" begin
     lake = Wflow.Lake{Float}(
-        Δt = 86400.0,
+        dt = 86400.0,
         lowerlake_ind = [0],
         area = [200_000_000],
         maxstorage = Wflow.maximum_storage(
diff --git a/test/run_sbm.jl b/test/run_sbm.jl
index e9dcd6ec3..62fdfba99 100644
--- a/test/run_sbm.jl
+++ b/test/run_sbm.jl
@@ -75,8 +75,8 @@ end
 
     @test model.clock.iteration == 1
 
-    @test sbm.θₛ[50063] ≈ 0.48755401372909546f0
-    @test sbm.θᵣ[50063] ≈ 0.15943120419979095f0
+    @test sbm.theta_s[50063] ≈ 0.48755401372909546f0
+    @test sbm.theta_r[50063] ≈ 0.15943120419979095f0
     @test mean(sbm.runoff) ≈ 0.04177459898728149f0
     @test mean(sbm.soilevap) ≈ 0.02122698830889417f0
     @test mean(sbm.actevap) ≈ 0.33545623834952154f0
@@ -92,8 +92,8 @@ model = Wflow.run_timestep(model)
 
 @testset "second timestep" begin
     sbm = model.vertical
-    @test sbm.θₛ[50063] ≈ 0.48755401372909546f0
-    @test sbm.θᵣ[50063] ≈ 0.15943120419979095f0
+    @test sbm.theta_s[50063] ≈ 0.48755401372909546f0
+    @test sbm.theta_r[50063] ≈ 0.15943120419979095f0
     @test mean(sbm.net_runoff) ≈ 0.23710522236947756f0
     @test mean(sbm.runoff) ≈ 0.23747368431143517f0
     @test mean(sbm.soilevap) ≈ 0.018750808322054897f0
@@ -304,7 +304,7 @@ model = Wflow.initialize_sbm_model(config)
 
 fp = model.lateral.river.floodplain.profile
 river = model.lateral.river
-Δh = diff(fp.depth)
+dh = diff(fp.depth)
 Δv = diff(fp.volume[:, 3])
 Δa = diff(fp.a[:, 3])
 
@@ -335,7 +335,7 @@ river = model.lateral.river
         297.8700179533214f0,
         463.35655296229805f0,
     ]
-    @test Δh .* fp.width[2:end, 3] * river.dl[3] ≈ Δv
+    @test dh .* fp.width[2:end, 3] * river.dl[3] ≈ Δv
     @test fp.a[:, 3] * river.dl[3] ≈ fp.volume[:, 3]
     # flood depth from flood volume (8000.0)
     flood_vol = 8000.0f0
@@ -449,7 +449,7 @@ Wflow.close_files(model, delete_output = false)
         @unpack vertical = model
         z = vertical.zi[i]
         kv_z = Wflow.hydraulic_conductivity_at_depth(vertical, z, i, 2, "exponential")
-        @test kv_z ≈ vertical.kvfrac[i][2] * vertical.kv₀[i] * exp(-vertical.f[i] * z)
+        @test kv_z ≈ vertical.kvfrac[i][2] * vertical.kv_0[i] * exp(-vertical.f[i] * z)
         @test vertical.z_exp == vertical.soilthickness
         @test_throws ErrorException Wflow.kh_layered_profile(
             vertical,
@@ -469,7 +469,7 @@ Wflow.close_files(model, delete_output = false)
         z = vertical.zi[i]
         kv_z =
             Wflow.hydraulic_conductivity_at_depth(vertical, z, i, 2, "exponential_constant")
-        @test kv_z ≈ vertical.kvfrac[i][2] * vertical.kv₀[i] * exp(-vertical.f[i] * z)
+        @test kv_z ≈ vertical.kvfrac[i][2] * vertical.kv_0[i] * exp(-vertical.f[i] * z)
         kv_400 = Wflow.hydraulic_conductivity_at_depth(
             vertical,
             400.0,
diff --git a/test/run_sbm_gwf.jl b/test/run_sbm_gwf.jl
index 9d2ba53ac..758ed4f0b 100644
--- a/test/run_sbm_gwf.jl
+++ b/test/run_sbm_gwf.jl
@@ -21,7 +21,7 @@ end
     sbm = model.vertical
 
     @test model.clock.iteration == 1
-    @test sbm.θₛ[1] ≈ 0.44999998807907104f0
+    @test sbm.theta_s[1] ≈ 0.44999998807907104f0
     @test sbm.runoff[1] == 0.0
     @test sbm.soilevap[1] == 0.0
     @test sbm.transpiration[1] ≈ 0.6117526566330049f0
@@ -32,7 +32,7 @@ model = Wflow.run_timestep(model)
 
 @testset "second timestep" begin
     sbm = model.vertical
-    @test sbm.θₛ[1] ≈ 0.44999998807907104f0
+    @test sbm.theta_s[1] ≈ 0.44999998807907104f0
     @test sbm.runoff[1] == 0.0
     @test sbm.soilevap[1] == 0.0
     @test sbm.transpiration[1] ≈ 1.0122634204681036f0
diff --git a/test/testing_utils.jl b/test/testing_utils.jl
index f343c0e37..06a660767 100644
--- a/test/testing_utils.jl
+++ b/test/testing_utils.jl
@@ -11,7 +11,7 @@
 #     https://github.com/stevengj/18S096-iap17/blob/master/pset3/pset3-solutions.ipynb
 
 # n coefficients of the Taylor series of E₁(z) + log(z), in type T:
-function E₁_taylor_coefficients(::Type{T}, n::Integer) where {T<:Number}
+function E1_taylor_coefficients(::Type{T}, n::Integer) where {T<:Number}
     n < 0 && throw(ArgumentError("$n ≥ 0 is required"))
     n == 0 && return T[]
     n == 1 && return T[-eulergamma]
@@ -26,8 +26,8 @@ function E₁_taylor_coefficients(::Type{T}, n::Integer) where {T<:Number}
 end
 
 # inline the Taylor expansion for a given order n, in double precision
-macro E₁_taylor64(z, n::Integer)
-    c = E₁_taylor_coefficients(Float64, n)
+macro E1_taylor64(z, n::Integer)
+    c = E1_taylor_coefficients(Float64, n)
     taylor = :(@evalpoly zz)
     append!(taylor.args, c)
     quote
@@ -38,9 +38,9 @@ macro E₁_taylor64(z, n::Integer)
 end
 
 # for numeric-literal coefficients: simplify to a ratio of two polynomials:
-# return (p,q): the polynomials p(x) / q(x) corresponding to E₁_cf(x, a...),
+# return (p,q): the polynomials p(x) / q(x) corresponding to E1_cf(x, a...),
 # but without the exp(-x) term
-function E₁_cfpoly(n::Integer, ::Type{T} = BigInt) where {T<:Real}
+function E1_cfpoly(n::Integer, ::Type{T} = BigInt) where {T<:Real}
     q = Polynomials.Polynomial(T[1])
     p = x = Polynomials.Polynomial(T[0, 1])
     for i = n:-1:1
@@ -51,8 +51,8 @@ function E₁_cfpoly(n::Integer, ::Type{T} = BigInt) where {T<:Real}
     return p, x * p + q
 end
 
-macro E₁_cf64(z, n::Integer)
-    p, q = E₁_cfpoly(n, BigInt)
+macro E1_cf64(z, n::Integer)
+    p, q = E1_cfpoly(n, BigInt)
     num_expr = :(@evalpoly zz)
     append!(num_expr.args, Float64.(Polynomials.coeffs(p)))
     den_expr = :(@evalpoly zz)
@@ -66,26 +66,26 @@ end
 
 # exponential integral function E₁(z)
 function expint(z::Union{Float64,Complex{Float64}})
-    x² = real(z)^2
-    y² = imag(z)^2
-    if real(z) > 0 && x² + 0.233 * y² ≥ 7.84 # use cf expansion, ≤ 30 terms
-        if (x² ≥ 546121) & (real(z) > 0) # underflow
+    xSq = real(z)^2
+    ySq = imag(z)^2
+    if real(z) > 0 && xSq + 0.233 * ySq ≥ 7.84 # use cf expansion, ≤ 30 terms
+        if (xSq ≥ 546121) & (real(z) > 0) # underflow
             return zero(z)
-        elseif x² + 0.401 * y² ≥ 58.0 # ≤ 15 terms
-            if x² + 0.649 * y² ≥ 540.0 # ≤ 8 terms
-                x² + y² ≥ 4e4 && return @E₁_cf64 z 4
-                return @E₁_cf64 z 8
+        elseif xSq + 0.401 * ySq ≥ 58.0 # ≤ 15 terms
+            if xSq + 0.649 * ySq ≥ 540.0 # ≤ 8 terms
+                xSq + ySq ≥ 4e4 && return @E1_cf64 z 4
+                return @E1_cf64 z 8
             end
-            return @E₁_cf64 z 15
+            return @E1_cf64 z 15
         end
-        return @E₁_cf64 z 30
+        return @E1_cf64 z 30
     else # use Taylor expansion, ≤ 37 terms
-        r² = x² + y²
-        return r² ≤ 0.36 ?
+        rSq = xSq + ySq
+        return rSq ≤ 0.36 ?
                (
-            r² ≤ 2.8e-3 ? (r² ≤ 2e-7 ? @E₁_taylor64(z, 4) : @E₁_taylor64(z, 8)) :
-            @E₁_taylor64(z, 15)
-        ) : @E₁_taylor64(z, 37)
+            rSq ≤ 2.8e-3 ? (rSq ≤ 2e-7 ? @E1_taylor64(z, 4) : @E1_taylor64(z, 8)) :
+            @E1_taylor64(z, 15)
+        ) : @E1_taylor64(z, 37)
     end
 end
 expint(z::Union{T,Complex{T},Rational{T},Complex{Rational{T}}}) where {T<:Integer} =
@@ -100,21 +100,21 @@ function expint(n::Integer, z)
     elseif n < 1
         # backwards recurrence from E₀ = e⁻ᶻ/z
         zinv = inv(z)
-        e⁻ᶻ = exp(-z)
-        Eᵢ = zinv * e⁻ᶻ
+        exp_minus_z = exp(-z)
+        Ei = zinv * exp_minus_z
         for i = 1:-n
-            Eᵢ = zinv * (e⁻ᶻ + i * Eᵢ)
+            Ei = zinv * (exp_minus_z + i * Ei)
         end
-        return Eᵢ
+        return Ei
     else
         # forwards recurrence from E₁
-        e⁻ᶻ = exp(-z)
-        Eᵢ = expint(z)
-        Eᵢ *= !isinf(Eᵢ)
+        exp_minus_z = exp(-z)
+        Ei = expint(z)
+        Ei *= !isinf(Ei)
         for i = 2:n
-            Eᵢ = (e⁻ᶻ - z * Eᵢ) / (i - 1)
+            Ei = (exp_minus_z - z * Ei) / (i - 1)
         end
-        return Eᵢ
+        return Ei
     end
 end