Changelog
All notable changes to this project will be documented in this file.
The format is based on Keep a Changelog, and this project adheres to Semantic Versioning.
Unreleased
Added
- Support for reinfiltration of surface overland flow water. Can be set through the
surface_water_infiltration
key in the[model]
section. This way, the water stored on the land surface is added to the water that is allowed to infiltrate into the soil. - Support for adding a water level flow-threshold for both kinematic wave and local inertial overland flow. This setting can be changed by changing the
lateral.land.h_thresh
key in the settings file. This can be added as a map or as a fixed value (through the use ofh_thresh.value
).
v0.8.1 - 2024-08-27
Fixed
- Reduce allocations in update of vertical
SBM
concept.
v0.8.0 - 2024-08-19
Fixed
- Added missing BMI function
get_grid_size
, it is used for unstructured grids, for example to get the length of arrays returned by BMI functionsget_grid_x
andget_grid_y
. - Added a check for the solution of the quadratic equation as part of the Modified Puls approach for lake outflow. Lower limit should be zero (very small negative values can occur).
- Limit lake evaporation (added variable
actevap
) and lake outflow to prevent negative lake storage. The variableactevap
has also been added to the reservoir module. - The
set_states
function for model typesbm
with local inertial routing for river and land component. - Inflow to reservoir and lake locations for local inertial routing: 1) with floodplain routing, the floodplain discharge was not added to the inflow of these locations, 2) the
to_river
variable from overland flow and lateral subsurface flow was not added to the inflow of these locations. - Close netCDF
NCDataset
with state variables in extended BMI functionsave_state
. - For the computation of Gash interception model parameter
e_r
multiply the precipitation input with the canopy fraction (this was only done for the potential evapotranspiration input). - Function
BMI.update_until
could throw anInexactError: Int64
caused by a not whole number, representing how many times (steps
) an update of a wflow model is required. This is fixed by usingdivrem
for the computation of the number ofsteps
in this function. An error is thrown when the absolute remainder ofdivrem
is larger thaneps()
, or when the number ofsteps
is negative. - Fixed internal and external broken links in docs.
- The internal time step of the local inertial model (
stable_timestep
function) can get zero whenLoopVectorization
is applied (@tturbo
) to the for loop of these functions. This issue occured on a virtual machine, Windows 10 Enterprise, with Intel(R) Xeon(R) Gold 6144 CPU (2 processors). This has been fixed by replacing@tturbo
withreduction
ofPolyester.jl
. - Fixed required states of the model type
sbm_gwf
: addedh_av
for the river and land domain. - For the
constanthead
boundary ofGroundwaterFlow
thehead
should not be changed (was set totop
elevation of theaquifer
ifhead
>top
), andexfiltwater
should be 0 for these boundary cells.
Changed
- Stop exposing scalar variables through BMI. The
BMI.get_value_ptr
function was not working correctly for scalar model variables (aview
was applied). Only a few scalar model parameters are defined, and it is not expected that exposing these variables is required (e.g. for model coupling) while code changes for these variables (including struct fields) are required. - The local inertial routing (constant) river boundary condition
depth
at a ghost node (downstream river outlet) was set at 0.0 in the code and can now be provided through the model parameter netCDF file (riverdepth_bc
). In addition, the boundary condition river lengthdl
at a ghost node should be set through the model parameter netCDF file (riverlength_bc
), providing this boundary condition through the[model]
settings in the TOML file (model.riverlength_bc
) has been deprecated. - As part of the vertical
SBM
concept: 1) add variablenet_runoff
(land runoff minus actual open water evaporationae_openw_l
) and 2) change the definition of variablerunoff
by removing the subtraction ofae_openw_l
(total land runoff). The variablenet_runoff
is now input to the kinematic wave for overland flow (replaces the originalrunoff
variable as input). The lower limit of the originalrunoff
variable was set at zero, which may result in water balance errors. This lower limit was implemented in the Python version of wflow_sbm, to avoid a very slow kinematic wave solution for surface flow computation (no distinction between a kinematic wave for overland and river flow). When in the Python version of wflow_sbm the kinematic wave for surface flow was split into a river and land component, the lower limit was removed for river runoff (net_runoff_river
), but not for land runoff. - 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 coefficientkc
(use ofet_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.
- Docs: 1) improved description of different model configurations in model-setup.md, also in relation to hydromt_wflow in docs, 2) citing info related to wflow_sbm publication in Geosci. Model Dev. (from in review to published).
- Root water uptake reduction (Feddes): 1) extend critical pressure head parameters with
h3_low
andh3_high
, 2) all critical pressure head parameters can be set (values forh1
,h2
andh3
were fixed) and 3) the root water uptake reduction coefficient $\alpha$ can be set at 0 (default is 1) at critical pressure headh1
(through the model parameteralpha_h1
). - For the actual transpiration computation in
SBM
, the potential transpiration is partitioned over the soil layers with roots according to the model parameterrootfraction
(fraction of the total root length in a soil layer). Previously, for each soil layer (from top to bottom) the actual transpiration was computed, and the remaining potential transpiration was used in the next soil layer.
Added
- Total water storage as an export variable for
SBM
concept. This is the total water stored per grid cell in millimeters. Excluded from this variable are the floodplain, lakes and reservoirs. - Checks to see if all states are covered in the .toml file. If not all states are covered, an error is thrown. If there are more states specified than required, these states are ignored (with a warning in the logging), and the simulation will continue.
- Support different vertical hydraulic conductivity profiles for the
SBM
concept. See also the following sections: The SBM soil water accounting scheme and Subsurface flow routing for a short description. - Local inertial routing to
sbm_gwf
model type. - Water demand and allocation computations for model types
sbm
andsbm_gwf
.
v0.7.3 - 2024-01-12
Fixed
- Documentation: add leakage term to the wflow_sbm figure, document external input parameter
ksathorfrac
and fix description of adding externalinflow
to the kinematic wave. - Fixed BMI functions (e.g.
BMI.get_value
) that deviated from BMI specifications (BasicModelInterface.jl), including function arguments, return types and the BMI specification that arrays are always flattened (this was not the case for variables stored as 2-dimensional arrays or as vector of SVectors). - Bump compat for NCDatasets to 0.13, 0.14.
- The solution for lake outflow as part of the Modified Puls Approach. The inflow and outflow variables are defined for period
Δt
, and not att1
andt2
(instantaneous) as in the original mass balance equation of the Modified Puls Approach. Because of this, the terms of the quadratic equation (and solution) were fixed. - Use
kvfrac
for the computation of vertical saturated hydraulic conductivity at the bottom of the soil layer, sincekvfrac
is also used for the computation of vertical unsaturated flow.
Changed
- For cyclic parameters different cyclic time inputs are supported (only one common cyclic time (for example daily or monthly) was allowed).
- BMI: 1) added grid information (type and location) and whether a variable can be exchanged to metadata Structs, 2) extend model grid functions for wflow components that store variables on
edges
(local inertial model) withget_grid_edge_count
andget_grid_edge_nodes
.
Added
- Functions for loading and saving states and getting the
starttime
in Unix time format. This is convenient for coupling with OpenDA (and other external) software. The set states functionality from the initialization function has been moved to a separateset_states
function for eachModel
type, to support loading states independent of initialization.
v0.7.2 - 2023-09-27
Fixed
- Water balance of modified Rutter interception model. The sum of the stemflow partitioning coefficient
pt
and free throughfall coefficientp
could get larger than 1, resulting in an overestimation of stemflow and throughfall amounts and negative net interception amounts. And the first drainage amountdd
, controlled by a change over time in canopy storage capacitycmax
, should not be subtracted from precipitation to compute net interception. - The
netinterception
(net interception) computed by the modified Rutter interception model was stored asinterception
inSBM
, while this should be theinterception
(interception loss by evaporation) output of the modified Rutter interception model. Theinterception
ofSBM
is used to compute the total actual evapotranspirationactevap
.
v0.7.1 - 2023-06-30
Fixed
- State initialization of 1D floodplain
volume
. In the initialization function the wrong field name of typeFloodPlainProfile
was used (area
instead ofa
).
v0.7.0 - 2023-06-12
Fixed
BMI.get_time_units
now gets called on the model rather than the type, like all other BMI functions, exceptBMI.initialize
. Also it returns "s" instead of "seconds since 1970-01-01T00:00:00", in line with the BMI specification.- Added the
interception
component to total actual evapotranspirationactevap
ofSBM
(was defined as the sum of soil evaporation, transpiration and open water evaporation).
Changed
- The time values returned in the BMI interface are no longer in seconds since 1970, but in seconds since the model start time. This is more in line with standard BMI practices.
- The
starttime
was defined one model timestepΔt
ahead of the actual model time (the initial conditions timestamp (state time)). As a consequence this was also the case for the current model time. To allow for an easier interpretation of wflow time handling, either through BMI or directly, thestarttime
is now equal to the state time, resulting in current model times without an offset. - Using more than 8 threads can result in too much overhead with
Threads.@threads
. After performance testing, this has been changed for kinematic wave routing and the verticalSBM
concept to spawning tasks withThreads@spawn
for number of threads <= 8, where each task iterates over a chunk of sizebasesize
. For more than 8 threads the low overhead threadingPolyester.@batch
(including theminbatch
argument) is used. For local inertial routing the use ofThreads.@threads
has been changed to threaded loop vectorization (river and 1D floodplain local inertial momentum equation) andPolyester.@batch
.
Added
- For (regulated) lakes with rating curve of type 1 (H-Q table), lake
storage
above themaximumstorage
(based on maximum water level from the H-Q table) is spilled instantaneously (overflow) from the lake. - Added support to use
sum
as a reducer function for csv and scalar output options.
v0.6.3 - 2023-03-01
Fixed
- Removed error when
_FillValue
is present in the time dimension of the forcing netCDF file. The simulation is allowed to continue with the attribute present, given that there are no missing values in the time dimension. This is checked by the code, and an error is thrown if this is the case. - Column index of daily lake rating curves. This was incorrectly based on
dayofyear
with a maximum of 365. The column index should be based on julian day (leap days are not counted).
Changed
NCDatasets
version. Reading thetime
dimension of multifile netCDF file became very slow sinceNCDatasets
v0.12.4, this issue has been solved in v0.12.11.- Store the
time
dimension of the forcing netCDF file as part of the structNCreader
instead of callingdataset["time"][:]
each time step when loading forcing data.
Added
- Show total duration of simulation in the log file (info), and show the current time at execution of each timestep (debug).
- Support for exponential decline in horizontal conductivity in the sbm_gwf concept. This can be enabled using the
conductivity_profile
setting (either "uniform" or "exponential"). If set to "exponential", it exponentially reduces thekh0
(orconductivity
) based on the value ofgwf_f
to the actual horizontal conductivity (k
). - An optional 1D floodplain schematization for the river flow inertial model, based on provided flood volumes as a function of flood depth per river cell. See also the following sections: SBM + Local inertial river and River and floodplain routing for a short description, and the following section for associated model parameters.
v0.6.2 - 2022-09-01
Fixed
- Two issues related to reservoir and lake locations as part of local inertial model: 1) added as boundary points to the update of overland flow, 2) fixed check reservoir and lake location in update river flow.
- Limit flow (
qx
,qy
andq
) in local inertial model when water is not available (set to zero). - In the grid output netCDFs, don't set the
_FillValue
attribute, since the CF conventions don't allow it. - Glacier parameter
g_sifrac
needs to be converted during initialization (time dependent). - The check that the sum of adaptive timesteps (
Δt
) of the local inertial model (1D and 2D) does not exceed the model timestep.
Changed
- Changed depth
h
for reservoir and lake locations as part of the river local inertial model frombankfull_depth
to zero.
Added
- External inflow to the SBM + Local inertial river (1D) and land (2D) model configuration.
v0.6.1 - 2022-04-26
Fixed
- Fixed an error with the log file, when writing to a folder that does not (yet) exists. Now, the folder is created prior to writing the log file.
- Fixed a MethodError for
read_dims
, thrown when reading netCDF data with NCDatasets.jl 0.12.3 or higher.
v0.6.0 - 2022-04-14
Added
- The FLEXTopo model.
- Set a (different) uniform value for each index of input parameters with an extra dimension. A list of values (instead of one uniform value) that should be equal to the length of the extra dimension can be provided in the TOML file. See also Modify parameters.
- Modify input parameters with an extra dimension at multiple indices with different
scale
andoffset
values, through the TOML file. See also Modify parameters. - Added support for netCDF compression for gridded model output, through the option
compressionlevel
in the[output]
section in the TOML file. The setting defaults to 0 (no compression), but all levels between 0 and 9 can be used.
Changed
- Re-organized the documentation: moved explanation of different model concepts to a model documentation section, added a user guide to explain setting up the model, added new figures to the description of wflow_sbm.
- The unit of lateral subsurface flow
ssf
ofLateralSSF
is now $m^3 d^{-1}$. The unit was $m^3 t^{-1}$, where $t$ is the model timestep. Other flow variables are already stored independently from $t$, this allows for easier interpretation and to use states independently of $t$. - Changed the reference level of water depth
h
andh_av
of 2D overland flow (ShallowWaterLand
) for cells containing a river from river bed elevationzb
to cell elevationz
.
Fixed
- Fixed calculation of average water depth
h_av
of 2D overland flow (ShallowWaterLand
) with the local inertial approach. The summation ofh
was not correct, resulting in too low values forh_av
. For river cells of 2D overland flowh_av
was only updated as part of the river domain (ShallowWaterRiver
), this value is now also updated as part of the land domain (ShallowWaterLand
). - Fixed the following two flow width issues for 2D overland flow of the local inertial model: 1) The flow widths
xwidth
andywidth
were mapped incorrectly (reversed) to the flow calculation inx
andy
direction, respectively. 2) For river cells the effective flow width for overland flow was not determined correctly: theriverwidth
vector supplied to the functionset_effective_flowwidth!
covered the complete model domain and should have covered the river domain, resulting in incorrect river widths and thus incorrect effective flow widths for river cells.
v0.5.2 - 2022-02-03
Changed
- Model types
sbm_gwf
andhbv
use the same approach for the calculation of the drain widthdw
as model typesbm
. - Renamed
h_bankfull
parameter tobankfull_depth
for consistency between kinematic-wave and local-inertial river routing. When using the old name under the[input.lateral.river]
section of the TOML file, it will work but it is suggested to update the name.
Added
- Additional log messages and log file as output, see also Logging.
- Option to use the local inertial model for river flow as part of the
sbm
model type. See also SBM + Local inertial river. - Option to use the local inertial model for 1D river flow combined with 2D overland flow as part of the
sbm
model type. See also SBM + Local inertial river (1D) and land (2D).
Fixed
- Model type
hbv
: the surface width for overland flow was not corrected with the river width. - Fixed use of absolute path for
path_forcing
in TOML file, which gave an error in wflow v0.5.1. - Fixed a crash when glacier processes are simulated as part of the
hbv
concept (Δt was not defined). - When the surface flow width for overland flow is zero, the water level
h
of the kinematic wave should not be calculated, otherwise this results inNaN
values. When the model is initialized from state files,q
andh
are set to zero for indices with a zero surface flow width. - Fixed how number of iterations
its
for kinematic wave river flow are calculated during initialization when using a fixed sub-timestep (specified in the TOML file). For a model timestep smaller than the fixed sub-timestep an InexactError was thrown. - Fixed providing a cyclic parameter when the netCDF variable is read during model initialization with
ncread
, this gave an error about the size of the netCDFtime
dimension. - Fixed CSV and netCDF scalar output of variables with dimension
layer
(SVector
).
v0.5.1 - 2021-11-24
Fixed
- Fixed calculation of
exfiltwater
as part of thesbm_gwf
model type. This was based directly on groundwater head above the surface level, without multiplying by thespecific_yield
, resulting in an overestimation ofexfiltwater
. This is required since the groundwater model estimates the head by dividing the volume by the specific yield or storativity of the aquifer. So, should the groundwater table rise above surface level, the head above surface level does not represent a water column one to one. (This also means the groundwater model (slightly) overestimates heads when the head rises above the surface level. However, this water is immediately removed, and the head will be set to surface level.)
Added
- Optional
dir_input
anddir_output
keys in the TOML, which can be used to quickly change the path for all input or output files that are given as a relative path.
v0.5.0 - 2021-11-12
Changed
- Scaling of potential capillary rise is replaced by a common approach found in literature, based on the water table depth
zi
, a maximum water depthcap_hmax
beyond which capillary rise ceases, and a coefficientcap_n
. See also Capillary rise. Multiplying the scaling factor with the ratio of model time step andbasetimestep
in the original approach resulted in (much) smaller capillary fluxes at sub-daily model time steps compared to daily model time steps, and is not used in the new approach. Parameterscap_hmax
andcap_n
can be set through the TOML file, parametercapscale
of the previous approach is not used anymore.
Fixed
- Conversion of
GroundwaterFlow
boundaries [$m^3 d^{-1}$] as part of model conceptsbm_gwf
to $m^3 s^{-1}$ for sub-daily model time steps. For the conversion thebasetimestep
(86400 s) should be used (and not the model time step).
v0.4.1 - 2021-11-04
Changed
- The $\alpha$ parameter of the kinematic wave has a fixed value now and is not updated because of changes in water height (this could result in large water balance errors). See also Surface routing.
- Cyclic input for other structs than vertical are also now supported (for example cyclic inflow to the river).
- Moved
update_forcing!
andupdate_cyclic!
functions to therun
function. It is now easier to implement a customrun
function with custom loading of input data (forcing and cyclic parameters).
Added
- Check if reservoirs and lakes have downstream nodes. Without downstream nodes is not supported and in that case an error message is thrown that is easier to understand than the previous one: "ArgumentError: Collection is empty, must contain exactly 1 element."
- For the
input.path_forcing
TOML setting we use Glob.jl to expand strings like"data/forcing-year-*.nc"
to a set of netCDF files that are split in time. - External water inflow (supply/abstractions) added to the kinematic wave
inflow
in m3/s. It is added/removed tosf.qlat[v]
before computing the newq[v]
with the kinematic wave equation. - Fixed values for forcing parameters are supported, see also Fixed forcing values.
Added
- Option to use the local inertial model for river flow as part of the SBM + Kinematic wave. See also SBM + Local inertial river.
Fixed
- River inflow for reservoirs and lakes in the kinematic wave. This inflow was based on
sf.q[v]
at the previous time step, and this has been fixed to the current time step.
v0.4.0 - 2021-09-02
Changed
- Changed length units for lateral subsurface flow component from millimeter to meter. This means that state netCDF files from previous versions can only be reused if
ssf
is divided by 10^9. - Add snow and glacier processes to wflow_sbm figure of the documentation.
Added
- Multi-threading of vertical SBM concept and lateral kinematic wave components (overland, river and subsurface flow) of wflow_sbm model SBM + Kinematic wave.
- Improved error message for CSV Reducer.
- The TOML keys
kv₀
,θᵣ
andθₛ
have been replaced with the ASCII versionskv_0
,theta_r
andtheta_s
, to avoid encoding issues with certain text editors. The old keys still work as well.
Fixed
- Calculation of volumetric water content of vertical SBM (soil layers and root zone).
- Update of
satwaterdepth
in SBM (evaporation was only subtracted from a local variable, and not fromsbm.satwaterdepth
). - Fixed the index
lowerlake_ind
ofNaturalLake
. Linked lakes were not simulated correctly (flows between upstream and downstream lake). - Interpolation function
interpolate_linear(x, xp, fp)
for CSV tables lakes. Whenx
was larger or smaller thanxp
, maximum(xp
) or minimum(xp
) was returned instead of maximum(fp
) or minimum(fp
). - Fixed model timestep of reservoirs (
SimpleReservoir
) and lakes (NaturalLake
) in relation to precipitation and evapotranspiration fluxes. This was set to the fixed wflowbasetimestep
of 86400 s, and should be set to the actual model time step from the TOML configuration file. - Add
flux
fromDrainage
(GroundwaterFlow
) in thesbm_gwf_model
to the overland flow component instead of the river component of the kinematic wave. - Fixed option
constanthead = false
(TOML file),constant_head
field ofGroundwaterFlow
was not defined in this case. Fixed this by initializing empty fields (Vector) for structConstantHead
. - Fixed return max(0, boundary.flux[index]) to return max(0, flux) the flux should be returned when cell is dry, no negative value allowed.
- Fixed path to initialize lake to: dirname(static_path)
- Fixed outflow = 0, when lake level is below lake threshold. Before a negative value could enter the log function and model would fail.
- Fixed the lake storage initialization. For continuation runs (
reinit = false
), this caused the lake to be reset to the initial conditions instead of the saved state.
v0.3.1 - 2021-05-19
Fixed
- Ignore extra dimensions in input netCDFs if they are size 1
v0.3.0 - 2021-05-10
Changed
- New deposition process for coarse sediment in the reservoirs with a new parameter
restrapefficiency
in the sediment model. - New variables added to the
LandSediment
andRiverSediment
structs in order to save more output from the sediment model. - Added variables
volume
andinwater
toSurfaceFlow
struct, this is convenient for the coupling with the water quality model Delwaq. - River water level (
h
) and discharge (q
) forced directly into theRiverSediment
struct (instead of using theOverlandFlowSediment
struct first). - Require Julia 1.6 or later.
Added
- Modify model parameters and forcing through the TOML file (see Modify parameters).
- Run wflow_sbm (SBM + kinematic wave) in two parts (until recharge and after subsurface flow) from BMI, including the option to switch off the lateral subsurface component of wflow_sbm.
- Support more netCDF dimension and axis order variants.
Fixed
- Corrected a bug in sediment deposition in the river (case when incoming sediment load is more than the river transport capacity).
- Fixed update of
snow
andglacierstore
fields of HBV and SBM concepts by theglacier_hbv
module.
v0.2.0 - 2021-03-26
Changed
- Removed dependency of the
f
model parameter of wflow_sbm on the parameters $\theta_{s}$, $\theta_{r}$ and $M$. This approach is used in Topog_SBM, but not applicable for wflow_sbm. Thef
parameter needs to be provided as part of the netCDF model parameter file. - Grid properties as cell length and elevation now stored as part of the
model.land.network
component and not as part of the vertical model components, as it is not used by these components.altitude
(elevation) should now be provided as part of the[input]
section of the TOML configuration file, and not as part of the[input.vertical]
section. - Removed parameter $\theta_{e}$ from SBM struct (not used in update). Parameters $\theta_{s}$ and $\theta_{r}$ included separately (instead of $\theta_{e}$) in
LateralSSF struct
, now directly linked to SBM parameters. - Improve error messages (netCDF and cyclic flow graph).
Added
- Export of netCDF scalar timeseries (separate netCDF file from gridded timeseries). This also allows for importing these timeseries by Delft-FEWS (General Adapter).
Fixed
- Model parameter Manning's
n
now used during the update of thestruct SurfaceFlow
, to change the related $\alpha$ parameter of the kinematic wave for channel flow.