.. _help-calibration: Calibration & validation ======================== A numerical model contains parameters that cannot be measured directly, most importantly the roughness of the riverbed. **Calibration** adjusts these parameters until the model reproduces field measurements. **Validation** then tests the calibrated model against measurements that were not used for the calibration. aXqua performs both with `HydroBayesCal `_. Two terms are used throughout this section. **Calibration parameters** are the uncertain inputs of the model that are adjusted, for example the roughness of a zone. **Calibration targets** are the measured quantities that the model is compared with, for example the water depth at surveyed points. .. _calibration-introduction: Introduction to Bayesian calibration ------------------------------------ A calibration by trial and error ends with one set of parameter values and no statement about how reliable these values are. Bayesian calibration answers a more useful question: which parameter values are probable, given the measurements and their errors? The user states for each calibration parameter the range in which its true value is expected. This range is the **prior** distribution. Bayes' theorem combines the prior with the measurements. Parameter values with which the model reproduces the measurements become more probable, and all others become less probable. The result is the **posterior** distribution of each parameter. A narrow posterior indicates that the measurements determine the parameter well. A posterior that is as wide as the prior indicates that the measurements contain no information about this parameter. The direct evaluation of Bayes' theorem requires many thousands of model runs, which is not feasible for a model that runs for hours. HydroBayesCal therefore replaces the model by a **surrogate model**, in the following steps: #. The numerical model is run for a limited number of parameter combinations that are spread over the prior ranges (30 by default). #. A Gaussian process emulator is trained with these runs. The emulator is a statistical model that predicts the result of the numerical model for any parameter combination within a fraction of a second, together with the uncertainty of its own prediction (Rasmussen and Williams 2006). #. The emulator is evaluated for thousands of parameter combinations to compute the posterior distribution. #. **Bayesian active learning** selects the parameter combination at which one additional run of the numerical model improves the emulator the most. The model is run for this combination, and the emulator is trained again (Oladyshkin et al. 2020). #. Steps 3 and 4 are repeated until the selected number of runs is reached (50 by default). In this way, a few dozen model runs are sufficient where a direct Bayesian calibration would require thousands. Applications to reservoir models are documented by Mouris et al. (2023) and Schwindt et al. (2023). .. _calibration-ground-truth: Ground truth data setup ----------------------- Field measurements are called ground truth. aXqua accepts them in a workbook with a fixed structure, which is called the calibration target template. Create it for a case with: .. code-block:: text axqua targets The command writes the workbook ``calibration-target-data.xlsx`` into ``user-sources/ground-truth/``. Its sheets are: .. list-table:: :header-rows: 1 :widths: 24 76 * - Sheet - Content * - ``hydraulics`` - One row per measurement point: velocity components, their fluctuations, water depth and bed elevation. The velocity magnitude and the turbulent kinetic energy are computed by formulas in the sheet. * - ``morphodynamics`` - Grain sizes (d16 to d90) and the fraction of fine sediment per sample. The measured bed change is filled in automatically from the DEM of difference (:ref:`case-morphodynamics`). * - ``parameters`` - The calibration parameters with their lower and upper limits, selected from a list of available TELEMAC and GAIA parameters. Each row carries an ID. The positions of the measurements are read from a point layer that contains the same IDs, so that coordinates do not have to be typed in. Enter the workbook and the point layers in the case file: .. code-block:: yaml ground_truth: targets: file: user-sources/ground-truth/calibration-target-data.xlsx hydraulics_positions: user-sources/geodata/measurement-points.gpkg sediment_positions: user-sources/geodata/sediment-samples.gpkg join_key: ID Then repeat the preprocessing, which compiles the measurements into the tables that the calibration reads. The IDs have to be unique within a sheet. Measurements of several survey days whose numbering restarts at zero therefore need distinct IDs. Two alternatives to the workbook exist. ``ground_truth.measurements`` names a table that the user compiled, with the coordinates x, y and z in the first three columns and one column per measured quantity. ``ground_truth.sources`` lists original instrument files together with the point layer of their positions, which aXqua compiles into such a table. **Measurements in one layer or table.** If the measured values and their positions already exist in one file, this file is sufficient. Enter it as a source of the kind ``points``: .. code-block:: yaml ground_truth: sources: - category: hydraulics kind: points positions: user-sources/ground-truth/velocity-measurements.gpkg The file is a point layer (GeoPackage, shapefile or GeoJSON) in any coordinate reference system, or a table (CSV file or Excel workbook) with the columns ``x`` and ``y`` in the coordinate reference system of the case. aXqua recognizes the following column names, in upper or lower case and with or without a unit in square brackets or parentheses, for example ``v(x) [m/s]`` or ``v_x (m/s)``: .. list-table:: :header-rows: 1 :widths: 44 56 * - Quantity - Column names * - Coordinates of a table - ``x`` and ``y``, or ``easting`` and ``northing`` * - Velocity components in m/s - ``u``, ``v``, ``w``, or ``v(x)``, ``v(y)``, ``v(z)``, or ``v_x``, ``v_y``, ``v_z``, or ``VelX``, ``VelY``, ``VelZ`` * - Total water depth in m - ``h``, ``depth``, ``water_depth`` or ``Total Depth`` * - Depth of a reading below the water surface in m - ``MeasD`` or ``Meas. Depth`` The velocity magnitude that is compared with the model is computed from the horizontal components. If a layer contains several readings at one position, aXqua uses the reading that is closest to 0.6 of the water depth below the surface, which approximates the depth-averaged velocity. Points outside the model outline are left out of the calibration. In the plugin, the files are entered in the case settings (*Open case settings...* on the *Case Setup* tab), page *Calibration and validation*, in the table *Measurement files* of the frame *Measurements*. *Add file...* adds a point layer or a table, and the column *Use for* assigns each file to the calibration or to the validation. A file for the validation needs the name of its flow situation and its discharge (:ref:`validation situation `). Velocity measurements with a SonTek FlowTracker2 are transferred into the ``hydraulics`` sheet with the script ``extract_flowtracker.py``, which aXqua places next to the workbook. Note that the velocity fluctuation is the standard deviation of the velocity samples. It is not the column ``VxErr`` of the instrument export, which is the standard error of the mean velocity. .. _calibration-setup: Bayesian calibration & validation setup --------------------------------------- The ``calibration`` block of the case file defines what is calibrated against which measurements: .. code-block:: yaml calibration: calibration_quantities: ["WATER DEPTH"] extraction_quantities: ["WATER DEPTH", "SCALAR VELOCITY"] measurement_error: 0.10 # used where no measured error is available init_runs: 30 # runs for the first training of the emulator max_runs: 50 # total number of runs parameters: - {name: zone1, min: 0.05, max: 0.50} # ks of roughness zone 1 [m] - {name: zone2, min: 0.10, max: 1.00} # ks of roughness zone 2 [m] ``calibration_quantities`` are the calibration targets. The names are those of the TELEMAC result variables: ``WATER DEPTH``, ``SCALAR VELOCITY``, ``TURBULENT ENERG.`` (requires the k-epsilon turbulence model) and ``CUMUL BED EVOL`` (requires GAIA). ``parameters`` lists the calibration parameters with the limits of their prior ranges. The name ``zone`` refers to the roughness of zone N. A name that starts with ``gaia`` refers to a GAIA keyword, for example ``gaiaCLASSES SHIELDS PARAMETERS 1`` for the critical Shields parameter of the first grain size class. Parameters that are selected in the ``parameters`` sheet of the workbook are added to this list. Select wide but physically plausible limits. If the posterior of a parameter accumulates at one of its limits, the true value probably lies outside the range, and the calibration has to be repeated with a wider range. .. _validation-setup: **Validation.** A validation examines whether the calibrated parameters describe the reach, or only the day on which the calibration data were measured. It therefore requires measurements that the calibration has not used and that stem from **another flow situation**: another discharge, measured in another survey. All measurements of the ``ground_truth`` block are calibration data. aXqua does not withhold a part of them for a validation. The measurements of one survey are taken within a few hours in one flow field, and neighboring points are therefore not independent of each other. A model that predicts some of these points from the others (leave-one-out cross-validation, k-fold cross-validation or a random share) demonstrates that it interpolates within that flow field. It does not demonstrate that its roughness values hold at another discharge, at which other parts of the bed are wetted and other roughness elements are submerged. The validation data are defined in the case file as a **situation**, with its measurements and its boundary values: .. code-block:: yaml calibration: validation: - name: September 2025 inflows: {1: 0.2, 3: 5.1} # discharge [m3/s] per inflow of the model sources: - category: hydraulics kind: points positions: user-sources/ground-truth/september.gpkg .. list-table:: :header-rows: 1 :widths: 26 74 * - Entry - Meaning * - ``name`` - A short name of the situation. It names the files of the validation. * - ``sources`` - The measurements, in the form of the sources of the ``ground_truth`` block. * - ``prescribed_flowrate`` - The total discharge of the situation in m³/s. For a model with several inflows, the discharge is distributed among them in the proportions of the calibrated case. * - ``inflows`` - The discharge of each inflow in m³/s, for a reach whose branches carried another share of the flow. The inflows are identified by the number of their liquid boundary. ``axqua validation `` lists the inflows of the model with the discharge that each carries in the calibrated case. * - ``prescribed_elevation`` - The water level at the outflow in m, for a case with ``outflow_condition: elevation``. With a rating curve, aXqua derives the water level from the discharge of the situation. * - ``duration`` - The simulated time in seconds, if it differs from that of the case. **In the plugin.** The tab *Calibration & validation* has the box *Validation with data of another flow situation*. Select the point layer with the validation data with *Browse...*, or with *From QGIS...* from the point layers that are open in QGIS. Enter the discharge of each inflow: the box lists the inflows of the model with the discharge of the calibrated case. *Save* writes the situation into the case file. The same situations appear in the case settings (*Open case settings...*, page *Calibration and validation*): in the table *Measurement files*, each file whose column *Use for* says *Validation* belongs to the situation named in the column *Situation*, with the discharge of that situation. The discharge belongs to the situation, so a change in one row changes it in every row of that situation, and a path in the columns *File* and *Positions* can be corrected in place. Settings of a situation that the table does not show, such as the discharge of each inflow, are kept. .. _calibration-run: Run HydroBayesCal ----------------- The calibration requires the additional package HydroBayesCal (:doc:`../installation/qgis-plugin`) and a built case with a completed steady simulation. Each run of the calibration repeats this steady simulation with another set of parameter values. In the plugin, click *Submit* on the *Calibration & validation* tab. The option *Prepare only* writes all input files without starting the runs, so that they can be inspected first. In a terminal: .. code-block:: text axqua submit --kind calibration --option prepare_only=true axqua submit --kind calibration .. figure:: ../img/gui-calibration.png :alt: The Calibration & validation tab of the plugin :width: 70% The tab *Calibration & validation* for a case with measurements. (1) The option *Prepare only*. (2) *Submit* starts the calibration. (3) The measurements of the validation situation, selected as a file or from the layers that are open in QGIS, and below them the discharge of each inflow in that situation. (4) *Save* stores the validation situation in the case file, and *Validate* runs it. The computing time is approximately the number of runs (``max_runs``) times the duration of one steady simulation. The job list shows the number of the current run, for example ``iter 5/12``. HydroBayesCal performs its runs in the folder of the built case (``axqua-case/simulation/``) and writes each tested set of parameter values into the friction table and the steering file there. aXqua restores both files when the calibration ends, so that a later simulation uses the case as it was built. For the same reason, jobs of one case run one after the other: a job that is submitted during a calibration waits until the calibration has ended. The folder ``cases/example-isar/`` of the repository contains a complete example with a run time of about half an hour, including a step-by-step guide for the plugin. .. _validation-run: **Validation run.** After the calibration, click *Validate* in the box *Validation with data of another flow situation*, or submit the job in a terminal: .. code-block:: text axqua submit --kind validation The validation runs the calibrated model with the boundary values of the validation situation. Its steering file is derived from that of the built case and differs only in the prescribed discharges and water levels, in the names of its own result files, and in the calibrated parameter values. The mesh and all other settings are those of the calibrated model, so that a difference between model and measurement cannot stem from a model that was built differently. The calibrated values are the most probable combination of the posterior (the joint posterior maximum) of the last calibration. The files of the built case and the result of the steady simulation are not changed. aXqua then reads the model at the validation points and writes three files into ``axqua-case/calibration-validation/validation/``: a table with the measured and the modeled value at each point (``validation-.csv``), a summary (``validation-.json``), and a figure of the modeled against the measured values (``validation-.png``). Measurement points outside the model are counted and excluded. **Several discharges.** Roughness values that were calibrated at one discharge are not necessarily valid at another one. The script ``run_Bayes_cal_multiflow.py`` in the case folder calibrates one common set of parameters against measurements at several discharges. Run it with ``--smoke`` first. This short test with three runs verifies the complete chain before days of computing time are spent. **OpenFOAM.** An OpenFOAM model can be calibrated against measured velocity components with the script ``run_Bayes_cal_openfoam.py``. Because a calibration requires dozens of runs, they are performed with a coarse mesh in the mode ``rigid-lid``. The script ``openfoam_verify_posterior.py`` then repeats the calibrated case at full resolution in the mode ``vof``. The roughness of an OpenFOAM model is calibrated as one value for the entire bed. It is therefore not comparable with the roughness values per zone of a TELEMAC calibration. .. _calibration-quality: Quality analysis ---------------- HydroBayesCal writes its results into the folder ``auto-saved-results-HydroBayesCal`` within ``axqua-case/calibration-validation/``. Assess a calibration in four steps: #. **Did the active learning converge?** HydroBayesCal records two indicators after every iteration. The Bayesian model evidence (BME) expresses how well the model reproduces the measurements on average over the parameter ranges. The relative entropy (RE) expresses how much information the measurements added to the prior. Both should level off toward the end of the calibration. If they still change strongly, increase ``max_runs``. #. **Is the posterior narrower than the prior?** Compare the two distributions of each parameter. A posterior that is clearly narrower than the prior and lies inside the range identifies the parameter. A posterior at a limit of the range calls for a wider range. A posterior that resembles the prior means that the measurements do not constrain this parameter. #. **Does the calibrated model reproduce the calibration targets?** The calibrated parameter set is the most probable combination of the posterior. HydroBayesCal compares the model results for this set with the measurements. Check whether the remaining differences are of the size of the measurement errors, and whether they are distributed randomly or concentrated in one part of the reach. A systematic deviation indicates an error in the model, for example in the terrain data or the boundary conditions, that no parameter value can compensate. #. **Does the calibrated model reproduce independent data?** Run the validation (:ref:`validation run `). The differences are typically somewhat larger than for the calibration targets. Much larger differences indicate that the calibration has compensated errors of the model structure through the parameters. The validation reports the following statistics for each quantity, for example the flow velocity and the water depth: .. list-table:: :header-rows: 1 :widths: 26 74 * - Statistic - Meaning * - Number of points - The validation points on the mesh of the model. Points outside the model, and points that are dry in the model although water was measured, are reported separately. * - Bias - The mean of the modeled minus the measured values. A negative bias of the velocity means that the model is too slow on average. The relative bias refers to the mean measured value. * - Root mean square error (RMSE) - The typical difference between model and measurement at one point. It contains the bias and the scatter. * - Mean absolute error - The mean magnitude of the differences, which is less sensitive to single large differences than the RMSE. * - Correlation - The correlation coefficient between modeled and measured values. A high correlation with a large bias indicates a correct spatial pattern at a wrong level. * - Share within the measurement error - The share of points at which model and measurement differ by no more than the error of the measurement. The error is that of the instrument combined with the relative ``measurement_error`` of the calibration. About 68 % are expected for a model without systematic deviation, and about 95 % within twice the error. The validation also states the remaining difference between inflow and outflow at the end of its run. If this difference exceeds the tolerance of the case, the run had not reached a steady state, and the ``duration`` of the validation situation has to be increased. The `HydroBayesCal documentation `_ describes the result files and figures in detail. References ~~~~~~~~~~ Mouris, K., Acuña Espinoza, E., Schwindt, S., Mohammadi, F., Haun, S., Wieprecht, S., and Oladyshkin, S. (2023). "Stability criteria for Bayesian calibration of reservoir sedimentation models." *Modeling Earth Systems and Environment*, 9, 3643-3661. `doi:10.1007/s40808-023-01712-7 `_ Oladyshkin, S., Mohammadi, F., Kroeker, I., and Nowak, W. (2020). "Bayesian3 active learning for the Gaussian process emulator using information theory." *Entropy*, 22(8), 890. `doi:10.3390/e22080890 `_ Rasmussen, C. E., and Williams, C. K. I. (2006). *Gaussian processes for machine learning*. MIT Press. Schwindt, S., Callau Medrano, S., Mouris, K., Beckers, F., Haun, S., Nowak, W., Wieprecht, S., and Oladyshkin, S. (2023). "Bayesian calibration points to misconceptions in three-dimensional hydrodynamic reservoir modeling." *Water Resources Research*, 59(3), e2022WR033660. `doi:10.1029/2022WR033660 `_