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 flow 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 wflow\_sediment, the new +SWAT's approach developed by Narasimhan et al. (2017). For wflow\_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, - , 2021. \ No newline at end of file + , 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/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/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 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