.. DO NOT EDIT. .. THIS FILE WAS AUTOMATICALLY GENERATED BY SPHINX-GALLERY. .. TO MAKE CHANGES, EDIT THE SOURCE PYTHON FILE: .. "auto_examples/Ex_ERT_workflow.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_ERT_workflow.py: Ex. ERT Workflow: From Hydrological Models to ERT responses and Inversion ==================================================================== This example demonstrates the complete workflow for integrating hydrological model outputs with ERT forward modeling and inversion using PyHydroGeophysX. The workflow includes: 1. Loading MODFLOW water contents data 2. Setting up 2D profile interpolation from 3D model data 3. Creating meshes with geological layer structure 4. Converting water content to resistivity using petrophysical relationships 5. Forward modeling to create synthetic ERT data 6. Performing ERT inversion to recover resistivity models This example serves as a comprehensive tutorial showing the integration of hydrological and geophysical modeling for watershed monitoring applications. .. GENERATED FROM PYTHON SOURCE LINES 22-55 .. code-block:: default import os import sys import numpy as np import matplotlib.pyplot as plt import pygimli as pg from pygimli.physics import ert from pygimli.physics import TravelTimeManager import pygimli.physics.traveltime as tt from mpl_toolkits.axes_grid1 import make_axes_locatable # Setup package path for development try: # For regular Python scripts current_dir = os.path.dirname(os.path.abspath(__file__)) except NameError: # For Jupyter notebooks current_dir = os.getcwd() # Add the parent directory to Python path parent_dir = os.path.dirname(current_dir) if parent_dir not in sys.path: sys.path.append(parent_dir) # Import PyHydroGeophysX modules from PyHydroGeophysX.core.interpolation import ProfileInterpolator, create_surface_lines from PyHydroGeophysX.core.mesh_utils import MeshCreator from PyHydroGeophysX.petrophysics.resistivity_models import WS_Model from PyHydroGeophysX.petrophysics.velocity_models import HertzMindlinModel, DEMModel from PyHydroGeophysX.forward.ert_forward import ERTForwardModeling from PyHydroGeophysX.inversion.ert_inversion import ERTInversion from PyHydroGeophysX.Hydro_modular import hydro_to_ert .. GENERATED FROM PYTHON SOURCE LINES 56-59 .. code-block:: default output_dir = "results/workflow_example" os.makedirs(output_dir, exist_ok=True) .. GENERATED FROM PYTHON SOURCE LINES 60-61 # Step by Step approach .. GENERATED FROM PYTHON SOURCE LINES 63-64 ## Loading domain information... .. GENERATED FROM PYTHON SOURCE LINES 66-67 These would be your actual data files .. GENERATED FROM PYTHON SOURCE LINES 67-76 .. code-block:: default data_dir = "data/" modflow_dir = os.path.join(data_dir, "modflow") # Load domain information from files # (Replace with your actual file paths) idomain = np.loadtxt(os.path.join(data_dir, "id.txt")) top = np.loadtxt(os.path.join(data_dir, "top.txt")) porosity = np.load(os.path.join(data_dir, "Porosity.npy")) .. GENERATED FROM PYTHON SOURCE LINES 77-78 ## Loading MODFLOW water content data.. .. GENERATED FROM PYTHON SOURCE LINES 80-81 Step 2: Exmaple of loading MODFLOW water content data .. GENERATED FROM PYTHON SOURCE LINES 81-94 .. code-block:: default # Note that to save the loading time, we only use a low resoluation model load for the example # In a real-world application, you would load the full resolution data # here we will load the npy file for the water content to save time # Load the water content from a .npy file for demonstration purposes Water_Content = np.load(os.path.join(data_dir, "Watercontent.npy")) water_content = Water_Content[5] # in the publication, we used the 50th time step data print(water_content.shape) .. GENERATED FROM PYTHON SOURCE LINES 95-96 ## Set up profile for 2D section .. GENERATED FROM PYTHON SOURCE LINES 98-99 Step 3: Set up profile for 2D section .. GENERATED FROM PYTHON SOURCE LINES 99-117 .. code-block:: default print("Step 3: Setting up profile...") # Define profile endpoints point1 = [115, 70] # [col, row] point2 = [95, 180] # [col, row] # Initialize profile interpolator interpolator = ProfileInterpolator( point1=point1, point2=point2, surface_data=top, origin_x=569156.0, origin_y=4842444.0, pixel_width=1.0, pixel_height=-1.0, num_points = 400 ) .. GENERATED FROM PYTHON SOURCE LINES 118-119 ## Interpolating data to profile .. GENERATED FROM PYTHON SOURCE LINES 121-122 Step 4: Interpolate data to profile .. GENERATED FROM PYTHON SOURCE LINES 122-129 .. code-block:: default # Interpolate water content to profile water_content_profile = interpolator.interpolate_3d_data(water_content) # Interpolate porosity to profile porosity_profile = interpolator.interpolate_3d_data(porosity) .. GENERATED FROM PYTHON SOURCE LINES 130-131 ## Creating mesh .. GENERATED FROM PYTHON SOURCE LINES 133-164 .. code-block:: default # Load structure layers bot = np.load(os.path.join(data_dir, "bot.npy")) # Process layers to get structure structure = interpolator.interpolate_layer_data([top] + bot.tolist()) # Create surface lines # Indicate the layer for the structure regolith, fractured bedrock and fresh bedrock top_idx=int(0) mid_idx=int(4) bot_idx=int(13) surface, line1, line2 = create_surface_lines( L_profile=interpolator.L_profile, structure=structure, top_idx=top_idx, mid_idx=mid_idx, bot_idx=bot_idx ) # Create mesh mesh_creator = MeshCreator(quality=32) mesh, geom = mesh_creator.create_from_layers( surface=surface, layers=[line1, line2], bottom_depth= np.min(line2[:,1])-10 #50.0 ) # Save mesh mesh.save(os.path.join(output_dir, "mesh.bms")) .. GENERATED FROM PYTHON SOURCE LINES 165-166 Visualize the result .. GENERATED FROM PYTHON SOURCE LINES 166-189 .. code-block:: default import matplotlib.pyplot as plt plt.figure(figsize=(15, 5)) top[idomain==0] = np.nan # Mask out the inactive cells in the top layer # Plot the surface and profile line plt.subplot(121) plt.imshow(top) plt.colorbar(label='Top Elevation (m)') plt.plot(point1[0], point1[1], 'ro', label='Start') plt.plot(point2[0], point2[1], 'bo', label='End') plt.plot([point1[0], point2[0]], [point1[1], point2[1]], 'r--') plt.legend() plt.title('Surface Elevation with Profile Line') # Plot the profile coordinates plt.subplot(122) plt.plot(surface[:,0], surface[:,1]) plt.title('Elevation Along Profile') plt.xlabel('Distance Along Profile') plt.ylabel('Elevation') plt.tight_layout() plt.show() .. GENERATED FROM PYTHON SOURCE LINES 190-203 Profile Setup and Topographic Analysis ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ Left panel shows the 3D MODFLOW domain with the selected 2D profile line (red dashed) connecting start (red) and end (blue) points. Right panel displays the extracted elevation profile that will guide mesh generation and geological structure definition. .. image:: /auto_examples/images/Ex_ERT_workflow_fig_01.png :align: center :width: 800px %% [markdown] ## Interpolating data to mesh .. GENERATED FROM PYTHON SOURCE LINES 205-206 Step 6: Interpolate data to mesh .. GENERATED FROM PYTHON SOURCE LINES 206-266 .. code-block:: default ID1 = porosity_profile.copy() ID1[:mid_idx] = 0 #regolith ID1[mid_idx:bot_idx] = 3 # fractured bedrock ID1[bot_idx:] = 2 # fresh bedrock # Get mesh centers and markers mesh_centers = np.array(mesh.cellCenters()) mesh_markers = np.array(mesh.cellMarkers()) # marker labels for the layers marker_labels = [0, 3, 2] # top. mid, bottom layers (example values) # assign the values for the regolith and fractured bedrock ID1 = porosity_profile.copy() ID1[:mid_idx] = marker_labels[0] #regolith ID1[mid_idx:] = marker_labels[1] # fractured bedrock # Get mesh centers and markers mesh_centers = np.array(mesh.cellCenters()) mesh_markers = np.array(mesh.cellMarkers()) porosity_mesh = np.zeros_like(mesh_markers, dtype=float) wc_mesh = np.zeros_like(mesh_markers, dtype=float) for i in range(len(marker_labels)): if i == 2: # assign the values for the fresh bedrock, since fresh bedrock is no flux boundary and assume the homogeneous as the last layer ID1[bot_idx] = marker_labels[2] # fractured bedrock # fresh bedrock # Interpolate porosity to mesh porosity_temp = interpolator.interpolate_to_mesh( property_values=porosity_profile, depth_values=structure, mesh_x=mesh_centers[:, 0], mesh_y=mesh_centers[:, 1], mesh_markers=mesh_markers, ID=ID1, # Use ID1 to indicate the layers for interpolation layer_markers = [marker_labels[i]], ) porosity_mesh[mesh_markers==marker_labels[i]] = porosity_temp[mesh_markers==marker_labels[i]] # Interpolate water content to mesh wc_temp = interpolator.interpolate_to_mesh( property_values=water_content_profile, depth_values=structure, mesh_x=mesh_centers[:, 0], mesh_y=mesh_centers[:, 1], mesh_markers=mesh_markers, ID=ID1, # Use ID1 to indicate the layers for interpolation layer_markers = [marker_labels[i]], ) wc_mesh[mesh_markers==marker_labels[i]] = wc_temp[mesh_markers==marker_labels[i]] saturation = wc_mesh / porosity_mesh .. GENERATED FROM PYTHON SOURCE LINES 267-338 .. code-block:: default import matplotlib.pyplot as plt import matplotlib as mpl import numpy as np import pygimli as pg from palettable.cartocolors.diverging import Earth_7 # Font settings for publication mpl.rcParams['font.family'] = 'Arial' mpl.rcParams['font.size'] = 12 mpl.rcParams['axes.labelsize'] = 14 mpl.rcParams['axes.titlesize'] = 14 mpl.rcParams['xtick.labelsize'] = 12 mpl.rcParams['ytick.labelsize'] = 12 mpl.rcParams['legend.fontsize'] = 12 mpl.rcParams['figure.dpi'] = 150 # Preprocessing top_masked = np.copy(top) top_masked[idomain == 0] = np.nan saturation = wc_mesh / porosity_mesh ctcolor = Earth_7.mpl_colormap # Create 2x2 figure fig, axs = plt.subplots(2, 2, figsize=(14, 10)) # --- Top Left: Surface elevation map --- im0 = axs[0, 0].imshow(top_masked, origin='lower', cmap='terrain') axs[0, 0].invert_yaxis() # Plot profile line and points axs[0, 0].plot(point1[0], point1[1], 'ro', label='Start') axs[0, 0].plot(point2[0], point2[1], 'bo', label='End') axs[0, 0].plot([point1[0], point2[0]], [point1[1], point2[1]], 'r--') # Remove ticks and axis borders axs[0, 0].set_xticks([]) axs[0, 0].set_yticks([]) for spine in axs[0, 0].spines.values(): spine.set_visible(False) # Title and colorbar cbar0 = fig.colorbar(im0, ax=axs[0, 0], orientation='vertical', shrink=0.8) cbar0.set_label('Elevation (m)') axs[0, 0].legend(loc='upper right') # --- Top Right: Elevation profile --- axs[0, 1].plot(surface[:, 0], surface[:, 1], color='darkgreen') axs[0, 1].set_xlabel('Distance (m)') axs[0, 1].set_ylabel('Elevation (m)') axs[0, 1].grid(True) # --- Bottom Left: Porosity mesh --- pg.show(mesh, porosity_mesh, ax=axs[1, 0], orientation="vertical", cMap=ctcolor, cMin=0.05, cMax=0.45, xlabel="Distance (m)", ylabel="Elevation (m)", label='Porosity (-)', showColorBar=True) # --- Bottom Right: Saturation mesh --- pg.show(mesh, saturation, ax=axs[1, 1], orientation="vertical", cMap='Blues', cMin=0, cMax=1, xlabel="Distance (m)", ylabel="Elevation (m)", label='Saturation (-)', showColorBar=True) # Layout adjustment plt.tight_layout(pad=3) plt.savefig(os.path.join(output_dir, "topography_and_properties.tiff"), dpi=300) .. GENERATED FROM PYTHON SOURCE LINES 339-351 Mesh Properties and Interpolation Results ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ Four-panel analysis showing: (1) surface elevation map with profile line, (2) elevation profile, (3) interpolated porosity field with three-layer structure, and (4) saturation distribution. The mesh preserves geological boundaries while providing appropriate resolution for ERT modeling. .. image:: /auto_examples/images/Ex_ERT_workflow_fig_02.png :align: center :width: 900px %% .. GENERATED FROM PYTHON SOURCE LINES 351-355 .. code-block:: default print("Water Content min/max:", np.min(wc_mesh), np.max(wc_mesh)) print("Saturation min/max:", np.min(saturation), np.max(saturation)) .. GENERATED FROM PYTHON SOURCE LINES 356-357 ## Calculating saturation .. GENERATED FROM PYTHON SOURCE LINES 359-360 Ensure porosity is not zero to avoid division by zero .. GENERATED FROM PYTHON SOURCE LINES 360-363 .. code-block:: default saturation = np.clip(wc_mesh / porosity_mesh, 0.0, 1.0) .. GENERATED FROM PYTHON SOURCE LINES 364-365 ## Converting to resistivity .. GENERATED FROM PYTHON SOURCE LINES 367-368 Convert to resistivity using petrophysical model .. GENERATED FROM PYTHON SOURCE LINES 368-413 .. code-block:: default marker_labels = [0, 3, 2] # top. mid, bottom layers (example values) sigma_w = [1/20,1/20,1/20] # Conductivity of water in S/m (example value) m = [1.2, 1.8, 1.9] # Cementation exponent for each layer (example values) n = [2.2, 1.8, 2.5] # Saturation exponent for each layer (example values) sigma_s = [1/100, 0, 0] # Saturated resistivity of the surface conductivity see Chen & Niu, (2022) for each layer (example values) # Convert water content back to resistivity res_models = np.zeros_like(wc_mesh) # Initialize an array for resistivity values mask = (mesh_markers == marker_labels[0]) top_res = WS_Model( saturation[mask], # Water content values for this layer porosity_mesh[mask], # Porosity values for this layer sigma_w[0], # Conductivity of water in S/m m[0], # Cementation exponent for this layer n[0], # Saturation exponent for this layer sigma_s[0] # Use a scalar value ) res_models[mask] = top_res mask = (mesh_markers == marker_labels[1]) mid_res = WS_Model( saturation[mask], # Water content values for this layer porosity_mesh[mask], # Porosity values for this layer sigma_w[1], # Conductivity of water in S/m m[1], # Cementation exponent for this layer n[1], # Saturation exponent for this layer sigma_s[1] # Use a scalar value ) res_models[mask] = mid_res mask = (mesh_markers == marker_labels[2]) bot_res = WS_Model( saturation[mask], # Water content values for this layer porosity_mesh[mask], # Porosity values for this layer sigma_w[2], # Conductivity of water in S/m m[2], # Cementation exponent for this layer n[2], # Saturation exponent for this layer sigma_s[2] # Use a scalar value ) res_models[mask] = bot_res .. GENERATED FROM PYTHON SOURCE LINES 414-488 .. code-block:: default import matplotlib.pyplot as plt import numpy as np from scipy.interpolate import griddata plt.figure(figsize=(12, 6)) # First subplot - narrower left side (spans both rows, 1 column) plt.subplot2grid((2, 3), (0, 0), rowspan=2, colspan=1) x1, y1 = np.mgrid[20:21:1j, 1568:1613:50j] res = griddata(mesh_centers[:,0:2], res_models, (x1, y1), method='linear') po = griddata(mesh_centers[:,0:2], porosity_mesh, (x1, y1), method='linear') sat = griddata(mesh_centers[:,0:2], saturation, (x1, y1), method='linear') wc = griddata(mesh_centers[:,0:2], wc_mesh, (x1, y1), method='linear') # Plotting the first interpolation plt.scatter(res.ravel(), y1.ravel(),c='k') plt.xlabel('Resistivity (Ω·m)') plt.ylabel('Elevation (m)') plt.xscale('log') plt.xlim([100,10000]) # Second subplot - larger top right position (spans 2 columns) plt.subplot2grid((2, 3), (0, 1), colspan=2) # First interpolation for the main plot x1, y1 = np.mgrid[20:21:1j, 1609.2:1613:8j] res = griddata(mesh_centers[:,0:2], res_models, (x1, y1), method='linear') po = griddata(mesh_centers[:,0:2], porosity_mesh, (x1, y1), method='linear') sat = griddata(mesh_centers[:,0:2], saturation, (x1, y1), method='linear') wc = griddata(mesh_centers[:,0:2], wc_mesh, (x1, y1), method='linear') # Plotting the first interpolation plt.scatter(po.ravel(), res.ravel(), c=sat, vmin=0.3, vmax=0.35, cmap='Blues',label='Regolith') plt.colorbar(label='Saturation (-)') plt.xlabel('Porosity (-)') plt.ylabel('Resistivity (Ω·m)') res0 = WS_Model( np.ones(100)*np.mean(sat), # Water content values for this layer np.linspace(np.nanmin(po.ravel()), np.nanmax(po.ravel()), 100), # Porosity values for this layer sigma_w[0], # Conductivity of water in S/m m[0], # Cementation exponent for this layer n[0], # Saturation exponent for this layer sigma_s[0] # Use a scalar value ) plt.plot(np.linspace(np.nanmin(po.ravel()), np.nanmax(po.ravel()), 100), res0, 'k--', lw=2,label = "Waxman-Smits model") plt.ylim([300,600]) plt.legend(frameon=False) # Third subplot - larger bottom right position (spans 2 columns) plt.subplot2grid((2, 3), (1, 1), colspan=2) # First interpolation for the main plot x1, y1 = np.mgrid[20:21:1j, 1578:1609:32j] res = griddata(mesh_centers[:,0:2], res_models, (x1, y1), method='linear') po = griddata(mesh_centers[:,0:2], porosity_mesh, (x1, y1), method='linear') sat = griddata(mesh_centers[:,0:2], saturation, (x1, y1), method='linear') wc = griddata(mesh_centers[:,0:2], wc_mesh, (x1, y1), method='linear') # Plotting the first interpolation plt.scatter(po.ravel(), res.ravel(), c=sat, cmap='Blues', vmin=0.4, vmax=1,label='Fractured Bedrock') plt.colorbar(label='Saturation (-)') plt.xlabel('Porosity (-)') plt.ylabel('Resistivity (Ω·m)') res0 = WS_Model( np.ones(100)*np.nanmean(sat), # Water content values for this layer np.linspace(np.nanmin(po.ravel()), np.nanmax(po.ravel()), 100), # Porosity values for this layer sigma_w[1], # Conductivity of water in S/m m[1], # Cementation exponent for this layer n[1], # Saturation exponent for this layer sigma_s[1] # Use a scalar value ) plt.plot(np.linspace(np.nanmin(po.ravel()), np.nanmax(po.ravel()), 100), res0, 'k', lw=2,label = "Archie's Law") plt.legend(frameon=False) plt.ylim([100,2500]) # plt.yscale('log') plt.tight_layout() plt.savefig(os.path.join(output_dir, "resistivity_porosity_saturation.tiff"), dpi=300) .. GENERATED FROM PYTHON SOURCE LINES 489-502 Petrophysical Relationship Validation ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ Multi-panel analysis validates rock physics models: vertical resistivity profile (left), Waxman-Smits model for regolith (top right), and Archie's law for fractured bedrock (bottom right). Scattered data points match theoretical curves, confirming realistic petrophysical transformations. .. image:: /auto_examples/images/Ex_ERT_workflow_fig_03.png :align: center :width: 800px %% [markdown] ## ERT forward modeling simulation .. GENERATED FROM PYTHON SOURCE LINES 504-532 .. code-block:: default xpos = np.linspace(15,15+72 - 1,72) ypos = np.interp(xpos,interpolator.L_profile,interpolator.surface_profile) pos = np.hstack((xpos.reshape(-1,1),ypos.reshape(-1,1))) schemeert = ert.createData(elecs=pos,schemeName='wa') # Step 10: Forward modeling to create synthetic ERT data mesh.setCellMarkers(np.ones(mesh.cellCount())*2) grid = pg.meshtools.appendTriangleBoundary(mesh, marker=1, xbound=100, ybound=100) fwd_operator = ERTForwardModeling(mesh=grid, data=schemeert) synth_data = schemeert.copy() fob = ert.ERTModelling() fob.setData(schemeert) fob.setMesh(grid) dr = fob.response(res_models) dr *= 1. + pg.randn(dr.size()) * 0.05 ert_manager = ert.ERTManager(synth_data) synth_data['rhoa'] = dr synth_data['err'] = ert_manager.estimateError(synth_data, absoluteUError=0.0, relativeError=0.05) #ert.showData(synth_data, logscale=True) synth_data.save(os.path.join(output_dir, "synthetic_data.dat")) .. GENERATED FROM PYTHON SOURCE LINES 533-560 .. code-block:: default from palettable.lightbartlein.diverging import BlueDarkRed18_18 import matplotlib.pyplot as plt fixed_cmap = BlueDarkRed18_18.mpl_colormap # --- Create figure with 1 row, 2 columns --- fig, axs = plt.subplots(1, 2, figsize=(12, 4)) # --- Left subplot (121): Resistivity with log scale --- pg.show(mesh, res_models, ax=axs[0], orientation="horizontal", cMap=fixed_cmap, logScale=True, showColorBar=True, xlabel="Distance (m)", ylabel="Elevation (m)", label='Resistivity (Ω·m)', cMin=100, cMax=2000) pg.viewer.mpl.drawSensors(axs[0], schemeert.sensors(), diam=1, facecolor='black', edgecolor='black') # --- Right subplot (122): ERT data --- ax_ert, cbar_default = ert.showData(synth_data, logscale=False, ax=axs[1], cMin=200, cMax=1000, cMap='jet', label='Apparent Resistivity (Ω·m)' ) axs[1].set_xlabel("Distance (m)") axs[1].spines['top'].set_visible(False) axs[1].spines['right'].set_visible(False) # Optional: Adjust spacing between subplots plt.tight_layout() plt.show() plt.savefig(os.path.join(output_dir, "res_model_and_synth_data.tiff"), dpi=300) .. GENERATED FROM PYTHON SOURCE LINES 561-574 Forward Modeling Results ~~~~~~~~~~~~~~~~~~~~~~~~ Left panel shows the true resistivity model with electrode positions (black circles) along the surface. Right panel displays synthetic ERT measurements showing characteristic pseudosection patterns that reflect the three-layer geological structure in the apparent resistivity data. .. image:: /auto_examples/images/Ex_ERT_workflow_fig_04.png :align: center :width: 800px %% Step 11: Run ERT inversion on synthetic data .. GENERATED FROM PYTHON SOURCE LINES 574-589 .. code-block:: default ## using my code to the inversion print("Step 11: Running ERT inversion...") # Create ERT inversion object inversion = ERTInversion( data_file=os.path.join(output_dir, "synthetic_data.dat"), lambda_val=10.0, method="cgls", use_gpu=True, max_iterations=10, lambda_rate= 1.0 ) inversion_result = inversion.run() .. GENERATED FROM PYTHON SOURCE LINES 590-591 # Using Pygimili default to the inversion .. GENERATED FROM PYTHON SOURCE LINES 591-594 .. code-block:: default mgr = ert.ERTManager(os.path.join(output_dir, "synthetic_data.dat")) inv = mgr.invert(lam=10, verbose=True,quality=34) .. GENERATED FROM PYTHON SOURCE LINES 595-615 .. code-block:: default fig, axes = plt.subplots(1, 3, figsize=(10, 12)) # True resistivity model ax1 = axes[0] cbar1 = pg.show(mesh, res_models, ax=ax1, cMap='jet', logScale=False, cMin=100, cMax=3000, label='Resistivity [Ohm-m]') ax1.set_title("True Resistivity Model") # Inverted model ax2 = axes[1] cbar2 = pg.show(inversion_result.mesh, inversion_result.final_model, ax=ax2, cMap='jet', logScale=False, cMin=100, cMax=3000, label='Resistivity [Ohm-m]',coverage=inversion_result.coverage>-1) ax2.set_title("Inverted Resistivity Model (Our Code)") ax3 = axes[2] cbar2 = pg.show(mgr.paraDomain, mgr.paraModel(), ax=ax3, cMap='jet', logScale=False, cMin=100, cMax=3000, label='Resistivity [Ohm-m]',coverage=mgr.coverage()>-1) ax3.set_title("Inverted Resistivity Model (Pygimli)") # Adjust layout plt.tight_layout() .. GENERATED FROM PYTHON SOURCE LINES 616-627 Inversion Results Comparison ~~~~~~~~~~~~~~~~~~~~~~~~~~~ Three-panel comparison showing: (1) true resistivity model from petrophysical conversion, (2) PyHydroGeophysX inversion result, and (3) PyGIMLi standard inversion. Both methods successfully recover the main geological features, validating the complete hydro-to-geophysical workflow. .. image:: /auto_examples/images/Ex_ERT_workflow_fig_05.png :align: center :width: 900px .. GENERATED FROM PYTHON SOURCE LINES 627-632 .. code-block:: default # The inversion results are almost same from this code and Pygimli default inversion. # the difference is that the chi2 value for stop inversion is not the same, we chose 1.5 while Pygimli is 1.0 .. GENERATED FROM PYTHON SOURCE LINES 633-634 # One step approach .. GENERATED FROM PYTHON SOURCE LINES 636-637 ## ERT one step from HM to GM .. GENERATED FROM PYTHON SOURCE LINES 639-640 Set up directories .. GENERATED FROM PYTHON SOURCE LINES 640-724 .. code-block:: default output_dir = "results/hydro_to_ert_example" os.makedirs(output_dir, exist_ok=True) # Load your data data_dir = "data/" idomain = np.loadtxt(os.path.join(data_dir, "id.txt")) top = np.loadtxt(os.path.join(data_dir, "top.txt")) porosity = np.load(os.path.join(data_dir, "Porosity.npy")) water_content = np.load(os.path.join(data_dir, "Watercontent.npy"))[5] # Time step 50 # Set up profile point1 = [115, 70] point2 = [95, 180] interpolator = ProfileInterpolator( point1=point1, point2=point2, surface_data=top, origin_x=569156.0, origin_y=4842444.0, pixel_width=1.0, pixel_height=-1.0, num_points=400 ) # Create mesh structure bot = np.load(os.path.join(data_dir, "bot.npy")) layer_idx = [0, 4, 12] # Example indices for top, middle, and bottom layers structure = interpolator.interpolate_layer_data([top] + bot.tolist()) surface, line1, line2 = create_surface_lines( L_profile=interpolator.L_profile, structure=structure, top_idx=layer_idx[0], mid_idx=layer_idx[1], bot_idx=layer_idx[2] ) # Create mesh mesh_creator = MeshCreator(quality=32) mesh, geom = mesh_creator.create_from_layers( surface=surface, layers=[line1, line2], bottom_depth=np.min(line2[:,1])-10 ) # Define layer markers marker_labels = [0, 3, 2] # top, middle, bottom layers # Define resistivity parameters for each layer rho_parameters = { 'rho_sat': [100, 500, 2400], # Saturated resistivity for each layer (Ohm-m) 'n': [2.2, 1.8, 2.5], # Cementation exponent for each layer 'sigma_s': [1/500, 0, 0] # Surface conductivity for each layer (S/m) } mesh_markers = np.array(mesh.cellMarkers()) # Generate ERT response directly synth_data, res_model = hydro_to_ert( water_content=water_content, porosity=porosity, mesh=mesh, mesh_markers = mesh_markers, profile_interpolator=interpolator, layer_idx=layer_idx, structure = structure, marker_labels=marker_labels, rho_parameters=rho_parameters, electrode_spacing=1.0, electrode_start=15, num_electrodes=72, scheme_name='wa', noise_level=0.05, abs_error=0.0, rel_error=0.05, save_path=os.path.join(output_dir, "synthetic_ert_data.dat"), verbose=True, seed=42, ) ert.showData(synth_data, logscale=True) .. GENERATED FROM PYTHON SOURCE LINES 725-736 One-Step Integrated Workflow Results ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ The streamlined hydro-to-ERT workflow produces synthetic apparent resistivity data directly from MODFLOW outputs. This integrated approach demonstrates automated processing from hydrological model results to geophysical measurements suitable for inversion and interpretation. .. image:: /auto_examples/images/Ex_ERT_workflow_fig_06.png :align: center :width: 700px .. GENERATED FROM PYTHON SOURCE LINES 738-744 Summary ~~~~~~~ This comprehensive workflow demonstrated the complete integration from hydrological models to ERT responses and inversion. Key achievements include realistic petrophysical transformations, successful forward modeling, and validated inversion results suitable for watershed monitoring applications. .. rst-class:: sphx-glr-timing **Total running time of the script:** (0 minutes 0.000 seconds) .. _sphx_glr_download_auto_examples_Ex_ERT_workflow.py: .. only:: html .. container:: sphx-glr-footer sphx-glr-footer-example .. container:: sphx-glr-download sphx-glr-download-python :download:`Download Python source code: Ex_ERT_workflow.py ` .. container:: sphx-glr-download sphx-glr-download-jupyter :download:`Download Jupyter notebook: Ex_ERT_workflow.ipynb `