Add the optional MODFLOW groundwater coupling - #357
Open
soaressgabriel wants to merge 10 commits into
Open
soaressgabriel wants to merge 10 commits into
soaressgabriel wants to merge 10 commits into
Conversation
Port of rubem/hydrological_processes/_modflow.py, rubem/configuration/modflow_configuration.py and the user-guide text from the branches feat_modflow_module and feat_calibrator_modflow at commit da41faa (commits ed24a30, 0954a20, 079b833, 299a51e, 63ebddb, c70370f, da41faa). The modules are not wired into the model yet; the configuration schema, validation, dynamic-model coupling, calibrator support, tests and documentation follow in the next commits of #356.
ruff format only, no code change; refs #356.
Rework the prototype MODFLOW schema into ModflowSettings: layers are listed from the top down and numbered like MODFLOW (layer 1 is the top), each with its own bottom map under a model top; pcraster_layer() and user_layer() convert to and from the bottom-up numbering of the PCRaster extension. RIV, GHB and DRN entries and the BCF wetting attach to explicit layer lists. The section is strict in both formats, so the prototype keys that no longer exist (wells, recharge, dis units, solver type, fail_on_non_convergence, baseflow_from_river_leakage, per-layer top, wetting multiplier and source layer) are refused by name instead of silently disabling a package. enabled is a boolean (0/1 still accepted). When the section is enabled, every rule is checked and reported at once: model top and layers, unique layer names, the river package (its leakage is the baseflow), layer numbers in range and once per package, storage per LAYCON in transient runs, LAYCON 1 only on the top layer, wetting only on LAYCON 1 or 3 layers and with its WETDRY map, a mask for a numeric river conductance, the water-table layer and the minimum root depth table of the root-depth coupling. The section is the MODFLOW key of the legacy file (null when absent) and the modflow key of format 1.0 (left out when absent); both anchor its paths on the configuration directory, the conversions carry it, and the migration rebases its paths like the others. ModelConfiguration exposes it with modflow_enabled and reports a present but disabled section as ignored. Refs #356
An enabled MODFLOW section is now checked when the configuration loads, so a coupled run or a calibration refuses bad inputs up front instead of failing inside the PCRaster MODFLOW extension, which aborts the process on some of them. check_modflow_inputs (rubem/validation/modflow_inputs.py) always reports, as blocking, a missing or empty file the run reads (files of disabled packages, wetting and root-depth coupling are not read) and a missing MODFLOW runtime. With validate_input it reads every raster on the clone and reports: rasters that cannot be read or do not share the clone geometry (PCRaster reads a map of another size silently), boundary values other than -1/0/1, a top or bottom missing or not strictly decreasing on columns active in a layer, heads, conductivities and storage missing on active cells, heads below the layer bottom (blocking in every cell, a warning in some), conductivity classes not covered by a positive table value, an enabled package layer without cells, package heads missing on package cells, riverbeds outside the elevation interval of their layer (non-blocking count) and a Dpz_min table that does not cover the soil classes or leaves (0, Zr]. rubem/_deps.py gains the groundwater tier: missing_groundwater_deps, groundwater_deps_message, resolve_mf2005 (PATH, then the interpreter prefix where conda installs mf2005) and ensure_mf2005_on_path. A configuration without the section runs no check and reports the same problems as before. Refs #356
Rework the ported groundwater module around the validated ModflowSettings. ModflowGroundwater(settings, cell_area_m2, run_directory, logger) reads every input once in initialize() with read_field, converts the top-down layer numbers of the configuration with pcraster_layer(), and gives the extension its packages in the order DIS, BAS, BCF, wetting, solver, then RIV, GHB and DRN. The stress packages are set once: the extension keeps them for every later run, so no map is read again per step. Conductance is zeroed outside the mask and the layer boundary, and the active cells are counted per package layer: a layer without cells is never set nor read, because the extension ends the process on a getter for a package layer without cells. Columns inactive in every layer keep the prototype's synthetic 1 m elevations. The kh lookup still reads a fresh copy of its table (PCRaster caches tables by name), now inside the run directory. run_step() updates PERLEN when transient, sets RCH (mm/step to m/day, highest active cell), runs mf2005 in the run directory so no pcrmf.* file lands in the working directory, raises RuntimeError on non-convergence (the next run would end the process), and returns the river exchange, the baseflow in mm and the heads, storage, drain flow and water-table head keyed by the configured layer numbers. The water table reads the cached boundaries and surfaces instead of re-reading maps each step. pcraster.initialise is imported only when a model starts, after ensure_mf2005_on_path(), and a missing extension or executable raises RuntimeError with the installation guidance. The prototype's own validation pass is gone (rubem.validation.modflow_inputs owns the rules); the module only refuses missing values on the cells it uses and layers that are not stacked, which the extension cannot take whatever validate_input says. Tests: a mocked extension for the call order, the numbering conversion, set-once stress data, skipped empty layers, unit conversions and the water table methods; the prototype's analytic 3x3 RIV, GHB, DRN and RCH balances against the real mf2005 (skipped when resolve_mf2005() finds none). Refs #356
When the MODFLOW section is enabled, initial() creates <output directory>/modflow, starts ModflowGroundwater there with the days of the first month and, with the root-depth coupling, keeps the DEM as the terrain surface and reads Dpz_min per soil class. The module is imported inside initial(), so a run without the section never loads it. Each step pushes the RUBEM recharge to MODFLOW for the days of the month and takes the aquifer-to-river RIV leakage of the result as bfw; the saturated-zone reservoir is then not updated. The root-depth coupling uses the water table of the previous step (one-step lag): the effective root depth is the water-table depth below the terrain bounded by Dpz_min and Zr, and only the vegetated-area water stress coefficient sees the storages of that root fraction; the bare soil keeps the whole-rootzone coefficient and the soil balance is unchanged. The enabled diagnostics go through the same two format branches as the output variables (mfh<n>, mfaq2rv, mfrv2aq, mfrvnet, mfst<n>, mfdrn<n>, and mfwt, mfgwd, mfzr, mfzfrac with output.root_depth), so an empty format set writes nothing. The run directory is removed after the last step's report and kept when a step raised, with pcrmf.lst for inspection. Without the section, or with it disabled, every output is identical to the run before this change. Refs #356
A configuration that enables MODFLOW is calibrated with the coupled model in every evaluation. Its MODFLOW parameters join the search only when the caller names them, so a calibration without MODFLOW names searches and reports exactly what it did before. rubem.calibration.modflow_parameters builds the catalog of a section: modflow.layers.<n>.specific_yield and .specific_storage (numbers the transient run reads for the layer's LAYCON), modflow.layers.<n>.kh.<class> (numeric rows of the layer's conductivity lookup table) and modflow.river.<i>.conductance (numeric entries). ModflowCatalog.apply writes a candidate into a format 1.0 document: numbers in place, named classes into a rewritten copy of the table whose path replaces the configured one. decision_space() takes the catalog. A MODFLOW name is searched only when bounded (finite, minimum below maximum, strictly positive, at most 1 for a specific yield), is appended after the eight hydrological parameters, and a fixed one keeps its value in every candidate; an unknown name lists the known ones, and a MODFLOW name without MODFLOW enabled is refused. The positions of w_1 and w_2, hence the weights constraint, are unchanged. The runner builds the catalog after the validated load, resolves the decision space there, starts from the configured MODFLOW values, puts mf2005 on the PATH the workers inherit, adds the MODFLOW names to evaluations.csv after x, reruns the best candidate with its MODFLOW values, and writes them into the calibrated configuration with its <name>-kh<n>.tbl tables. A dead worker of a coupled calibration also names non-convergence before the last time step of a period (dis.nstp > 1), which ends the worker process instead of failing the evaluation. The --bound and --fix help texts name the MODFLOW parameters. Refs #356
Add the Groundwater Coupling page (doc/source/groundwater.rst, in the toctree after the calibration page): what the coupling exchanges (recharge to RCH in m/day, the aquifer-to-river RIV leakage as bfw in mm, the saturated-zone reservoir no longer updated, GHB/DRN/wetting shaping the heads only), the runtime it needs, the top-down layer numbering, every key of the section with its type, default and unit, the rules the configuration and the validation enforce, the raster diagnostics, the run directory, the experimental root-depth coupling and the limitations, among them the solver failure before the last time step of a period with dis.nstp above 1, which ends the process instead of raising. Replace the prototype text appended to the user guide, which numbered the layers from the bottom up and described keys that no longer exist, with a MODFLOW section listing the keys in the style of the other entries, add the section (disabled) to the configuration file template, and state what the coupling changes for bfw, for -s and in the calibrate help. Add a MODFLOW parameters section to the calibration page (names, the values they may take, when a name exists, the mandatory bound, an example) and update the sentences that the MODFLOW names make incomplete: what a calibration changes, the budget, the columns of evaluations.csv, result.json, the calibrated configuration and its kh tables, the process model and the dead-worker message. Add the changelog entry and a README feature bullet, and a unit test that parses every JSON example of the new page and validates its MODFLOW section, so the page cannot drift from the schema. Refs #356
Default dis.nstp to 1, one MODFLOW time step per stress period. The
PCRaster MODFLOW extension reads the heads of a period only when its last
time step converged: with several time steps a solver failure at an earlier
one leaves no head file, and the extension ends the whole process inside
run() ("Can not open head value result file") before converged() can be
asked. With one time step a failure is always at the last one, and the run
raises the RuntimeError that names the period and the listing, which the
calibrator records as a failed evaluation ranked last, as decided. The
scientist's value, 5, stays available as an explicit choice with the
documented limitation.
Warn in rubem calibrate when the configuration enables MODFLOW with
dis.nstp above 1: a candidate that fails before the last time step then
ends its worker process, and the pool cannot replace it, so the whole
search ends on one bad candidate. The value is not refused because the
discretization of a calibration should be the one of the production run.
Tests: the child-process non-convergence test now runs the default section
and asserts the raised error, and the former characterization test pins the
opt-in nstp 5 case; the calibrator warning has a positive and a complement
case; the defaults and the mocked DIS calls read 1. Documentation: the dis
table, the Limitations entry, the user guide entry, the calibration
warning and the changelog state the default and the warning.
Refs #356
An initial head equal to -888, -999, -999.9, -999.99 or -9999 in an active cell is a blocking problem that names the layer, the count and the first cell. The extension reads such a marker as a head: the cell starts dry and the marker reaches the head outputs, which the Batalha exercise showed (5041 cells of one map at -888 came out as heads of -1393 m). The below-bottom warning no longer counts those cells; a marker on an inactive cell stays accepted. Documented in the blocking rules of the groundwater page and in the changelog. Refs #356
Codecov Report❌ Patch coverage is
Additional details and impacted files@@ Coverage Diff @@
## main #357 +/- ##
==========================================
+ Coverage 93.36% 94.19% +0.83%
==========================================
Files 68 72 +4
Lines 5260 6394 +1134
Branches 675 888 +213
==========================================
+ Hits 4911 6023 +1112
- Misses 274 285 +11
- Partials 75 86 +11 ☔ View full report in Codecov by Harness. 🚀 New features to boost your workflow:
|
This was referenced Sep 24, 2026
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment
2 participants
Add this suggestion to a batch that can be applied as a single commit.This suggestion is invalid because no changes were made to the code.Suggestions cannot be applied while the pull request is closed.Suggestions cannot be applied while viewing a subset of changes.Only one suggestion per line can be applied in a batch.Add this suggestion to a batch that can be applied as a single commit.Applying suggestions on deleted lines is not supported.You must change the existing code in this line in order to create a valid suggestion.Outdated suggestions cannot be applied.This suggestion has been applied or marked resolved.Suggestions cannot be applied from pending reviews.Suggestions cannot be applied on multi-line comments.Suggestions cannot be applied while the pull request is queued to merge.Suggestion cannot be applied right now. Please check back later.
Checklist
Description
Optional coupling of RUBEM with MODFLOW-2005 through the PCRaster MODFLOW extension (
pcraster.initialiseand themf2005executable, both shipped by the conda-forgepcrasterpackage on linux-64 and win-64). The first commit is LINA OSORIO's prototype (feat_modflow_module/feat_calibrator_modflow, commits ed24a30..da41faa) ported as she wrote it, with her authorship; the next commits rework it by concern. Without aMODFLOWsection nothing changes: the default path is proven byte-identical by a test that runs the synthetic dataset with the section absent, present-but-disabled and{"enabled": false}and compares every raster and CSV, and by theexactjob.rubem/configuration/modflow_configuration.py, both file formats,rubem config schemaandrubem config migrate): sectionMODFLOW(legacy) /modflow(1.0),extra="forbid"in both. Layers are listed from the top down and numbered like MODFLOW (layer 1 = top); the code converts to the extension's bottom-up numbers (ModflowSettings.pcraster_layer). Each layer has its ownbottom,initial_head,boundary,laytype,horizontal_conductivity(map, number or{map, table}class lookup),vertical_conductivity,specific_storage,specific_yield; a modeltopmap;dis,solver(PCG),wetting,river/ghb/drainwithentriesthat list their target layers,coupling.dynamic_root_depth,output. Rules enforced by the model:river.enabledwhen the section is enabled, LAYCON 1 only on the top layer, wetting only on LAYCON 1/3 layers, storage keys per LAYCON in transient runs, every layer reference in range and once per package, at most 9 layers (themfh<n>output names),output.storageonly in transient runs, numeric river conductance needsmask; every violation is reported at once. A present but disabled section yields the non-blocking "MODFLOW section is ignored." problem.rubem/validation/modflow_inputs.py, fourth tier ofrubem/_deps.py): always, the files the run reads exist and the runtime is available (pcraster._pcraster_modflowimportable,mf2005on PATH or under the interpreter prefix, prepended to PATH for the extension and the calibration workers), both blocking even with-s; with validation on, the clone geometry of every map, boundary values in {-1, 0, 1},topabove the bottoms and bottoms decreasing downward on active columns, finite heads/conductivities/storage on active cells, heads below the layer bottom (blocking when every active cell, warning with the count otherwise), kh tables covering the classes of the active cells with positive values, at least one cell per enabled package layer (the extension ends the process on a getter for an empty package), the riverbed-versus-layer diagnostic ("N of M river cells of layer k have their bed outside the layer's elevation interval: B below the layer bottom, A above the layer top", non-blocking), and theDpz_mintable rules.rubem/hydrological_processes/_modflow.py,ModflowGroundwater): inputs read once with the repository readers; DIS -> BAS -> BCF -> wetting -> solver, stress packages set once (the extension keeps them across periods); one RUBEM step = one stress period of the days of the month (updateDISParameter); recharge in mm/step to m/day, RCH option 3;mf.run(<output>/modflow)so nothing lands in the working directory; getters only for package layers with cells; non-convergence raisesRuntimeErrornaming the period andpcrmf.lst; heads, storage, drain flow and river exchange keyed by the user's layer numbers. The extension records LAYCON by call order, so BCF is set bottom-up (proved with the realmf2005on a three-layer grid).bfwis the aquifer-to-river RIV leakage in mm (the saturated-zone reservoir is no longer updated), diagnostics go through the standard raster writer (mfh<n>,mfst<n>,mfdrn<n>,mfaq2rv,mfrv2aq,mfrvnet,mfwt,mfgwd,mfzr,mfzfrac), so the formats apply and the calibrator's raster-free evaluation writes none; the run directory<output>/modflowis created at initialization (refused when it already exists) and removed after the last step, kept when a step raised. The optional root-depth coupling (prototype formulation, one-step lag, vegetation water stress on the root fraction) is off by default.--bound/--fix:modflow.layers.<n>.specific_yield,modflow.layers.<n>.specific_storage(only the storage the layer's LAYCON reads, transient runs),modflow.layers.<n>.kh.<class>(numeric rows of the layer's lookup table),modflow.river.<i>.conductance(numeric entries); bounds are mandatory, finite, positive,specific_yieldwithin (0, 1]. Per-evaluation kh tables are written in the evaluation directory; the calibrated configuration and<config>-calibrated-kh<n>.tblcarry the best values;evaluations.csv,result.jsonandbest_<variable>.csvinclude them.rubem calibratewarns whendis.nstpis above 1 (see below).doc/source/groundwater.rst(exchange equations, layer convention, every key with units, rules, outputs, run directory, root-depth coupling, limitations, references), theMODFLOWkeys and template in the user guide, "MODFLOW parameters" incalibration.rst, changelog entry, README bullet. A test validates every JSON example of the groundwater page against the schema.Resolution of the five review observations on the prototype (the extension numbers layers bottom-up, the observations were written top-down, hence the new convention):
initial_head; the exercise configuration putshead3.mapon layer 1 (top) andhead1.mapon layer 3 (base) as asked, and the validation reports heads below their layer bottom.top_model > botton3 > botton2 > bottonis enforced cell by cell on active columns; layer 1 bottom =botton3, layer 3 bottom =botton.layers: [1]and the diagnostic reports the 383 of 2277 Batalha river cells whose bed lies outside layer 1 (315 below the layer bottom, 284 of them in layer 2's interval and 31 in layer 3's, 68 above the model top), the "83 %" figure.layers: [1, 2, 3]and the same maps; the total boundary conductance per column is then three times the map value, documented, not split.laytype3 on the lower ones; the configuration refuses a listed layer whose LAYCON is not 1 or 3.Deviations from the approved plan, each recorded in the ledger with its reason:
dis.nstpdefaults to 1 (the plan said 5, the scientist's value, still available explicitly). With several time steps per period a solver failure at an earlier time step makes the extension end the whole process insiderun()("Can not open head value result file"), beforeconverged()can be asked; with one time step the failure raises theRuntimeErrorthe calibrator records as a failed evaluation ranked last.rubem calibratewarns onnstp > 1instead of refusing it.<output>/modflowdirectory is refused instead of reused (the run removes the directory at the end).map(the prototype'ssource_boundary_layer/multiplieralternative is gone);laytypeis required; numeric properties are strict floats; the legacyto_dict()writesMODFLOW: nulllikeTABLES.lai_max.[3,5]) are not calibratable (no name); the catalog names only the storage the layer's LAYCON reads..codespellrcignores "botton" (the scientist's file names used in tests and docs), so a real misspelling of "bottom" inrubem/ortests/now passes the blocking spelling step.-888,-999,-999.9,-999.99,-9999) in an active cell are a blocking problem (added after the Batalha exercise showed the marker propagating into the head outputs).Related Issue
Motivation and context
main, predate the Pydantic configuration layer,rubem.apiandrubem calibrate, write MODFLOW files into the working directory, bypass the output formats and import the extension unconditionally.How has this been tested
python -m pytest --ignore=tests/integration/doc -n 2 -q -p no:cacheprovider: 2238 passed, 2 skipped (the byte-exact golden test, CI-only, and the ipojuca dataset smoke) on this branch (Linux, Python 3.13, PCRaster 4.4.2, GDAL 3.11). New tests:test_modflow_settings.py(148),test_modflow_inputs.py(60),test_module_modflow.py(66, mocked extension),test_dynamic_model_modflow.py(22),test_modflow_coupling.py(44 with the realmf2005: the prototype's analytic RIV/GHB/DRN/RCH cases, a three-layer run asserting the BCF LAYCON line, full coupled runs, default-path identity, non-convergence in a child process),test_modflow_parameters.py(34),test_deps.py(+21),test_modflow_documentation.py(4), plus additions to the calibration and CLI tests (a real tiny search withmodflow.layers.1.specific_yieldsearched andmodflow.layers.1.kh.2fixed).uvx ruff@0.16.4 check .,uvx ruff@0.16.4 format --check .,codespell rubem tests: clean. Sphinx was not run locally (not installed); the CI docs job is the check.calibracao/HKs): validation reports the riverbed diagnostic above and heads below the layer bottom (5136 cells ofhead3.mapon layer 1, 5041 of them the-888sentinel, 277 on layer 2, 18 on layer 3). With the sentinel rule added at the end of this branch the 5041 sentinels are a blocking problem, so the exercise runs with--allow-blocking-problems. The run as configured (wetting withwet.map= 1 everywhere on layer 1,nstp5) does not converge in the first period (1334 iterations, 803609 dry and 797756 wet cell conversions, budget discrepancy 43 %). Sensitivity probes (two steps,nstp1):nstp1,--allow-blocking-problems)HKs,KXvertical, wetting on layer 1)KXas horizontal andKYas vertical (the scientist's base JSON)head1,HK1,KX1on layer 1)bfwup to 1228 mm (river cells), aquifer-to-river 1.19e4 m3/day, heads 448-686 m on layers 1-2KXas horizontal andKYas verticalKXhorizontal)multiplier: -1.0, rewetting from below only)Every variant with the BCF wetting enabled fails in the first period and every variant without it converges in about 20 s, whatever the suffix or conductivity convention: the wetting configuration together with the 5041 top-layer cells that start dry is what stops the solver, and the layer conventions do not decide convergence. The exercise can proceed today with wetting off, as a modelling choice for the scientist to confirm; the other questions remain.
The calibration exercise and the
datasetsmoke test are therefore not in this PR; they follow once the scientist settles the questions below.Left for the reviewer and for the scientist (LINAMARIAOSORIO):
botton= base), butKY1..3(KH 0.5/0.1) and theHK1..3class maps put the Marília-rich map at suffix 1 and Adamantina only at suffix 3, andhead1 > head2 > head3on average, which reads top-down.KX<n>.mapcarry the VANIS values ofvariacoes_parametrosMODFLOW.xlsx(4.81, 5.19): anisotropy ratios (vertical K = KH / VANIS) or vertical conductivities in m/d?-888values inside the active domain (5041 cells ofhead3.map, 277 ofhead2.map, 18 ofhead1.map; the validation now blocks them as no-data markers): dry cells, absent layer or no-data?top_model(they differ by -143 to +93 m)?wet.mapis 1 on every active cell (positive WETDRY: cells rewet from the side and from below with a 1 m threshold); the prototype'smultiplier: -1.0was never applied when a map was given. The probes above show the run converging without wetting.hclose 5.0andrclose 3.0are the prototype's values.Screenshots