Gravity Inversion Study Reveals Subsurface Density Architecture at Mexico’s Los Humeros Geothermal Field

Geophysical researchers evaluating the volcanic subsurface of the Los Humeros Geothermal Field in central Mexico have established that gravity anomalies primarily reflect lithological structure and rock compaction rather than direct reservoir temperature or fluid signatures. Utilizing gravimetric observations alongside petrophysical data and a three-dimensional geological framework, the investigation determined that standard forward gravity modeling exhibits systematic negative discrepancies reaching several milligals when compared directly to field observations. The researchers demonstrated that stochastic inverse modeling accounting for depth-dependent compaction and fluid-saturated matrix densities successfully resolves this deficit, bringing computed gravity anomalies into close alignment with observed records to within a misfit of less than one milligal. This finding clarifies how structural density variations govern gravity signals within high-temperature volcanic settings, providing an evidence-based pathway to refine three-dimensional geological representations without relying on speculative temperature-driven gravity interpretations.

The study focused on the Los Humeros Geothermal Field, situated in the northeastern sector of the Trans-Mexican Volcanic Belt, which stands as Mexico’s third-largest producing geothermal reservoir. The complex volcanic succession across the field comprises Pleistocene to Holocene basaltic to rhyolitic sequences deposited above a folded Mesozoic sedimentary basement cut by granitic and syenitic plutons. Stratigraphically, the terrain is grouped into pre-caldera, caldera-forming, and post-caldera units spanning nine distinct geological units dissected by twenty-one modeled tectonic faults. The exploited reservoir is largely hosted within pre-caldera lavas, notably the andesites and basaltic andesites of the Teziutlán Formation dating from 2.6 to 1.4 million years ago, alongside older lavas of the Cuyoaco Formation. Overlying these lavas are extensive explosive volcanic deposits, prominently the 164-thousand-year-old Xáltipan ignimbrite, which reaches thicknesses up to 780 meters and local porosities up to 70 percent, the 70-thousand-year-old Faby Tuff, and the 69-thousand-year-old Zaragoza ignimbrite associated with the nested Los Potreros caldera collapse structure.

Geothermal fluid circulation within the field occurs under both liquid-dominated and superheated vapor-dominated regimes, with deeper reservoir temperatures surpassing 300 degrees Celsius across superheated sectors. Despite these thermal extremes, reservoir fluids maintain an exceptionally low dissolved chemical concentration, with total dissolved solids typically hovering around 1 gram per liter and rarely exceeding 4 grams per liter. This dilute, predominantly meteoric composition ensures that chemical salinity yields a negligible effect on subsurface fluid density calculations. However, interpreting regional and local gravity anomalies remains inherently ambiguous due to the non-uniqueness of potential field solutions, where multiple subsurface mass configurations can produce identical gravimetric measurements at the surface. To resolve these ambiguities, geophysicists must establish rigorous parameterization protocols that test laboratory-derived dry rock densities against realistic in-situ conditions incorporating water saturation, overburden compaction, and complex lithological proportions across three-dimensional volumes.

Field gravity data collection across the northern and central regions of the geothermal concession encompassed 263 gravity stations distributed across ten east-west profiles, spaced 500 meters apart with 200-meter station intervals. Observations were collected using a Scintrex CG-5 Autograv meter providing an instrumental precision of one microgal, with spatial positions recorded via differential GPS achieving five-millimeter vertical accuracy. Following rigorous quality filtering that discarded unreliable stations and applied standard deviation thresholds of 0.05 milligals for gravimetry and 0.2 meters for positioning, 247 high-quality stations were retained for structural modeling. Data processing incorporated tidal corrections, instrumental drift tracking via repeated daily base-station observations, latitude adjustments, and detailed complete Bouguer corrections using a reduction density of 2,670 kilograms per cubic meter. Terrain corrections integrated digital elevation models ranging in resolution from 1 meter near measurement stations to 167 kilometers regionally, yielding a cumulative mean measurement uncertainty of approximately 0.7 milligals.

To evaluate how density parameterization affects gravity predictions, the researchers formulated three distinct starting models within a voxelized domain derived from an existing local geological framework. The domain spans a grid of 17,388 voxels with cell dimensions of 200 meters horizontally and 100 meters vertically down to sea level. The first baseline parameterization assigned single dry bulk density values to each geological unit based on laboratory core plug measurements and dry outcrop analogues, supplemented by bounded least-squares density optimization. The second parameterization implemented a single saturated bulk density per unit, corrected for depth-dependent compaction using Athy’s exponential porosity relationship parameterized with unit-specific compaction coefficients derived from seismic compressional and shear wave velocities. The third parameterization applied a multi-modal saturated density scheme assigning normal density distributions across up to three lithological sub-fractions within each geological unit to capture internal stratigraphic proportions.

Forward modeling calculations executing the gravimetric response equations of Plouff, Okabe, and Holstein revealed that all three initial parameterizations produced substantial discrepancies when compared against observed residual gravity anomalies ranging between minus 12 and positive 2 milligals. Computed gravity signatures exhibited a pronounced, systematic negative bias, spanning values from minus 50 to minus 16 milligals across the various configurations. While the compaction-corrected saturated density configuration generated a misfit histogram that shifted from left-skewed toward a more Gaussian distribution, forward models consistently underestimated total subsurface mass by several milligals. This persistent underestimation confirmed that forward modeling alone, even when incorporating optimized rock matrix densities and compaction corrections, cannot capture the full three-dimensional mass distribution without localized stochastic adjustments to density values and structural interfaces.

To overcome forward modeling deficiencies, the authors performed stochastic inverse modeling using a Markov Chain Monte Carlo algorithm within GeoModeller software, systematically perturbing model voxels across 173,880 iterations to ensure each grid cell underwent at least ten perturbations. Inversion calculations systematically relaxed voxel properties subject to commonality, shape ratio, and stratigraphic vertical order constraints. When density and structural geometries were allowed to adjust with equal fifty-percent probabilities, the compaction-corrected saturated parameterization achieved dramatic misfit reductions, narrowing residual gravity discrepancies to a tight Gaussian distribution with a half-width of plus or minus 0.2 milligals. In contrast, applying the complex multi-modal density distribution failed to improve the inversion fit and introduced substantial positive skewness, an outcome attributed to uncertainties in estimating precise volumetric fractions from sparse borehole sampling across heterogeneous volcanic units.

To determine whether model convergence relied more heavily on modifying lithological boundaries or updating internal rock densities, the team executed sensitivity inversions by varying the relative probabilities of structural versus density adjustments in twenty-percent increments across forty perturbations per voxel cell. Evaluating the resulting misfit distributions revealed that prioritizing density modifications generated significantly superior convergence compared to altering geometry. A configuration assigning a seventy-percent probability to density adjustments and a thirty-percent probability to structural alterations achieved the optimal result, characterized by the lowest standard deviation below 0.3 milligals, minimal skewness, and low kurtosis. This outcome demonstrates that the initial three-dimensional geological model possesses robust geometric fidelity, indicating that misfits between predicted and observed fields stem primarily from petrophysical property variability within units rather than misplaced structural contacts or incorrect fault dip angles.

Detailed examination of the recovered inversion densities revealed that initial input values had systematically underestimated optimal in-situ rock densities across the geothermal field by 60 to 130 kilograms per cubic meter. The largest post-inversion density increases occurred within the uppermost volcanic sequences, specifically Units 1 and 4, which exhibited shifts of up to 130 kilograms per cubic meter, whereas deeper pre-caldera lavas and basement rocks displayed only minor property adjustments. The resulting density histograms showed distinct log-normal distributions within the upper volcanic units, reflecting high structural heterogeneity, localized welding variations within ignimbrites, and hydrothermal alteration. Conversely, the deeper pre-caldera formations and Mesozoic basement rocks exhibited symmetrical Gaussian distributions with skewness values below 0.1, confirming that deep volcanic and metamorphic basements maintain substantially higher lithological homogeneity across the reservoir footprint.

The inverted three-dimensional density model further revealed localized vertical columns of elevated density aligned directly along major fault corridors, matching regions of residual misfit observed near primary structural conduits. Depth sensitivity tests extending the inversion domain downward to 7,000 meters below sea level confirmed that calculated gravity responses remained virtually unchanged below roughly 3,000 meters depth, establishing that gravimetric anomalies at Los Humeros are predominantly controlled by density contrasts in shallow to intermediate formations above sea level. Crucially, the gravity signatures did not differentiate between conventional production sectors and superheated, vapor-dominated zones of the field. The investigators concluded that gravity variations across this active caldera complex reflect volcanic lithology, structural faulting, and mechanical compaction rather than in-situ thermal anomalies or reservoir fluid distributions, providing a dependable quantitative framework for validating complex volcanic structural models.

Subject of Research: Density parameterization and 3-D gravity inversion in the Los Humeros geothermal field

Article Title: In-situ density and litho-constrained modeling of gravity data in the Los Humeros geothermal field, Mexico

Article References: Cornejo-Triviño, N., Leandra, M., Weydt, Sass, I., & Schill, E. (2026). In-situ density and litho-constrained modeling of gravity data in the Los Humeros geothermal field, Mexico. Geothermal Energy, 14(1), Article 19. https://doi.org/10.1186/s40517-026-00395-4

Image Credits: AI Generated

DOI: 10.1186/s40517-026-00395-4

Keywords: geothermal energy, gravity inversion, Los Humeros, volcanic reservoir, density modeling, petrophysics, geophysics, structural geology, In-situ, density, litho-constrained, modeling