8. Known issues¶
Known bugs, risks and open scientific questions in ASYNCH, found by building and running the code and by reading it (September 2026). Any fix should be verified with the regression harness (09_reproducibility.md).
How each item was established:
Confirmed: reproduced by running the code. The evidence is quoted: a crash, an AddressSanitizer report, a comparison.
Code reading: visible in the source, but not (yet) triggered by an example.
Open question: a scientific or design question, not necessarily a bug. These need a hydrologist’s judgement, not just a programmer’s.
Severity scale: Critical memory corruption, wrong results or a crash in normal use; High a crash or wrong result in a plausible configuration; Medium; Low cosmetic, or only in unusual situations. This is the technical severity. Chapter 7 ranks the same issues by their effect on a user of the model, which can differ (a crash is serious technically, but less dangerous than a silent error).
Line numbers refer to commit 84da43a (the state of master at the time of writing).
Summary table¶
ID |
Severity |
Status |
Evidence |
One-line description |
|---|---|---|---|---|
Critical |
fixed |
confirmed |
Model 254 uses model 256’s snapshot filter: heap buffer overflow, crashes clearcreek on 1 process |
|
High |
fixed |
code reading |
Snapshot values are filtered only for links owned by MPI rank 0: output depends on process count |
|
High |
fixed |
confirmed |
Output file closed twice at shutdown: every debug build aborts at the end of a run |
|
High |
fixed |
confirmed |
Solver index 3 or 4 (advertised as “implicit”) segfaults; the index is never validated |
|
Medium |
fixed |
code reading |
|
|
Medium |
fixed |
confirmed |
|
|
Medium |
fixed |
code reading |
|
|
Medium |
fixed |
confirmed (UB sanitizer) |
Misaligned |
|
Low |
fixed |
code reading |
Model 402 dam check prints a debug line on every call |
|
Low |
fixed (except riversys.c indentation) |
compiler |
Missing prototype for |
|
Low |
open |
code reading |
~75 |
|
Medium |
fixed |
code reading |
|
|
Critical |
fixed |
confirmed (ASan) |
Solver methods 0 and 1 used Butcher coefficients from freed stack memory: random results or endless runs |
|
High |
fixed |
confirmed |
Reading an |
|
High |
fixed |
confirmed |
A missing output folder loses all results, yet the run ends with a success exit code |
|
High |
fixed |
confirmed (ASan) |
Links with more than 8 parents overflowed memory; reader and solver disagreed on the limit |
|
Medium |
fixed |
confirmed |
Library functions that crashed with built-in models or were missing (global parameters, duration, init file) |
|
Medium |
fixed |
code reading |
|
|
Low |
fixed |
code reading |
Model 190 read a third forcing value that does not exist (unused, no effect on results) |
|
High |
fixed |
confirmed |
Custom outputs of non-interpolated states were written as 0 / memory contents |
|
Low |
fixed |
code reading |
Custom models inherited the snapshot filter of the built-in model with the same number |
|
Low |
fixed |
unit test |
Wrong constant in the Dormand-Prince dense-output derivative (no effect: multiplied by 0) |
|
High |
fixed |
unit test |
Models 263 and 601-603 wrote/read one or two parameters past the per-link array |
|
Medium |
fixed |
unit test |
7 model numbers without equations crashed; RK tables shared by all solvers of a program |
|
High |
fixed |
test |
Binary rain files: file past the range read, last file lasted 0.0001 min, overflow, crash on a missing file |
|
High |
fixed |
valgrind, unit test |
Models 105 and 263 left derivatives unset: states changed by values left in memory |
|
Low |
fixed |
valgrind |
Consistency check read the parents’ non-dense states uninitialised (no effect on results) |
|
High |
fixed |
unit test, examples |
Models 257, 258, 259, 261, 262 evaporated ponded water 1000 times too fast (since 2015) |
|
High |
fixed |
unit test |
Model 255 turned channel storage into 8 to 15 times too much discharge (since 2022) |
|
Medium |
fixed |
unit test |
Models 258, 259: the baseflow equation read the accumulated evaporation (state 6) |
|
Medium |
fixed |
code reading |
Model 249: baseflow mixed m³/min and m³/s; its reservoir version printed at every evaluation |
|
Low |
fixed |
code reading |
Model 257: the accumulated evaporation output was 720 times too large |
|
Low |
fixed |
code reading |
Models 0–6, 105, 200, 2000: the area column of |
|
Medium |
fixed |
unit test |
Models 225, 601–609: NaN when every hillslope storage is empty (division by zero) |
|
Low |
fixed |
unit test |
Model 606 left derivative 5 unset when the tile storage is empty |
|
Medium |
open |
code reading |
Model 249 with reservoirs: the reservoir function returns derivatives where the solver expects states |
|
High |
fixed |
confirmed |
No automated regression tests; only one unit test ( |
|
Medium |
resolved |
confirmed |
clearcreek references (2015) differ: another configuration, and a 2021 change to model 254 |
|
Medium |
open |
confirmed |
Model 259 benchmark cannot be reproduced from the files in the repository |
|
Medium |
fixed |
confirmed |
Examples 258/259 pointed to a file on the original developers’ cluster |
|
n/a |
information |
confirmed |
Results change at noise level with the number of MPI processes |
|
High |
resolved |
confirmed |
Old Python API broken beyond repair; replaced by the |
|
Medium |
fixed |
confirmed |
The Python package made the symbols of libasynch’s HDF5 global: h5py could not be imported after it |
|
Low |
resolved |
confirmed |
~6 500 lines (15 %) of C were never compiled; removed |
|
Medium |
open |
code reading |
A model is defined in 7 different places; duplicated unreachable code |
|
Low |
partly fixed |
confirmed |
CI (Travis) was dead: replaced by GitHub Actions; |
|
Low |
fixed |
confirmed |
CLI slept 1 s during initialisation |
|
Medium |
open |
code reading |
Snapshots gather every link through rank 0 one message at a time |
|
n/a |
open |
hypothesis |
Scheduler, barriers and step-size resets in |
|
Medium |
resolved |
confirmed |
Explicit solvers are limited by stiffness: steps far smaller than accuracy needs; stiff solver added (index 4) |
|
Low |
open |
confirmed |
Explicit solvers record peaks only at the end of a step: peak values slightly low |
|
n/a |
open question |
code reading |
Parent states indexed with |
|
High |
resolved |
confirmed |
Model 254 baseflow floor |
|
n/a |
open question |
code reading |
Potential evaporation assumes a 30-day month |
|
n/a |
open question |
code reading |
Models 400–405: |
|
n/a |
open question |
code reading |
Snapshot filter rewrites cumulative states with |
|
n/a |
open question |
code reading |
Model 254 evaporation always runs at the full potential rate; clamping then creates water |
|
n/a |
open question |
code reading |
Models 105 and 263: states without equations; what was intended? |
|
n/a |
open question |
code reading |
Model 603: the lateral flow of the deepest soil layer never reaches the channel |
|
n/a |
open question |
code reading |
Model 605: the third intercept of the subsurface runoff curve counts the first one twice |
|
n/a |
open question |
code reading |
Models 1–5: hillslope loss and channel gain differ by a factor h_b^(2/3) |
|
n/a |
open question |
code reading |
Models 400–405: storages and rates mixed (1-minute assumption); 402/403 dams drop local runoff |
|
n/a |
open question |
code reading |
Models 604, 606, 30, 21: unexplained or inconsistent constants and units |
|
Low |
fixed |
code reading |
|
|
Medium |
fixed |
code reading |
Units in |
Bugs¶
B-01¶
Model 254 uses the snapshot filter of model 256: heap buffer overflow. Critical, confirmed.
Fixed (2026-09-25): the missing break; was added. clearcreek now runs on 1 process; with
2 processes all 27 output files are bit-identical to the original code (see CHANGELOG).
SetOutputConstraints in src/models/definitions.c:885-893 has no break; after
case 254:. Execution falls through into case 256: (a classic C pitfall: a
switch keeps executing the following cases until it meets a break). Model 254
therefore gets OutputConstraints_Model256_Hdf5, which reads and writes states[7].
Model 254 only has 7 states (states[0] … states[6]).
The filter runs on a packed buffer (src/processdata.c:1892-1910): each record is a
4-byte link id followed by dim doubles. states[7] is the first 8 bytes of the
next link’s record. For the last link it is 8 bytes past the end of the malloc
(a heap buffer overflow).
Evidence:
==7192==ERROR: AddressSanitizer: heap-buffer-overflow ... READ of size 8
#0 OutputConstraints_Model256_Hdf5 src/models/output_constraints.c:60
#1 DumpStateH5 src/processdata.c:1910
#2 Advance src/advance.c:108
0x... is located 0 bytes after 381540-byte region (allocated at processdata.c:1892)
With a normal (non-sanitizer) build, examples/clearcreek.gbl crashes on 1 MPI process:
Fatal glibc error: malloc.c:2599 (sysmalloc): assertion failed ...
It “works” with mpirun -n 4 (as in the README) only because the last link then
belongs to another process, so the overflowing write never happens (see B-02).
Adding the missing break; was verified (temporarily) to remove the crash.
Proposed fix: add break;. Additionally make every filter receive dim and
never index beyond it.
B-02¶
Snapshot values are filtered only on rank 0. High, code reading (consistent with B-01 crashing only on 1 process). Fixed (2026-09-25): the filter is applied to every link. In a 2-process clearcreek snapshot, the original code left 4 values in (0, 1e-12) unfiltered; the fixed code leaves none.
In DumpStateH5 (src/processdata.c:1903-1911) the output filter
(OutputConstrainsHdf5) is applied only to links that live on process 0. Values
received from other processes (MPI_Recv branch) are written unfiltered. The same
simulation therefore writes a different snapshot file depending on the number of
MPI processes. That is a reproducibility problem, because snapshots are used as
initial conditions for the next run.
Proposed fix: apply the filter after MPI_Recv too (or on the sending side).
B-03¶
outputfile is closed twice. High, confirmed.
Fixed (2026-09-25): the pointer is set to NULL after each fclose. Debug and sanitizer
builds now run every example to the end with exit code 0.
Asynch_Delete_Temporary_Files (src/asynch_interface.c:1077-1078) calls
fclose(asynch->outputfile) but does not set it to NULL. Asynch_Free
(src/asynch_interface.c:442-443) then closes it again. Closing a FILE* twice is
undefined behaviour. glibc detects it and aborts:
free(): double free detected in tcache 2
... exited on signal 6 (Aborted)
Asynch_Free is only called when NDEBUG is not defined (src/asynch_cli.c:351-353),
so this only hits debug builds. That is probably why the README tells you to compile with
-DNDEBUG. Every example currently exits with code 134 in a debug build, even though
the results were written correctly.
Proposed fix: asynch->outputfile = NULL; after the fclose, and call Asynch_Free
in all builds.
B-04¶
Solver index 3/4 crashes. High, confirmed.
Fixed (2026-09-25): the index is validated (in the .gbl file and in .rkd files). An invalid
value stops the run with Error: numerical solver index 4 in the global file is not valid. Use 0 (RK 3(2)), 1 (RK 4(3)) or 2 (Dormand-Prince 5(4)). The comments in the example .gbl files and in the docs were corrected.
The global file comment says %Numerical solver index (0-3 explicit, 4 implicit).
In reality (src/riversys.c:686-690) the table contains:
index |
method |
type |
|---|---|---|
0 |
RK3(2) dense (“RKDense3_2”) |
explicit |
1 |
RK4(3) dense (“TheRKDense4_3”) |
explicit |
2 |
Dormand–Prince 5(4) dense (“DOPRI5_dense”) |
explicit |
3 |
Radau IIA order 3 |
implicit, but it has no stepper (the old |
The index read in src/config_gbl.c:628-630 is never validated. Running examples/test.gbl with index 3,
4 or 7 gives a segmentation fault (exit 139). With index 3 you first see
Warning: No solver selected for link ID 1.
Proposed fix: validate the index (0–2) with a clear error message; fix the comments in all .gbl files and docs.
B-05¶
Destroy_ErrorData frees the wrong addresses. Medium, code reading + compiler warning.
Fixed (2026-09-25) together with B-14: it frees the arrays and the per-link structure.
src/system.c:204-207: free(&error->abstol) frees the address of the field
(inside a struct) instead of the memory the field points to. It should be
free(error->abstol). It is reached only when tolerances come from an .rkd file and
the solver is freed (debug builds).
B-06¶
Asynch_Get_Num_Links truncates. Medium, code reading.
src/asynch_interface.c:500 returns unsigned short (max 65 535). State-wide
networks (e.g. Iowa, ~400 000 links) are silently truncated. The CLI does not use it,
but any external program would.
Fixed (2026-09-25): returns unsigned int. On a synthetic network of 70 000 links it returns 70 000; the old
type would have returned 4 464.
B-07¶
DumpStateH5 edge cases. Medium, code reading.
Fixed (2026-09-25): the search is bounded and the buffers are freed.
src/processdata.c:1861-1863: while (assignments[i] != my_rank) i++; runs past the
end of the array if process 0 owns no link (possible with many processes and a small
network). The buffer data_storage allocated at line 1892 is never freed (a memory
leak at every recurrent snapshot).
B-08¶
Misaligned double access. Medium, confirmed by UndefinedBehaviorSanitizer.
Fixed (2026-09-25): states are filtered in an aligned array before being packed. The sanitizer
reports for the examples went from 21 to 0.
Snapshot records are packed as uint32 + doubles, so the doubles sit at addresses
that are not multiples of 8, and the filters access them through double*
(src/models/output_constraints.c). This is undefined behaviour in C. x86 tolerates it,
but other architectures (and optimisers) may not. Fix: filter a properly aligned
copy before packing it.
B-09¶
Debug print in model 402. Low. src/models/check_state.c:48-55 has int debug = 1;
and prints found dam_check_qvs_402 on every call. That floods the output and slows runs with dams.
Model 403 did the same. Fixed (2026-09-25): debugging is off by default.
B-10¶
Compiler-detected mistakes. Low.
src/forcings.c:176callsCreate_Rain_Data_Par_IBinwithout a prototype (implicit declaration). The call works only by luck of the calling convention.src/models/check_state.c:84prints anunsigned intwith%f.src/riversys.c:361-363: misleading indentation. The secondprintfruns on all processes.
B-11¶
Unchecked input reading. Low → medium for users. About 75 calls to fscanf/fread
(mostly src/riversys.c, src/processdata.c) ignore the return value. A truncated or
malformed .rvr/.prm file is not reported and leaves variables uninitialised.
B-12¶
.uini reader does not detect missing values. Medium, code reading.
src/riversys.c:1163-1170 reads no_ini_start values and checks
if (fscanf(...) == 0). At end of file fscanf returns EOF (−1), not 0, so a file
with too few values is accepted, and the missing states keep whatever was in memory
(malloc does not zero). The header’s model id (%*i) is ignored as well.
examples/clearcreek.uini is such a file: its header says model 252 and it gives
4 values, while model 254 reads 7 (no_ini_start = dim). It works only because
ReadInitData for model 254 happens to overwrite exactly the 3 missing states (4, 5, 6).
Fixed (2026-09-25): missing values are detected. ASYNCH warns and sets them to 0 (the buffer is
zero-initialised); a value that is not a number is an error; a model number that differs from the
.gbl gives a warning. The check revealed that examples/more/common/test.uini (3 values) is also
short for models 258 and 259, which read 4. Their subsurface storage started from uninitialised memory
in the original code, which happened to be 0. It is now 0 by design, and results are bit-identical.
examples/clearcreek.uini now has the right model number and all 7 values.
B-13¶
RK 3(2) and RK 4(3) read their coefficients from freed memory. Critical for users of
methods 0 and 1, confirmed. Found while testing the fix of B-04.
RKDense3_2 and TheRKDense4_3 (src/solvers/rk3_2_dense.c, src/solvers/rk4_3_dense.c)
stored in the RKMethod struct pointers to their local coefficient arrays (A, b, c, d,
e). Local arrays live on the stack and disappear when the function returns, so every time step
then read whatever else happened to be on the stack. AddressSanitizer:
stack-use-after-return ... in ExplicitRKSolver src/steppers/explicit.c:120, pointing at
array c of TheRKDense4_3. With method 1, examples/test.gbl ran forever in most runs (with 1
or 2 processes, in the original code too). Any result computed with solver index 0 or 1 by
earlier versions is unreliable. Method 2 (Dormand–Prince), used by all examples, stores its
tables as static and was never affected.
Fixed (2026-09-25): the tables are static. Methods 0 and 1 now finish, are free of sanitizer
reports and bit-reproducible with 1 process. They agree with method 2 within the expected accuracy
(clearcreek, 6 359 links: peak discharge differs by at most 7.7e-4 m³/s, median relative
difference 0.1 %).
B-14¶
.rkd files (tolerances and method per link) could not be used. High for users of this
option, confirmed. Found while fixing B-04 and B-05. With a solver flag of 1 in the .gbl, the
original code never got past “Reading dam and reservoir data…”. Build_RKData
(src/riversys.c) had five defects:
the loops reading the tolerances incremented
iinstead ofj(for (j = 0; j < num_states; i++)), an endless loop running off the arrays;the array of method indices had
num_statesentries instead of one per link;the broadcast to the other processes sent
methods(the table of RK methods) instead of the method indices, overwriting memory;the per-link
ErrorDatastructure was never allocated before being written to;the link ids in the file were read but not used: rows were assumed to be in network order.
Fixed (2026-09-25): the reader was rewritten. Rows are matched to links by id; every value is
checked, and a bad file stops the run with a clear message (unknown or repeated id, missing link or
value, invalid method, too few tolerances for the model). The format is now documented in
docs/input_output.rst. examples/test_rkd.gbl + examples/test.rkd repeat the settings of
examples/test.gbl for every link. They give bit-identical results with 1 process, and are part of
the regression tests.
B-15¶
Unwritable outputs are not fatal. High, confirmed (while testing the documentation). If a
folder named in the .gbl for hydrographs, peaks or snapshots does not exist, every write prints
Error: could not open h5 file ... / Error: Cannot open peakflow file ... (27 lines for test.gbl), but the
simulation runs to the end and asynch exits with code 0, the code for success, having written no results.
In a script or an operational chain, the failure goes unnoticed. The simulation time is also wasted.
Fixed (2026-09-25): Read_Global_Data (src/config_gbl.c) checks that the folders of the hydrograph, peak,
snapshot and temporary files exist and are writable, and stops before computing otherwise. main
(src/asynch_cli.c) now checks the return values of the final output functions, and exits with EXIT_FAILURE if
one failed.
B-16¶
Links with many parents overflowed memory. High, confirmed (while testing B-06 on a synthetic network).
The time-step routines (src/steppers/*.c) keep the parents’ data in an array of ASYNCH_LINK_MAX_PARENTS = 8
entries, but the network reader (Create_River_Network, src/riversys.c) accepted up to 10 parents, and stored
each parent before checking the limit. A link with 9 or 10 parents made every time step write past the array
(AddressSanitizer: stack-buffer-overflow in ExplicitRKSolver, the run aborted). With 11 or more parents the reader
itself wrote past its buffer.
Fixed (2026-09-25): one limit, ASYNCH_LINK_MAX_PARENTS = 16, is used by the reader and the solvers; both readers
(file and database) check it before storing. A network with 70 000 links and 10 parents per main-channel link now
runs cleanly (release and sanitizer builds, 1 and 2 processes). One with 21 parents stops with
Error: link 1 has 21 parents; ASYNCH supports at most 16 (ASYNCH_LINK_MAX_PARENTS in src/constants.h).
B-17¶
Several library functions crashed or did not exist. Medium, confirmed (while writing the Python package).
They only concern programs that use ASYNCH as a library; the asynch program does not call them.
Asynch_Get_Size_Global_ParametersandAsynch_Set_Global_Parametersreadasynch->model, which only exists for custom models: with a built-in model they dereferenced a null pointer.Asynch_Set_Global_Parametersalso did not update the count kept inglobals, soAsynch_Get_Global_Parameterscopied the old number of values.Asynch_Set_Total_Simulation_Durationwas declared inasynch_interface.hbut never written: a program calling it did not link.Asynch_Set_System_Statepassed the discontinuity state wherecheck_stateexpects the dam flag.Asynch_Custom_Modelallocated a model that it immediately replaced (a leak), and a model created byAsynch_Custom_Partitioningalone (no equations) made the global-file reader call a null function.Asynch_Set_Init_Fileread before the start of names shorter than 4 characters, did not accept.h5initial states, and wrote into a null pointer when the global file took the initial states from a database.
Fixed (2026-09-25) in src/asynch_interface.c, src/config_gbl.c and src/riversys.c: the counts come from
globals, the missing function exists, the flags are passed in the right place, a custom model is recognised by
its callbacks (so a partitioning-only “model” keeps the built-in equations), and the file type is taken from the
extension (.ini, .uini, .rec, .h5).
B-18¶
Initial-state readers stored the discontinuity state on the wrong link. Medium, confirmed by reading.
In Load_Initial_Conditions (src/riversys.c), the two .ini branches computed the state of the link at
location loc but stored it in system[i], where i is only a loop counter. The .rec, .dbc and .h5 branches
passed the state where check_state expects the dam flag. Only models with discontinuity states (the dam models) are
affected: their first step could start in the wrong regime. The shipped examples have no dams and do not change.
Fixed (2026-09-25).
B-19¶
Model 190 read a forcing that does not exist. Low, confirmed by reading. LinearHillslope_MonthlyEvap
(src/models/equations.c) read forcing_values[2] into an unused variable, but model 190 has two forcings (the
array holds two values). The optimiser removes the unused read, so results were never affected; the line, added in
2017 together with model 195, is removed.
B-20¶
Custom time series outputs of states that are not interpolated were wrong. High, confirmed. A program can add
its own outputs (Asynch_Set_Output_Int/Double/Float) and say which states they use. Those states must be
interpolated (“dense output”) at the print times. The code added the states that were already interpolated, and
skipped the others; those were then written as whatever was in memory. Test: an output returning state 1 of model 190,
with State1 not otherwise printed, wrote 0 at every time instead of values around 0.001. Also, when called at the
documented moment (right after reading the global file), the links did not exist yet and nothing was added at all.
Fixed (2026-09-25): the used states are recorded before the model is initialised and added afterwards; the output
now equals State1 to its print precision (1, 2 and 3 processes). The same off-by-one (> instead of >= when checking
that a state exists) was fixed in Initialize_Model.
B-21¶
A custom model inherited the snapshot filter of the built-in model with the same number. Low, confirmed by reading.
SetOutputConstraints chooses a filter for .h5 snapshots from the model number of the global file, even for a
custom model, whose states may mean something else. Custom models now get no filter.
B-22¶
Wrong constant in the derivative of the Dormand-Prince dense output. Low, confirmed by a unit test.
DOPRI5_bderiv (src/solvers/dopri5_dense.c) had 1144640195640 where 32805 x 3489224 = 114463993320 belongs
(the fifth coefficient). The derivative is only evaluated at theta = 1 (src/steppers/explicit_index1_dam.c), where
this term is multiplied by 0, so no result was affected. Fixed (2026-09-25).
B-23¶
Models 263, 601, 602 and 603 wrote and read past the parameter array of every link. High for users of these
models, confirmed by a unit test. SetParamSizes (src/models/definitions.c) declared fewer parameters per link
(num_params) than it reads from disk (num_disk_params): 14 < 15 (263 and 601), 16 < 17 (602), 20 < 21 (603). The
reader stores every value read into an array of num_params doubles, so the last one was written past it. The
precalculations of 601-603 then read it (v_0, used for invtau) from there, and the equations of 263 read two values
(v_B, k_tl) past the array: results depended on whatever memory followed. Fixed (2026-09-25): the arrays have
room for every value read (15, 17, 21; 16 for model 263, whose equations use 16 values, all read from disk: its
parameter files or database queries must give 16 values per link).
B-24¶
Seven model numbers crashed at the first step; the Runge-Kutta tables were shared by all solvers. Medium, confirmed by unit tests.
Models 200, 260, 300, 301, 315, 607 and 2000 have sizes in
SetParamSizesbut no equations (200 is meant for another program; 260’s equation is commented out; the others have noInitRoutinesbranch). A run called a null function.Initialize_Modelnow stops withError: model N cannot be integrated by ASYNCH ... (no equations).Build_RKData(src/riversys.c) built the tables into astaticarray, shared by every solver of the program. With two solvers (easy from Python), creating the second rebuilt the tables of the first. Each solver now owns its array, andDestroy_RKMethodfrees what the constructors allocate (it freed nothing; RK 4(3) pointedbto a static table, so it could not be freed consistently).
B-25¶
Rain from binary files: a file past the range was read, the last file lasted 0.0001 min, memory overflow, crash on
a missing file. High, confirmed by a test (the same rain as a .str file and as binary files). Binary forcing
files (global file flags 2 and 6: one file per time step, used for radar rainfall) are read in passes of chunk size
files. In src/forcings.c the last file of a pass was capped at last + 1, a file outside the declared range, and
src/forcings_io.c:
read it (flag 2: from a NULL file pointer if it did not exist, which crashed; flag 6: stopped);
gave the value of the last file for 0.0001 min only, then 0 (it was only right when the file after the range existed);
allocated
number of files + 1values per link, but wrotechunk size + 1: a last pass with fewer files than the chunk size wrote past the array.
Fixed (2026-09-25): files first to last are read, the last one applies for a full time step, then 0 (as
documented); the arrays are large enough; a missing file stops the run with its name. Test
(tests/python/test_forcings.py): rain different at every link, written as .str, binary, gzipped binary and
irregular binary files, gives identical results. Runs whose files covered one step more than the declared range differ
only by the 0.0001 min of that extra file.
B-26¶
Models 105 and 263 did not set every derivative: states took values left in memory. High for users of these
models, confirmed with valgrind and a unit test. A model function must write the derivative of each of its states
into ans. river_rainfall_summary (model 105, 2 states) writes only the first; model263 (8 states) adds the
parents’ state 4 to ans[4] without setting it first, and never writes ans[5] to ans[7]. The solver passes a
work array (Create_Workspace, allocated with malloc) that holds values of an earlier stage or of another link, so
these states changed by amounts that depended on memory contents. It was found because one run of make check from
the source archive (built with -O2) did not finish in the loop over all models, right after model 263; valgrind then
showed the uninitialised values, and the models were checked one by one.
Fixed (2026-09-26): states without an equation keep their initial values (ans[k] = 0); the sum of model 263
starts at 0. The intended equations of these states are an open question (S-07). The unit test that evaluates every
model’s equations now fills ans with NaN first and fails if any derivative is left unset: it reports exactly models
105 and 263. The results of these two models change (they were undefined before); no example uses them.
B-27¶
The consistency check of the steppers read the parents’ non-dense states uninitialised. Low, confirmed with
valgrind; no effect on results. When a link steps, the states of its parents are interpolated at the stage times, but
only the dense ones (for model 254: states 0 and 6); ExplicitRKSolver (src/steppers/explicit.c:88) then applies
the consistency check to the whole state vector of each parent, reading entries that were never written. The equations
never use those entries. Fixed (2026-09-26): the work arrays are allocated with calloc (src/system.c). All
examples give bit-identical results.
B-28¶
Models 257, 258, 259, 261 and 262 evaporated ponded water 1000 times too fast. High for users of these models,
confirmed with a unit test and the examples of models 258 and 259. The Top Layer models share the potential
evaporation e_pot [m/min] between the ponded water, the top soil and the deeper soil, so that the three parts add up
to e_pot. In these models the ponded part was e_p = s_p * 1e3 * e_pot / Corr: no conversion factor belongs there
(the storages are in m, e_pot in m/min). Whenever water was ponded it evaporated up to 1000 times too fast, total
evaporation could exceed the potential evaporation of the input file, and surface runoff and flood peaks were too low.
History: the factor was added on 2015-06-22 (commit fb21cb1, “fixed evaporation bug in models”) to every Top Layer
model, model 254 included, and spread to later models copied from them. It was then removed one model at a time:
263 (2019, 547f40c), 254 (2020-04-23, 14054bf, “I erase the 1000 from the ponded evaporation”), 252 (2020),
255 (2022). Every release from 1.0.0 to 1.4.3 has it in model 254; the old branches assim and json still do. The
clearcreek reference results (uploaded 2015-05-21) predate it.
Fixed (1.6.0): e_p = s_p * e_pot / Corr in all five models. A unit test checks, for 15 Top Layer models, that the
evaporation taken from the storages equals the potential evaporation. The results of models 257–262 change; the
examples of models 258 and 259 now differ from their 2018 benchmarks, which were produced with the factor: outlet
peaks rise by up to 23 % (the benchmark files are kept unchanged; see R-03).
B-29¶
Model 255 turned channel storage into 8 to 15 times too much discharge. High for users of model 255, confirmed with
a unit test. Model 255 computes the discharge from the channel storage S (dam_model255). Since 2022-07-24 (commit
82fcfc1, titled “changes to model 402”) it used invtau/60 * S^(1/(1-λ₁)), which is not in m³/s; the correct
relation, ((1-λ₁) invtau/60 · S)^(1/(1-λ₁)), had been written in 2021 (1d793c1) and was left commented out on the
next line. It is also the relation used by models 261 and 262 and by the initial conditions of model 255 itself, so the
first step of every run jumped. For v₀ = 0.33 m/s, λ₁ = 0.2, λ₂ = −0.1 the error is ×7.8 at 1 km², ×9.8 at 10 km²,
×14.5 at 1000 km²: flood waves travelled too fast and too sharp.
Fixed (1.6.0): the correct relation is restored. A unit test checks, for models 255, 261 and 262, that the discharge computed from the storage gives back the discharge the storage was computed from.
B-30¶
Models 258 and 259: the baseflow equation read the accumulated evaporation. Medium, confirmed with a unit test.
The baseflow is state 7 (ReadInitData initialises it, the equation writes ans[7]); state 6 is the accumulated
evaporation. The equation read q_b = y_i[6], added the parents’ state 6, and passed state 6 downstream
(dense_indices). The baseflow output was meaningless; the discharge (state 0) does not use it and was not affected.
The mistake came with the models in 2018. Fixed (1.6.0): the baseflow is read from state 7 at the link and at its
parents, and state 7 is passed downstream.
B-31¶
Model 249: the baseflow mixed m³/min and m³/s; its reservoir version printed at every evaluation. Medium, code
reading. The baseflow q_b is routed like the discharge, but the local terms were in m³/min
(q_sl · A_h − 60 q_b) and the parents’ baseflow in m³/s, a factor 60 between them. The reservoir version printed two
lines each time it was evaluated. Fixed (1.6.0): every term in m³/s (q_sl · A_h / 60 − q_b + Σ q_b,parents);
the prints are removed. See also B-36.
B-32¶
Model 257: the accumulated evaporation output was 720 times too large. Low (output only), code reading. State 5
accumulated forcing_values[1] * c_1: the evaporation forcing (mm/month) converted with the rain factor (mm/h → m/min).
Fixed (1.6.0): it accumulates the potential evaporation e_pot [m/min], as models 258 and 259 do.
B-33¶
Models 0–6, 105, 200 and 2000: wrong area in .pea files. Low (output only), code reading. These models declared
convertarea_flag = 1 (areas converted to m²) but never convert the upstream area, which stays in km². The peak-flow
file multiplies the area by 10⁻⁶ for such models, so it showed the area in km² × 10⁻⁶. Fixed (1.6.0): the flag is
0 for these models; the area is written in km², as for every other model.
B-34¶
Models 225 and 601–609 gave NaN when every hillslope storage was empty. Medium, confirmed with a unit test. The
evaporation is shared with weights 1/(C_p + C_l + C_s), without the check model 254 has: with every storage at 0
(a dry start, or storages emptied and reset to 0) this is 1/0 and the derivatives became NaN, which stops or corrupts a
run. Fixed (1.6.0): no evaporation when the sum is exactly 0; results are unchanged otherwise. A unit test
evaluates every model with empty storages.
B-35¶
Model 606 left derivative 5 unset when the tile storage was empty (the same kind of bug as B-26). Low, confirmed with a unit test. Fixed (1.6.0): the derivative is set in every case (0 when the storage is empty).
B-36¶
Model 249 with reservoirs does not work. Medium, code reading, open. At links with a reservoir,
ForcedSolutionSolver uses the output of model249_reservoirs as the new states, but the function returns
derivatives for states 1 to 5, leaves state 0 unset when the reservoir forcing is 0 or less, and routes the total
discharge as the open-loop discharge. It looks like an unfinished data-assimilation experiment (2021). What was
intended needs the model’s author; model 249 without reservoirs is not affected.
Reproducibility¶
R-01¶
No regression testing. High. make check runs a single unit test (days_in_month).
Nothing checks that the model still produces the same hydrographs. Addressed by
tests/regression/run_examples.py (see 09_reproducibility.md). Since 2026-09-25
make check runs 23 C unit tests (tests/check_asynch.c), 73 tests of the Python package (tests/python) and the
9 example comparisons; the tests found B-22 to B-25. Line coverage: 66.5 % (chapter 9).
R-02¶
Clearcreek reference is from another configuration, and model 254 changed in 2021. Medium, confirmed.
examples/results/clearcreek.pea (with .dat and .rec) was committed in May 2015 (commit
b73fc2d, the first import). For the outlet (link 2527) the reference peak is at 3001 min,
while clearcreek.gbl only simulates 1440 min (one day), so the reference cannot validate
today’s example as it stands.
Explained (2026-09-25). The 2015 global file (Global254.gbl in b73fc2d) simulated 6000
minutes starting 2014-05-01. All input files are byte-for-byte the same today. That configuration
is now examples/clearcreek_2015.gbl. Running it:
with today’s code, the references are not reproduced (hydrograph differences up to 0.085 m³/s; zeros where the reference has small baseflow values);
with today’s code and one line of
model254restored to its 2015 form (double q_b = y_i[6];instead ofmax(0.001, y_i[6])), hydrographs and final states agree with the references within the solver tolerance (largest difference 9.4e-5, solver abs tolerance 1e-4), and peak discharges within 1.7e-4 m³/s.
So the difference is model 254 itself: the floor max(0.001, q_b) was added in January 2021, in commit
93241a3, whose message is “added model 194”. That change of the operational model is not
mentioned anywhere. See S-02. The reference files are kept unchanged.
Resolved (2026-09-25): the 2015 equation was restored (S-02). examples/clearcreek_2015.gbl now
reproduces all three reference files: hydrographs within 9.5e-5, final states within 1.5e-5, and peak values
within 1.7e-4 m³/s (only 9 of 6 359 links above 1e-4; their peak times moved by 3–7 minutes, because peaks
are recorded at solver steps). examples/clearcreek.gbl itself still simulates a different period (one day
from 2017-01-01), so its comparison with the 2015 file stays a known mismatch by design.
R-03¶
Model 259 benchmark cannot be reproduced. Medium, confirmed. The 2018 commit that
added the benchmark (cba763b) was built and run: it produces output bit-identical to
today’s code, and both differ from the benchmark (outlet peak 0.696 vs 0.755 m³/s).
So the code has not changed. The benchmark was produced with an input that is not
in the repository, most likely model 259’s own evap.mon on the original cluster
(/Dedicated/IFC/.../mdl259a/evap.mon). With zero evaporation the peak is 0.845, so the
original file lies between the two. The benchmark files are kept unchanged; the case stays
a known mismatch unless the original evaporation file is found.
Update (1.6.0): the benchmark was also computed with the evaporation of ponded water multiplied by 1000 (B-28) and the baseflow read from the wrong state (B-30). Both are fixed, so today’s results differ from the benchmark for these reasons as well; the benchmark stays as it is.
R-04¶
Examples 258/259 referenced a cluster path. Fixed. Their .gbl
files pointed to /Dedicated/IFC/projects/asynch_1_4_3b/tests/mdl25Xa/evap.mon. They now
use ../common/evap.mon, like the other examples. With this change, model 258 reproduces its
benchmark bit for bit. Update (1.6.0): no longer, on purpose: the benchmark was computed with the evaporation
error B-28 and the baseflow error B-30, which are fixed. The regression harness declares the difference as intended.
R-05¶
Results depend (slightly) on the number of processes. Information. ASYNCH is asynchronous: each link has its own adaptive time step, and with several processes the order of computation and the data available from upstream change. Observed differences: relative ≤ 3·10⁻⁴ on small values, ≤ 10⁻⁵ on peaks. That is below the solver tolerances, so it is not a bug, but it means results should be compared with a tolerance, never byte by byte.
Python API¶
A-01¶
The Python API does not work and cannot be repaired incrementally. High, confirmed.
py/is not part of the build (SUBDIRS = src tests toolsinMakefile.am).py/asynch_interface_py.cdoes not compile against the current headers (unknown type name 'MPI_Comm',incomplete typedef 'AsynchSolver', …).py/asynch_interface.pyandasynchdist.pyuse Python 2 syntax (print 'x').The library path is hard-coded:
ASYNCH_LIBRARY_LOCATION = '/home/ssma/NewAsynchVersion/libs/libasynch_py.so'.14 of the 60 C functions it calls no longer exist (for example
Asynch_Get_Number_Links,Asynch_Set_Output,Asynch_Get_Total_Simulation_Time).It re-declares C structs (
UnivVars, …) inctypeswith the old field layout. Even if it loaded, it would read and write memory at the wrong offsets.
Resolved (2026-09-25): replaced by the package in python/ (chapter 10), a ctypes binding over the
shared library libasynch.so and a small C interface (src/asynch_api.h) that exposes only numbers, arrays and
opaque handles, never structure layouts. py/, asynchdist.py and asynchdist_custom.py were removed; the custom
model of asynchdist_custom.py is ported in examples/python/custom_model.py and reproduces model 191 exactly.
A-02¶
The Python package made HDF5’s symbols global: h5py could not be imported after asynch. Medium, confirmed.
asynch/_lib.py loaded libasynch.so with RTLD_GLOBAL (for Open MPI, whose plugins need the MPI symbols). That
also made the symbols of every library it depends on global, HDF5 among them. h5py, which carries its own HDF5, then
failed to import in the same program (ValueError: Not a datatype), and a library linked with another MPI could be
handed the wrong MPI functions. Found while testing the ready-made wheel. Fixed (2026-09-26): the library is
loaded with RTLD_LOCAL, and only Open MPI’s library is made global (RTLD_NOLOAD), as mpi4py does. Tested: Open MPI
build (make check, 2 and 3 processes, mpi4py) and the MPICH wheel (h5py imported after asynch).
Maintainability¶
M-01¶
Dead code. These files are in the repository but not compiled (src/Makefile.am):
file |
lines |
what it is |
|---|---|---|
|
3132 |
old version of solvers/steppers (old |
|
1215 |
old forcing reader |
|
666 |
old custom-model example |
|
477 |
stepper for discontinuous models |
|
464 |
Radau implicit stepper (see B-04) |
|
279 |
data-assimilation stepper |
|
239 |
old output routines |
|
47 |
That is about 6 500 lines, 15 % of the C code. Also in src/models/definitions.c,
ReadInitData, the branches for models 200, 254, 255, 256 and 257 at lines 3629-3704
are unreachable duplicates, because the same if/else chain already matched those models above.
Resolved (2026-09-26): the eight files were deleted (git keeps their history), together with other unused
files: the Visual Studio projects (ide/, which referred to files that no longer existed), an old conda recipe, the
Travis CI and Read the Docs settings (replaced by GitHub Actions and GitHub Pages), cluster job scripts, an editor
workspace, a personal build script, the old LaTeX manual (replaced by the .rst reference) and a generated
Makefile.in. The unreachable branches of ReadInitData are still there.
M-02¶
A model is spread over 7 places. Adding or reading model N means finding its
case/if in SetParamSizes, SetOutputConstraints, ConvertParams,
InitRoutines, Precalculations, ReadInitData (all in definitions.c, 4 063
lines), plus its equations in equations.c (6 293 lines). Mistakes like B-01 are a
direct consequence. Possible improvement: one descriptor per model (struct
with sizes + function pointers) in one file per model family.
M-03¶
Tooling. .travis.yml targeted travis-ci.org, which shut down in 2021, so there
was no working CI. Fixed (2026-09-26): .travis.yml was removed; GitHub Actions now build ASYNCH and run
make check on every push (.github/workflows/tests.yml), and build and publish this documentation
(.github/workflows/docs.yml).
The root .gitignore also lists examples (and *.rvr, *.str, …), although those files
are tracked. git add examples/... therefore refuses to stage changes to the examples, and
git add -u (or -f) is needed. Recommendation: ignore only generated outputs (examples/**/results/).
Performance¶
These are hypotheses from reading the code. They must be measured (profiling a large network) before anything is changed.
P-01¶
src/asynch_cli.c:320 sleeps for 1 s (ASYNCH_SLEEP(1)) before the computation.
Negligible for large runs, but it makes up 99 % of the runtime of the small examples.
(On Windows the same macro sleeps 1 ms, because Sleep takes milliseconds.)
The pause came with the first upload (2015), right after each process prints “good to go” and before a barrier:
most likely to let those lines reach the screen before the next message. It has no numerical role: without it every
example gives bit-identical results with 1 process, and differences with 2 and 4 processes are of the same size as
between two runs of the unchanged program (R-05). Fixed (1.7.0): the pause is removed; output is flushed before
the barrier. Initialisation of the 6 359-link test network: 1.3 s → 0.27 s. The pauses left in
src/asynch_interface.c and src/processdata.c are on error paths (they let process 0 print its message before the
program aborts) and stay.
P-02¶
Recurrent HDF5 snapshots (DumpStateH5, src/processdata.c:1893-1936) send one
synchronous MPI message per link to process 0. For 400 000 links that is 400 000
round-trips per snapshot. A single MPI_Gatherv would do the same work in one collective call.
P-03¶
In Advance (src/advance.c):
the next link to compute is found with a linear scan (
aroundloop),O(my_N);two
MPI_Barriers plus a duplicatedTransfer_Data_Finishrun at every forcing period (the code itself says “This is sloppy”);InitialStepSizeis recomputed for every link at every forcing period.
P-04¶
The explicit solvers are limited by stiffness. Medium, confirmed (profiling, model 254, 6 359 links, 100 h). Some
states react much faster than others (ponded water on small hillslopes drains within minutes), so the explicit
methods must take steps of the order of the fastest reaction. Making the tolerances 100 times looser removed only 9 % of
the steps of Dormand–Prince, and 21 % of its steps were rejected. By instruction count, 67 % of the run is the solver’s
own work and 28 % the model equations (pow alone 17 %); input/output and MPI are below 2 %. Resolved by an option
(1.7.0): numerical solver index 4, Rodas5P, a Rosenbrock method whose steps are limited only by accuracy. With the same
tolerances it is 16 times faster; with tolerances 100 times smaller, 6.6 times faster with similar accuracy
(chapter 4, section 4.7). Index 2 stays the default.
P-05¶
The explicit solvers record peaks only at the end of a step. Low, confirmed. The peak flow of a link is updated
with the state at the end of each accepted step (ExplicitRKSolver), so a crest between two step ends is missed. The
2015 reference peaks are lower than those of a run at tolerance 10⁻⁸ at 10 of 11 links (test) and 5 270 of 6 359 links
(larger example), by up to 4.3·10⁻⁴ m³/s. The stiff solver (index 4) searches the crest inside each step with its dense
output. The explicit solvers are left unchanged, so that their results stay identical to the reference results.
Scientific review items (open questions)¶
These are not necessarily bugs. They are places where a modelling choice is hidden in the code and should be written down, and possibly revisited, by a hydrologist.
S-01¶
In the model equations, the parents’ states are read as y_p[i * dim]
(e.g. model254, src/models/equations.c:1667-1684), but the solver stores them with
stride max_dim (src/steppers/explicit.c, stages_parents_approx + i * max_dim).
Both are equal as long as every link has the same number of states, which is true today.
It becomes silently wrong if links with different dimensions are ever mixed (for example
reservoirs or coupled models). Recommendation: use max_dim (it is already passed to every function).
S-02¶
Model 254’s baseflow equation uses q_b = max(0.001, y[6]) in its outflow term
(equations.c:1638). When the baseflow is below 0.001 m³/s, the channel still loses baseflow as if
it were 0.001 m³/s, so it drains faster than the linear reservoir would and is then clamped to 0 by
check_consistency. This changes the water balance at low flow.
This line was not part of the original model: it was added in January 2021 (commit 93241a3,
“added model 194”) and changed model 254’s results (see R-02). With the 2015 form (q_b = y[6]),
today’s code reproduces the original repository’s clearcreek references within the solver tolerance.
Resolved (2026-09-25, owner’s decision): the 2015 form q_b = y[6] is restored. The original
reference results for clearcreek (examples/results/clearcreek.dat, .pea, .rec) are now reproduced
within the solver tolerance.
S-03¶
Potential evapotranspiration is converted from mm/month to m/min with a fixed 30-day
month (equations.c:1621). February and 31-day months are off by up to 7 %.
S-04¶
Models 400–405 test if (temperature != 0 & temperature < temp_thres)
(equations.c:2514, 2676, 2834, 3029, 3296, 3437). A temperature of exactly 0 °C
(probably used as “no data”) is never treated as snow. & (bitwise) is used instead
of && (logical); the result is the same here, but it looks accidental.
S-05¶
The snapshot filters (src/models/output_constraints.c) set values below 10⁻¹² to 0
and replace values above 10²⁰⁰ with fmod(x, 1e200). Snapshots are used as initial
conditions, so this silently changes the state of the next run. That is fine for
round-off negatives, but the fmod behaviour should be justified or removed.
S-06¶
In model 254 (equations.c:1642-1654) evaporation is split between the storages in
proportion to their relative fullness, so e_p + e_t + e_s = e_pot whenever any
storage holds water, however little. Near-empty storages are then over-drawn,
become negative, and are clamped back to 0 by CheckConsistency_Nonzero_AllStates_q,
which silently adds water. The water balance does not close in dry periods. A
common remedy is to scale actual ET by availability (e.g. min(e_pot, storage/Δt) or
a smooth factor s/(s+ε)), but this changes the model and must be a deliberate,
documented decision.
S-07¶
Models 105 and 263 have states without equations (B-26): the storage of model 105, and states 5 to 7 of model 263,
whose equations are commented out in src/models/equations.c. Since 1.5.0 these states keep their initial values.
What the authors intended (for model 263, state 7 is even passed to the downstream links) needs the model’s author or a
hydrologist.
S-08¶
Model 603 (VariableTriLayer): the lateral flow of the deepest soil layer, qh4 = v4 · s3^a4, is removed from that
layer but not added to the channel (ans[0] adds qh1 + qh2 + qh3). Water disappears. It may be a deliberate deep loss,
or a forgotten term; the author must say.
S-09¶
Model 605: the intercepts that keep the piecewise-linear subsurface runoff continuous are I1 = v_s1 (S2 − S1),
I2 = v_s2 (S3 − S2) + I1, I3 = v_s3 (S4 − S3) + I2 + I1 (definitions.c, Precalculations). I2 already contains
I1, so the runoff jumps by I1 when the storage reaches S4. Probably + I2 only was meant.
S-10¶
Models 1–5 (2011, “September 18, 2011 document”): the hillslope loses c₄ s^(5/3) [m/min] but the channel receives
c₁ s^(5/3), and c₁ contains an extra factor h_b^(2/3) compared with c₄ · A_h/60/Q_r. Water is not conserved
unless the state s has a meaning that the code does not show. Needs the 2011 document.
S-11¶
Models 400–405: several lines add or compare a storage [m] and a rate [m/min] (x2 = max(0, x1 + h1 − Hu),
min(e_pot, h1), min(h5, melt rate)), which is only meaningful because the time unit is 1 minute (a storage above
capacity is drained in one minute). At links with a dam, models 402 and 403 do not add the runoff of the link’s own
hillslope to the reservoir (model 405 does); model 402 uses one hard-coded storage-discharge curve for every dam.
S-12¶
Constants whose unit or origin is not explained:
model 604 uses the runoff speed column
v_ras read, models 605–609 convert the same column to 1/min (v_r · L/A_h · 60);model 606 multiplies tile outflow by 2500;
model 30: two thresholds on
deriv_a_Ilack the factor 10⁶ that the author’s own commented line has (“Should use this”);model 21: the groundwater fraction
F_etis 0.05 on links without a dam and 1.0 on links with a dam.
Documentation errors¶
D-01 to D-03¶
docs/builtin_models.rst, section Top Layer Hydrological Model (model 254):
D-01 1/τ shows
L · 10⁻³with L in km. It should beL · 10³. The code is correct.D-02 the baseflow equation shows
+ q_b,in(t). The code has+ 60·q_b,in(dimensionally consistent).D-03 V_r is described as m³/s. It is an accumulated depth in m.
Fixed (1.6.0) in docs/builtin_models.rst.
Details: 05_model_254_explained.md §5.7.
D-04 to D-07¶
Found in the units check of all models (2026-09-27):
D-04 the same
L · 10⁻³as D-01 in the sections of models 190 and 21.D-05 model 191: potential evaporation given in mm/hour; the code reads mm/month. A user following the page got evaporation 720 times too small.
D-06 model 21: 1/τ defined like the other models; the code uses
60 · (v_r (A/A_r)^λ₂ / L)^(1/(1−λ₁)), the form that goes with its storage-discharge relation.D-07 models 400–405: code comments gave the melt factor in mm/hour/degree in one place, mm/day/degree in another; the equations use mm/day/degree.
Fixed (1.6.0).