TELEMAC backend

axqua.solvers.telemac builds, runs and evaluates TELEMAC models. The modules appear in the order of the preprocessing steps, followed by the evaluation of runs and the extensions.

Build sequence

End-to-end orchestration of the TELEMAC case-build pipeline.

class axqua.solvers.telemac.pipeline.Artifacts(geometry_slf: 'Path | None' = None, boundary_cli: 'Path | None' = None, friction_tbl: 'Path | None' = None, initial_conditions: 'Path | None' = None, cas_file: 'Path | None' = None, unsteady_cas: 'Path | None' = None, liquid_boundaries: 'Path | None' = None, liquid_boundaries_meta: 'Path | None' = None, gaia_cas: 'Path | None' = None, source_regions: 'Path | None' = None, user_fortran: 'Path | None' = None, ground_truth: 'Path | None' = None, calibration_csv: 'Path | None' = None, hbc_config: 'Path | None' = None, rasters: 'dict[str, Path]'=<factory>, surfaces: 'list[Path]' = <factory>)

Bases: object

Parameters:

• geometry_slf (Path | None)

• boundary_cli (Path | None)

• friction_tbl (Path | None)

• initial_conditions (Path | None)

• cas_file (Path | None)

• unsteady_cas (Path | None)

• liquid_boundaries (Path | None)

• liquid_boundaries_meta (Path | None)

• gaia_cas (Path | None)

• source_regions (Path | None)

• user_fortran (Path | None)

• ground_truth (Path | None)

• calibration_csv (Path | None)

• hbc_config (Path | None)

• rasters (dict[str, Path])

• surfaces (list[Path])

geometry_slf: Path | None = None

boundary_cli: Path | None = None

friction_tbl: Path | None = None

initial_conditions: Path | None = None

cas_file: Path | None = None

unsteady_cas: Path | None = None

liquid_boundaries: Path | None = None

liquid_boundaries_meta: Path | None = None

gaia_cas: Path | None = None

source_regions: Path | None = None

user_fortran: Path | None = None

ground_truth: Path | None = None

calibration_csv: Path | None = None

hbc_config: Path | None = None

rasters: dict[str, Path]

surfaces: list[Path]

axqua.solvers.telemac.pipeline.run(cfg, *, validate_env=True, dry_run=False, log_to_file=True)

Build the full case. Returns the produced Artifacts.

Parameters:

• validate_env (check that the TELEMAC environment can be sourced first.)

• dry_run (after building, launch the solver once to confirm the case is valid.)

• log_to_file (attach this build's ``<model_dir>/axqua.log` (set False`) – when a caller already routes everything to its own compound log, e.g. the mesh-convergence study logging the per-level builds to postprocessing/).

• cfg (Config)

Return type:

Artifacts

Mesh and mesh quality

Stage 2 - mesh generation, bathymetry, and SELAFIN geometry.

Builds a triangular mesh with gmsh from the ROI boundary polygon and the breaklines, with per-region (MATID) size refinement, then:

• extracts nodes + triangles,

• walks the mesh boundary to build the ordered contour (drives both TELEMAC’s IPOBO array and the .cli row order - they must agree),

• interpolates the clipped DEM onto the nodes (bathymetry),

• assigns each element a MATID (friction zone),

• writes the geometry .slf via axqua.core.selafin.

The boundary contour and the per-node/element classification are returned so the boundary-condition stage can tag inflow/outflow nodes consistently.

class axqua.solvers.telemac.mesh.Mesh(x: 'np.ndarray', y: 'np.ndarray', triangles: 'np.ndarray', bottom: 'np.ndarray', ipobo: 'np.ndarray', boundary_nodes: 'np.ndarray', element_matid: 'np.ndarray', node_matid: 'np.ndarray', boundary_loops: 'np.ndarray | None' = None, roughness: 'np.ndarray | None' = None, quality: 'object | None' = None)

Bases: object

Parameters:

• x (np.ndarray)

• y (np.ndarray)

• triangles (np.ndarray)

• bottom (np.ndarray)

• ipobo (np.ndarray)

• boundary_nodes (np.ndarray)

• element_matid (np.ndarray)

• node_matid (np.ndarray)

• boundary_loops (np.ndarray | None)

• roughness (np.ndarray | None)

• quality (object | None)

x: np.ndarray

y: np.ndarray

triangles: np.ndarray

bottom: np.ndarray

ipobo: np.ndarray

boundary_nodes: np.ndarray

element_matid: np.ndarray

node_matid: np.ndarray

boundary_loops: np.ndarray | None = None

roughness: np.ndarray | None = None

quality: object | None = None

property npoin: int

property nelem: int

axqua.solvers.telemac.mesh.roi_polygon(cfg)

Return the ROI as a shapely Polygon (polygonising lines if needed).

Under structures.solid_mode: cut the solid footprints are removed from the domain here, which is the whole of that feature: the geometry code already turns every interior ring into a gmsh curve loop and hands it to addPlaneSurface, so a wall subtracted from the ROI becomes a meshed hole with no change to the mesher at all. Its edges are then mesh edges, the opening beside it is exact, and none of the blanking machinery applies - there is nothing to blank, because the elements were never created.

Parameters:

cfg (Config)

axqua.solvers.telemac.mesh.nominal_channel_size(cfg)

Effective (scaled) target channel edge length [m]. See axqua.core.geodata.nominal_channel_size() - it moved to the core because it reads the shared mesh zones, and both meshers size themselves by it.

Parameters:

cfg (Config)

Return type:

float

axqua.solvers.telemac.mesh.channel_node_mask(cfg, mesh)

Boolean (NPOIN,) flag of mesh nodes on/inside the *channel* mesh-zones.

Uses the same geodata.mesh_zones polygons (Zone Name contains channel) that drive the anisotropic mesh, so the pre-wetting region coincides with the meshed channel. Raises if no channel zones are configured.

The footprint is buffered by ~one cell before the point-in-polygon test so that boundary nodes lying on the channel’s outer edge are included. The inflow/outflow liquid boundaries coincide with that edge, and strict containment would drop those nodes - leaving the prescribed-discharge inflow cross-section dry, which makes TELEMAC’s DEBIMP abort at t=0 (“PROBLEM ON BOUNDARY NUMBER … CHECK THE WATER DEPTHS”). The seeded depth is still max(water_level - bed, 0), so the buffer only wets nodes actually below the hotstart surface; the dry banks stay dry.

Parameters:

• cfg (Config)

• mesh (Mesh)

Return type:

numpy.ndarray

axqua.solvers.telemac.mesh.build_mesh(cfg, dem_initial_roi=None)

Generate the mesh (geometry + MATID). When dem_initial_roi is given, its elevations are interpolated onto the nodes; otherwise the bottom is left at zero and can be filled later with interpolate_elevations().

Parameters:

• cfg (Config)

• dem_initial_roi (Path | None)

Return type:

Mesh

axqua.solvers.telemac.mesh.interpolate_elevations(mesh, dem_path, *, decimals=4, cfg=None)

Interpolate DEM elevations onto the mesh nodes (in place) and return it.

Elevations are rounded to decimals digits after the decimal point (4 by default). NaNs (nodata / just outside the grid) are nearest-neighbour filled.

With cfg, the case’s structures (axqua.core.structures) are burnt into the bed afterwards. A depth-averaged model has no vertical wall to remove, so both structure modes end up as terrain here: an overflow structure raises the bed to its crest, and a solid one raises it to the crest plus structures.solid_freeboard_2d so it cannot be overtopped. That is the standard way a floodwall is represented in TELEMAC, and it is the one place where the two solvers necessarily differ - OpenFOAM removes the footprint and gets a real vertical face, which is reported in the build notes rather than left implicit.

Parameters:

• mesh (Mesh)

• dem_path (Path)

• decimals (int)

• cfg (Config | None)

Return type:

Mesh

axqua.solvers.telemac.mesh.read_roughness_table(path)

Read the zone-roughness CSV -> {zone_id: nominal roughness}. See axqua.core.geodata.read_roughness_table().

Parameters:

path (Path)

Return type:

dict[int, float]

axqua.solvers.telemac.mesh.read_roughness_zones(path)

Read the zone-roughness CSV -> {zone_id: ZoneRoughness} (nominal ks plus the calibration prior bounds and flag). See axqua.core.geodata.read_roughness_zones().

Parameters:

path (Path)

Return type:

dict[int, ZoneRoughness]

axqua.solvers.telemac.mesh.interpolate_roughness(cfg, mesh)

Map the roughness zones onto the mesh: friction zone id + roughness value.

Each node (and element, by its centroid) is tagged with the roughness_zone_field id of the polygon it falls in (nearest polygon for points just outside any). That id becomes the per-node FRIC_ID written to the geometry - overriding the MATID-derived ids so the friction zonation comes from the roughness zones - via mesh.node_matid / mesh.element_matid. The matching geodata.roughness_table value (e.g. a Nikuradse k_s) is stored in mesh.roughness (BOTTOM FRICTION); HydroBayesCal later perturbs it.

Parameters:

• cfg (Config)

• mesh (Mesh)

Return type:

Mesh

axqua.solvers.telemac.mesh.write_mesh(mesh, path, *, title='mesh')

Write a Mesh to a TELEMAC geometry SELAFIN file.

Includes the per-node roughness as BOTTOM FRICTION when it has been set (see interpolate_roughness()).

Parameters:

• mesh (Mesh)

• path (Path)

• title (str)

Return type:

Path

axqua.solvers.telemac.mesh.run(cfg, dem_initial_roi)

Build the mesh and write the geometry SELAFIN; returns (mesh, slf path).

When roughness zones are configured, their Zone ID drives the per-node FRIC_ID and the table ks is written as BOTTOM FRICTION (so the geometry is consistent with the friction .tbl derived from the same table).

Parameters:

• cfg (Config)

• dem_initial_roi (Path)

Return type:

tuple[Mesh, Path]

Mesh-quality assessment and validity checks for the generated triangular mesh.

Reports the shape/size metrics that matter for a stable TELEMAC-2D run while respecting the intended anisotropy of the channel (triangles deliberately elongated along the flow). Ideals follow common hydrodynamic-meshing guidance:

• internal angles ~60 deg, ideally within [30, 90] (avoid obtuse / sliver cells);

• aspect ratio (longest/shortest edge) near 1 — but channel cells are meant to be stretched, so channel and floodplain are reported separately (and the shape warnings are relaxed for the channel) when a channel mask is supplied;

• smooth size transitions: adjacent triangles should not jump in area by >20%;

• the SHORTEST EDGE is surfaced prominently — it drives the CFL-limited adaptive time step.

Validity (fatal) checks: zero-area / degenerate triangles and inverted (mixed orientation) elements raise; duplicate / orphan nodes, non-manifold edges, extra boundary loops (hanging gaps) and IPOBO inconsistency (wrong boundary tags) are flagged. All metrics are vectorized with numpy so they scale to large meshes.

class axqua.solvers.telemac.mesh_quality.RegionQuality(label, relaxed, n_elem, min_angle, max_angle, frac_below_min, frac_obtuse, aspect_median, aspect_p95, aspect_max, skew_median, skew_max)

Bases: object

Shape statistics for one subset of elements (whole mesh, channel, …).

Parameters:

• label (str)

• relaxed (bool)

• n_elem (int)

• min_angle (float)

• max_angle (float)

• frac_below_min (float)

• frac_obtuse (float)

• aspect_median (float)

• aspect_p95 (float)

• aspect_max (float)

• skew_median (float)

• skew_max (float)

label: str

relaxed: bool

n_elem: int

min_angle: float

max_angle: float

frac_below_min: float

frac_obtuse: float

aspect_median: float

aspect_p95: float

aspect_max: float

skew_median: float

skew_max: float

class axqua.solvers.telemac.mesh_quality.QualityReport(n_elem: 'int', n_nodes: 'int', min_edge: 'float', median_edge: 'float', max_edge: 'float', min_area: 'float', min_angle_deg: 'float', max_angle_deg: 'float', area_jump_threshold: 'float', frac_area_jump_over: 'float', max_area_jump_ratio: 'float', regions: 'list[RegionQuality]' = <factory>, n_zero_area: 'int' = 0, n_inverted: 'int' = 0, n_duplicate_nodes: 'int' = 0, n_orphan_nodes: 'int' = 0, n_nonmanifold_edges: 'int' = 0, n_boundary_loops: 'int' = 0, ipobo_consistent: 'bool' = True)

Bases: object

Parameters:

• n_elem (int)

• n_nodes (int)

• min_edge (float)

• median_edge (float)

• max_edge (float)

• min_area (float)

• min_angle_deg (float)

• max_angle_deg (float)

• area_jump_threshold (float)

• frac_area_jump_over (float)

• max_area_jump_ratio (float)

• regions (list[RegionQuality])

• n_zero_area (int)

• n_inverted (int)

• n_duplicate_nodes (int)

• n_orphan_nodes (int)

• n_nonmanifold_edges (int)

• n_boundary_loops (int)

• ipobo_consistent (bool)

n_elem: int

n_nodes: int

min_edge: float

median_edge: float

max_edge: float

min_area: float

min_angle_deg: float

max_angle_deg: float

area_jump_threshold: float

frac_area_jump_over: float

max_area_jump_ratio: float

regions: list[RegionQuality]

n_zero_area: int = 0

n_inverted: int = 0

n_duplicate_nodes: int = 0

n_orphan_nodes: int = 0

n_nonmanifold_edges: int = 0

n_boundary_loops: int = 0

ipobo_consistent: bool = True

property is_fatal: bool

Defects that would make TELEMAC reject or mis-read the geometry.

axqua.solvers.telemac.mesh_quality.assess_quality(mesh, *, channel_mask=None, min_angle_deg=30.0, max_angle_deg=90.0, area_jump_threshold=0.2)

Compute the full QualityReport for mesh.

channel_mask (optional, (NELEM,) bool) flags the deliberately anisotropic channel elements so they are reported in their own, shape-warning-relaxed region; the rest are reported as ‘floodplain’. Without a mask the whole mesh is one region.

Parameters:

• mesh (Mesh)

• min_angle_deg (float)

• max_angle_deg (float)

• area_jump_threshold (float)

Return type:

QualityReport

axqua.solvers.telemac.mesh_quality.log_report(report, logger=<Logger axqua (WARNING)>)

Log the quality report: size summary, per-region shape, and warnings.

Parameters:

• report (QualityReport)

• logger (Logger)

Return type:

None

Boundary conditions

Stage 3 - boundary conditions (.cli).

Classifies each contour node (in IPOBO order) against the user’s liquid-boundary lines and writes the TELEMAC boundary-conditions file. Codes:

• solid wall -> 2 2 2 (every non-liquid outer bound)

• inflow (prescribed flowrate) -> 4 5 5 (depth free, velocity from Q)

• outflow, prescribed elevation -> 5 4 4 (outflow_condition: stage_discharge

  [default] or elevation)

• outflow, free / Neumann -> 4 4 4 (outflow_condition: free)

The liquid-boundary lines come from boundaries.liquid_boundaries (a line layer whose Type (inflow/outflow) field tags each line inflow or outflow; there may be several of each). They must coincide with the outer bounds of the mesh zones so that contour nodes fall on them. Liquid boundaries are numbered exactly the way TELEMAC’s FRONT2 does - per contour loop, starting at the south-westernmost node and walking the contour with the domain on its left, merging a run that wraps the loop start (see _front2_runs()) - so the PRESCRIBED FLOWRATES / PRESCRIBED ELEVATIONS line up with the boundary numbering TELEMAC derives from the geometry.

For numerical stability the total inflow-node count should be within ~10% of the total outflow-node count; otherwise a stability-risk warning is logged (raise the mesh resolution near the boundaries or adjust the line lengths to rebalance).

class axqua.solvers.telemac.boundary.LiquidBoundary(index: 'int', kind: 'str', n_nodes: 'int', discharge: 'float | None' = None)

Bases: object

Parameters:

• index (int)

• kind (str)

• n_nodes (int)

• discharge (float | None)

index: int

kind: str

n_nodes: int

discharge: float | None = None

class axqua.solvers.telemac.boundary.InternalSourceRegion(name, discharge, polygon, area, n_nodes=0, porous_depth=None)

Bases: object

An internal source/sink REGION for a losing-gaining reach.

A 2D depth-averaged model has no subsurface, so hyporheic / underflow exchange where the surface flow loses water in one place and gains it back downstream is represented as TELEMAC source regions (SOURCE REGIONS DATA FILE): a withdrawal (discharge < 0) over the losing polygon and an injection (discharge > 0) over the gaining polygon. TELEMAC spreads each region’s Q uniformly over the mesh nodes inside the polygon (Q/area as a depth rate, prosou.f), so the exchange acts as a gentle distributed flux instead of hammering the few nodes under the line - a sink concentrated on single nodes dries them, spikes velocities and collapses the CFL-adaptive time step (TELEMAC has no depth guard on negative sources).

Parameters:

• name (str)

• discharge (float)

• polygon (object)

• area (float)

• n_nodes (int)

• porous_depth (float | None)

name: str

discharge: float

polygon: object

area: float

n_nodes: int = 0

porous_depth: float | None = None

axqua.solvers.telemac.boundary.dump_liquid_boundaries(liquids, path)

Serialize the FRONT2-ordered liquid boundaries to JSON.

Written during the build (axqua.solvers.telemac.pipeline) so the standalone unsteady / 3D scripts recover the exact TELEMAC boundary numbering (index/kind/node count, and any per-line discharge) without rebuilding the mesh - all load_liquid_boundaries() needs to write the hydrograph liquid-boundaries file and the prescribed-value arrays.

Parameters:

• liquids (list[LiquidBoundary])

• path (str | Path)

Return type:

Path

axqua.solvers.telemac.boundary.load_liquid_boundaries(path)

Load the liquid boundaries dumped by dump_liquid_boundaries().

Parameters:

path (str | Path)

Return type:

list[LiquidBoundary]

axqua.solvers.telemac.boundary.exchange_regions(cfg, mesh=None)

The internal exchange a build applies: the int-* lines, if the case asks.

The lines are a property of the boundary layer, and a layer is shared between cases. They used to be applied whenever the layer held them, so a case that said nothing about an exchange got one anyway, in its least stable form: a fixed-rate source region on a thin strip. On the isar-2025 reach that ends a run at its first time step. An exchange now needs the case to ask for it, with a gain_lose block or with boundaries.internal_sources: true for the fixed-rate regions. Without either, the lines are reported and left out.

Parameters:

• cfg (Config)

• mesh (Mesh | None)

Return type:

list[InternalSourceRegion]

axqua.solvers.telemac.boundary.load_internal_source_regions(cfg, mesh=None)

Internal losing/gaining lines -> TELEMAC source regions (see InternalSourceRegion). A line whose type tag starts with ‘int’ becomes a withdrawal (lose) or injection (gain) of its flow-column discharge, spread over a polygon region. Empty when the layer has no internal lines.

The region polygon is the line buffered to a strip of boundaries.internal_source_region_width - except in percolation.mode: region, where a losing line whose line intersects a percolation patch uses the whole patch polygon instead (same -Q over a far larger area, so the per-node depth rate drops by orders of magnitude).

With a mesh the nodes inside each region are counted; a region containing no node raises (TELEMAC aborts on it at run time).

Parameters:

• cfg (Config)

• mesh (Mesh | None)

Return type:

list[InternalSourceRegion]

axqua.solvers.telemac.boundary.classify_nodes(cfg, mesh)

Classify contour nodes; return per-node kind list and liquid boundaries.

kinds stays in mesh.boundary_nodes order (so it drives the .cli row order, which must match the geometry’s IPOBO). The returned LiquidBoundary list is numbered the way TELEMAC’s FRONT2 does (see _front2_runs()), so the steering file’s per-boundary prescribed values line up with the boundary numbering TELEMAC derives from the geometry.

Parameters:

• cfg (Config)

• mesh (Mesh)

Return type:

tuple[list[str], list[LiquidBoundary]]

axqua.solvers.telemac.boundary.write_cli(cfg, mesh)

Parameters:

• cfg (Config)

• mesh (Mesh)

Return type:

tuple[Path, list[LiquidBoundary]]

Steering and friction files

Stage 4 - steering files: friction .tbl and the TELEMAC .cas (+ GAIA stub).

The friction table lists one row per MATID zone (<id> <LAW> <coef> NULL, terminated by END); node-to-zone linkage is carried by the FRIC_ID variable in the geometry SELAFIN (written in stage 2). HydroBayesCal perturbs the coefficient column of this table during calibration.

The .cas is assembled from a small set of building blocks plus user overrides so the case is reproducible and the prescribed boundary values stay consistent with the liquid-boundary numbering from stage 3.

axqua.solvers.telemac.steering.eddy_viscosity_estimate(cfg)

Depth-averaged turbulent eddy viscosity [m2/s] from the velocity guess.

Elder-type closure nu_t = alpha * u_star * h with the shear velocity u_star ~ 0.1 * U (a skin-friction proxy C_f ~ 0.01) and the transverse-mixing coefficient alpha = 0.6; h is the turbulence length scale (flow depth). Used as the constant VELOCITY DIFFUSIVITY (model 1) and reported for context.

Parameters:

cfg (Config)

Return type:

float

axqua.solvers.telemac.steering.select_turbulence_model(cfg, mesh=None)

Pick the TELEMAC turbulence model from the mesh resolution and velocity guess.

An explicit hydrodynamics.turbulence_model (int or name) is honored as-is. With "auto" the choice follows the channel cell size dx relative to the turbulence length scale L (the flow depth) - a mesh-size/velocity criterion:

• Smagorinski LES (4) when the mesh resolves >=80% of the TKE (dx/L <= 0.2**1.5, ~11+ cells per L): the sub-grid model then handles <20% of the dissipation, the LES validity requirement.

• k-epsilon (3) at moderate resolution (>= ~4 cells per L but below the LES gate): a full two-equation RANS closure resolving the mean shear.

• Spalart-Allmaras (6) on coarse meshes (< ~4 cells per L): a robust, economical one-equation RANS model, tolerant of the stretched channel cells.

The velocity guess sets the eddy-viscosity scale (and the steady initial flow field); for a fixed L it cancels from the resolution ratio, so the cell size drives the choice. Returns (model_number, rationale).

Parameters:

cfg (Config)

Return type:

tuple[int, str]

axqua.solvers.telemac.steering.turbulence_pick_for_dx(cfg, dx)

The mesh-resolution auto-selection of select_turbulence_model() for an arbitrary channel cell size dx [m], ignoring any explicit hydrodynamics.turbulence_model override.

Used by the mesh-convergence validity check to ask “what would auto-selection pick at this refinement?” while the study runs with the closure pinned to the baseline choice. Returns (model_number, rationale).

Parameters:

• cfg (Config)

• dx (float)

Return type:

tuple[int, str]

axqua.solvers.telemac.steering.write_friction_tbl(cfg)

Write the FRICTION DATA FILE (.tbl), one row per friction zone.

Zones come from friction.zones (MATID) or, failing that, from the roughness zones (Zone ID -> ks via roughness_table) so the rows match the geometry’s per-node FRIC_ID.

Parameters:

cfg (Config)

Return type:

Path

axqua.solvers.telemac.steering.spill_elevations(mesh, seed_mask)

Per-node spill elevation relative to the seed nodes: the lowest water level at which the node is hydraulically connected to them.

This is Barnes et al.’s priority flood run on the mesh graph: starting from the seed nodes, the elevation of every node is raised to the highest bed it must cross on the cheapest path to a seed, i.e.

\[S_i = \min_{\text{paths } i \to \text{seed}} \; \max_{j \in \text{path}} z_j\]

So S_i == z_i on ground that drains freely to a seed node, and S_i > z_i behind a rim, where S_i is the elevation of the lowest saddle that must be overtopped. Seeded from the outflow it answers “can water here get out?”; seeded from the inflow, “could the flow ever have got here?” - the pair defines the through-flowing corridor used by write_initial_conditions().

axqua.solvers.telemac.steering.PREWET_SLOPE_LIMITS = (0.002, 0.03)

bounds on the per-transect bed slope used for the normal-depth seed [-]

axqua.solvers.telemac.steering.write_initial_conditions(cfg, mesh)

Write the hotstart SELAFIN that pre-wets the channel for a stable start.

Seeds WATER DEPTH on the *channel* mesh-zone nodes (zero velocity, dry elsewhere) on the geometry’s own mesh; TELEMAC reads it as the PREVIOUS COMPUTATION FILE (see write_cas()) so the run continues from a pre-wetted channel instead of advancing the wetting front from a dry bed.

The seeded water surface comes from initialization.prewet_mode:

• normal-depth (default): the normal-flow stage of the real cross-sections at the case discharge, scaled by prewet_fill (_normal_depth_prewet_depth()). The seeded surface then tracks the surface the run converges to, so no water is laid down above it.

• constant: the older longitudinally smoothed “low bed percentile + prewet_depth” surface (_longitudinal_prewet_depth()). Also the automatic fallback when there is no centerline or no discharge to size the normal depth with - a plain constant depth would make the seeded surface as jagged as the bed and leave a depth “dam” at every channel/floodplain edge, which diverges at t=0 on steep terrain.

Whatever the mode, nodes that would receive less than initialization.prewet_min_depth are left dry: a feathered seed margin carries no flow worth hotstarting and is exactly what turns into stagnant film.

Drainable seeding (initialization.drainable_prewet, on by default): the depth is measured from the node’s spill elevation (spill_elevations()) rather than from the bed, so no water is ever seeded below the rim of a closed depression. Seeding a bowl the flow would never fill leaves a pond that cannot drain - TELEMAC-2D has no infiltration or evaporation - and it survives to the end of the run as a stagnant, wrongly wetted alcove or side channel. Referencing the spill elevation removes exactly that water and nothing else: on ground that drains freely the spill elevation is the bed, so the open channel is seeded unchanged.

Parameters:

cfg (Config)

Return type:

Path

axqua.solvers.telemac.steering.write_dry_start_conditions(cfg, mesh)

Write the DRY-START initial-conditions SELAFIN: a thin water plug only on the nodes near the inflow line(s), dry everywhere else.

A fully dry bed (ZERO DEPTH) makes TELEMAC’s DEBIMP abort at a prescribed-Q inflow (“PROBLEM ON BOUNDARY … CHECK THE WATER DEPTHS”): with no water at the boundary the discharge cannot be distributed. So instead of pre-wetting the whole channel, we wet only a strip within dry_start_extent of the inflow line(s) to dry_start_depth (deep enough to keep the inflow subcritical) and leave the rest of the domain dry. TELEMAC continues from this file and the flow wets the reach from the inflow. Returns the path, or None when there is no inflow line (the caller then uses the analytical dry initial condition).

Parameters:

cfg (Config)

Return type:

Path | None

axqua.solvers.telemac.steering.write_liquid_boundaries(cfg, liquids, inflow, outflow_wse_fn=None)

Write the TELEMAC liquid-boundaries (hydrograph) file for the unsteady run.

One column per liquid boundary, in boundary-index order: Q(i) for inflows (the total reach discharge at each time split across the inflow boundaries by node share, exactly as the steady prescribed flowrates) and SL(i) for prescribed-elevation outflows (outflow_wse_fn(Q) - the rating-curve stage at that discharge). A free (Neumann) outflow prescribes nothing and gets no column. The first column is the time T in seconds.

Parameters:

• cfg (Config)

• liquids (list[LiquidBoundary])

Return type:

Path

axqua.solvers.telemac.steering.write_source_regions(cfg, regions)

Write the TELEMAC SOURCE REGIONS DATA FILE (source-regions.txt).

Format (read_source_data.f): # comment lines; per region a header line X(i)   Y(i) followed by one x y vertex pair per line; a # comment line terminates each region block. The region order defines TELEMAC’s region numbering and must match the order of WATER DISCHARGE OF SOURCES in the .cas (both come from the same regions list here). TELEMAC assigns every mesh node inside a polygon to its region and spreads the region discharge uniformly over the enclosed area (telemac2d_init.F / prosou.f).

No blank lines, and no line may start with a space or tab. read_source_data.f skips leading whitespace with a GO TO 2 that jumps back onto 2 CONTINUE - which re-runs IDEB=1. Any whitespace-only line reaching that block therefore loops forever: the solver spins at 100% CPU right after RESCUE : SPALART ALLMARAS, with no error, no time step and no end (only MAXIMUM NUMBER OF TIME STEPS, unused here, would ever stop it). Hence the # block terminators rather than blank separator lines.

Parameters:

• cfg (Config)

• regions (list)

Return type:

Path

axqua.solvers.telemac.steering.write_cas(cfg, liquids, inflow_q, outflow_wse=None, gaia_cas=None, previous_computation=None, turbulence_model=None, unsteady=False, liquid_boundaries_file=None, duration=None, out_name=None, results_name=None, n_distributive_corrections=None, sections_input=None, sections_output=None, hotstart_note=None, source_regions=None)

Write a TELEMAC-2D steering (.cas) file.

The default (unsteady=False, out_name=None) writes the steady initial run to cfg.cas_file (steady2d.cas). With unsteady=True and a liquid_boundaries_file it writes a hydrograph-driven run (Q(t)/SL(t) read from that file, total time from duration) - the pipeline emits it as unsteady2d.cas when the inflow carries a varying series, and axqua.solvers.telemac.unsteady writes the recommended one hotstarted from the steady result.

turbulence_model is the resolved TELEMAC model number (see select_turbulence_model()); when None it is selected from the config here.

When previous_computation is given the case is continued from that SELAFIN - the flow field is read from it instead of the analytical INITIAL CONDITIONS, with the clock reset to zero. This is either the pre-wetted/dry-start hotstart from write_initial_conditions() (the steady run) or the converged steady result the unsteady run continues from; hotstart_note overrides the comment line.

results_name overrides the RESULTS FILE (so an unsteady run continued from the steady r2d.slf writes to a distinct file rather than clobbering its own input). n_distributive_corrections emits NUMBER OF CORRECTIONS OF DISTRIBUTIVE SCHEMES (the developers’ >=2 recommendation for quasi-steady runs). sections_input / sections_output wire the CONTROL SECTIONS keywords so the run reports the flux across each open boundary (see axqua.unsteady()). source_regions (a axqua.boundary.InternalSourceRegion list) adds the internal losing/gaining exchange - as SOURCE REGIONS keywords referencing the file written by write_source_regions(), or, in percolation.mode: fortran, as the FORTRAN FILE + RAIN keywords for the generated USER_RAIN routine (see axqua.solvers.telemac.fortran).

Parameters:

• cfg (Config)

• liquids (list[LiquidBoundary])

• inflow_q (float)

• outflow_wse (float | None)

• gaia_cas (str | None)

• previous_computation (str | None)

• turbulence_model (int | None)

• unsteady (bool)

• liquid_boundaries_file (str | None)

• duration (float | None)

• out_name (str | None)

• results_name (str | None)

• n_distributive_corrections (int | None)

• sections_input (str | None)

• sections_output (str | None)

• hotstart_note (str | None)

• source_regions (list | None)

Return type:

Path

axqua.solvers.telemac.steering.write_hotstart_cas(cfg, duration, out_name='hotstart2d.cas')

Derive a hotstart steering file from the built steady case.

Rewrites cfg.cas_file (steady2d.cas) as a continuation run: the initial conditions block is switched from the pre-wetted/dry-start seed to PREVIOUS COMPUTATION FILE : <results_slf> (the steady run’s own result), the results file is renamed so the run does not clobber its hotstart input, and DURATION is capped at duration - the simulated time the initial run needed to reach a sustained boundary-flux balance (flux_convergence.find_steady_window). Everything else - notably the constant PRESCRIBED FLOWRATES / PRESCRIBED ELEVATIONS - is carried over unchanged, so the hotstart drives the same steady Q and downstream H as the initial run.

Parameters:

• cfg (Config)

• duration (float)

• out_name (str)

Return type:

Path

axqua.solvers.telemac.steering.write_gaia_cas(cfg)

Write the GAIA steering file (bedload and/or suspended load + bed processes).

Emits a structurally valid GAIA .cas - referencing the same GEOMETRY FILE and BOUNDARY CONDITIONS FILE as the coupled hydrodynamic run (GAIA needs both even in coupled mode) - then enables the transport modes and bed-process capacities set on cfg.morphodynamics:

• bedload -> BED LOAD FOR ALL SANDS : YES + BED-LOAD TRANSPORT FORMULA FOR ALL SANDS (bedload_formula, 1 = MPM);

• suspended_load -> SUSPENSION FOR ALL SANDS : YES (GAIA transports the classes as TELEMAC tracers through the coupling);

• MORPHOLOGICAL FACTOR accelerates bed evolution over a hydrograph;

• SLOPE EFFECT (+ formula + friction angle) steers bedload down transverse bed slopes; SECONDARY CURRENTS (+ alpha) adds the spiral-flow deviation in bends; NUMBER OF LAYERS FOR INITIAL STRATIFICATION / ACTIVE LAYER THICKNESS set the bed stratigraphy;

• PRESCRIBED SOLID DISCHARGES imposes a sediment supply per liquid boundary (else GAIA uses its equilibrium inflow).

The same file is shared by the steady case and by both hydrograph runs (unsteady2d.cas / unsteady3d.cas) - see axqua.config.Morphodynamics. Sediment-transport tuning is otherwise left to the user / calibration (e.g. the Shields parameters or class diameters perturbed by HydroBayesCal via gaia*).

Parameters:

cfg (Config)

Return type:

Path | None

Evaluation of runs

Reader for the TELEMAC .sortie listing (the solver’s printed output).

A steady TELEMAC-2D/3D run prints a water-volume balance block at every listing printout:

                  BALANCE OF WATER VOLUME
VOLUME IN THE DOMAIN :    6317.532     M3
FLUX BOUNDARY    1:    0.8000116     M3/S  ( >0 : ENTERING  <0 : EXITING )
FLUX BOUNDARY    2:   -7.3518840     M3/S  ( >0 : ENTERING  <0 : EXITING )
RELATIVE ERROR IN VOLUME AT T =        28.57     S :   -0.1008281E-14

which is everything needed to judge whether a steady run has converged: at steady state the boundary fluxes balance (total inflow = total outflow) and the domain volume stops changing. This module turns that listing into arrays.

TELEMAC-3D prints a different block (telemac3d/bil3d.f), and it has to be read too - a 3D pre-run seeds OpenFOAM exactly as a 2D one does, and until this was handled its flux balance could not be judged at all:

                MASS BALANCE

  WATER
VOLUME AT THE PREVIOUS TIME STEP              :     105.4575
VOLUME AT THE PRESENT TIME STEP               :     105.4497
VOLUME LEAVING DOMAIN DURING THIS TIME STEP   :    0.7793448E-02
ERROR ON THE VOLUME DURING THIS TIME STEP     :   -0.1874828E-12
  BOUNDARY FLUXES FOR WATER IN M3/S ( >0 : ENTERING )
FLUX BOUNDARY      1                          :    0.1350000
FLUX BOUNDARY      2                          :   -0.4133374

Three things differ and each is handled explicitly. (1) There is no VOLUME IN THE DOMAIN line; the domain volume is VOLUME AT THE PRESENT TIME STEP. (2) The block carries no time. TELEMAC-3D prints it only in the ITERATION ... TIME ... ( N S) header that immediately precedes the block, so that header’s time is attached to the block that follows it - still a single-pass association made while reading, not two sequences zipped afterwards, so it cannot slip. (3) The error line is the absolute volume error of one time step, where TELEMAC-2D prints a relative one; Sortie.dimension says which listing a result came from, and volume_error must be read accordingly. The fluxes - the only thing convergence is actually judged on - mean the same in both.

Why axqua parses this itself. The obvious alternative is TELEMAC’s own postel.parser_output (and the pythomac package that adapts it), but both are GPL-3 while axqua is BSD-3, so neither can be vendored here. The file format is not the licensed part, so this is an independent reader written against the listing itself - which also lets it be stricter and simpler than a general-purpose one, because axqua needs exactly one thing from the listing:

• Each balance block is parsed as a unit, and its simulation time is taken from the RELATIVE ERROR IN VOLUME AT T = ... line inside the block. The general parser instead collects times from the separate ITERATION ... TIME: lines and zips the two sequences afterwards, which needs a synthetic zero-flux sample prepended to make the lengths match - so the first sample is fabricated and every later sample can slip by one if the listing has an extra ITERATION line (PRINTING CUMULATED FLOWRATES and the listing printout period do not have to agree). Reading the time from inside the block makes misalignment impossible.

• Nothing is fabricated: what the run printed is what is returned.

• Blocks are only accepted when complete (volume + every boundary flux + the time), so a listing truncated mid-block by a crash or a still-running job yields the blocks written so far instead of raising.

Parallel runs write one main listing plus a *_p0000N.sortie per processor; only the main one carries the merged balance. latest_sortie() returns it.

class axqua.solvers.telemac.sortie.Sortie(path, study, exec_seconds, iteration, time, volume, fluxes, volume_error, dimension='2d')

Bases: object

The water-volume balance history of one TELEMAC listing.

All arrays are one entry per listing printout and share the same ordering, so time[i], volume[i] and fluxes[i] describe the same instant.

Parameters:

• path (Path)

• study (str)

• exec_seconds (float)

• iteration (numpy.ndarray)

• time (numpy.ndarray)

• volume (numpy.ndarray)

• fluxes (numpy.ndarray)

• volume_error (numpy.ndarray)

• dimension (str)

path: Path

study: str

exec_seconds: float

iteration: numpy.ndarray

time: numpy.ndarray

volume: numpy.ndarray

fluxes: numpy.ndarray

volume_error: numpy.ndarray

dimension: str = '2d'

"2d" (TELEMAC-2D) or "3d" (TELEMAC-3D). It matters for volume_error alone, which TELEMAC-2D prints relative and TELEMAC-3D prints as an absolute per-step volume; the value is passed through as printed rather than converted, so the two are not silently made to look alike. Everything else means the same in both.

Type:

which listing shape this came from

property n_boundaries: int

property gross_in: numpy.ndarray

the sum of the positive boundary fluxes.

Summing by sign rather than by boundary index means any number of inflow / outflow boundaries works, and a boundary that momentarily reverses is counted on the side it is actually flowing.

Type:

Total inflow [m3/s]

property gross_out: numpy.ndarray

the sum of the negative fluxes, positive.

Type:

Total outflow magnitude [m3/s]

property printout_interval: float

Median simulated time between listing printouts [s].

Taken from the times actually printed, because a CFL-adaptive run has no fixed time step: LISTING PRINTOUT PERIOD x TIME STEP is meaningless there.

to_frame()

The balance history as a pandas DataFrame indexed by simulated time.

axqua.solvers.telemac.sortie.read_sortie(path)

Parse the water-volume balance history out of a .sortie listing.

Reads both the TELEMAC-2D and the TELEMAC-3D balance block (see the module docstring); which one a listing carried is reported as Sortie.dimension.

Raises ValueError when the listing holds no complete balance block - normally because the run was launched without -s or without PRINTING CUMULATED FLOWRATES : YES.

Parameters:

path (str | Path)

Return type:

Sortie

axqua.solvers.telemac.sortie.latest_sortie(model_dir, cas_name)

Newest main listing for cas_name in model_dir.

TELEMAC names it <cas>_<YYYY-MM-DD-HHhMMminSSs>.sortie; a parallel run adds a _p0000N copy per processor, which carries only that subdomain’s balance and is skipped here.

Parameters:

• model_dir (str | Path)

• cas_name (str)

Return type:

Path | None

axqua.solvers.telemac.sortie.processor_sorties(model_dir, cas_name)

The per-processor *_p0000N.sortie copies of a parallel run.

Parameters:

• model_dir (str | Path)

• cas_name (str)

Return type:

list[Path]

axqua.solvers.telemac.sortie.LIMIT_STOP = 'LIMIT VALUES TRESPASSED'

What TELEMAC-2D prints when CONTROL OF LIMITS ends a run.

telemac2d.F calls this a “neat (programmed) stop”: the solver returns from its time loop, prints CORRECT END OF RUN and exits with code 0. The exit code can therefore not tell a run that blew up at the first time step from one that finished, and a job that trusted it reported the first as completed.

class axqua.solvers.telemac.sortie.LimitStop(bound, variable, x, y, value, limit)

Bases: object

One node at which a value left the range LIMIT VALUES allows.

Parameters:

• bound (str)

• variable (str)

• x (float)

• y (float)

• value (float)

• limit (float)

bound: str

variable: str

x: float

y: float

value: float

limit: float

describe()

Return type:

str

axqua.solvers.telemac.sortie.read_limit_stops(text)

The nodes a listing names as having left their limits, in listing order.

The node number is not kept: in a parallel run it counts within one subdomain and means nothing on the mesh the user sees. The coordinates do.

Parameters:

text (str)

Return type:

list[LimitStop]

axqua.solvers.telemac.sortie.find_limit_stops(model_dir, cas_name)

Where the newest run of cas_name left its limits; empty when it did not.

In a parallel run only the subdomain that owns the node reports it, into its own _p0000N listing, so those are read as well as the main one.

Parameters:

• model_dir (str | Path)

• cas_name (str)

Return type:

list[LimitStop]

axqua.solvers.telemac.sortie.sediment_mass_profile(path)

GAIA sediment mass balance per class over the run.

Returns {class_number: {quantity: array}} with one entry per listing printout; class 0 holds the “OVER ALL CLASSES” summary block. Quantities are the keys of _SEDIMENT_FIELDS (total/lost/cumulated mass, erosion and deposition fluxes, bed evolutions, the relative mass errors, …). An empty dict means the run was not coupled to GAIA.

The mass balance is the morphodynamic counterpart of the water-volume balance: relative_error_active_layer and cumulated_lost_mass are what tell you a GAIA run stayed mass-conservative.

Parameters:

path (str | Path)

Return type:

dict[int, dict[str, numpy.ndarray]]

axqua.solvers.telemac.sortie.tracer_mass_profile(path)

Tracer mass balance per tracer over the run.

Returns {tracer_number: {"total": array, "error": array}}, one entry per listing printout. GAIA carries suspended sediment as TELEMAC tracers, so this is how a suspended-load run is mass-checked. The listing’s WATER balance block uses the same wording and is skipped.

Parameters:

path (str | Path)

Return type:

dict[int, dict[str, numpy.ndarray]]

axqua.solvers.telemac.sortie.find_lines(path, pattern, *, with_line_numbers=False)

Every listing line matching pattern (a regular expression).

The escape hatch for anything this module does not model - a user-Fortran printout, a solver warning, a keyword echo. USER_RAIN percolation diagnostics in the isar-2025 case are read this way, for instance.

Parameters:

• path (str | Path)

• pattern (str)

• with_line_numbers (bool)

Temporal (flux / mass-balance) convergence analysis of the steady test run.

The initial_run.py test run serves two purposes: confirm the built case does not crash, and – because that same steady result is the hotstart seed for the HydroBayesCal calibration – find the simulation time at which the model has reached a steady, mass-conservative state. Calibrating against a not-yet-converged hotstart would bake a transient bias into every perturbed run.

Convergence criterion (after mass-balance-bc.md / convergence.md in the hyhome numerics notes): a steady Telemac2d run is converged when the liquid boundaries are in mass balance, i.e. total inflow equals total outflow. We track the dimensionless relative flux imbalance per listing printout:

eps_t = | ``|Q_in,t|`` - ``|Q_out,t|`` | / ``|Q_in,t|``

(magnitudes because Telemac reports inflow positive, outflow negative; normalizing by the inflow makes it dimensionless), and take the run as converged at the first printout beyond which eps_t stays permanently below a target tolerance. For a steady run this mass-flux balance is the governing criterion: once it holds, the depth and velocity fields have stopped evolving as well. (Pointwise depth / velocity grid-convergence at the measurement probes is a different question, handled by the mesh-convergence study.)

The tolerance depends on what the result is for, and the imbalance being a relative measure is what makes that judgment possible: the discharge, depth and velocity read off a steady run inherit it directly. 1e-3 (0.1%) is therefore the default – a tenth of a percent moves Q, h or U by far less than the uncertainty of the field data they are compared against, and by far less than the 5% the mesh-convergence study accepts as grid-independent. The stricter 1e-4 is reserved for a hotstart seed, where a residual transient is not averaged out but inherited by every one of the dozens of perturbed calibration runs restarted from it; 1e-6 is a validation grade whose slow asymptote is rarely worth the simulated time. Set it per case with hydrodynamics.flux_tolerance, or pass tolerance= to analyze_flux_convergence().

This module is self-contained: the listing is parsed by axqua.solvers.telemac.sortie and the convergence maths and plots are below. It previously delegated to the pythomac package, whose analysis this reproduces (and whose four output files it still writes, with the same names, so existing workflows and documentation are unaffected). pythomac is GPL-3 and derives its listing parser from TELEMAC’s own GPL-3 postel code, which cannot be vendored into BSD-3-licensed axqua; and having the maths here removes a set of adapters that existed only to work around the package boundary – the printout spacing had to be back-calculated because a CFL-adaptive run has no fixed time step, the rate plot had to be replaced whenever the rate was still undefined, and convergence-rate.csv had to be written separately because only the .png was produced.

class axqua.solvers.telemac.flux_convergence.FluxConvergence(converged, tolerance, final_imbalance, converged_time_steps, converged_seconds, listing_printout_period, time_step, fluxes_csv, flux_plot, rate_csv, rate_plot, steady_seconds=None, steady_abs_tolerance=0.001, steady_window=10, steady_strict=True, hotstart_cas=None)

Bases: object

Outcome of the flux/mass-balance convergence analysis.

Parameters:

• converged (bool)

• tolerance (float)

• final_imbalance (float)

• converged_time_steps (int | None)

• converged_seconds (float | None)

• listing_printout_period (int)

• time_step (float)

• fluxes_csv (Path | None)

• flux_plot (Path | None)

• rate_csv (Path | None)

• rate_plot (Path | None)

• steady_seconds (float | None)

• steady_abs_tolerance (float)

• steady_window (int)

• steady_strict (bool)

• hotstart_cas (Path | None)

converged: bool

tolerance: float

final_imbalance: float

converged_time_steps: int | None

converged_seconds: float | None

listing_printout_period: int

time_step: float

fluxes_csv: Path | None

flux_plot: Path | None

rate_csv: Path | None

rate_plot: Path | None

steady_seconds: float | None = None

steady_abs_tolerance: float = 0.001

steady_window: int = 10

steady_strict: bool = True

hotstart_cas: Path | None = None

axqua.solvers.telemac.flux_convergence.relative_imbalance(gross_in, gross_out)

Relative flux imbalance eps_t = ||Q_in| - |Q_out|| / |Q_in| per printout.

Telemac reports boundary fluxes signed (inflow positive, outflow negative), so the magnitudes are differenced: balance means |Q_in| == |Q_out|, hence eps_t -> 0 at steady state. Using magnitudes also keeps the metric robust to a boundary that momentarily reverses.

Return type:

numpy.ndarray

axqua.solvers.telemac.flux_convergence.convergence_rate(imbalance)

Convergence rate iota_t = log_{eps_t}(eps_{t+1}), aligned with eps[1:].

A rate of 1 means the imbalance is shrinking at a constant order; above 1 it is accelerating. It is undefined wherever the imbalance is flat or not yet decaying (log of a base near 1, or of a non-positive value), which is the normal state during the filling/transient phase - those entries come back as NaN rather than raising, and the plots below simply show gaps there.

Return type:

numpy.ndarray

axqua.solvers.telemac.flux_convergence.convergence_index(imbalance, tolerance=0.0001)

First index beyond which imbalance stays permanently below tolerance.

“Permanently” matters: a transient dip below the tolerance during the filling phase is not convergence. Returns None when the tolerance is never reached, or is reached but later lost again.

Parameters:

tolerance (float)

Return type:

int | None

axqua.solvers.telemac.flux_convergence.find_steady_window(times, gross_in, gross_out, *, abs_tolerance=0.001, window=10)

Simulated time [s] at which the boundary fluxes reach a sustained balance.

Returns the time of the first listing printout whose absolute flux imbalance ||Q_in| - |Q_out|| is below abs_tolerance (m3/s) and stays below it for *window* consecutive printouts (so a momentary crossing of the oscillating imbalance does not count), or None if no such window exists. This time is the end time (DURATION) for the generated hotstart case: continuing from the steady result, the flow re-equilibrates well within the span the initial run needed to balance its fluxes.

Parameters:

• abs_tolerance (float)

• window (int)

Return type:

float | None

axqua.solvers.telemac.flux_convergence.plot_boundary_fluxes(sortie, path)

flux-convergence.png: the signed flux at every liquid boundary over time, with the gross inflow/outflow totals the convergence metric is built from.

Parameters:

• sortie (Sortie)

• path (Path)

Return type:

Path

axqua.solvers.telemac.flux_convergence.plot_convergence(times, imbalance, rate, path, *, tolerance=0.0001)

convergence-rate.png: the relative flux imbalance (log scale, with the tolerance marked) and, where it is defined, the convergence rate beneath it.

Drawing the imbalance unconditionally is the point: during the filling phase the rate is undefined everywhere, and a rate-only figure would then be empty exactly when one most wants to see how far from balance the run still is.

Parameters:

• path (Path)

• tolerance (float)

Return type:

Path

axqua.solvers.telemac.flux_convergence.analyze_flux_convergence(cfg, *, tolerance=None, abs_tolerance=None, steady_window=10, write_hotstart=True, pythomac_dir=None)

Analyze boundary-flux convergence of the steady run in cfg.model_dir.

Reads the latest <cas>_<timestamp>.sortie listing (so the run must have used the -s flag and PRINTING CUMULATED FLOWRATES : YES) and writes four files into model_dir - extracted-fluxes.csv + flux-convergence.png (the per-boundary fluxes) and convergence-rate.csv + convergence-rate.png (the relative imbalance and convergence rate). Returns the converged NUMBER OF TIME STEPS (the time at which the relative flux imbalance drops permanently below tolerance).

On top of the relative-imbalance criterion, the absolute steady window is detected (find_steady_window: ||Q_in|-|Q_out|| < abs_tolerance m3/s over steady_window consecutive printouts); when found and write_hotstart is on, a hotstart2d.cas continuing from the steady results file with that time as DURATION is written next to the steady case (steering.write_hotstart_cas). The per-processor *_p0000N.sortie copies of a parallel run are deleted after a successful extraction (only the merged main listing is analyzed and kept).

pythomac_dir is accepted and ignored; the analysis no longer shells out to that package (see the module docstring).

Parameters:

• cfg (Config)

• tolerance (float | None)

• abs_tolerance (float | None)

• steady_window (int)

• write_hotstart (bool)

• pythomac_dir (str | Path | None)

Return type:

FluxConvergence

Wetted-extent diagnostics for a converged 2D result.

A steady run can balance its boundary fluxes perfectly and still be wrong about where the water is. The two failure modes this module measures are:

Stagnant film. A 2D model has neither infiltration nor evaporation, so any water that ends up somewhere the flow cannot carry it away from stays there for the rest of the run. That happens when the initial condition seeds a water surface above the one the run converges to: the excess drains only where it can flow, and on flat ground over coarse gravel it stops as a one-to-five-centimetre immobile film. The film does not shrink with a longer run - on isar-2025 it fell by 1.8 % over the last 3100 s - so it has to be measured, not waited out. wetting_report() splits the wetted area into what is active (deep and moving), what is film (wet but at a standstill) and what is an isolated puddle (wet but not connected to the main water body), and reports how much of each the initial condition put there.

A mis-set outflow stage. A prescribed-elevation outlet that sits above the reach’s own level backs water up over ground that should be dry; one that sits below pulls the flow into a drawdown and can drive it supercritical. Either shows up as an anomaly in the free-surface slope over the last few meters before the boundary, which outlet_profile() measures against the slope of the reach above it.

Both write a CSV next to the result so successive runs can be compared directly.

axqua.solvers.telemac.wetting.WET_DEPTH = 0.01

below this depth [m] a node counts as dry

axqua.solvers.telemac.wetting.ACTIVE_DEPTH = 0.05

an active node is at least this deep [m] …

axqua.solvers.telemac.wetting.ACTIVE_VELOCITY = 0.15

… and moving at least this fast [m/s]

axqua.solvers.telemac.wetting.FILM_VELOCITY = 0.05

a film node is wet but slower than this [m/s]

axqua.solvers.telemac.wetting.OUTLET_BANDS = ((0.0, 3.0), (3.0, 6.0), (6.0, 10.0), (10.0, 20.0), (20.0, 40.0), (40.0, 70.0), (70.0, 120.0))

distance bands from the outflow line [m] used by outlet_profile()

axqua.solvers.telemac.wetting.OUTLET_SLOPE_FACTOR = 1.5

a near-boundary surface slope this many times the reach slope is an anomaly

class axqua.solvers.telemac.wetting.WettingReport(time, wet_area, wet_volume, active_area, active_volume, film_area, film_volume, isolated_area, isolated_volume, isolated_count, supported_area=0.0, supported_volume=0.0, seeded_area=None, seeded_volume=None, film_seeded_area=None, film_history=<factory>)

Bases: object

Where the water sits at the end of a run (see wetting_report()).

Parameters:

• time (float)

• wet_area (float)

• wet_volume (float)

• active_area (float)

• active_volume (float)

• film_area (float)

• film_volume (float)

• isolated_area (float)

• isolated_volume (float)

• isolated_count (int)

• supported_area (float)

• supported_volume (float)

• seeded_area (float | None)

• seeded_volume (float | None)

• film_seeded_area (float | None)

• film_history (list[tuple[float, float]])

time: float

wet_area: float

wet_volume: float

active_area: float

active_volume: float

film_area: float

film_volume: float

isolated_area: float

isolated_volume: float

isolated_count: int

supported_area: float = 0.0

water an external source holds in place (a groundwater-fed pool on a bar): legitimately wet, so it is reported on its own and kept out of film/puddles

supported_volume: float = 0.0

seeded_area: float | None = None

seeded_volume: float | None = None

film_seeded_area: float | None = None

film_history: list[tuple[float, float]]

(time, film area) of the sampled frames - a flat tail means it will not drain

property film_fraction: float

Film as a share of the wetted area.

property film_seeded_fraction: float | None

property film_trend: float | None

Relative change of the film area over the sampled history (negative = still draining). None when fewer than two frames were sampled.

summary()

User-facing report lines.

Return type:

list[str]

axqua.solvers.telemac.wetting.wetting_report(results, *, geometry=None, initial_conditions=None, supported=None, wet_depth=0.01, active_depth=0.05, active_velocity=0.15, film_velocity=0.05, history_frames=6, out=None)

Split the wetted area of a result into active flow, film and puddles.

results is the run’s SELAFIN, geometry the build’s double-precision geometry (its coordinates are used, as in axqua.solvers.telemac.sections), and initial_conditions the pre-wet SELAFIN, which lets the report attribute the film to the seed rather than to the flow. history_frames earlier frames are re-read to show whether the film is still shrinking; pass 0 to skip that (it is the expensive part on a many-frame result).

supported is an optional boolean node mask (or a depth array) marking water an external source holds in place - a pool fed by the water table under a gravel bar (axqua.solvers.telemac.watertable). Such water is standing and disconnected by construction, so without the mask it would be counted as film and as an isolated puddle, i.e. reported as a defect when it is exactly what the model intends.

Writes wetting-report.csv to out (a directory or an explicit file path) when given.

Parameters:

• wet_depth (float)

• active_depth (float)

• active_velocity (float)

• film_velocity (float)

• history_frames (int)

• out (Path | None)

Return type:

WettingReport

class axqua.solvers.telemac.wetting.OutletProfile(bands, reach_slope, outlet_slope, verdict)

Bases: object

Free-surface profile approaching the outflow boundary.

Parameters:

• bands (list[dict])

• reach_slope (float)

• outlet_slope (float)

• verdict (str)

bands: list[dict]

reach_slope: float

outlet_slope: float

verdict: str

summary()

Return type:

list[str]

axqua.solvers.telemac.wetting.outlet_profile(cfg, results, *, geometry=None, bands=((0.0, 3.0), (3.0, 6.0), (6.0, 10.0), (10.0, 20.0), (20.0, 40.0), (40.0, 70.0), (70.0, 120.0)), out=None)

Measure the free-surface profile approaching the outflow boundary.

Bins the actively flowing nodes by their distance from the outflow liquid boundary and reports, per band, the discharge-weighted free surface plus the mean depth, speed and Froude number. Comparing the surface slope over the last few meters with the slope of the reach above it distinguishes

• backwater - the boundary holds the level above the reach’s own, the surface flattens or rises into it, and the extra depth floods ground near the outlet;

• drawdown - the boundary sits at or below the reach’s level and the surface steepens into it (mild is normal, especially where the section narrows);

• neutral - the boundary continues the reach profile.

Writes outlet-profile.csv to out when given.

Parameters:

out (Path | None)

Return type:

OutletProfile

Discharge across internal cross-section (“baffle”) lines.

TELEMAC’s own CONTROL SECTIONS keyword reports the flux across sections defined by node pairs in the steering file, which means the section geometry has to be pinned to mesh node numbers before the run. This module does the same job afterwards, from the result file, so the sections can be drawn freely in GIS as a line layer and changed without re-running the solver.

For each line the unit-discharge vector q = (H*U, H*V) is interpolated from the mesh onto dense sample points along the line and integrated against the line normal:

\[Q = \int_L (q \cdot n) \, ds\]

The normal is the right-hand normal of the digitized line; because the digitizing direction of a GIS layer is arbitrary, the sign is normalized so that a section reports the discharge positive in the direction of its own net flow (the sign convention is reported in orientation). A braided reach can therefore be split into its threads and each thread’s share read off directly.

Accuracy notes

• Interpolation is linear within the triangle containing each sample point (matplotlib.tri.LinearTriInterpolator), i.e. the same P1 basis TELEMAC itself uses for the depth and velocity fields - so the integral is consistent with the discretization rather than a nearest-node approximation.

• Sample spacing defaults to a quarter of the local mesh edge length, which keeps the trapezoidal integration error well below the solver’s own mass imbalance.

• Coordinates come from the geometry SELAFIN (written in double precision) when available: the result file is single precision, and on UTM easting ~677 000 m the float32 spacing is ~0.06 m - a fifth of the 0.3 m channel cell size, which would visibly distort the section geometry.

axqua.solvers.telemac.sections.SAMPLE_FRACTION = 0.25

sample spacing along a section, as a fraction of the mean mesh edge length

axqua.solvers.telemac.sections.MAX_SAMPLE_SPACING = 0.5

never sample coarser than this (m), whatever the mesh

axqua.solvers.telemac.sections.line_discharges(results, lines, *, geometry=None, name_field=None, crs_epsg=None)

Integrate the discharge across each line of the lines layer.

lines is a path to a line layer (any CRS - reprojected to crs_epsg) or an already-loaded GeoDataFrame. Returns a DataFrame with one row per line: name, discharge (m3/s, positive along the section’s net flow), wetted_width, mean_depth, max_depth, mean_velocity, max_velocity, orientation. The maxima are over the wet part of the section and exist for comparing against laboratory data reported per cross-section rather than per point (axqua.labdata).

Parameters:

• results (Path)

• geometry (Path | None)

• name_field (str | None)

• crs_epsg (int | None)

axqua.solvers.telemac.sections.write_line_discharges(results, lines, out, *, geometry=None, name_field=None, crs_epsg=None, name_column='Baffle Name', discharge_column='discharge (m3/s)')

Extract the section discharges and write the two-column CSV to out.

The CSV carries exactly the section name and its discharge in m3/s; the full table (wetted width, mean depth, mean velocity) is returned for logging.

Parameters:

• results (Path)

• out (Path)

• geometry (Path | None)

• name_field (str | None)

• crs_epsg (int | None)

• name_column (str)

• discharge_column (str)

TELEMAC-3D

Extrude the calibrated 2D case into a non-hydrostatic TELEMAC-3D run.

add3d.build_3d_cas turns the converged 2D steady result (the initial_run output r2d.slf) into a <case-name>3d.cas for telemac3d.py:

• Hotstart from the 2D run (TELEMAC v9+) - by default FILE FOR 2D CONTINUATION reads the final 2D state and initializes the 3D free surface + horizontal velocities from it (vertical velocity 0; horizontal velocity copied up each column). In v9 this keyword alone triggers the continuation (the old 2D CONTINUATION flag is gone). If a thin/supercritical 2D inflow makes telemac3d abort with DEBIMP_3D on the inflow boundary, swap in the commented-out cold-start CONSTANT DEPTH seed (median wetted 2D depth) - the hotstart="constant_depth" mode the vertical convergence study uses, which over-wets the inflow so the prescribed-Q boundary stays subcritical. INITIAL GUESS FOR DEPTH reuses the previous depth and a MINIMAL VALUE FOR DEPTH floor guards the sigma transform against dry (zero-depth) columns.

• Vertical discretization - TELEMAC has no Delft3D-style constant-dz Z-layer option; the vertical mesh is sigma planes (NUMBER OF HORIZONTAL LEVELS + MESH TRANSFORMATION : 1), so each water column is split into equal sigma slices and the physical dz is proportional to the local depth. We therefore infer the level count from the 2D results so a representative dz lands near the horizontal cell size: dz ~ dx/2 (TARGET_ASPECT) while never letting dz get more than 4x finer than dx (MAX_ASPECT, the user’s “max 4 times less” rule), clamped to a sane level count for a workable 3D mesh.

• Turbulence - the same quality-checked selection as 2D (axqua.steering.select_turbulence_model()), mapped onto the TELEMAC-3D numbering (k-epsilon 3, Smagorinski 4, Spalart-Allmaras 5 - note 3D numbers S-A 5, not the 2D 6). A vertical closure is paired with it (mixing length, or k-epsilon when the horizontal model is k-epsilon).

• Non-hydrostatic - NON-HYDROSTATIC VERSION : YES (with the pressure-Poisson solver), since a depth-resolved river model wants the non-hydrostatic pressure. non_hydrostatic=False writes the cheaper hydrostatic solver instead (no PPE keywords) - used for the long steady flux-convergence check add3d.py emits.

• Courant - TELEMAC-3D has no DESIRED COURANT NUMBER / variable time step, so we size the fixed TIME STEP to hit a target Courant of TARGET_COURANT (0.6) from the 2D velocity and the cell sizes.

TELEMAC-3D cannot run with finite volumes, so the 3D case is always the finite-element kernel (the 2D finite_volumes flag is irrelevant here). Zonal friction is not carried over - TELEMAC-3D has no FRICTION DATA FILE - so a single representative bottom-friction law/coefficient is written instead.

class axqua.solvers.telemac.threed.VerticalDiscretization(n_levels: 'int', dz: 'float', dx: 'float', depth: 'float', aspect: 'float', reason: 'str')

Bases: object

Parameters:

• n_levels (int)

• dz (float)

• dx (float)

• depth (float)

• aspect (float)

• reason (str)

n_levels: int

dz: float

dx: float

depth: float

aspect: float

reason: str

class axqua.solvers.telemac.threed.ThreeDSetup(cas: 'Path', n_levels: 'int', dz: 'float', dx: 'float', depth: 'float', time_step: 'float', n_time_steps: 'int', h_turbulence: 'int', v_turbulence: 'int', courant: 'float', turbulence_reason: 'str', layers_reason: 'str', non_hydrostatic: 'bool' = True)

Bases: object

Parameters:

• cas (Path)

• n_levels (int)

• dz (float)

• dx (float)

• depth (float)

• time_step (float)

• n_time_steps (int)

• h_turbulence (int)

• v_turbulence (int)

• courant (float)

• turbulence_reason (str)

• layers_reason (str)

• non_hydrostatic (bool)

cas: Path

n_levels: int

dz: float

dx: float

depth: float

time_step: float

n_time_steps: int

h_turbulence: int

v_turbulence: int

courant: float

turbulence_reason: str

layers_reason: str

non_hydrostatic: bool = True

axqua.solvers.telemac.threed.representative_scales(data)

(median wetted depth h, representative wetted horizontal cell size dx) [m].

Parameters:

data (dict)

Return type:

tuple[float, float]

axqua.solvers.telemac.threed.characteristic_velocity(cfg, data)

A high-percentile wetted flow speed [m/s] (for the Courant time step), floored at the configured velocity guess.

Parameters:

• cfg (Config)

• data (dict)

Return type:

float

axqua.solvers.telemac.threed.infer_vertical_layers(cfg, data)

Infer NUMBER OF HORIZONTAL LEVELS from the 2D results.

The sigma layer thickness at the representative depth h is dz = h/(n-1). We aim for dz ~ dx/2 (TARGET_ASPECT) but never thinner than dx/4 (MAX_ASPECT - the user’s “max 4x finer than dx/dy” rule), and keep the level count within [N_LEVELS_MIN, N_LEVELS_MAX] so the 3D mesh stays workable.

Parameters:

• cfg (Config)

• data (dict)

Return type:

VerticalDiscretization

axqua.solvers.telemac.threed.select_3d_turbulence(cfg, mesh)

(horizontal model, vertical model, rationale) for TELEMAC-3D.

Reuses the 2D mesh-resolution/velocity selection and maps it onto the 3D numbering; pairs a vertical closure (k-epsilon when horizontal is k-epsilon, else a mixing-length model, which is robust for the vertical log profile).

Parameters:

cfg (Config)

Return type:

tuple[int, int, str]

axqua.solvers.telemac.threed.courant_time_step(dx, dz, u_char, courant=0.6)

Fixed time step [s] giving an advective Courant ~ courant on the limiting (smallest) cell dimension. TELEMAC-3D has no DESIRED COURANT NUMBER keyword, so this is how the requested Courant is honored.

Parameters:

• dx (float)

• dz (float)

• u_char (float)

• courant (float)

Return type:

float

axqua.solvers.telemac.threed.cas3d_name(cfg, suffix='')

<case-name>3d<suffix>.cas with a filesystem-safe case name.

suffix distinguishes variants written next to each other (e.g. "-unsteady" -> <case-name>3d-unsteady.cas).

Parameters:

• cfg (Config)

• suffix (str)

Return type:

str

axqua.solvers.telemac.threed.build_3d_cas(cfg, results_2d=None, cas_2d=None, total_time_factor=4.0, n_levels=None, hotstart='continuation', unsteady=False, liquid_boundaries_file=None, duration=None, gaia_cas=None, out_name=None, results3d_name=None, results2d_name=None, non_hydrostatic=True, n_time_steps=None, listing_period=None, graphic_period=None, min_depth=None, data=None)

Write <case-name>3d.cas next to the 2D case and return the setup summary.

results_2d defaults to <model_dir>/<results_slf> (the initial_run output); cas_2d to the steady 2D steering, whose prescribed boundary values are reused. total_time_factor sets the simulated duration as that many reach flow-through times (length / characteristic velocity). n_levels overrides the inferred NUMBER OF HORIZONTAL LEVELS (used by the vertical convergence study; the time step then adapts to the resulting layer thickness).

With unsteady and a liquid_boundaries_file the run is hydrograph-driven: the LIQUID BOUNDARIES FILE supplies Q(t)/SL(t) and the step count spans duration (TELEMAC-3D has no DURATION keyword, so NUMBER OF TIME STEPS = duration / dt); the hotstart still lifts the initial 3D field from the 2D result. gaia_cas couples GAIA morphodynamics. out_name / results3d_name / results2d_name override the output filenames so an unsteady variant does not clobber the steady 3D run’s files.

non_hydrostatic toggles NON-HYDROSTATIC VERSION (False writes the hydrostatic solver and drops the pressure-Poisson keywords - the cheap variant for a long boundary-flux convergence run). n_time_steps forces the step count (overriding the flow-through/duration sizing); listing_period / graphic_period override the printout spacing (e.g. a short listing period so the flux-convergence analysis gets enough printouts). data is an optional pre-read 2D result (selafin.read_slf) so several variants can be built from one read of a large results file.

hotstart selects how the 3D run is initialized from the 2D result:

• "continuation" (default) - the TELEMAC v9+ true 2D->3D hotstart (FILE FOR 2D CONTINUATION): the 3D free surface + horizontal velocities are lifted from the final 2D state. The most faithful start, but a thin/supercritical 2D inflow can make telemac3d abort with DEBIMP_3D: PROBLEM ON BOUNDARY; the CONSTANT DEPTH fallback is then emitted commented-out to swap in.

• "constant_depth" - a cold start at the median wetted 2D depth (uniform), so the prescribed-Q inflow is deeply wet and robust. Used by the vertical convergence study, where every per-layer run must reach steady state.

Parameters:

• cfg (Config)

• results_2d (str | Path | None)

• cas_2d (str | Path | None)

• total_time_factor (float)

• n_levels (int | None)

• hotstart (str)

• unsteady (bool)

• liquid_boundaries_file (str | None)

• duration (float | None)

• gaia_cas (str | None)

• out_name (str | None)

• results3d_name (str | None)

• results2d_name (str | None)

• non_hydrostatic (bool)

• n_time_steps (int | None)

• listing_period (int | None)

• graphic_period (int | None)

• min_depth (float | None)

• data (dict | None)

Return type:

ThreeDSetup

axqua.solvers.telemac.threed.build_3d_cases(cfg, hydrostatic_steps=None, hydrostatic_listing=None)

Write the three 3D steering files add3d.py provides and return their setups as {"hydrostatic": ThreeDSetup, "hydrodyn": ThreeDSetup, "unsteady": Unsteady3DSetup | None}.

1. hotstart3d_hydrostatic.cas - NON-HYDROSTATIC VERSION : NO (the cheap hydrostatic solver), constant in-file Q/H, marching hydrostatic_steps (default HYDROSTATIC_N_STEPS) fixed steps with a short listing period (hydrostatic_listing, default HYDROSTATIC_LISTING_PERIOD) so the sortie carries enough FLUX BOUNDARY printouts to check when the boundary fluxes converge in a steady telemac3d run.

2. hotstart3d_hydrodyn.cas - the steady non-hydrostatic run with the same in-file prescribed Q (inflow) and H (outflow).

3. unsteady3d.cas - non-hydrostatic, hydrograph-driven: Q(t) inflow and the stage-discharge outflow SL(t) come from the same LIQUID BOUNDARIES FILE as unsteady2d.cas. Requires a varying boundaries.inflow series; None (with a log notice) otherwise.

All three hotstart from the 2D steady result (v9 2D->3D continuation) and share the sigma-level count, turbulence closure and Courant-sized time step inferred from one read of r2d.slf.

Parameters:

• cfg (Config)

• hydrostatic_steps (int | None)

• hydrostatic_listing (int | None)

Return type:

dict

Varying discharge

Unsteady (quasi-steady hydrograph) TELEMAC runs, built on the converged 2D case.

The steady initial run (preprocessing.py + initial_run.py) establishes a well-converged flow field and writes it to r2d.slf. This module turns that into a hydrograph-driven run - a natural flood wave discretized as a sequence of steady discharges (see telemac2d-unsteady.md) - hotstarted from that steady result, for 2D (telemac2d) or, with the same forcing, 3D (telemac3d).

Design decisions (deliberately different from the manual-Notepad tutorial):

• Boundaries are NOT re-defined. axqua’s generated boundaries.cli is already correct for the unsteady case - the inflow is 4 5 5 (prescribed Q, depth free) and the outflow 5 4 4 (prescribed elevation). So there is no 5 5 5 -> 4 5 5 find-and-replace and no separate boundaries-unsteady.cli: the inflow discharge is driven by the LIQUID BOUNDARIES FILE (Q(t)) and the downstream water level by its prescribed elevation (constant, or the rating stage SL(t) at each Q when a stage-discharge rating exists). The tutorial’s STAGE-DISCHARGE CURVES route is the alternative; we take the liquid-file route because it needs no .cli edit and reuses the steady boundary numbering.

• Hotstart from the steady result (PREVIOUS COMPUTATION FILE : r2d.slf, clock reset to zero) instead of re-seeding a dry/pre-wetted bed - the flood wave then propagates from an already-converged field.

• Quasi-steady numerics: NUMBER OF CORRECTIONS OF DISTRIBUTIVE SCHEMES : 2 (the developers’ recommendation for predictor-corrector schemes on quasi-steady flows), and optional CONTROL SECTIONS so the run reports the flux across each open boundary (verifying Q_in(t) vs Q_out(t)).

• GAIA (bedload / suspended load) is coupled when morphodynamics.enabled.

The routines recover the FRONT2 liquid-boundary numbering from the sidecar liquid-boundaries.json the build wrote (axqua.boundary.dump_liquid_boundaries()), so no mesh rebuild is needed.

class axqua.solvers.telemac.unsteady.UnsteadySetup(cas: 'Path', liquid_boundaries: 'Path', control_sections: 'Path | None', gaia_cas: 'Path | None', duration: 'float', results: 'Path', previous_computation: 'str', turbulence_model: 'int', n_inflows: 'int', n_outflows: 'int')

Bases: object

Parameters:

• cas (Path)

• liquid_boundaries (Path)

• control_sections (Path | None)

• gaia_cas (Path | None)

• duration (float)

• results (Path)

• previous_computation (str)

• turbulence_model (int)

• n_inflows (int)

• n_outflows (int)

cas: Path

liquid_boundaries: Path

control_sections: Path | None

gaia_cas: Path | None

duration: float

results: Path

previous_computation: str

turbulence_model: int

n_inflows: int

n_outflows: int

class axqua.solvers.telemac.unsteady.Unsteady3DSetup(cas: 'Path', liquid_boundaries: 'Path', gaia_cas: 'Path | None', duration: 'float', n_levels: 'int', dz: 'float', time_step: 'float', n_time_steps: 'int', results3d: 'Path')

Bases: object

Parameters:

• cas (Path)

• liquid_boundaries (Path)

• gaia_cas (Path | None)

• duration (float)

• n_levels (int)

• dz (float)

• time_step (float)

• n_time_steps (int)

• results3d (Path)

cas: Path

liquid_boundaries: Path

gaia_cas: Path | None

duration: float

n_levels: int

dz: float

time_step: float

n_time_steps: int

results3d: Path

axqua.solvers.telemac.unsteady.load_hydrograph(cfg)

Load the inflow hydrograph Q(t) from boundaries.inflow (the full series).

Raises with an actionable message when no series is configured or when it is a single constant value (a genuine hydrograph is required for an unsteady run - a constant discharge is the steady case).

Parameters:

cfg (Config)

Return type:

Inflow

axqua.solvers.telemac.unsteady.write_control_sections(cfg)

Write a CONTROL SECTIONS file: one cross-section per open-boundary segment.

Parses the generated boundaries.cli and, for each contiguous run of inflow / outflow nodes (in file / IPOBO order), emits a section defined by that run’s first and last global node numbers (node-id mode, -1). TELEMAC then reports the cumulated flux across each section, so the unsteady run can be checked for mass balance (Q_in(t) vs Q_out(t)). Returns None if the .cli has no liquid nodes. (A boundary that wraps the .cli start/end is reported as two half sections; sum them when reading the output.)

Parameters:

cfg (Config)

Return type:

Path | None

axqua.solvers.telemac.unsteady.build_unsteady_case(cfg, *, control_sections=True)

Write the unsteady telemac2d case, hotstarted from the steady r2d.slf.

Reuses boundaries.cli unchanged, drives the inflow with the hydrograph liquid file and the outflow with its prescribed elevation / rating, adds the quasi-steady correction keyword and (optionally) control sections, and couples GAIA when morphodynamics.enabled. Returns the UnsteadySetup.

Parameters:

• cfg (Config)

• control_sections (bool)

Return type:

UnsteadySetup

axqua.solvers.telemac.unsteady.build_unsteady_3d_case(cfg, *, total_time_factor=4.0, out_name=None, data=None)

Write the unsteady telemac3d case: the same hydrograph forcing on the 3D, non-hydrostatic mesh, hotstarted from the 2D result.

Writes the hydrograph liquid file (as for 2D) and delegates to axqua.threed.build_3d_cas() with the unsteady forcing (the 3D case has no DURATION keyword, so the run spans the hydrograph via its step count). Couples GAIA when enabled. The 3D horizontal mesh should already be validated (2D mesh-convergence) before running this. out_name overrides the steering filename (default <case-name>3d-unsteady.cas); data is an optional pre-read 2D result so several 3D variants share one read of a large results file.

Parameters:

• cfg (Config)

• total_time_factor (float)

• out_name (str | None)

• data (dict | None)

Return type:

Unsteady3DSetup

Exchange with a porous body

Where a gain-lose reach exchanges water, derived from the water table.

A river that loses flow into a porous body and regains it downstream has no sharp exchange faces: it is genuinely fuzzy where infiltration begins and ends. Asking the user to draw the two faces therefore asks for a number they do not have.

The water table supplies it instead. The body’s saturated zone is bounded by the two channel levels it exchanges with, so once the phreatic surface is known (axqua.solvers.telemac.watertable) each node’s role follows from a head comparison:

• it loses where it is wet and its free surface stands above the table - there is head pushing water down into the body;

• it gains where the table stands above the bed - the saturated zone cuts the ground surface and groundwater emerges.

Both conditions are evaluated inside the zone polygon, optionally buffered (gain_lose.zone_buffer): a bar polygon usually outlines the bar itself, while the water it takes in is in the channel beside it, so the losing face has to be allowed to reach the adjacent wetted edge.

This module is the build-time view - it reports the faces and the exchange they imply so both are visible before a run. The generated USER_RAIN routine (axqua.solvers.telemac.fortran) applies the same two rules at run time, so the faces move with the stage instead of being frozen at build time.

class axqua.solvers.telemac.gainlose.ExchangeFaces(losing, gaining, losing_area, gaining_area, losing_head, gaining_head, plane)

Bases: object

The losing and gaining faces of a gain-lose reach at one water level.

Parameters:

• losing (object)

• gaining (object)

• losing_area (float)

• gaining_area (float)

• losing_head (float)

• gaining_head (float)

• plane (object)

losing: object

gaining: object

losing_area: float

gaining_area: float

losing_head: float

gaining_head: float

plane: object

summary()

Return type:

list[str]

axqua.solvers.telemac.gainlose.derive_faces(cfg, mesh, plane, surface, *, tol=0.01)

Split the porous body into its losing and gaining faces at surface.

surface is a per-node free-surface elevation (bed + depth). A node is assigned to at most one face: standing water pushing head into the body wins over emergence, so a node whose surface is above the table loses even if the table is also above its bed.

Raises when either face is empty - a gain-lose reach with nowhere to lose or nowhere to gain is a geometry error the user needs to see, not something to paper over with a zero exchange.

Parameters:

tol (float)

Return type:

ExchangeFaces

axqua.solvers.telemac.gainlose.estimate_exchange(cfg, mesh, faces, surface)

Discharge [m3/s] the losing face implies at the configured conductivity.

Green-Ampt’s saturated limit per node, f = kf (h + Lz + hf) / Lz with h the head over the table, integrated over the losing face. Logged at build time so the modeled exchange is visible before a run - it is easy to be an order of magnitude out on kf, and this is the number to compare against a measured one.

Parameters:

faces (ExchangeFaces)

Return type:

float

axqua.solvers.telemac.gainlose.zone_mask(cfg, mesh, *, buffer=None)

Nodes inside the porous body, buffered by gain_lose.zone_buffer.

The buffer defaults to one channel cell: without it the losing face of a bar polygon is only the bar’s own toe, while the water actually being lost is in the channel beside it.

Parameters:

buffer (float | None)

The water table under a porous gravel bar, as a phreatic plane.

A losing-gaining gravel bar exchanges with the channel at two places: it takes water in at the losing line and returns it at the gaining line (see axqua.solvers.telemac.fortran). The bar’s saturated zone is therefore bounded by two known water levels, and the surface joining them is the water table of that through-flow - its gradient is what drives the exchange in the first place. Over a bar of a hundred meters or so that surface is planar to within a few centimeters, so it is represented here as

\[z_{wt}(x, y) = c_0 + c_x (x - x_0) + c_y (y - y_0)\]

fitted by least squares through sample points along the internal lines, each carrying its own water level.

Two consumers need it, and they need it to agree:

• the pre-wet (axqua.steering.write_initial_conditions()) seeds any bar ground lying below the table, because a closed depression on a bar above the channel can never fill from the surface - no flow path reaches it - and the drainable-seed filter deliberately refuses to seed it. Without this a real pool stays dry for the whole run;

• the patch drain in the generated USER_RAIN routine tapers to zero at the table rather than at an absolute depth, so it clears standing water off the bar top but cannot empty a pool that cuts below the water table.

Because a plane is five numbers, the Fortran side needs no per-node array: the coefficients are baked into the generated routine.

The plane is only meaningful inside the patch. Extended across a whole reach it sits above any ground lower than the bar - on isar-2025, 22 598 m² of it. Every use therefore goes through water_table_depth(), which is the single place the clip-to-the-patch rule lives.

class axqua.solvers.telemac.watertable.PhreaticPlane(c0, cx, cy, x0, y0, levels, residual=0.0)

Bases: object

A planar water table z = c0 + cx (x - x0) + cy (y - y0) [m a.s.l.].

Parameters:

• c0 (float)

• cx (float)

• cy (float)

• x0 (float)

• y0 (float)

• levels (dict)

• residual (float)

c0: float

cx: float

cy: float

x0: float

y0: float

levels: dict

residual: float = 0.0

property gradient: float

Slope of the table [-] (its steepest descent).

elevation(x, y)

Water-table elevation at x, y (arrays or scalars).

axqua.solvers.telemac.watertable.describe(plane, mesh, mask, area=None)

One-line summary of what a table would wet (for the build log).

Parameters:

plane (PhreaticPlane | None)

Return type:

str

axqua.solvers.telemac.watertable.fit_phreatic_plane(cfg, mesh, *, surface=None, faces=None)

Fit the water table of the porous body.

Two ways to find the two bounding levels, chosen by faces (default gain_lose.faces):

• water-table - from the zone polygon alone: the channel water surface at each end of the zone’s reach extent (_samples_from_zone());

• lines - from the int-* lines the user drew (_samples_from_lines()).

gain_lose.water_table_levels overrides either. Returns None (with a warning) when fewer than two levels can be established, since a plane through one level is not determined.

Return type:

PhreaticPlane | None

axqua.solvers.telemac.watertable.patch_node_mask(cfg, mesh)

Boolean (NPOIN,) flag of the nodes inside the percolation patch(es).

axqua.solvers.telemac.watertable.water_table_depth(plane, mesh, mask)

Water depth [m] the table supports: max(table - bed, 0) inside mask.

Zero everywhere outside mask. This is the only place the clip is applied: a plane fitted to a gravel bar is meaningless across the rest of a reach, where it would stand above every bed lower than the bar.

Parameters:

plane (PhreaticPlane)

Generated TELEMAC user-Fortran for percolation (percolation.mode: fortran).

Emits a USER_RAIN routine (the TELEMAC-2D hook for spatialised rain / infiltration, called from prosou.f whenever RAIN OR EVAPORATION : YES) that models the losing-gaining exchange of a porous patch depth-limited:

• over each losing region (the percolation patch polygon) it withdraws the target discharge as a uniform depth rate Q/area, tapered linearly to zero as the local depth approaches percolation.min_depth and additionally capped at half the water available above that floor per time step - so the sink can never dry a cell (TELEMAC’s own source terms have no such guard, which is what collapses the CFL-adaptive time step when a sink node dries);

• the discharge actually extracted (which may be less than the target while the patch is still shallow) is summed - in parallel via P_SUM, using the same partial-VOLU2D idiom as telemac2d_init.F - and reinjected exactly over the gaining region(s), so the exchange is mass-exact by construction;

• PLUIE is assigned, not accumulated, so the double CALL USER_RAIN in prosou.f (once inside the RAIN block, once unconditionally) is harmless.

The generated file lands in <model_dir>/<user_fortran_dir>/user_rain.f and is referenced by FORTRAN FILE in the .cas (see steering.write_cas()); TELEMAC compiles it into the executable at run time.

axqua.solvers.telemac.fortran.KIND_PRESCRIBED = 1

region kinds dispatched inside the generated routine

axqua.solvers.telemac.fortran.patch_drain_regions(cfg, used)

Percolation patches to drain, as extra InternalSourceRegion regions.

Standing surface water on a porous gravel bar infiltrates instead of ponding, so it is an artifact of a 2D model that has no subsurface. With percolation.losing_region: line the prescribed exchange is taken from the channel at the losing line and nothing removes water sitting on the bar; these regions do (see percolation.patch_drain).

Returns an empty list when the drain is off, or when the patch is already the losing region (losing_region: patch), where the prescribed withdrawal already covers it.

Parameters:

• cfg (Config)

• used (list)

Return type:

list

axqua.solvers.telemac.fortran.write_user_fortran(cfg, regions, *, plane=None)

Generate user_fortran/user_rain.f for the given source regions.

regions is the axqua.boundary.load_internal_source_regions() list (signed discharges; in fortran mode the losing region polygon is the percolation patch). plane is an optional axqua.watertable.PhreaticPlane: when given, the patch drain tapers to zero at the water table instead of at percolation.min_depth, so it cannot empty a pool that cuts below the bar’s saturated zone. A plane is five numbers, so it is baked in as coefficients - no per-node array is needed.

Returns the user-fortran directory (the FORTRAN FILE target).

Parameters:

• cfg (Config)

• regions (list)

Return type:

Path

axqua.solvers.telemac.fortran.write_water_table_fortran(cfg, plane)

Generate user_rain.f for a gain-lose reach whose faces come from the water table (gain_lose.faces: water-table).

Unlike the lines route, nothing about where the exchange happens is baked in. Only the zone polygon and the five plane coefficients are; each call then classifies every zone node from its own head:

• free surface above the table -> losing (head pushes water into the body);

• table above the bed -> gaining (the saturated zone cuts the ground).

So the faces move with the stage - a rising river widens the losing face, which is what actually happens and what a build-time mask cannot represent.

Magnitude is Green-Ampt at conductivity, or discharge normalized over the taper-weighted wet losing area when a measured total is prescribed. Whatever is withdrawn is reinjected over the gaining face in the same step, so the routine is mass-exact by construction and opens no net sink in the boundary budget.

Parameters:

cfg (Config)

Return type:

Path

Backend adapter

The TELEMAC backend: what the job system calls, adapting what already exists.

Deliberately thin. Every piece of work here is a function the case-template scripts have been calling for months - pipeline.run, run_solver_streaming, analyze_flux_convergence, run_mesh_convergence, run_single_flow_calibration. The backend adds no hydraulics and no numerics; it maps a verb plus a capability onto the right one of them and passes the job’s event sink and cancel token down.

That thinness is the point. cases/*/initial_run.py and a submitted steady job must drive the same code, because the standalone scripts remain the supported path (plan §29) and a second implementation would drift the moment one of them was fixed.

This module is loaded only through BackendSpec.load(). It imports gmsh, numpy and the rest transitively, so it must never be imported by spec.py or by the package __init__ - tests/test_capabilities.py asserts exactly that from a subprocess.

class axqua.solvers.telemac.backend.TelemacBackend

Bases: BaseBackend

TELEMAC-2D / TELEMAC-3D, with GAIA where the case asks for it.

name: str = 'telemac'

check_environment(cfg)

True when the TELEMAC environment can actually be entered.

Must not raise (the protocol says so), because this is also called from axqua status --check-env where a broken install is information, not a crash.

Return type:

bool

build(cfg, capability=None, *, ctx=None, **kwargs)

Parameters:

• capability (Capability | None)

• ctx (Any)

run(cfg, capability=None, *, ctx=None, **kwargs)

Parameters:

• capability (Capability | None)

• ctx (Any)

study(cfg, capability=None, *, ctx=None, **kwargs)

Parameters:

• capability (Capability | None)

• ctx (Any)

postprocess(cfg, ctx=None)

Flux convergence, wetting and section reports - what initial_run.py does.

Only for the kinds where those reports mean something: running a flux analysis after a mesh-convergence study would read the last level’s listing and report it as if it were the case’s own.

Parameters:

ctx (Any)

extract_objective(cfg, ctx=None)

The number this job was minimizing, when there is one.

A steady run’s is the relative flux imbalance - the thing its own convergence is judged on - which makes “how good was this run” answerable from the job record without opening anything.

Parameters:

ctx (Any)

Return type:

float | None

export_qgis_results(cfg, ctx=None)

Describe what QGIS should open. No conversion, and no PyQGIS.

TELEMAC results are SELAFIN and QGIS reads SELAFIN natively through MDAL, so converting them would cost a full copy of a multi-hundred-megabyte file to gain nothing (plan §26). The export is therefore a manifest: paths plus a hint about how to style each one.

Parameters:

ctx (Any)