Mud flow dynamics at gas seeps ‑ Nirano Salse , Italy ‑

The Nirano Salse , known since the Roman Times, are one of the most beautiful and scenic mud volcanoes areas of Italy with thousands of visitors every year. In this work, we apply novel (for the context) hydrogeological techniques to characterize mud levels in the Salse by means of GPS-RTK positioning and continuous level logging within mud conduits. This is important to quantify the gas–liquid ratio in the conduits and evaluate the potential for dangerous abrupt mud eruptions. The results presented suggest that different mud levels in mud volcanoes clusters are due to the different gas–liquid ratio in the conduits and not necessarily exclude interconnection at depth, a hypothesis, on the other hand, that seems strengthened by mud level time series correlations. The presence of shallow aquifers at a depth of 5 to 30 m is also supported by our field data and allows us to delineate the boundaries of the shallow mud reservoir—pipes system and its overall shape. The shallow aquifers may provide a temporary storage for the ascending gas and when fluid pressure in these aquifers exceeds the tensional strength of the sedimentary rock, leakage of fluids to the surface would occur. In this case, if the gas–liquid ratio is high, mud volcanoes develop into tall gryphons and tend to have a discontinuous activity with sudden eruptions of mud after long periods of quiescence. This, together with the knowledge of shallow conduits localization has an important impli-cation for site safety in proximity to the mud volcanoes. Our inferences based on mud level relationships to mud extrusion dynamics can be applied to lower risk in other mud volcanoes areas of the world with high geo-tourist visits, such as those of Trinidad, Azerbaijan, and Colombia.


Introduction
Mud volcanoes are broadly distributed throughout the globe, both on land and below the seas (Milkov 2000(Milkov , 2005. In the last few decades mud volcanism has been widely investigated from geological, geophysical, and geochemical perspectives with important implications in energy resource exploration, seismicity, geo-hazard, and greenhouse gas emissions (Mazzini and Etiope 2017). Mud volcanoes form through a combination of gravitative instability, due to the overall low density of clay-bearing strata on surrounding units, and fluid overpressures, which can develop both in the same strata or in surrounding sedimentary rocks, shales, and reservoir rocks (Kopf 2002;Revil 2002). Hydrofracturing, due to the increase of fluid pressure or tectonic stresses, fault reactivation and seismicity, allows for the pressurised gas-water-sediment motion towards the surface and the brecciation of sedimentary units (Mazzini and Etiope 2017). The manifestation to the surface can happen via progressive and slow release of mud and fluids, or in violent and explosive forms (Mazzini et al. 2021). The sudden violent forms of mud eruptions may have serious consequences for the life of nearby visitors.
Mud volcanoes are peculiar natural features that have always attracted visitors in different contexts, such as cultural (geo-tourists), recreational (mud bathers), and religious (worshippers) settings. Geo-tourists interested in visiting mud volcanoes are numerous in Italy (Salse di Nirano, Macalube di Aragona, Salse di Montegibbio, etc.), Azerbaijan, Georgia, Turkmenistan, and Indonesia with visitations in the hundreds of thousands. Recreational mud bathing exploiting the allegedly beneficial properties of mud for the skin is popular in Trinidad, Colombia (Cartagena, El Totumo), and Azerbaijan (Alal). Religious significance is associated with mud volcanoes in many parts of the world: Hindus pilgrimage to the Chandrakup mud volcano in Pakistan and the Balka Devi in Trinidad. Famous is the temple of fire worshippers Ateshgah, which is located 30 km from Baku and was revered at different times by Zoroastrians, Hindus, and Sikhs. This territory is known for burning natural gas outlets associated with mud volcanoes (Gamkrelidze et al. 2021). Site safety, where frequentation is high, requires knowledge of the potential for sudden gas/mud eruptions and the shallow structure of the pipes-conduits system in proximity to the mud volcanoes, which may also be prone to sudden cave in and soil liquefaction to a distance of a few meters from the mud pool.
In the literature, however, not much attention has been given to the shallow structure of a mud volcano apparatus. Summit calderas have been described by several authors such as Fowler et al. (2000) in the Shah Deniz of the Caspian region and Evans et al. (2008) in the South Caspian Basin. A comprehensive model for a "mud volcano system evolution" has been presented by Stewart and Davies (2006). In their model, by means of seismic data analysis, they describe the growth and set of structures associated with a large mud volcano in the South Caspian Basin, which they assimilate to volcanic complex with associated radial and cylindrical dikes, circular faults, and vertical pipes creating a densely intruded zone where the volcano will develop. Wall rock erosion and compaction of the intruded zone leads to the collapse of a downward-tapering cone enveloping the cylindrical zone, capped by ring faults that define a km-scale caldera that downthrows the overlying mud volcano (Stewart and Davies 2006). While the model of Stewart and Davies (2006) is at the spatial scale of the "mud volcano system" (~ 100-1000 m), given it has been acquired by seismic data interpretation, our model is focused on a spatial scale which is about two order of magnitude smaller (~ 1-10 m). Accordingly, the novelty of our model is that we can study mud flow dynamics and conduits-shallow reservoirs plumbing system at small-scale and with more detail, being the system mostly concerned with 2-10 m wide gryphons and mud pools.
The mud volcanoes of Italy are clustered in four main geographical zones: in northern Apennines (mainly in the Emilia-Romagna Region); in central Apennines (Marche and Abruzzo Regions); in southern Apennines (in Basilicata, Calabria, and Campania Regions), and in Sicily (Gattuso et al. 2020). Italian mud volcanoes are typically small (area < 500 m 2 ), low-lying (only 5% exceed 2 m in height), and characterized by the continuous, but relatively quiet expulsion of gas, water, and mud. Martinelli and Judd (2004) ascribed most of their formation to tectonic compression in areas with thick sedimentary sequences. The fluids rise from considerable depth within these sedimentary sequences as consequence of overpressure and tectonic stress, possibly aided by the loss of density associated with gas generation (Martinelli and Judd 2004). Specifically, in the northern Apennine foothills near Fiorano Modenese a 4-km-long mud volcanoes belt is currently active and located in the Regional Natural Reserve of Nirano Salse (Modena, Emilia-Romagna Region) (Capozzi and Picotti 2002;Bonini 2012;Borgatti et al. 2019). This area has been widely studied since the Reserve was established in 1982. Several studies focused on the structural geological setting of the area to explain relationships between mud volcanism and active tectonics (Capozzi and Picotti 2002;Bonini 2008), as well as between mud volcanism, fluid venting, and fault/fracture zones (Bonini 2007). Seismic and geo-electrical investigations were carried out to detect the sub-vertical structures of the superficial outlet of the volcanic conduits and chimneys (Accaino et al. 2007;Lupi et al. 2016). Micropaleontological studies on cones, level-pool mud vents, and mud flows were used to identify the fluid source and reservoirs depths (Papazzoni 2017). Several geochemical studies were carried out to characterize the fluid venting, deep fluid compositions, origin of these fluids, and verify the possible contribution of meteoric water (Boschetti et al., 2011;Conti et al. 2000;Martinelli et al. 2012;Mattavelli et al. 1983;Mattavelli and Novelli 1988;Minissale et al. 2000;Oppo et al. 2017Oppo et al. , 2013; extensive geochemical soil gas survey and exhalation fluxes (CO 2 and CH 4 ) were also carried out in the mud volcano field (Sciarra et al. 2019) to show geochemical anomalies that would indicate the presence of high permeability conduits and preferential leakage pathways for the gas migration.
To the best of the authors' knowledge a detailed hydrogeological characterization of the Nirano area has not been attempted with our novel approach, yet. Mud pool levels, as well as conduits connections have been explained so far in terms of overpressured fluids that are periodically expelled from the main deep reservoir by the creation or the reactivation of pre-existing fractures/faults or pipes, without considering the hydrogeological aspects. Bonini (2007) stated that, most likely, the persistent mud extrusion could be due to short-lived escape of overpressured fluids along permeable structures connecting the reservoir(s) to surface. In this context, our paper aims at understanding the fluid dynamics and the connection of conduits and shallow aquifer-reservoirs to assess eventual elements of risk. First, we focus on the review of previous integrated geophysical, morphological, and structural geological studies carried out in the study area (Sect. 2). Then, insights on hydraulic connectivity and the effect of a variable mix between gas and mud on fluid dynamics are achieved by continuous monitoring of mud levels and chemical-physical parameters in several conduits, measurements of gas bubbles frequency, as well as mud grain size analysis and density determinations (Sect. 3). Results are then discussed, and a conceptual model related to the presence of shallow small fluid-saturated aquifers in which gas is temporarily stored during its ascent from the deepest reservoirs, before the final emission, is presented (Sects. 4 and 5). We believe that these latter results are very important for assessing the geosite safety in proximity to the mud volcanoes and the potential for abrupt eruptions not only at Nirano, but also in other mud volcano areas of the world (Gamkrelidze et al. 2021).

Morphology
The Regional Nature Reserve of the Nirano Salse is located at the edge of the Northern Apennines (Modena, Italy), an area that is characterized by several active fluid vents (Fig. 1).
With a surface of approximately 75,000 m 2 , Nirano is one of the largest mud volcano areas of Italy and Europe, as well as one of the most visited (more than 30,000 visitors each year). It is situated at the bottom of a caldera-like structure within a hilly landscape. The oval depression (caldera) is interpreted to be the result of the gravitational collapse of a mud diapir top, which broke to the surface within the Plio-Pleistocene succession, due to degassing, emissions of mud and fluid, and emptying of mud chambers (Bonini 2008).
There are several individual active cones, gryphons, and subcircular pools distributed at the base of the caldera within an area of recently extruded mud that has not been vegetated yet. The specific number and location of the cones vary over time as the area constantly evolves. Volcano morphology at the site depends mostly on the characteristics of extrusion style (abrupt Strombolian-type eruptions forming gryphons), persistence of degassing produced by the rising to the surface of salty and muddy water mixed with gas (mostly CH 4 ), and connection with groundwater aquifers (mud pools) (Onderdonk et al. 2011).
The largest cones have a basal diameter of about 4-5 m, a rim with a diameter of 0.2-2 m and, a height of 2-3 m above ground; they intermittently emit gas bubbles and muddy water in variable flow rates. Generally, there is a continuous supply of small bubbles (about a few centimetres in diameter), but their frequency may vary and every few minutes larger bubbles break the surface. The main Nirano mud volcano emits an estimated 100 to 300 m 3 of gas per day (Martinelli and Rabbi 1998). According to Oppo et al. (2013), the volume of emitted fluids varies significantly from well-developed mud volcanoes to gryphons, and mud pools. In the first case, the slow dense mud generated by abrupt eruptions runs down the sides of the cones, dries out, and increases the steep slopes of the structure; in the second case, mud pools seem to contain less viscous mud and have a more abundant and continuous emission of saline water and gas.

Mud characterization
Our observations and measurements with an electrical conductivity (EC) meter and gas accumulation chamber connected with a spectrometer show that fluid venting at the surface consists of a mixture of cold clay mud, brackish water, gas mixtures dominated by methane and sometimes CO 2 , liquid hydrocarbons, and peat fragments from levels encountered during fluid's ascent (data under publication). The mud extrusion is driven by the adiabatic expansion of methane during the ascent and the consequent decrease in the fluid density compared to the surrounding rocks (Mattavelli et al. 1983;Mattavelli and Novelli 1988;Martinelli et al. 2012), and the increase in overpressure where gas accumulates (bottlenecks in the conduits and reservoirs).
The saline composition, δD and δ 18 O signatures indicate that part of the waters originates from marine connate pore waters entrapped in the Miocene and Plio-Pleistocene sediments during marine deposition, with no contamination from recent meteoric water (Conti et al. 2000;Minissale et al. 2000;Martinelli and Judd 2004;Heller et al. 2011;Oppo et al. 2013). Geochemical analyses, of both water and gas samples, show that the gas contained in these reservoirs is mainly a mixture of primary and secondary thermogenic gases (due to thermal cracking of oil), and secondary biogenic methane (due to biodegradation of oils), with minor condensates and oil (Tassi et al. 2012;Oppo et al. 2017). Several authors (Mattavelli et al. 1983(Mattavelli et al. , 1993Riva et al. 1986;Martinelli et al. 2012) recognize that in this region of the Northern Apennines, the hydrocarbons (gas and oil) originate in the turbiditic sequences of the Tertiary Marnosa-Arenacea Fm. Micropaleontological analysis on cones, level-pool mud vents, and mud flows, however, indicated  Bonini (2007) homogeneous fossil assemblages, all consistent with the Plio-Pleistocene age of the Argille Azzurre Formation (FAA) (Papazzoni 2017) (Fig. 1b).
More recent data suggest that a common source rock responsible for the generation of these hydrocarbons could be located at greater depth than the Tertiary reservoir units Picotti 2002, 2010;Picotti et al. 2007;Bonini 2009;Oppo et al. 2013).

Structural geology
The Nirano area has a typical badland morphology due to the outcrops of the Argille Azzurre Fm (Fig. 1b), composed of fine-grained sediments of Pliocene to Pleistocene age. The study site is 4 km south of the main seismogenic Pedeapennine thrust (Fig. 1a) and Bonini (2008Bonini ( , 2009) maps a NW-SE oriented anticline axis dividing in half the geomorphological bowl (caldera) where the gas emissions are active today (Fig. 1b). Current seismicity in the area is due mainly to thrust (Magnitudes 3.0 to 4.0) and strike-slip faulting (Magnitudes 3.5-5).
The Nirano mud volcanoes are interpreted to be just above a NW-SE blind thrust anticline (Bonini 2008(Bonini , 2009(Bonini , 2012. The gas emissions would be the surface expression of fluids escaping from a deep leaky reservoir (about 1.5 km depth) located in the Marnoso-Arenacea Fm, and eventually from shallower (400-500 m from the topographic surface) reservoirs located in permeable Epiligurian units of Eocene-Miocene Age (Bonini 2007(Bonini , 2008. Any gas accumulation within a subsurface reservoir would increase the pore pressure due to buoyancy (Dasgupta and Mukherjee 2020); if the sealing unit is of poor quality, the gas may rise towards the surface either following faults and fracture zones or fingering through poorly consolidated mud sediments.
The mud volcanoes are roughly oriented along a NE-SW trend ( Fig. 1a and Fig. 2), which correlates with the direction of maximum horizontal in-situ stress (S H ) due to Apennine shortening. Bonini (2009Bonini ( , 2008 interprets both the alignment of the mud volcanoes and the shape of the geomorphological bowl as the result of the in-situ stress conditions present in the Nirano area.

Source of fluids
The source of fluids (gas and water) for mud volcanoes, in particular at the Nirano Salse site, is a matter of debate. According to Kopf (2002), the possible fluid sources for overpressuring and mud extrusion can be summarized as: (1) pore fluids from compaction; (2) biogenic methane from degradation of organic matter; (3) fluid migration along deep-seated thrusts; (4) thermogenic methane; (5) fluids from mineral dehydration; (6) hydrothermal fluids, and alteration of crustal rock; (7) fluid expulsion from internal deformation within the diapiric intrusion. By reviewing the abovementioned mechanisms for the Nirano Salse, we can argue that: (1) fluid expulsion during compaction has already occurred at Nirano, because the FAA already went through primary compaction; (2) biogenic methane from degradation of organic matter is not supported by the geochemical isotopic characteristics of the gas analyzed (Martinelli et al. 2012); (3) (Martinelli et al. 2012); (7) fluid expulsion from internal deformation within the diapiric intrusion does not fit well Nirano Salse dynamics and formation given that the gas source is below the FAA. Furthermore, mud diapirs are interpreted today as areas of seismic attenuation caused by intrusion of mud dike swarms formed by hydrofracturing due to overpressures (Fig. 1c) and not by the density driven phenomena of viscous flow that are typical of salt diapirs (Tingay et al. 2009). Bonini (2008) presents a geological model where pressurized fluids move up through discontinuities in the Ligurian Units, and accumulate in shallower reservoirs controlled by the lithological boundary between the impermeable claystone FAA and the underlying, more permeable, Epi-Ligurian Units and Colombacci Fm deposits. Bonini (2007) stated that anticlines in combination with the brittle structures associated with folding provide an efficient system for trapping and transferring fluids to the surface, where overpressured fluids are periodically expelled from the main reservoir (Marnoso-Arenacea Fm) through the creation or the reactivation of pre-existing fractures/faults (Fig. 1c). Permeability contrasts, tectonic loading, and gas generation likely represent the main factors triggering fluid overpressures (Bonini 2007). Vannoli et al. (2021) suggest that a sealbypass system, such as extrados fracture corridors developed along the fold crest, is needed to allow fluids to reach the surface. Furthermore, they claim that mud volcanoes can persist if they are connected by a network of deep long-lived structures possibly associated with background seismicity.

Mud monitoring and analysis
Two mud pools and two small gryphons were equipped with level-loggers LTC Solinst to monitor mud hydrostatic pressure, temperature (T), and electrical conductivity (EC) at different time intervals (every second for a total of 4.5-h monitoring, and every 5 min) (Fig. 2). The probes were installed at different depths as reported in Table 1. The time series collected ranges from July 7, 2020, to January 7, 2021. The hydrostatic pressure data measured by the divers were converted into mud levels using the following Eq. 1 and atmospheric pressure data recorded by the weather stations of Vignola (Modena) and available from the Regional Agency for Prevention, Environment, and Energy of Emilia-Romagna (ARPAE website; https:// simc. arpae. it/ dext3r/) (Fig. 3).
where MC is the mud column above the diver (m), P diver is the pressure exerted by the mud column and recorded by the diver (m H 2 0), P baro is the atmospheric pressure (m H 2 O), g (1) MC = 9806.65 P diver − P baro g , Shown are also logger depths (m below the ground-m b.g.) and altimetry (m a.s.l.). Gray cells mean that data are available in that period; white cells mean that data are not available. Refer to Fig. 2 for the mud pools locations is the acceleration of gravity (m/s 2 ), and ρ is the mud density (kg/m 3 ). The mud level (ML) in relation to the vertical reference datum (mean sea level, m.s.l.) is as follows: where TOC is the elevation of the mud pool rim in which the diver was installed (m a.s.l., Table 1), and CL is the cable length (m, Table 1).
Electrical conductivity values (mS/cm at 25 °C) were converted into salinity (g/L) using the conversion formula of Lewis and Perkin (1981).
To highlight hydraulic connectivity between volcanoes, variance of daily mud level (ML) elevation with respect to the mean mud level within the period of data analysis, as well as correlation coefficients between the monitored parameters (such as ML-Patm; ML-EC; ML-T; Patm-EC) and between mud volcanoes were calculated.
Mud level measurements were calculated as hourly averages to evaluate the barometric efficiency and describe how mud levels in the volcanoes fluctuate in response to atmospheric pressure. We used the slope method with mud level change (ΔML) on the y-axis and barometric-pressure change (ΔB) on the x-axis (Hare and Morse 1999; Gonthier 2007). Measurable ΔML is the mud level (ML) in the volcano at time (t + 1) minus the mud level in the volcano at time (t): (2) ML = TOC − CL + MC, while ΔB that causes the mud-level change (ΔML) is measurable as the barometric-pressure at time (t) minus the barometric pressure at time (t + 1): The order of time (t) and time (t + 1) in (3) is inverted with respect to (4) to follow the convention that ΔB is negative for an increase in barometric pressure during a time interval.

Fluid characterization
To characterize the fluid, a total of 22 mud samples (500 ml) were collected in N1, N2, and N3 (Fig. 2) and mud densities were measured in the laboratory using an electronic scale and graduated beaker.
Mud samples were also analysed for grain size. The samples were wet sieved with a 63 μm clear sieve to separate the muddy from the coarser fraction. The fraction retained by the sieve was dried in a natural convection oven at a temperature of 105 °C. The coarse fraction was subjected to analysis by sieving in the dimensional range between 8000 μm (-3 φ, Krumbein φ scale; Krumbein, W. C. 1934) and 63 μm (4 φ) using a battery of American Standard Test Sieve Series (ASTM) sieves with a particle size range of 1/2 φ and a vibro-tilting sieve. The passer-through (muddy fraction) was measured by aerometry, which correlates the size and percentage of particles suspended in the liquid with their sedimentation rate, based on specific weights of the liquid (4) ΔB = B(t) − B(t + 1). and particles themselves. The particle size distribution was reconstructed by software Gradistat v. 8.0.
Topographic surveys were carried out during the monitoring campaigns using a real-time kinematic differential global positioning system (RTK-DGPS, Fig. 2) with the aim of referring the rims and the mud level of the main gryphons and mud pools, as wells as the volcano field, to the mean sea level (m a.s.l.). Moreover, the conduit depths of the main mud pools and gryphons were measured by lowering a stainless-steel cylindrical weight down the conduit until progression was stopped.
The frequency and diameter of gas bubbles in eighteen mud pools and gryphons were visually counted and measured with a stopwatch, by photographs and video recordings; the average gas flow rates and hydraulic heads were compared with the aim of identifying the presence of gas bearing aquifers at different depth.
Finally, all data were collated with regard to the geological information available for the area, such as penetrometer, cores, and well logs data available from the archive of the ViDEPI project database (https:// www. videpi. com/ videpi/ videpi. asp) concerning Italian oil exploration within the National Mining Service for Hydrocarbons and geothermal energy of the Italian Ministry for Economic Development (UNMIG). The well logs considered for the present study are: Levizzano 1, Levizzano 2, and Maranello 1 (Fig. 1a), which were drilled in the 1960s by former Agip S.p.A. (now Eni S.p.A.), the Italian multinational energy company. Figure 4 reports the time series measured at 5-min intervals at locations N1 (a), N2 (b), N3 (c), and N8 (d) (Fig. 2). For N2, only the data acquired in the largest time window are shown (see Table 1).

Results
To better understand the behavior of the signals (Fig. 4) on long (weeks-months) time scale, we removed the high frequency noise by applying a Gaussian moving average filter on the raw time series of ML, T, and EC. To have a uniform dataset, the time series of atmospheric pressure (1 data point per hour) were re-sampled, with a linear interpolation, at the same time-step of the measured time-series (1 data point per 5 min). The noise filtered out from Fig. 4 is shown in Fig. 5. The EC sensor in the divers takes maximum 2 days to adjust to an EC value with a precision of 0.1 mS/ cm as apparent in Fig. 4c. A sudden drop in mud level (about 0.2 m) occurs at N3 (Fig. 4c) around 10/27/2020 without any change in P atm , T, and EC. There is another sudden drop of mud level at N8 (Fig. 4d) at the beginning of December, where mud level and EC drop (−0.2 m and 7 mS/cm, respectively); no correlation is shown with N3.
The high frequency noise (Fig. 5) is due to bubble activity (the bubble flowing next to the sensor causes a change in pressure; the more bubbles, the more activity) and it is the largest in N1 (amplitude of about 0.3 m), which has always been one of the most active (in terms of gas bubble emissions) gryphons during the monitored interval. In general, the noise is low at N1, N2, and N3 until September and it is high in the period September-October (Fig. 5).
The correlation coefficients of mud level, EC, T, Patm of monitored volcanoes are listed in Table 2. Within the 1-s dataset, the highest value is between N1 and N2 (0.8), while the lowest is for N2 and N3 (0.54), which are the closest pools. Correlation values decrease if the 5-min dataset is considered, with no correlation for N2-N3 (0.06) and good correlation for N1-N3 (0.5).
Correlation coefficients calculated on 1-s mud level dataset (Table 2) show high values due also to the short interval (about 0.5-1 s) between sampling frequency and gas bubble frequency in the conduit. The correlation coefficients for the 5-min mud level dataset show good correlation between N1 and N3, while it is lower between N2 and N3.
Mud in all volcanoes shows no correlation between ML and EC except in N1, which has good correlation (0.58). In the whole database, EC ranges between 9 and 20 mS/cm, with average values of 14.2 mS/cm.
The daily mud level variance within the periods of data analysis in N1, N2, and N3 is shown in Fig. 6. The graphs show variance peaks in different dates for each monitored mud pool. In N1, the maximum variance values are on 09/21 and 10/23, with ML variance values of 0.01 and 0.005 m 2 , respectively (6a). In N2, the maximum variance values are 0.0015 and 0.0005 m 2 on 07/21 and 07/31, respectively (Fig. 6b). In N3, the maximum variance values are recorded on 08/10 (0.0005 m 2 , Fig. 6c), and on 10/27 (0.009 m 2 , Fig. 6d). Figure 7 shows the barometric efficiency calculated for N1, N2, N3, and N8. Figure 8 and Fig. 9 show gas flow rates and mud level in all gryphons and mud pools presenting gas emission. In Fig. 8, mud level and gas flow are grouped based on their location (refer to Fig. 2). The highest flow rate is recorded in N14 with 0.0013 m 3 /sec, followed by N1 and N2 (0.0012 and 0.0010 m 3 /sec, respectively), while N4 and N5 present no gas bubbles ( Table 1 in Supplementary Material). In Fig. 9, mud level and gas flow data highlight two main groups: N1-N14 on the Eastern, and N15-N18 on the Western side of the Nirano Salse.
Regarding fluid characterization, the results of mud densities analysis are listed in Table 3. Mud density ranges from 1145 to 1350 g/l. To convert the hydrostatic pressure data measured by level-loggers into mud levels (1), we used the mud density value (1265.3 g/l), which was measured in N2, and it is consistent with an average value measured in N1, N2, N3, and N4.  Table 1). The time series of ML (black line), atmospheric pressure (magenta line), measured temperature T (green line) and electrical conductivity, EC (orange line) are reported in each panel  The particle size distributions of mud samples are shown in Fig. 10. All samples show similar grain size distribution and consist of clayey silt, with average values of 55.4 and 36.4% of silt and clay, respectively.

Discussion
Based on the results above, a conceptual novel model of flow dynamics is shown in Fig. 11 and justified in the following discussion section.
Our model, to the best of our knowledge, is the first to explain mud flow dynamics in the shallow subsurface (from the surface to several tens of meters) by seeking a connection between gas/fluid migration, mud levels in conduits, and shallow aquifers. From our data, we cannot focus on the deep migration routes and reservoir; for that we refer to the work of Bonini (2007Bonini ( , 2008. Nevertheless, our novel structural/aquifer connection model has some implications also for what proposed by Bonini (2007Bonini ( , 2008. Our data, in fact, suggest the presence of small shallow aquifers between 4-20 m depth where rising gas, from deep structures in the Epi-Ligurian units (Capozzi and Picotti 2010), is trapped and temporarily stored, so that deep fluid ascent is buffered in the shallow subsurface. In the present paper this process is recognized for the first time. The impermeable FAA (in which the shallow aquifers are encased) allows pore fluid pressure to build up till this overcomes the tensional strength of the rock at which point gas can escape along fractures or fractures-aligned pipe conduits, reaching the surface and dragging and fluidifying solid material along the way up (Fig. 11c). This flow dynamics allows both clay particles, gas, and eventually some silt (from the shallow aquifers) to rise to the surface with formation of gryphons and mud pools. As most of the mud pools are located around the gryphons, it is suggested that the overburden of the tall gryphons causes collapse and fracturing through which the fluids migrate, mixing with shallow meteoric waters similar to what proposed elsewhere by Mazzini (2009) and Mazzini and Etiope (2017). Most of the gas is rising from deep reservoirs (gas dominated fluid) and, perhaps, some deep connate water is carried along with the gas. The fluids reaching the surface may also mix with brackish water, which is present in the shallow aquifers, and this justifies the EC values that we measured in the field.
The presence of local shallow aquifers is supported by particle size distribution of mud samples collected in the field (Fig. 10) showing coarser granulometry than clay sediments. The presence of coarse sediments (sand and silt) within FFA is also confirmed by the AGIP S.p.A. core logs near the study area (Maranello 001, Levizzano 001 and 002, refer to Fig. 1a for location) and the explanatory notes of the Geological Maps of Italy (Scale 1:50,000, Sassuolo sheet; RER, 1999). Here, the FAA is described as consisting of silty and slightly marly clays, with a thin to medium stratification marked by intercalations of fine sands in flat, isolated, or connected lenses. On the right of the Secchia River, where Nirano is located, the sandy levels become thicker, laterally continuous and the stratification more powerful. Core logs show clays with frequent sandy and silty intercalations within the first 100 m of the FAA sequence. Sciarra et al. (2019), who carried out extensive geochemical soil gas and exhalation fluxes (CO 2 and CH 4 ) surveys, indicate the presence of high permeability areas that act as preferential leakage pathways for gas migration. These areas positively correlate with the dome-shaped conductive anomalies mapped by Lupi et al. (2016) at 20 m depth. We suggest that these permeable areas are none other than shallow aquifers with variable size and thickness possibly leaking to the surface along circular faults formed during the collapse of the area and the formation of the caldera-like morphology (bowl). The salinity of the muds is low (around 7 g/l) and correlates with the salinity of formation waters recovered during DST tests in the Levizzano 1 and 2 wells, as well as Maranello 1 from sandy layers within the FAA. The mud reaching the surface also does not contain any microfossil older that the age of the FAA, suggesting no direct mass transport from below them (Papazzoni 2017).
The connection between shallow aquifers is variable in time and depends on the gas flow activity. The measurement of conduits depths indicates that they range varies from 0.5 to 5 m with a mode around 1.5 m; at N2, however, we measured a depth of 8 m and 15 m on two distinct occasions. Our observations suggest that conduit opening is variable through time, some conduits may close and then reactivate; the whole conduit network appears to be in a state of continuous change and individual conduits are temporary features. Furthermore, the observations of Kopf (2002) that mud pools conduits are larger than those of gryphons is also confirmed by our observations. The N2 mud pool, in fact, is Correlation coefficients of mud levels variations between mud volcanoes depend on the conduit geometry, type, as well as their connection (i.e., direct connection and degree of tortuosity). High correlation coefficients between ML vents (Table 2 and Fig. 4) could be explained by the connection to the same source (shallow reservoir aquifer) whereas low correlation coefficients may represent separation of the conduits feeding the different pools and mud vents.
In the Nirano system, ML and P atm are not correlated (Table 2, Fig. 4) and this is further confirmed by the barometric efficiency analysis in Fig. 7. The influence of barometric pressure on a groundwater surface can follow several scenarios. In confined aquifers, the change in water level is caused by a change in the force applied to the Earth's surface by the atmospheric pressure-higher atmospheric pressure causes a greater load, which is transferred to water at depth, causing the water pressure within the aquifer to rise (Rasmussen 2005). The unconfined aquifers show delayed responses due to the delay in transmitting the atmospheric pressure signal through the vadose zone to the water table surface (Spane 2002). If the system were affected by barometric-pressure change, the slopes of the fitted linear curves in Fig. 7 would be positive. As it can be seen, the coefficients of determination (R-squared values) are all small, indicating no significant barometric efficiency (Fig. 7). Gas storage and overpressure in small shallow aquifers and flow from those to the surface seems more important than the barometric-pressure change.
The mud level variance (Fig. 6) is an indicator of gas activity in the monitored mud pools. Peaks appear following periods of quiescence during which gas accumulates and overpressure increases in the small aquifers system. The trend of these values indicates a chaotic system characterized by non-constant gas flow rates, with an alternation of quiescence and extrusion activity periods ( Fig. 4 and Fig. 6). This is further confirmed by some drawdown tests carried out on site (by emptying the N3 mud pool), during which no linear recovery of mud level was observed, testifying the presence of a non-continuous and constant flow in the conduits.
The different mud levels measured at the vents seem to be controlled by the variable gas-liquid ratio in the mud conduit ( Fig. 8 and Fig. 9). Low mud levels correspond to the vents that do not allow gas accumulation (fast continuous degassing) in their conducts and their pressure head is dominated by mud density (low gas-liquid ratio) whereas volcanoes whose conduits are occupied by gas bubbles that move slowly (high gas-liquid ratio) have a higher mud level dominated by the gas volume. The latter are more likely to erupt mud suddenly and unexpectedly than the former type. In Fig. 9, mud level and gas flow data highlight two main groups (N1-N14 on the Eastern side, and N15-N18 on the Western side of the Nirano Salse), which could indicate two different shallow aquifers at different depths (Fig. 11c). Neighbouring volcanoes have comparable but not equal mud level values. The mud viscosity is variable and heterogeneous, depending on climatic conditions, depth of measurements along the conduits, and, more importantly, gas flow. Also mud density is heterogenous and as gas flow increases, density will decrease along with viscosity. However, as shown in Table 2, gas flow is not the only cause. Dilution with rainwater and evaporation during hot and dry periods could contribute to density variability.   By considering the possible fluid sources that we discussed in Sect. 2.4 and our observations, we confirm that the reason for fluid release at the Nirano Salse is leakage of a hydrocarbon seal (Fig. 11c) as already proposed by other authors (Capozzi and Picotti 2002;Bonini 2007Bonini , 2008Oppo et al. 2013). Abnormal overpressure in a deep reservoir (Marnoso-Arenacea Fm) would be generated by gas accumulation following secondary migration. Fluids could escape from a seal broken by a fault, or gas may leak from the spill-point of a faulted reservoir layer. The seal could also have a valve behavior and fail when the overpressure in the gas reservoir increases due to the continuous gas migration from below. The novelty of our model is that gas following upward migration routes, such as faults and fractures (Fig. 11c), would accumulate in shallow aquifers confined within the FAA and then be released when the fluids overpressure would exceed the tensional strength of the seal (Gibson, 1994). In this way, the conduits forming the mud volcanoes would start at the depth of this shallow aquifers (5 to 30 m from the surface) as also measured by our soundings, which has not been recognized in the past. Conduits may have different shapes (cylindrical to fracturelike) and may be interconnected where mud volcanoes are in proximity as also suggested by our mud level correlations. The system of fractures and conduits feeding from the deep source into the shallow aquifers cannot be assessed by our work and is better addressed by geophysical methodologies.
The implications of recognizing the role of shallow aquifers in mediating gas/mud release at mud volcanoes is important as aquifer dynamics (recharge, depletion, etc.) is critical for the extrusion of mud at the surface. We observed that mud levels are low in conduits with continuous degassing, whereas they are high in conduits with discontinuous activity. These latter are associated with abrupt "Strombolian-type" eruptions that are a risk for visitors in proximity to the mud volcanoes, as demonstrated by the casualties at the Macalube di Aragona (Agrigento) (https:// www. ansa. it/ sicil ia/ notiz ie/ 2018/ 01/ 30/ esplo sione-vulca nelli-6-anni-a-fonta na_ 51897 811-841d-4e0a-835b-1727c db322 a2. html; consulted 15/05/202; https:// www. agi. it/ archi vio/ stori co/ small_ volca no_ erupts_ in_ sicily_ killi ng_ two_ child ren-15620/ news/ 2014-09-27/; consulted 17/05/2022). Mud levels dynamics mediated by the presence of shallow aquifers could be used as an indicator for geosite safety, also in other mud volcanoes of the world with high tourist frequentations.

Conclusions
This work has presented a novel conceptual model for mud flow dynamics in the shallow subsurface below mud volcanoes by highlighting the important role that surface aquifers have in controlling extrusion frequency and type. The model, despite being developed locally at the Salse of Nirano, has implications for other mud volcanoes areas in Italy (i.e., Macalube di Aragona-Agrigento) and of the world, such as in the Caucasus, Trinidad, and Colombia. In the specific, the model can help forecast activity of the mud volcanoes by means of mud level monitoring and contribute to improve the geosite safety for visitors.
The important points and conclusions that we can draw from our study are the following: 1. Mud level in the different mud pools and mud volcanoes depends on the gas-liquid ratio within the individual conduits, and on the geometric characteristics of the conduits. For this reason, different mud levels at mud volcanoes clusters do not exclude a connection of the conduits at depth. Mud levels time series in a cluster of mud volcanoes, in fact, suggests that they are somehow connected at the groundwater of shallow aquifers. The morphology of the mud volcanoes is also related to the frequency of degassing: mud volcanoes with low degassing form gryphons associated with abrupt mud emissions, whereas fast and continuous degassing is associated with mud pools or gentler slope volcanoes, whose activity consists of fluid mud extrusions. 2. Geologic data indicate the presence of shallow aquifers located in the first 30 m below the ground. These units could be distal turbidite lobes with coarser grain size than the surrounding clay, and act as small local temporary storage reservoirs where gas, raising from the main deep reservoir in the Marnoso-Arenacea Fm and Epi-Ligurian Units, accumulates generating fluid overpressures. The gas alimentation from the depth is not constant over time, and as the pressure in the aquifers increases and overcomes the tensional strength of the seal, it comes to the surface following new or pre-set conduits depending on the pressure in the aquifer. The gas liquefies and drags solid material from the conduit walls carrying along to the surface clay, silt, and peat fragments. Shallow aquifers recharge from the surface and from the fluid conduits at depth may contribute to mud flow extrusions frequency in time.
The presented conceptual model integrates and expands the geological and structural models already existing for mud volcano areas and completes the understanding of processes of mud extrusions and gas seepage by highlighting the links with shallow aquifers. More work, however, is required to understand at what stage of development is a mud volcanoes field (in terms of amount of fluids ascending from a deep reservoirs) and the connection between shallow aquifers processes and ascending deep fluids processes. Important issues to be addressed in the future research are: (1) the evolutionary stage of mud volcanism at Nirano: Is this phenomenon fading in time or is it dormant and preparing for new temperamental activity phases? (2) The use of the degassing level as an indicator for an imminent abrupt extrusion. (3) The change of mud level in conduits as an indicator to forecast important extrusion phases. These are important research questions that we can in part address with our observations and that are important issues, in particular, for areas where mud volcanoes attract geo-tourists around the world.
provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http:// creat iveco mmons. org/ licen ses/ by/4. 0/.