.. DO NOT EDIT. .. THIS FILE WAS AUTOMATICALLY GENERATED BY SPHINX-GALLERY. .. TO MAKE CHANGES, EDIT THE SOURCE PYTHON FILE: .. "auto_examples/Ex_posterior_uncertainty.py" .. LINE NUMBERS ARE GIVEN BELOW. .. only:: html .. note:: :class: sphx-glr-download-link-note :ref:`Go to the end ` to download the full example code. .. rst-class:: sphx-glr-example-title .. _sphx_glr_auto_examples_Ex_posterior_uncertainty.py: Ex. Posterior Uncertainty of Water Content Recovered from ERT ============================================================== Estimate how uncertain the water content read from an ERT inversion is, then check that estimate against the MODFLOW state that produced the survey. As in Ex_sensitivity_analysis, the water content of the Treeline catchment MODFLOW model on the Ex_ERT_workflow transect becomes a 72-electrode Wenner survey through ``hydro_to_ert``, and ``ERTInversion`` inverts it. The inversion's objective is then read as a Gaussian posterior of log resistivity: its data weights give the data covariance, and its regularization, smoothness with a weak pull toward the reference model, gives the prior covariance. ``linearized_posterior`` returns the posterior covariance at the recovered model, and ``propagate_petro_uncertainty`` samples it through the Waxman-Smits model of each unit into water-content intervals. The inversion runs at two strengths of that prior, to show how much the intervals depend on it. Unlike Ex_MC_Hydro, the petrophysical parameters stay at the values that produced the survey, so the intervals carry the survey's uncertainty alone. Requires PyHydroGeophysX with pyGIMLi. Run time: about 2 minutes. .. GENERATED FROM PYTHON SOURCE LINES 21-22 .. code-block:: Python :dedent: 1 .. GENERATED FROM PYTHON SOURCE LINES 24-45 .. code-block:: Python import os import numpy as np import matplotlib.pyplot as plt import pygimli as pg try: current_dir = os.path.dirname(os.path.abspath(__file__)) except NameError: current_dir = os.getcwd() if os.path.isdir(os.path.join(current_dir, "examples")): current_dir = os.path.join(current_dir, "examples") from PyHydroGeophysX.core.interpolation import ProfileInterpolator, create_surface_lines from PyHydroGeophysX.core.mesh_utils import MeshCreator from PyHydroGeophysX.forward.ert_forward import ertforandjac2 from PyHydroGeophysX.Hydro_modular import hydro_to_ert from PyHydroGeophysX.inversion.ert_inversion import ERTInversion from PyHydroGeophysX.inversion.ert_mesh import build_inversion_mesh from PyHydroGeophysX.petrophysics import resistivity_to_water_content from PyHydroGeophysX.uncertainty import linearized_posterior, propagate_petro_uncertainty .. GENERATED FROM PYTHON SOURCE LINES 46-52 MODFLOW state, transect and predicted survey -------------------------------------------- The same steps as Ex_sensitivity_analysis: water content, porosity and layer bottoms from the MODFLOW arrays, the transect of Ex_ERT_workflow, and the survey this state predicts, with 5% noise. Layers 1-4 are regolith, 5-12 fractured bedrock and 13-14 fresh bedrock. .. GENERATED FROM PYTHON SOURCE LINES 54-78 .. code-block:: Python data_dir = os.path.join(current_dir, "data") top = np.loadtxt(os.path.join(data_dir, "top.txt")) bottoms = np.load(os.path.join(data_dir, "bot.npy")) porosity = np.load(os.path.join(data_dir, "Porosity.npy")) water_content = np.asarray(np.load(os.path.join(data_dir, "Watercontent.npy"), mmap_mode="r")[5]) interpolator = ProfileInterpolator(point1=[115, 70], point2=[95, 180], surface_data=top, origin_x=569156.0, origin_y=4842444.0) structure = interpolator.interpolate_layer_data([top] + bottoms.tolist()) layer_idx = [0, 4, 12] surface, regolith_base, fractured_base = create_surface_lines( interpolator.L_profile, structure, *layer_idx) mesh, _ = MeshCreator(quality=32).create_from_layers( surface=surface, layers=[regolith_base, fractured_base], bottom_elevation=fractured_base[:, 1].min() - 10.) units = [0, 3, 2] # mesh markers of regolith, fractured bedrock and fresh bedrock unit_names = {0: "Regolith", 3: "Fractured bedrock", 2: "Fresh bedrock"} petrophysics = {"rho_sat": [100., 500., 2400.], "n": [2.2, 1.8, 2.5], "sigma_s": [1 / 500, 0., 0.]} survey, _ = hydro_to_ert( water_content=water_content, porosity=porosity, mesh=mesh, profile_interpolator=interpolator, layer_idx=layer_idx, structure=structure, marker_labels=units, rho_parameters=petrophysics, electrode_spacing=1., electrode_start=15., num_electrodes=72, scheme_name="wa", noise_level=.05, rel_error=.05, seed=2026) .. GENERATED FROM PYTHON SOURCE LINES 79-87 Invert, then read each inversion as a posterior ----------------------------------------------- With a smallness weight alpha the regularization is a proper Gaussian prior on log resistivity, with precision lambda (Wm^T Wm + alpha I) for the smoothness matrix Wm, so lambda sets how strong the prior is. The inversion's data weights 1/log(1+e), for relative error e, are inverse standard deviations of log apparent resistivity; the data covariance holds the squares of their reciprocals. lambda=3 fits the survey to its noise; lambda=10 smooths more. .. GENERATED FROM PYTHON SOURCE LINES 89-108 .. code-block:: Python inversion_mesh = build_inversion_mesh(survey, mesh_quality=34) def posterior_for(lam, alpha=.01): inversion = ERTInversion(survey, mesh=inversion_mesh, lambda_val=lam, reference_weight=alpha, max_iterations=10, verbose=False) result = inversion.run() log_resistivity = np.log(result.final_model) _, J = ertforandjac2(inversion.fwd_operator, log_resistivity, inversion.mesh) Wm = inversion.Wm_r.toarray() prior_cov = np.linalg.inv(lam * (Wm.T @ Wm + alpha * np.eye(Wm.shape[1]))) return {"cells": inversion.fwd_operator.paraDomain, "chi2": result.meta["chi2"], "log_resistivity": log_resistivity, "prior_std": np.sqrt(np.diag(prior_cov)), "posterior_cov": linearized_posterior(J, 1. / inversion.Wdert_diag ** 2, prior_cov)} posteriors = {lam: posterior_for(lam) for lam in (3., 10.)} cells = posteriors[3.]["cells"] .. GENERATED FROM PYTHON SOURCE LINES 109-113 MODFLOW water content and porosity on the inversion cells --------------------------------------------------------- Each inversion cell belongs to the MODFLOW unit its centre lies in, and takes the profile's water content and porosity interpolated within that unit. .. GENERATED FROM PYTHON SOURCE LINES 115-126 .. code-block:: Python x, z = np.array(cells.cellCenters())[:, :2].T unit = np.where(z > np.interp(x, *regolith_base.T), 0, np.where(z > np.interp(x, *fractured_base.T), 3, 2)) profile_water_content = interpolator.interpolate_3d_data(water_content) layer_unit = np.zeros_like(profile_water_content) layer_unit[layer_idx[1]:layer_idx[2]], layer_unit[layer_idx[2]:] = 3, 2 modflow_theta, cell_porosity = ( interpolator.interpolate_to_mesh(interpolator.interpolate_3d_data(values), structure, x, z, unit, layer_unit, units) for values in (water_content, porosity)) .. GENERATED FROM PYTHON SOURCE LINES 127-132 Propagate each posterior into water content ------------------------------------------- Each sample of log resistivity is converted with the Waxman-Smits parameters of its cell's unit. Saturation is clipped to [0, 1], so a sample beyond full saturation returns the porosity. .. GENERATED FROM PYTHON SOURCE LINES 134-151 .. code-block:: Python def to_water_content(log_rho): theta = np.empty_like(log_rho) for i, marker in enumerate(units): sel = unit == marker theta[sel] = resistivity_to_water_content( np.exp(log_rho[sel]), petrophysics["rho_sat"][i], petrophysics["n"][i], cell_porosity[sel], petrophysics["sigma_s"][i]) return theta for posterior in posteriors.values(): estimate = propagate_petro_uncertainty(posterior["log_resistivity"], posterior["posterior_cov"], to_water_content, n_samples=1000, seed=2026) posterior["mean"], posterior["std"] = estimate["mean"], estimate["std"] posterior["low"], posterior["high"] = np.quantile(estimate["samples"], [.025, .975], axis=0) .. GENERATED FROM PYTHON SOURCE LINES 152-157 Check the intervals against MODFLOW ----------------------------------- Only a synthetic test can do this. Coverage is the share of each unit, by area, where the MODFLOW water content lies inside the 95% interval; the standard-deviation ratio is posterior over prior, in log resistivity. .. GENERATED FROM PYTHON SOURCE LINES 159-175 .. code-block:: Python area = np.array([cell.size() for cell in cells.cells()]) print(f"{'lambda':>6s} {'Unit':18s}{'chi2':>6s}{'Std ratio':>11s}{'RMSE':>8s}" f"{'95% width':>11s}{'Coverage':>10s}") for lam, posterior in posteriors.items(): posterior_std = np.sqrt(np.diag(posterior["posterior_cov"])) inside = (modflow_theta >= posterior["low"]) & (modflow_theta <= posterior["high"]) for marker in units: sel = unit == marker if sel.any(): weights = area[sel] / area[sel].sum() rmse = np.sqrt(weights @ (posterior["mean"][sel] - modflow_theta[sel]) ** 2) print(f"{lam:6.0f} {unit_names[marker]:18s}{posterior['chi2']:6.2f}" f"{np.median(posterior_std[sel] / posterior['prior_std'][sel]):11.2f}{rmse:8.3f}" f"{np.median((posterior['high'] - posterior['low'])[sel]):11.3f}" f"{weights @ inside[sel]:10.0%}") .. GENERATED FROM PYTHON SOURCE LINES 176-180 Save the posterior and plot the interpretation ---------------------------------------------- The maps are for lambda=3. The column on the right follows the dotted line at 50 m: the MODFLOW layers, and the 95% intervals at both prior strengths. .. GENERATED FROM PYTHON SOURCE LINES 182-233 .. code-block:: Python output_dir = os.path.join(current_dir, "results", "posterior_uncertainty") os.makedirs(output_dir, exist_ok=True) chosen = posteriors[3.] np.savez_compressed(os.path.join(output_dir, "posterior.npz"), cell_centers=np.array(cells.cellCenters()), unit=unit, modflow_water_content=modflow_theta, posterior_cov=chosen["posterior_cov"], **{f"{key}_lambda{lam:g}": value for lam, posterior in posteriors.items() for key, value in posterior.items() if key not in ("cells", "posterior_cov")}) chosen_std = np.sqrt(np.diag(chosen["posterior_cov"])) fig = plt.figure(figsize=(16, 9)) grid = fig.add_gridspec(2, 3, width_ratios=[1, 1, .55], left=.05, right=.98, top=.96, bottom=.04, wspace=.3, hspace=.3) maps = [ (grid[0, 0], modflow_theta, dict(cMap="Blues", cMin=.05, cMax=.35, label="m³/m³"), "MODFLOW water content"), (grid[0, 1], chosen["mean"], dict(cMap="Blues", cMin=.05, cMax=.35, label="m³/m³"), "Posterior mean water content"), (grid[1, 0], chosen_std / chosen["prior_std"], dict(cMap="viridis", cMin=.3, cMax=1., label="Posterior / prior"), "Standard deviation of log resistivity"), (grid[1, 1], (chosen["mean"] - modflow_theta) / chosen["std"], dict(cMap="RdBu", cMin=-3, cMax=3, label="(mean - MODFLOW) / std"), "Error in posterior standard deviations"), ] for position, values, style, title in maps: ax = fig.add_subplot(position) pg.show(cells, values, ax=ax, orientation="horizontal", **style) for line in (regolith_base, fractured_base): ax.plot(*line.T, "k--", lw=.8) ax.axvline(50., color="k", lw=.6, ls=":") ax.set(title=title, xlim=(10, 92), ylim=(1585, 1626), xlabel="Distance (m)", ylabel="Elevation (m)") column = np.argmin(np.abs(interpolator.L_profile - 50.)) elevation = np.linspace(structure[0, column] - .05, z.min(), 300) hits = [cells.findCell(pg.Pos(50., e)) for e in elevation] keep = np.array([cell is not None for cell in hits]) ids = np.array([cell.id() for cell in hits if cell is not None]) ax = fig.add_subplot(grid[:, 2]) ax.stairs(profile_water_content[::-1, column], structure[::-1, column], orientation="horizontal", baseline=None, color="k", lw=1.5, label="MODFLOW layers") for lam, color in [(3., "C0"), (10., "C3")]: posterior = posteriors[lam] ax.fill_betweenx(elevation[keep], posterior["low"][ids], posterior["high"][ids], color=color, alpha=.25, label=f"95% interval, λ = {lam:g}") ax.plot(posterior["mean"][ids], elevation[keep], color=color, lw=1) ax.set(ylim=(z.min(), structure[0, column] + .5), xlabel="Water content (m³/m³)", ylabel="Elevation (m)", title="Column at 50 m") ax.legend(loc="lower right", fontsize=9) fig.savefig(os.path.join(output_dir, "posterior_uncertainty.png"), dpi=150, bbox_inches="tight") plt.show() .. GENERATED FROM PYTHON SOURCE LINES 234-252 The survey reduces the standard deviation of log resistivity most near the surface (bottom left): to a median of 0.55 of its prior value in the regolith, against 0.74 in the fractured bedrock. Through the fixed petrophysics, the lambda=3 posterior gives 95% intervals about 0.09 m³/m³ wide in the regolith and 0.14 in the fractured bedrock, and they contain the MODFLOW water content in 99% of both units. The errors are not random, though (bottom centre): the smooth section underestimates the water content of the saturated bedrock below the water table and errs both ways along the base of the regolith, where it blurs the MODFLOW interfaces. The stronger prior (lambda=10) fits the survey a little worse (chi-squared 1.09 against 0.93) and returns intervals about 40% narrower, which contain the MODFLOW state in only 80-84% of the section (right). Linearized intervals are only as reliable as the regularization they inherit, so choose lambda by the data misfit, and read the intervals as conditional on the petrophysics and the unit boundaries. .. image:: /auto_examples/images/Ex_posterior_uncertainty_fig_01.png :align: center :width: 900px .. _sphx_glr_download_auto_examples_Ex_posterior_uncertainty.py: .. only:: html .. container:: sphx-glr-footer sphx-glr-footer-example .. container:: sphx-glr-download sphx-glr-download-jupyter :download:`Download Jupyter notebook: Ex_posterior_uncertainty.ipynb ` .. container:: sphx-glr-download sphx-glr-download-python :download:`Download Python source code: Ex_posterior_uncertainty.py ` .. container:: sphx-glr-download sphx-glr-download-zip :download:`Download zipped: Ex_posterior_uncertainty.zip `