Gravity Data Allow to Image the Shallow‐Medium Subsurface Below Mud Volcanoes

The debate about the conceptual model of mud volcanoes functioning is still alive in the literature. A large part of the literature focuses on the characterization of the deep reservoir where expelled fluids are expected to originate. Another part of literature is focused on the study of the shallow system of mud volcanoes, which could influence the short‐term variations in mud volcanoes activity. We present and analyze a new data set of micro‐gravimetric data to study the area of the Nirano Salse, Italy. Unlike what is commonly assumed for the study area, our results suggest that the geomorphology of the Nirano Salse is not related to a caldera collapse above a shallow mud chamber, but to the surface expression of slip distribution of a fault termination along which the fluids ascended to the surface. We believe that gravimetric data can significantly improve the study of hydrocarbon seeps and mud volcanism.

• A microgravimetric survey coupled with laboratory density measurements allows to image the shallow to medium subsurface below mud volcanoes • Gravimetric data allow to distinguish between morphologies due to caldera collapse and the surface expression of slip along a fault • For the first time the attenuation volume caused by mud intrusion is imaged by gravimetry along a fault

Supporting Information:
Supporting Information may be found in the online version of this article.
is also favored by the presence of faults, which can act as preferential pathways for the ascent of these deep fluids (Mazzini, 2009).
In Italy, mud volcanoes are characterized by the almost continuous, but relatively quiet, expulsion of material (Giambastiani et al., 2022).In the present work we focus on the mud volcanoes located in the northern Apennines Italy (e.g., Borgatti et al., 2019;Capozzi & Picotti, 2002), near the city of Fiorano Modenese (Emilia-Romagna region).Such volcanoes, belonging to the Regional Natural Reserve of Nirano Salse are well known since ancient times and attract thousands of visitors every year.To date, the physical conceptual model of the Nirano Salse ascribes the eruptions to the presence of overpressurized fluids that are expelled from a leaky deep reservoir (Bonini, 2008).
In general, most studies devoted to the characterization of the Nirano Salse focus on the deep reservoir where expelled fluids originate (e.g., Martinelli et al., 2012;Oppo et al., 2017;Sciarra et al., 2019).However, much more important to characterize short term variations in gas-mud emissions activity might be the presence of shallow depth secondary reservoirs.Within them, gas and connate water coming from depth may accumulate and remain at rest until overpressure building up opens the conduits to the surface.According to Giambastiani et al. (2022) and ref. therein, the reservoirs are put into communication with the surface due to the episodical reactivation of pre-existing faults or pipes.In this view, Lupi et al. (2015) and Giambastiani et al. (2022) also underline the importance of the shallow aquifers (tens of meters of depth), which may act as buffers where rising gas (methane) is trapped and temporarily stored.This configuration could lead to the localization of small groups of gryphons, mud pools, and mud cones (Figure 1a) just above the shallow reservoirs.
This view suggests that small scale geophysical surveys could be determinant for the identification of these shallow reservoirs or of a large mud caldera at shallow surface below the whole Nirano area as suggested by Bonini (2008Bonini ( , 2009Bonini ( , 2020)).The implications for unraveling the subsurface structure are important to evaluate potential hazards at different spatial scales connected to the fluid emissions.Short scale resistivity surveys were carried out in the Nirano area (e.g., Accaino et al., 2007;Lupi et al., 2015;Romano et al., 2023) also integrated with geochemical methods (e.g., Martinelli et al., 2012;Oppo et al., 2017;Sciarra et al., 2019).To better understand the conceptual model of fluids, rise and storage as well as the subsurface structural setting of Nirano a new data set of gravimetric data is here presented.The gravimetric data can be valuable in distinguishing the volumes characterized by a high gas/fluid saturation (i.e., low density) from the hosting materials with higher density.Gravity surveys performed in mud volcanoes zones (e.g., Kadirov et al., 2005;Mauri, Husein, Mazzini, Irawan, et al., 2018, Mauri, Husein, Mazzini, Karyono, et al., 2018;Napoli et al., 2020;Odonne et al., 2021;Riffo et al., 2021;Sayyadul, 2005) are useful to constrain the actual location of gas saturated volumes.We show that they can also give important information about the mechanisms that influenced the actual Nirano morphology and provide guidance about hazards of mud-fluid emissions in the area.

Geological Setting
The Nirano Salse are located at the foothills of the Northern Apennines (Figure 1) in an area where Plio-Pleistocene transgressive clays (Argille Azzurre Fm, FAA) blanket the surface down to a depth of 50-150 m.Below this clay unit, a series of Miocene to Oligocene marly and sandy deposits (Termina Fm-sand, Pantano Fm-sand/marl turbidites, Antognola Fm-marls, and Ranzano Fm-sand; Gasperi et al., 2005) form a stratigraphic package of 1,000-1,500 m in thickness.The lithostratigraphic series continues with Cretaceous-Eocene Ligurian clay units constituting the basement of the sequence.The Miocene-Oligocene units are crosscut by a series of high angle faults and fracture zones with a NW-SE and NE-SW orientations (Gasperi et al., 2005).No brittle structures, however, are observed in the upper FAA.The mud volcanic features in the area consist of mud pools, gryphons, and dry seeps (Figures 1a and 1b) (Etiope, 2015), which are located within a geomorphological depression.Bonini (2008Bonini ( , 2009) ) argued a collapse of a shallow mud filled caldera along bordering ring faults and the presence of an anticline just above the depression (Bonini, 2007).Field evidence and the published geologic map of Gasperi et al. (2005), however, do not show an anticline but the hinge of a monocline (Figure 1).More recently, Bonini (2020) and Maestrelli et al. (2019) stressed the importance of far structures and near structures in controlling fluid ascent in the compressional setting of the Northern Apennines defining two separate structural and reservoir systems (deep source and shallow source reservoirs), which would be present also in the Nirano area.
Field mapping, structural planes attitudes measurements, and sampling were also performed during the geophysical field study to integrate the existing information.In particular, rock saturated density was determined on 18  S1 in Supporting Information S1).

Material and Methods
The data set was acquired during a field campaign performed in 2022.The measures were obtained with a LaCoste-Romberg Model D-149 relative gravimeter in 61 different stations (Figure 1c) distributed over an area of about 4 km 2 and centered on the mud volcanoes zone.The reference point of the data set is reported in Figure 1c.
The altitudes of gravity measurement stations were obtained with leveling in the central zone (Leica Sprinter 250M) whereas the farthest points were measured with GPS-RTK (Leica CS15).
The regional Digital Elevation Model (DEM) (Figure S1 in Supporting Information S1) has a resolution of 5 m × 5 m (https://geoportale.regione.emilia-romagna.it)and it was used to compute the terrain correction.The mean value of gravity, g i obs , at each observation station, i, was first corrected in order to account for to the tides contribution and the drift rate of the instrument (e.g., McCubbine et al., 2018;Reilly, 1970).For each location, i, we also consider the latitude correction g i n (λ i ), the free-air correction δg i FA and the Bouguer correction δg i B (see Supporting Information S1 for further details).For the latter, we assume a density of the plate of ρ = 1,850 kg/m 3 , that is the average density of the shallow lithologic units making up the stratigraphic sequence of the area (Table S1 in Supporting Information S1).The free air ∆g i FA and the Bouguer ∆g i B anomalies can be then computed respectively, as following: Due to the complex topography of the study area, the terrain correction for each location, ∆g i T , was also computed.Its computation was performed with the prism method (e.g., Almeida et al., 2018;Nagy, 1966) considering the regional DEM shown in Figure S1 in Supporting Information S1 and compared with another independent method (Hammer, 1939).The complete Bouguer anomaly can be then computed as Relative free air, ∆g i RFA , Bouguer and ∆g i RB anomalies, as well as the relative terrain correction ∆g i

RT
, with respect to the reference location REF (Figure 1c), can be computed respectively as following (5) The resulting ∆g i RFA , ∆g i RT , and ∆g i RB are plotted in Figures 2c-2e, respectively.In the same way, we can compute the relative, complete Bouguer anomaly (Figure 2f) as (7) As we can see from Figure 2d, the terrain correction has a minimum close to the mud volcanoes due to the flat topography of this area.The complete Bouguer anomaly will be considered in the following data inversion process.

Inversion of Gravity Data
The 3D inversion for density contrasts was performed with the GROWTH 3.0 software (Camacho et al., 2021), which allows to determine the subsurface density distribution starting from discrete gravity anomaly observations.The inversion was performed in a 3D domain (32,171 cells), which encompasses the distribution of all measurement stations and a depth of 403 m (a.m.s.l.).We used as input data the complete Bouguer anomaly, ∆g i RCB shown in Figure 2f.The offset O f was estimated by the inversion as O f = −210 μGal.
The resulting complete Bouguer anomaly, after the removal of the offset, is shown in Figure 3a.Anomalies range from about −1.5 to 1.5 mGal in the whole study area.The strongest negative anomaly of about 1.5 mGal is located just below the four mud volcanoes and in the southeast part of the study area.Positive anomalies up to 1.5 mGal are located to the southwest of the mud volcanoes.The trade-off curve used to select the smoothing parameter (Camacho et al., 2021) is shown in Figure S3 of Supporting Information S1.The model sensitivity (Figure S4 in Supporting Information S1) is higher in the vicinity of the area hosting the mud volcanoes, where the distribution of stations is finer.The resolution tests (Figure S5 in Supporting Information S1) indicate that our distribution of the stations is suitable to resolve anomalous density volumes with a size (width and depth) of the order of some hundreds of meters.
The subsurface density anomalies retrieved by the inversion are reported in Figures 3b-3d.Results suggest the presence of a volume (∼2 × 10 8 m 3 ) with an anomalous negative density contrast up to about −350 kg/m 3 , located below the mud volcanoes.Such a volume extends from a few meters below the surface down to about 400 m depth.From the map section plotted at z = 50 m (Figure 3b), the N-S vertical section AAʹ (Figure 3c), and the BBʹ vertical section (Figure 3d) we identify some volumes with positive density contrast up to about +350 kg/m 3 .The inferred positive density contrasts extend from a few meters below the surface (about 200 a.m.s.l.) down to about 0 m (a.m.s.l.) elevation.The residuals between measured and modeled Bouguer gravity anomalies are generally low except for a few stations (Figure S1 in Supporting Information S1).

Discussion
The Bouguer anomaly map (Figure 3a) shows a low-density anomaly (about −1,500 μGal) striking in a NW-SE direction that terminates to the NW where a high-density anomaly (about +1,500 μGal) is observed.The tip of the anomaly termination corresponds to the region where three different groups of gryphons and mud pools are observed at the surface (Figure 3a).Published geological maps and interpretations (Bonini, 2008;Gasperi et al., 2005) report the axis of an anticline at the same location of the low-density anomaly.Field mapping  faults.The latter, that present just south of the Nirano area in the exposed Miocene sequence covered by the transgressive FAA clays, have offsets ranging from a few meters to a few hundred meters (Gasperi et al., 2005).
Even if the resolution of the model used for the inversion of the gravity data does not allow to see the shallow, very small-scale (tens of meters) details of the study area, it has allowed us to reveal the main structures that characterize the study area on larger spatial scales (hundreds of meters).The gravity inversion results (Figures 3b-3d) define a NW-SE oriented low-density zone with ρ ≈ 1,500-1,600 kg/m 3 and a density contrast ∆ρ ≈ −350 kg/m 3 with respect to the average density of the surface clays (1,850 kg/m 3 ) used in the inversion.Such a low-density zone has a width of about 200 m and a length of about 1,300 m (Figure 3b).It extends from beneath the topographic surface (about 200 m.a.s.l.) down to about −300 m.a.s.l.(Figures 3c and 3d).Considering that the average density of the Miocene rocks from 100 m down to a depth of 1,000 m from the surface is around 2,000-2,200 kg/m 3 , the low-density zone may be interpreted as the result of gas and/or mud saturation in a rock with a large amount of matrix or fracture porosity (e.g., Mauri, Husein, Mazzini, Irawan, et al., 2018;Mauri, Husein, Mazzini, Karyono, et al., 2018), or by intrusion of viscous mud-gas slurry dykes.Given that in the Nirano area the geothermal gradient is normal (25 K/km, Viganò et al., 2012), the methane down to a depth of 700-800 m from the surface would be in gaseous form and, eventually, in supercritical state only at the bottom of the low-density zone (McCain, 1990).In these conditions, the low-density zone is compatible with a rock of the Miocene sequence having a gas-saturated porosity in the order 30%-40%.This would be possible in the damage zone of a fault with large matrix and fracture porosity or where some fluid filled mud/gas cracks and feeders intruded the rock.Considering the densities of both mud (ρ = 1,200 kg/m 3 , Giambastiani et al., 2022) and rocks, and the gas content within them, we can easily obtain the inferred density anomalies assuming volumetric ratios (volume of gas/volume of mud and/or rock) in the range 0.55-0.70.By examining the isosurfaces computed for density contrasts |Δρ| = 300 kg/m 3 in the inverted gravity model (Figure 4a), we can infer that: (a) there is a low-density zone (Δρ = −300 kg/m 3 ), which is laying on a vertical NW-SE striking plane and (b) the top of the high density zone (Δρ = +300 kg/m 3 ), which forms a flexure, is structurally higher in the South-SW than in the North-NE area.Although the NE area is constrained by a lower number of stations with respect to the SW area, the presence of a strong positive density anomaly in the SW area is coherent with the high Bouguer anomaly observed there (Figure 2f).
We interpret the top of the high-density zone as the top of the Miocene marly formations that have densities compatible with those obtained in the inversion.These features could be the result of gas/mud intrusion along a NW-SE oriented vertical fault zone whose displacement is diminishing in the NW direction.This conceptual model, which is represented in Figure 4b, is coherent with the structural frame of the area (Gasperi et al., 2005), the density values measured in the rocks of the sedimentary sequence, the surface geology and morphology.The offset of the fault, based on the thickness of the sedimentary sequence, is in the order of 100-150 m with a down-to the NE-throw.As the low-density zone terminates at a depth from the surface of about 800 m, we suppose that the gas reservoir is located within the Ranzano Fm, as also suggested by Capozzi and Picotti (2002).This reservoir would correspond to the shallow reservoir buffering deep ascending gas, which was proposed in Bonini (2020) model.At the Nirano reservoir, gas likely escaped updip at the termination or in the linkage zone of a fault above a poorly sealed reservoir (Knipe, 1997).The gas leaked into the FAA probably in association with over-pressured mud dykes (Tingay et al., 2009) and/or directly into the fractured porosity of the fault damage zone.Similar situations of mud diapirism and intrusion along fault zones have been described also by Bonini ( 2020 Wang et al. (2022), andZhong et al. (2021).The gravity data that we present, however, are not in agreement with the seismic line interpretation presented in Maestrelli et al. (2019) who at 0.5-0.6 s depth still interpret the reflectors in the sequence as belonging to the FAA; we believe that those reflectors are within the Miocene sequence where there is enough significant contrast in acoustic impedance.The homogeneous FAA, on the other hand, would not cause any reflection.With reflectors within the Miocene, our gravity data interpretation would be coherent with the one of the seismic data (Maestrelli et al., 2019).The conceptual structural model that we present, however, concerns only the upper zone above a thrust-folded area, which was investigated by the gravity survey and is not in disagreement with the deep plumbing system proposed by Maestrelli et al. (2019).The novelty of our work is to show the importance of faults in controlling mud gas emissions at the surface of the Nirano Salse.
Our conceptual model for the Nirano Salse based on the gravity data is valid down to a depth of about 1,000 m from the surface.Below that depth, we cannot exclude the existence of other deeper gas reservoirs (DR in Bonini ( 2020)) associated with the thrusted folds of the Apennines, whose axes have the same NW-SE direction.
The thrust faults proposed by Bonini (2020) could crop out in the Argille Azzurre Fm to the NE of the investigated area, in which case the fault that we describe could have developed from extrados fractures of a deep anticline connected to fault propagation folds (Bonini, 2020).
The geophysical results of this study also point out some important characteristics of the Nirano mud volcanism, which have implications for the safety of the reserve visitors and the local communities.First, there is no evidence of a shallow extensive mud caldera as proposed by past studies (Bonini, 2008;Sciarra et al., 2019) as the low-density zone appears to lay on a NW-SE-oriented vertical plane.Furthermore, there is no evidence of a diapir (sensu strictu) below the Nirano area as suggested by Oppo (2012), Bonini (2008), andSciarra et al. (2019).The "bowl" morphology (Figure 1) so reminiscent of a caldera at Nirano does not seem to be caused by collapse along ring faults above a mud diapir or a mud chamber, but it is simply the surface expression of fault offset with a strong slip gradient toward its tip termination (Figure 4b).
The implications are that potential fluid escape structures be located along the fault zone or fracture zones (NE-SW) normal to the fault as also shown by the small isolated low-density zones in Figure 4a.These gryphons form above shallow (down to a few tens meters) mud/gas accumulations just below the vents whose plumbing system is connected to the low-density zone along the fault.Such a scenario was also suggested by hydrogeological and georesistivity surveys (Giambastiani et al., 2022;Romano et al., 2023).This information could be used to better plan access routes to the reserve and to fence off potential areas prone to sinking and subsidence.

Conclusions
The gravity inversion results indicate the existence of a low-density zone (1,200-1,500 m long, 100-200 m wide, 800 m deep) laying on a NW-SE-oriented vertical plane in the structural trend typical of the Northern Apennines.This zone likely represents the intrusion of fluids (mud and gas) in the damage zone of a subvertical fault, which feeds shallow fluid reservoirs (meters to a few tens of meters deep) just below the mud volcanoes both in the footwall and hanging wall of the fault.For the first time, this process is clearly documented along a fault zone.
The caldera like morphology of the Nirano Salse summital region is not related to the presence of a mud diapir s.s., or to a caldera collapse above a shallow extensive mud chamber, but to the surface expression of slip distribution at the fault termination along which the mud/gas fluids ascended to the surface.This has implications to assess better the safety of the site and the location of access points.The structure imaged by the gravity inversion has implications for the understanding of fluid flow migration along fault zones, especially where fault seal traps are involved and the combination of slip gradients, bed dips, and sealing units contribute to define the sealing risk of a geofluids reservoir.Gravity in our case was essential in testing between the caldera and fault intrusion hypotheses.This methodology, therefore, has shown to still have potential in subsurface characterization within complex structural areas.

Figure 1 .
Figure 1.Illustration of the study area.(a) Sketch of the main features of the study area.(b) Aerial image view of the study area (Image Landsat-Copernicus, Google).Magenta arrows indicate the four gryphons.(c) Zoom on the study area: the yellow circles represent the gravity stations.The biggest one represents the reference base station.The black line indicates the trace of the monocline as inferred by the present work.
during our work, pointed out only the presence of a NW-SE oriented monocline axis marking a steepening of the strata north of the mud volcanoes in the ductile clays of the FAA.The direction of the Bouguer low-density anomaly corresponds to the strike of the monocline axis and of the NW-SE oriented subvertical

Figure 2 .
Figure 2. Gravity data analysis.Circles represent the location of the gravity measures (stations).(a) Locations of the stations.The greatest circle is the reference station.(b) Altitude from DEM; (c) Relative free air anomaly, ∆g i RFA ; (d) Terrain Correction, ∆g i RT ; (e) Bouguer Anomaly ∆g i RB ; (f) complete Bouguer anomaly ∆g i RCB .

Figure 3 .
Figure 3. (a) Relative, complete Bouguer anomaly without offset estimated by the inversion procedure.(b-d) Results of the inversion for density contrasts performed with GROWTH 3.0 (Camacho et al., 2021).Plot of ∆ρ on a horizontal section at z = 50 m, the black line includes the area of mud volcanoes, the dashed gray lines indicate the traces of the AAʹ and BBʹ sections (b); vertical N-S (c) and NW-SE section (d).

Figure 4 .
Figure 4. (a) Isosurfaces computed for density contrasts |Δρ| = 300 kg/m 3 .The orange plane represents the fault.(b) Conceptual model showing the different geological structures: the fault termination with a steep slip gradient, the monocline, the mud/gas intrusion in the fault zone, and some ancillary faults are also shown.BAI are clay breccias of Ligurian origin (Eocene), RAN is the Ranzano Fm (Oligocene-Miocene), ANT is the Antognola Fm (Miocene), PAT is the Pantano Fm (Miocene), TER is the Termina Fm (Miocene), FAA are the Argille Azzurre Fm (Pliocene-Pleistocene); this latter formation (light yellow) has been removed in the 3D conceptual model to show the top of the deformed Miocene.The inset of panel (b) shows a cross section.Light green is gas, light blue is the zone of mud/gas intrusion into the fault zone.Density values are also reported.