Characterising the heterogeneous nature of tufa mounds by integrating petrographic, petrophysical, acoustic and electromagnetic measurements

Determination of the physical properties of subsurface geological bodies is essential for georesource management and geotechnical applications. In the absence of direct measurements, this usually passes via geophysical methods such as seismic and ground‐penetrating radar. These require conversion to physical properties, and measurements at different scales to test for consistency. This approach is non‐trivial in geobodies with heterogeneous patterns of properties. Tufa mounds—in‐situ terrestrial carbonate buildups precipitating from geothermal waters—are characterised by high contrasts in facies and petrophysical properties from microscale to macroscale, and are therefore ideally suited to test the ability of non‐invasive geophysical methods to estimate such contrasts, and to develop petrophysical models based on geophysical properties. Here, a laboratory‐based study of a Pleistocene tufa mound in Spain is presented that combines (1) petrography, (2) digital 2D pore network analysis, (3) gas porosity and permeability measurements, (4) acoustic velocity measurements and (5) electromagnetic wave velocity and porosity determination from ground‐penetrating radar, to develop empirical petrophysical models. These results show the consistency of petrophysical properties determined with different methods across various observational scales. Electromagnetically derived porosity positively correlates with gas porosity. Petrophysical properties depend on measurable rock fabric parameters and the degree of cementation, which provide predictive tools for subsurface geobodies. Strongly cemented peloidal‐thrombolitic fabrics with intergranular and intercrystalline pores, and a dominance of small complex pores best transmit acoustic waves. Weak cementation and a significant fraction of large simple pores (framework, vegetation moulds) increase porosity and permeability of shrubby fabrics, while causing lower acoustic velocity. This study demonstrates that ground‐penetrating radar models can be used in combination with direct measurements of physical subsurface properties to capture highly contrasting physical properties associated with different sedimentary facies that would not be achievable with other methods, thus improving the understanding of formational processes.


| INTRODUCTION
Geophysical methods have been routinely applied for decades to infer petrophysical properties of the subsurface.Early studies include the characterisation of geological reservoirs using electrical methods (Archie, 1942;Schön, 2011;Waxman & Smits, 1968), or empirical models to describe the permittivity behaviour of natural liquids, solids or mixtures of both (Cole & Cole, 1941).Other more recent studies include efforts to link seismic interpretation to petrophysical properties of the subsurface via data inversion (Bosch, 2004), amplitude versus offset (Yardley et al., 1991) or attribute analysis (Possato et al., 1984).Advantages of these methods over other approaches such as coring include their minimally invasive nature, the wide range of sampling volumes or the potential for high-resolution measurements (Koesoemadinata & McMechan, 2003).However, there is a need for calibration in order to convert geophysical signatures into physical properties.For example, correct interpretation of seismic data depends on properly correlating the seismic response to the rock lithological and petrophysical properties, such as porosity and density, for which more than one solution may be possible (Kearey et al., 2002).This is non-trivial in carbonate rocks where impedance variations depend less on compositional changes, but-since all carbonate minerals have similar elastic properties-on complex and interrelated variations in fabric, porosity, pore type, pore geometry and diagenesis (Anselmetti & Eberli, 1993;Weger et al., 2009).An increasing body of literature has examined the controls on carbonate seismic response (Anselmetti & Eberli, 1993;Eberli et al., 2003;Weger et al., 2009), including the use of micro-computed tomography (microCT), nuclear magnetic resonance and pore network analysis to quantify controlling factors (Archilha et al., 2016;Bailly et al., 2022;Reijmer et al., 2022;Ronchi & Cruciani, 2015;Soete et al., 2015;Vasquez et al., 2019;Weger et al., 2009).Analogue data sets from buried carbonate systems as well as outcrops provide 2D and 3D constraints on geobodies and their seismic response.
Similar to seismic reflection methods, ground-penetrating radar (GPR) has been used extensively to link electromagnetic (EM) parameters, such as bulk dielectric permittivity, to subsurface petrophysical properties (Baker et al., 2001).This occurs via calibration of different mixing models that express the bulk dielectric permittivity in a system of solid, liquid and gas phases (Birchak et al., 1974;Huisman et al., 2003).Despite a large body of both field and laboratory-based studies using GPR over the last four decades to estimate petrophysical properties in carbonates (Conti et al., 2019;Cunningham, 2004;Harbi & McMechan, 2011;Mount & Comas, 2014;Mount et al., 2014), few studies have combined methodologies to specifically test the correspondence between physical parameters estimated with different methods (i.e.seismic versus GPR; Al-Shuhail & Adetunji, 2016;Ghose & Slob, 2006;Koesoemadinata & McMechan, 2003).The combination of these complementary techniques can more reliably determine parameters such as porosity and water saturation, and provide input to numerical modelling, particularly in geobodies characterised by markedly heterogenous architecture of facies and physical properties that may offer information on formational processes.A good example of such geobodies are continental carbonate tufas where a complex interplay of physico-chemical and biological processes, and secondary diagenetic overprint result in a heterogeneous network of primary and secondary pore types at a range of scales (Armenteros, 2010;De Boever et al., 2017;Della Porta, 2015).Semi-lithified tufas are ideally suited to reconstruct 3D geobodies using EM methods (Linares et al., 2010;Pedley et al., 2000;Pellicer et al., 2014).
This study investigates the Isona Tufa Mound Complex (ITMC, North-East Iberian Peninsula), which displays a wide range of subaerial tufa facies.The study expands on previous geophysical surveys in the area (Pellicer et al., 2014), by applying a combination of petrophysical, acoustic and EM methods constrained with direct petrographic and mineralogical measurements at the laboratory scale.The approach is used to infer changes in the physical properties (i.e.porosity, permeability or dielectric permittivity) for different distinctive facies of the ITMC.The results show the complementary nature of these measurements, allowing further definition of the spatial variability of physical properties in these unique environments.

STUDY SITE
Tufas of the ITMC are located in the Tremp Basin, an E-W trending structural piggy-back basin deformed acoustic properties, electromagnetic properties, ground-penetrating radar, porosity, tufa mounds and carried southwards by the Montsec thrust during the Pyrenean orogeny (Figure 1A; Vergés et al., 2002).Local topography reflects the underlying geology that developed during thrusting: the basin is bound to the north by the Boixols thrust and San Corneli hangingwall anticline, and to the south by the Montsec range, which developed above the Montsec thrust (Figure 1A; Puigdefabregas et al., 1992;Vergés et al., 2002).Upper Cretaceous-Eocene syn-orogenic sediments fill the Tremp Basin (Vergés et al., 2002).Marine conditions lasted to the Uppermost Cretaceous, when shallow-water carbonates, deep-water marls and the deltaic sandstones of the Aren Formation were deposited (Nagtegaal et al., 1983).Progressive basin fill established continental conditions in the Uppermost Maastrichtian; these are recorded by red claystones, coarse siliciclastics and carbonates of the Tremp Group (Pujalte & Schmitz, 2005;Rosell et al., 2001).
The ITMC is located in a continuation of the Isona anticline, a secondary fold structure plunging west (Figure 1A).The mound complex developed over an artesian spring system (Linares et al., 2010).The groundwater recharge area includes Aren Formation sandstones and the updip area of the Isona anticline (Linares et al., 2010) outcropping north and east of the ITMC.In addition, the Upper Cretaceous karstified limestone bedrock along the Montsec range to the south is hydraulically connected with the Aren sandstone aquifer (Figure 1A).This aquifer carries HCO 3 − -rich water into the Tremp Basin, where the overlying Tremp Formation acts as an aquitard (Linares et al., 2010;Rosell et al., 1994).Fracturing and erosional thinning of the aquitard, as well as denudation of topography to below the piezometric surface, allows groundwater discharge, and the establishment of an artesian spring system along the buried crest of the Isona anticline (Linares et al., 2010).
The ITMC forms a mesa caprock about 9 km 2 in area, and contains several discharge outlets with tufa deposits aligned along E-W and N-S trending faults (Figure 1B).Deposits represent subaerial carbonate precipitation from an artesian karstic groundwater system with water temperatures of 14.8 to -18.8°C (Linares et al., 2010).Tufa deposits are mapped as three distinct morpho-stratigraphic units (Figure 1B; Pellicer et al., 2014): 1. Three inactive outlets located on the 'Mont de Conques' mesa (Figure 1B) comprise 47 m thick tufa deposits with ages older than 350 ka (Linares et al., 2010).2. A series of fossil mounds formed between >350 and 214 ± 11 ka (Linares et al., 2010), occur at lower elevation west and north of the oldest mound.3. The youngest mounds with tufa deposits reaching up to 10 m thickness and ages ranging from 6 to 103 ka (Linares et al., 2010) are located in the northern area of the ITMC.These are associated with active groundwater outlets forming circular lakes (Basturs lakes).

| MATERIALS AND METHODS
A combination of petrographic, petrophysical, acoustic and EM measurements were carried out in the laboratory to characterise the physical properties of selected carbonate samples.

| Sampling and petrography
Samples for this study cover the whole age range of the ITMC (Figure 1B), and represent lithified carbonates from five facies (Table 1; details in Section 4.1).Blocks were plugged both parallel and vertical to depositional features such as lamination, and plugs trimmed for thin section preparation.Petrographic analysis of 11 thin sections was conducted under transmitted light to characterise depositional fabrics, cementation and pore geometries.Cementation was estimated on thin sections.Plugs were carefully trimmed for a cylindrical shape (diameter 2.5 cm, length 4.0-8.8cm) without fractures or chipped surfaces, which would disturb petrophysical measurements.Anisotropy in the plugs could lead to poor wave propagation through the plugs causing poor wave picks and unrepresentative results (Singh, 2007).Plugs were used for determination of porosity, permeability and acoustic properties.
A second set of limestone blocks (approximate dimensions 40 × 20 × 20 cm) were trimmed by saw in the field.These were used both for laboratory-based EM measurements, as well as determination of porosity and permeability on plugs cut from the blocks. 3.2 | Digital image analysis of pore space Thirteen representative thin section photographs from all five facies were analysed for 2D pore shape parameters using ImageJ© software (Table 2).The software allowed manual segmentation of high-resolution photographs into pore space and rock matrix.The resulting binary images were edited with the Despeckle and Fill Holes filters in ImageJ© to remove individual pixels and holes.Each image was manually checked to remove pixels along image edges, as well as erroneous data such as air bubbles that falsified segmented pore shapes, and non-porosity that was segmented together with porosity.At the selected image resolution, individual pixels had a size of about 3 × 3 μm, and so pore shape parameters were calculated on pixels with an area ≥10 μm 2 to further filter out individual pixels.The following 2D shape parameters were measured, following standard practice in earlier studies (Weger et al., 2009): • Area (mm 2 ) of segmented pores, median pore area (mm 2 ) and total thin section porosity (as a fraction of photograph area).As measured pore areas span a large range from 11 μm 2 to 2.5 mm 2 , they were grouped into 20 logarithmic classes analogous to the phi grain size scale.• Dominant pore size (DomSize) calculated as the upper boundary of pore areas making up at least 50% of thin section porosity.It is given as the radius (in μm) of a sphere with the same area.This parameter indicates the pore size dominating the sample.• Perimeter measures the length of the outside boundary of a pore (mm).• Perimeter over Area (PoA), calculated as the ratio of perimeter and pore area (mm −1 ).High values correspond to more complex pore shapes.• Circularity, calculated as with a value of 1 indicating a perfect circle, and increasingly elongate shapes as the value approaches 0.
• Aspect ratio (AR) is the ratio of major and minor axes of an ellipse fitted to the pore.

| Porosity and permeability measurement
Petrophysical data were acquired on 18 plugs (Table 3).Porosity was calculated using measurements of bulk volume (V b ) and grain volume (V g ).Plugs were dried in an oven to remove bound water as this could affect the bulk volume.Each plug was heated to 60°C and weighed every (1) 10 min until the drop-off in mass became negligible.After cooling, the length and diameter of each plug was measured using a calliper.Each plug was measured six times in different places to calculate an average value of cylinder volume V b .Equilibrium pressure was measured using a ResLab™ DHP-100 digital helium porosimeter.After calibration using samples of known volumes, equilibrium pressure was measured by letting helium gas expand into the sample chamber of known volume.These measurements were repeated three times per sample.Samples were left for 30 min between each run to reduce the chance of residual gas affecting subsequent results.These measurements were used to calculate grain volume V g for each sample using Boyle's law.Once grain volume was calculated, porosity was determined by: Porosity values were averaged using the three measurements of V g .
A ResLab™ DGP-200 digital gas permeameter calibrated for nitrogen gas was used to measure permeability using a steady state method.After completely confining each plug in a rubber jacket-lined core holder under ca 1.8 MPa confining pressure, nitrogen was allowed to flow through each plug.This flow was adjusted until the differential pressure in the jacket equilibrated.Three runs were performed per sample and an average value calculated to increase result accuracy.Each sample was left for 30 min to allow gas to filter out between each run.The differential pressure could then be used to calculate permeability using Darcy's law.Viscosity is temperature dependent and was calculated based on the temperature displayed by the permeameter.

| Laboratory-based acoustic velocity analysis
Acoustic velocities of p-waves (V p ) and s-waves (V s ) were determined on 12 plugs with benchtop equipment at ambient conditions (Table 3).Due to the fragile nature of samples, no measurements under confining pressure were carried out.Measurements determined the travel time of an ultrasonic pulse through the core plug.The signal was generated by a Tektronix AFG 2021 signal generator and fed through a Falco System amplifier to produce a wave of 18 ns and 100 ν.The sample was held in a vice between two transducers smeared with ultrasound gel to aid wave propagation.The vice generates a small (<0.5 MPa) pressure on the samples.The wave was picked up by a LeCroy Wave-ace 1002 60 MHz oscilloscope.Dead time of the transducers was also calculated and subtracted from the measured values.A different set of transducers was used for both P and S waves.The velocity was then calculated based on the distance-time ratio.Acoustic velocities served to calculate Poisson's ratio according to equation: The acoustic impedance Z is given by the equation: where ρ is density and V p the p-wave velocity.
Time average equations relate porosity and acoustic velocity based on theoretical and empirical considerations.The two equations applied here are the Wyllie Time Average (WTA) and Raymer-Hunt-Gardner (RHG) equations (Raymer et al., 1980;Wyllie et al., 1958): where V p min and V p fl are acoustic velocity in the mineral matrix, and pore fluid respectively (m/s).

| Laboratory-based EM velocity analysis
A Mala ProEx GPR system paired with two 1200 MHz GPR antennas were used in the laboratory to estimate EM wave velocity along seven sample blocks in transmission mode (Table 4).In this mode, only EM waves travelling one-way from the transmitter to receiver are used (instead of reflections).Measurements were collected for a total of seven samples and under two different conditions, (1) fully dry; and (2) fully saturated.In each case, an EM wave velocity is calculated through: where v is velocity (m/ns), c a constant (0.3 m/ns) and ε r(b) the bulk dielectric permittivity of the material.The distance that the EM wave will travel through the material allows calculation of v and thus the calculation of . This permittivity is then used to apply a petrophysical model, the Complex Refractive Index Model (CRIM).The CRIM is a three-phase dielectric mixing model that can be used to express the bulk dielectric permittivity of a solid, liquid and gas in a system (Huisman et al., 2003;Robinson et al., 2003).The CRIM expresses the bulk relative dielectric permittivity ( r(b) ) as: where ε r(s) is the relative dielectric permittivity of limestone (i.e.solid phase with specific values), ε r(a) and ε r(w) the relative dielectric permittivity of air (1) and water (79.5) based on water temperature in the sample measured at 22°C (Buchner et al., 1999), ϕ is the porosity, and S w the water saturation with values between 0 and 1, with 0 and 1 representing fully dry and fully saturated conditions, respectively.The factor α accounts for the orientation of the electrical field with respect to the geometry of the limestone, with values between −1 and 1, with 0.5 being used here as explained below.By obtaining both dry and wet EM wave velocity, a system of two equations can be generated from Equation (4).These equations are characterised by a unique ε r(b) , and a S w equal to 1 or zero for completely saturated and dry conditions, respectively, allowing for the isolation of both n and ε r(s) .This approach follows the laboratory setup described in Mount and Comas (2014).

| X-ray diffraction
Mineralogy of eight samples from all facies was determined by x-ray diffraction.Analyses were conducted using a Bruker D8 Advance Diffractometer (Cu-Kα x-ray source) at the Williamson Research Centre (University of Manchester).Samples were scanned from 5° to 70° 2θ, using a step size of 0.02° and a counting time of 0.2 s per step.

| Facies
Inactive and modern outlets and deposits can be subdivided into five facies representing distinct tufa sub-environments (Figure 2; Table 1).The central part of most inactive mounds and the modern lakes forms sub-circular pool depressions 30 to 140 m in diameter (Figure 2; Pellicer et al., 2014).Apart from the modern lakes, these depressions are filled by basal unlithified palustrine lime muds with oncoids, coated grains and stromatolites (Pellicer et al., 2014), overlain by detrital sediments reworked from the rim of the enclosed depression following lake desiccation (Pellicer et al., 2016).Locally in the 'La Cassola' mound, vents that acted as fluid conduits, are preserved within the pool (Pellicer et al., 2014).Vents form largely circular geobodies, are several metres in diameter, and rise about 1 m above the ground.The pool is separated from the outwardfacing slopes by a rim composed of a limestone ridge that rises 0.5 to 1 m above the pool ground (Figure 2; Pellicer et al., 2014).Slopes have angles of 5° to 7° in each direction (Figure 2) and are characterised by a stepped morphology of terracettes, microgours and channels (Table 1; Pellicer et al., 2014).Beds are subparallel to the morphological surface.North of the Basturs lakes, cascades drape the topography along the Abella River left bank and near active spring outlets (Pellicer et al., 2014).These consist of parallel beds with dips of up to 40°.Cascades generally are the youngest deposits from ca 6 ka to sub-recent (Linares et al., 2010).

| Petrography
All samples were crystalline and consisted of >99.5% calcite.A matrix density of 2.71 g/cm 3 was therefore assumed for all further analyses.Petrographically, the samples represent five lithotypes (Table 1).The peloidal lithotype is dominated by peloids embedded in microspar and 30 to 40% equigranular sparite (Figure 3A,B).This creates a granular-crystalline fabric with intergranular, intercrystalline, vuggy pores and fractures (Figure 3B).The outline of the intercrystalline pores is largely determined by the crystals surrounding them, and thus tends to be irregular (Figure 3B).Some of the pores have remains of organic matter, probably from plants.
A framework of sparry shrubs and peloids constitutes the sparry shrub lithotype.Shrubs are locally arranged in beds, and cemented by small sparite crystals (Figure 3E,F).The sparry shrubs consist of sparite crystals forming around hollow or micrite-filled upright and branching tubes (Figure 3G).They are equivalent to the spar-rhomb shrubs described by Guo andRiding (1994, 1998), and the sparry filaments of Gradzinski (2010).Cements make up 5 to 15% of the thin sections, but can increase to 20 to 25% where sparry shrubs are better developed.Pore types include intercrystalline pores between the small sparite crystals (Figure 3F), growth framework and vegetation moulds, possibly after bryophytes (Figure 3H; Melón & Alonso-Zarza, 2018).
The peloidal shrub lithotype (Figure 4A,B) shares the abundance of peloids with peloidal and thrombolitic fabrics.Peloids can be arranged into peloidal shrubs (Figure 4B), which are analogous to the micritic shrubs of Guo and Riding (1994) or the micritic dendrites of Della Porta (2015).Pore types include intercrystalline, framework, vuggy and vegetation moulds (Figure 4A,B).Equigranular and acicular cements make up 25 to 35% of the thin sections.
In summary, vent facies are associated with the peloidal lithotype (Table 1).Samples from the older 'La Cassola' rim   (sample RM-09) have the thrombolitic lithotype, whereas their younger counterparts at the 'Basturs' rim (sample RM-07) consist of the sparry shrub lithotype (Table 1).The cascade facies is associated with the sparry shrub lithotype, and the pool lithified stromatolites have the oncoidal lithotype (Table 1).A gradation of lithotypes is observed in the slope samples.Slopes from the older 'La Cassola' mound are dominated by the peloidal lithotype (sample SL-10).

| Petrophysical and acoustic properties
Petrophysical measurements resulted in gas porosity and permeability from 18.5 to 43.0%, and 2 to 14,000 mD, respectively (Table 3), showing a statistically significant positive linear correlation between the two variables (R 2 = 0.5438, R < 0.005; Figure 5A).The lowest porosity and permeability values occur in vent and older ('La Cassola' mound) rimstone samples (peloidal and thrombolite lithotypes), whereas the highest values are in the sparry shrub lithotype of the younger ('Basturs') cascade and rim facies (Table 3, Figure 5A).In slope samples, an increase of sparry shrubs at the expense of peloids (sample SL-10 to SL-02) parallels an increase in porosity and, to a lesser extent, in permeability (Figure 5A).No porosity and permeability data were obtained for the pool facies (oncoidal lithotype), as samples were damaged during preparation.Given the high porosity and permeability of cascade samples, it was difficult to maintain laminar flow during analysis, and the values obtained may overestimate actual porosity and permeability.Acoustic measurements resulted in V p values from 2055 to 6197 m/s, while V s ranges from 901 to 2859 m/s (Table 3, Figure 5B).In both cases, a statistically significant inverse linear correlation exists between porosity and acoustic velocity.The highest velocities occur in vent and older rim facies (peloidal and thrombolitic lithotypes) (V p = 5670-6197 m/s; V s = 2573-2859 m/s), whereas the lowest values correspond to cascade samples (sparry shrub lithotype) (Table 3, Figure 5B).Similarly, a strong inverse correlation also exists between gas porosity and acoustic impedance (Figure S1), with values ranging between 2.3 Ns/m 3 (cascade) and 12.4 Ns/m 3 (rim).In porosity-acoustic velocity space, most data significantly positively from time average curves (Figure 6), so that time average curves will underestimate V p values for a given porosity.The divergence is strongest for peloidal and thrombolitic lithotypes at the high-V p end of this study.The V p /V s ratios for ITMC samples fall between 1.5 and 2.3 (Table 3).Acoustic velocity and Poisson's ratio covary (Figure 7).Peloidal and thrombolitic lithotypes cluster at high velocities and Poisson's ratios, whereas the sparry lithotype dominates at the lower end.Two low-V p samples (CS-01 and SL-10) fall outside of the central data trend.If these are excluded, the correlation is good and statistically relevant (R 2 = 0.788, p < 0.005); otherwise it is poor and statistically insignificant (R 2 = 0.259, p = 0.091).

| Pore shape analysis
Although thin section porosity calculated from digital image analysis (DIA) consistently underestimates gas porosity measurements, there is a positive correlation between both (Figure S2A), with a correlation coefficient R 2 = 0.765 and p < 0.005 (Table 2).The DIA correctly replicates the low porosity of thrombolitic and peloidal lithotypes, and the high porosity of the sparry shrub lithotype (Figure S2A).In a diagram of cumulative pore area distributions (Figure 8), high-V p samples (VT-08 and SL-10; peloidal lithotype) tend to have normal distributions with a 50 percentile in the range 7.3-14.5 × 10 5 mm 2 .As sample porosity increases and Vp decreases (cascade and high-porosity rim and slope samples CS-01, SL-02 and RM-07, with sparry shrub lithotype), pore-size distributions develop a marked tail of larger pore sizes, and the 50 percentile increases to 12.0-44.0× 10 5 mm 2 (Figure 8).The exception here is sample RM-09 (thrombolitic lithotype), which despite having the lowest measured porosity, has an intermediate 50 percentile (11.0-19.0× 10 5 mm 2 ) and a distinct tail of larger pores (Figure 8).This is consistent with the thin section observation of mainly larger framework pores between mesoclots.
F I G U R E 5 (A) Gas porosity versus permeability for the ITMC data set, plotted according to facies (symbols) and lithotypes (colours).Note that no porositypermeability pair was obtained for pool facies/oncoidal lithotype.Slope samples are highlighted that show a progressive increase in porosity, which is related to changes in lithotype from more peloidal to sparry shrub.(B) Gas porosity versus V p and V s .Data are plotted according to facies (symbols) and lithotype (colours).

(A) (B)
The relationship between DomSize and gas or thin section porosity is weak (DomSize of low-porosity samples 9.1-223.7 μm, or 9.1-18.4μm if sample RM-09 is excluded; DomSize of high-porosity samples 120.9-180.6 μm; R 2 = 0.146, p < 0.005; Figure S2B).The high DomSize sample  reflects the dominance of framework pores between mesoclots.The sample aspect ratios and circularity do not show a relationship amongst themselves or with pore size parameters.On the other hand, there is an inverse relationship between median PoA and F I G U R E 6 Gas porosity versus V p for the current data set compared against various published data sets of marine (Eberli et al., 2012;Weger et al., 2009)

F I G U R E 7
Poisson's ratio plotted as a function of V p .Data are plotted according to facies (symbols) and lithotype (colours).Two outlier samples are indicated, which could be analytical artefacts of acoustic measurements, or realistic data.The linear regression shown does not take these samples into account.median pore size or DomSize (Figure 9, Figure S3).Samples with lower porosity and higher V p (samples VT-08, SL-10; peloidal and peloidal shrub lithotypes) cluster at the high PoA end (451.5-527.8mm −1 ), whereas samples with higher porosity and lower V p (CS-01, SL-02 and RM-07, sparry shrub lithotype) preferentially occur at the low PoA end of the spectrum (229.5-382.5 mm −1 ).
The DIA demonstrates that low-porosity, high-V p samples tend to have the peloidal lithotype (Figures 5B, 8 and 9), which is dominated by intergranular and intercrystalline pore types.These pores tend to be small and have a more irregular 2D outline due to polygonal crystalline margins (Figure 3B).As porosity increases (and V p decreases), larger pores become more important in peloidal shrub and sparry shrub lithotypes (Figures 3F,4B,5B,8 and 9).These are mainly mouldic pores after vegetation, as well as framework and vuggy pores (Figure 3F,H).The outline of these pores is simpler and smoother, reflected in their lower median PoA.Samples RM-09 F I G U R E 8 Cumulative pore area curves derived from image analysis.Colours represent V p measured on corresponding plugs.The lithotype determined on thin sections is given as well.Image analysis on the two pool stromatolite samples was considered unreliable, and these samples are excluded from the graph.and SL-05 form exceptions, with intermediate to very high PoA values (314.6-557.1 mm −1 ) despite having high and intermediate velocities respectively.The rock fabric of sample SL-05, although a peloidal shrub lithotype, has a high proportion of smaller intercrystalline pores, and so is similar to the peloidal lithotype.This is combined with larger framework pores that contribute to elevated porosity and lower V p in this sample.
Although framework pores dominate in sample RM-09 and create the tail in the pore area distribution (Figure 8), pores are lined or completely filled by well-developed cements, which contribute to lower porosity and higher V p .

| EM properties
Electromagnetic measurements gave values between 0.185 and 0.101 m/ns for dry conditions, and from 0.059 to 0.999 m/ns for fully saturated conditions (Table 4; Figure 10).Solving Equation ( 8) provided GPR estimated solid dielectric permittivity (ε) values between 3.4 and 5.7, and porosity values ranging from 10 to 46% (Table 4), with a linear inverse correlation between the two (Figure 11).Slope and cascade samples CS-01 and SL-05 were too brittle during analysis, and therefore no measurements under saturated conditions were possible.

| Controls on petrophysical and acoustic properties
The petrophysical and acoustic data (porosity, permeability and velocity) show clear relationships amongst each other.Whereas a positive correlation exists between porosity and permeability (Figure 5A), data span a large range, and permeability varies for a given porosity.This reflects the inherent heterogeneity of carbonate pore systems, in particular in continental carbonates (De Boever et al., 2017;Soete et al., 2015).The inverse relationship of porosity and V p (Figure 6) is well known for carbonates (Anselmetti & Eberli, 1993).Acoustic velocity is a function of bulk density (Kearey et al., 2002).Carbonate systems generally exclude large variations in matrix density, and carbonate minerals have very similar elastic properties (Anselmetti & Eberli, 1993;Verwer et al., 2008).Consequently, bulk density in almost pure carbonate systems such as the ITMC will depend strongly on porosity, hence the observed inverse relationship (Anselmetti & Eberli, 1993;Verwer et al., 2008).
Porosity and velocity data from this study are consistent with previously published data from marine and continental carbonates, and specifically continental tufas and travertines (Figure 6; Bailly et al., 2022;Eberli et al., 2012;Regnet et al., 2019;Reijmer et al., 2022;Soete et al., 2015;Vasquez et al., 2019;Weger et al., 2009).Datasets were variably obtained at ambient and confined pressures; the only tufa and travertine dataset shown that was measured under confined pressure is from Soete et al. (2015).It may be shifted to higher velocities due to faster wave transmission at confining pressure.The ITMC dataset has somewhat higher porosity and/or velocity values relative to the other tufa and travertine datasets.Despite the large natural heterogeneities in these lithologies, tufa and travertine datasets remain, however, broadly consistent.
All datasets diverge to higher V p values from time average curves, and this divergence is most prominent in continental tufas and travertines (Figure 6; Bailly et al., 2022;Soete et al., 2015).In the ITMC, mineralogy is predominantly calcite, and thus the large velocity and impedance variation cannot relate to mineralogical variations, such as clay minerals, which are commonly invoked for marine and lacustrine carbonates (Bailly et al., 2022 et al., 2019;Reijmer et al., 2022;Weger et 2009).Instead, velocity variations are probably a function of changes in rock fabric, pore type and pore geometry (Anselmetti & Eberli, 1993;Weger et al., 2009).The positive correlation between V p and Poisson's ratio in the ITMC (Figure 7) reflects a crystalline rock fabric, which has different acoustic behaviour than purely granular fabrics (Kenter et al., 2007).The development of an indurated framework at deposition is further supported by the observed V p /V s ratios for ITMC samples (1.5-2.3,Table 2) that largely fall within the range considered normal for indurated carbonates (1.8-2.2;Anselmetti & Eberli, 1993).The crystalline rock framework of travertines and tufas consists of framestones and cementstones with considerable compressive and shear strength (Soete et al., 2015).They often contain mechanically stiff mouldic and vuggy pores (Xu & Payne, 2009) that explain the positive V p divergence.The cascade sample has the highest V p /V s (2.3) and highest Poisson's ratio (0.38).This anomalous behaviour can relate to V s values dropping off relative to V p values in highly porous fabrics (Anselmetti & Eberli, 1993;Pickett, 1963), during closure of microcracks (Reijmer et al., 2022), or can be an analytical artefact.
In the ITMC, the relationship between rock fabric parameters and acoustic-petrophysical properties is approximated by variations between two fabric end members.Samples with high velocity and low porosity are dominated by peloidal and thrombolitic lithotypes with intergranular and intercrystalline pores (Figure 5B).The original granular fabric is strongly overprinted by microspar and equigranular cement (Figure 3B), the latter forming extensive bridges between grains (Figure 3B,D), together creating a very crystalline fabric.Interlocking crystals and the high degree of cementation create numerous contact points and longer contact lengths, which enhance acoustic wave propagation (Bailly et al., 2022;Eberli et al., 2003;Kenter et al., 2007;Storvoll & Bjørlykke, 2004).The more irregular pore shapes reduce permeability, but equally contribute to wave transmission (Archilha et al., 2016).Therefore, at one end of the spectrum in the ITMC, irregular and small pores, higher degree of cementation and interlocking crystalline fabric are probably key factors for the observed acoustic and petrophysical properties.
At the other end of the spectrum, the sparry shrub lithotype contains a tail of large macropores with relatively smoother outlines (Figures 8 and 9).To a large degree these are framework pores and vegetation moulds that were only lightly coated by cement crystals (Figure 3F,H).These macropores are visually very well connected, and the relatively low degree of cementation largely preserves open pore throats, favouring high permeabilities (Figure 5A).Intercrystalline porosity occurs between the crystals, adding to the high porosity of this fabric (Figure 3F).This very open crystalline framework with a combination of larger and smaller pores provides a lower number and area of framework contacts, but will equally promote scattering of acoustic waves (Soete et al., 2015).Both factors will reduce wave velocities.Finally, the acoustic velocities of some of these highly porous (>35%) samples are lower than predicted by the Wyllie equation (Figure 6; Anselmetti & Eberli, 1993;Pickett, 1963).The elastic rigidity of the hollow framework characteristic for the sparry shrub lithotype is low compared to marine frameworks such as corals, causing slower sonic velocities relative to predictions (Soete et al., 2015).Presence of micropores and microcracks (both original as well as artefacts from preparation of the fragile sample material) can also reduce sonic velocities (Bailly et al., 2022;García del Cura et al., 2012;Reijmer et al., 2022;Xu & Payne, 2009).
Between these two endmembers, the shrub lithotype combines peloids, peloidal shrubs and sparry shrubs in a largely crystalline fabric.Pore types include intercrystalline, vegetation moulds, and elongate framework pores.As the degree of cementation is comparable to peloidal and thrombolitic lithotypes, a key factor for reducing V p in the peloidal shrub lithotype is the development of more open rock frameworks via mouldic and framework pores that are associated with sparry shrubs.This is clearly visible in the increase in porosity and the decrease in V p that accompany the increase of sparry shrubs at the expense of peloids in the slope facies (sample SL-10 to SL-02; Figure 5A).Presence of vegetation moulds with a more peloidal-micritic framework tends to stiffen the rock (Bailly et al., 2022;Soete et al., 2015) and cause positive deviation from time average curves.
Acoustic and petrophysical properties in the ITMC are thus a function of development of a crystalline fabric, presence of framework and mouldic pores versus intergranular and intercrystalline pores, and the degree of cementation.As all older ('La Cassola') samples are fairly strongly cemented, and the sparry shrub lithotype has only been observed in the younger ('Basturs') samples, the age of the deposit is a secondary control.

| Relationship between porosity and dielectric permittivity
Laboratory-based EM wave velocity estimates under fully dry and saturated conditions are consistent with EM wave velocities estimated in the field for each ITMC morphological element by Pellicer et al. (2014;Figure 10).At low porosity, the bulk dielectric permittivity is buffered by the dielectric permittivity of the rock matrix, leading to only small variations between saturated and dry conditions (Figure 10; Mount & Comas, 2014).With increasing porosity, the relative contributions of the contrasting dielectric permittivities of water (79.5)versus air (1) will increase (Equation 8), causing a decrease or increase of bulk permittivity respectively (Figure 10;Neal, 2004;Topp et al., 1980).
Porosity and solid dielectric permittivity estimated from EM measurements show an inverse relationship, however given the limitation in datasets, the correlation coefficient is low (R 2 = 0.5; Figure 11) and not statistically significant (p > 0.05).As both were unknowns in solving Equation ( 5), coupled variability is to be expected (Mount & Comas, 2014).The determined permittivity data are within the range for carbonate rocks (4-8; Davis & Annan, 1989). 5.2.2 | Comparison between measured gas porosity and estimated porosity The advantage of the approach presented here is that porosity estimates from different analytical methods can be compared.This allows porosity calculation where only individual methods are available.For example, where V p is available (e.g. from wireline logs), it can be used to estimate porosity via the RHG equation (Equation 6).In the current study, comparing porosity calculated in this way against measured gas porosity, results in a good and statistically significant correlation (R 2 = 0.823 and p < 0.05; Figure 12).This calculation strongly underestimates measured gas porosity (by a factor of ≤4) at low porosity, but produces similar values at high porosity (Figure 12).As discussed above, due to the dominance of stiff rock frameworks and stiff pores in continental carbonates, the RHG equation would predict lower V p for a given porosity (Figure 6), so this divergence is to be expected.At highest porosity, lower stiffness of hollow rock frameworks and the effects of microcracks reduce V p , so the RHG equation better predicts V p values (Figure 6).
However, given that time average curves do not account for fabric and pore network variation, in particular in crystalline carbonates (Kenter et al., 2007), it may be more appropriate deriving porosity directly from measured velocities on travertine and tufa samples (regression line in Figure 6).For these samples, including those of the ITMC, velocity relates to porosity as: This results in a good and statistically significant correlation (R 2 = 0.825 and p < 0.05; Figure 12), which produces similar values as time average curves at low porosities, but diverges to overestimation at high porosities (cf.Soete et al., 2015).
Measured gas porosity and EM porosity determined from the CRIM model show a good linear correlation that is statistically significant (R 2 = 0.812 and p < 0.05; Figure 12).Porosities relate as: With the exception of one sample, gas porosities are higher (up to a factor of 2.3) than corresponding EM porosities (Figure 12).This shift may relate to incomplete saturation achieved during the EM measurements.If certain pore spaces were not fully saturated, remaining air (permittivity = 1) would cause a higher EM wave velocity, and consequently a lower permittivity and a lower porosity estimate (Mount & Comas, 2014).The helium gas used in petrophysical analysis on the other hand would have been able to enter micropores, and thus gives a more realistic porosity value.Nevertheless, the trends of gas and EM porosity suggests that the two methods can be combined for estimates of physical properties of subsurface materials.observed differences may reflect the different behaviour of EM waves and helium gas during measurements, the scale of observation volume, and rock heterogeneities captured within this volume (cf.Conti et al., 2019;Mount & Comas, 2014;Mount et al., 2014).The latter point is particularly true for continental carbonates.The dependency of porosity measurements on the reference volume being investigated has been recognised in reservoir petrophysics (Bailly et al., 2022;Fitch et al., 2015;Frykman & Deutsch, 2002) and GPR studies (Mount & Comas, 2014;Mount et al., 2014).Combination of analytical techniques at several scales of observation is therefore important for coherent characterisation of subsurface properties.Mount and Comas (2014) have demonstrated consistency in porosity determination between outcrop GPR (metre-scale) and laboratory-based GPR (decimetre-scale) measurements.The present study shows a consistency between laboratory-based GPR and plug petrophysics (millimetre/centimetre-scale).

| Implications for formational processes
The consistency between different methods at predicting physical properties can be related to formational processes in tufa mounds.Methods were consistent at identifying: (1) lowest porosity and permeability, highest velocity and well-cemented crystalline rock frameworks correspond to older vent and rim facies; (2) highest porosity and permeability, lowest velocity and less cemented rock frameworks correspond to cascade facies; and (3) intermediate values correspond to slope facies.Previous studies describing evolutionary models for the ITMC tufa mounds (i.e.Linares et al., 2010;Pellicer et al., 2014) used geophysical methods like GPR or electrical resistivity imaging (ERI) to image subsurface contrasts in physical properties (like dielectric permittivity or electrical resistivity), which were then used to interpret a concentric increase in porosity from the centre of the mound (where mound conduits were infilled by chemical precipitates as supersaturated groundwater emerged) to the edges (where cascade and slope facies accumulated below the rimstones).For example, electrical resistivity decreases from the vents towards the outer zones (Linares et al., 2010).In a similar manner, Pellicer et al. (2014) used 2D GPR profiles to define the interfaces between stratigraphic facies based on differences in EM wave amplitude and reflector signature (or radar facies).Given the limitations imposed in this study by the spatial distribution of samples, the limited number of analyses due to the fragile nature of samples, and the inherent geological heterogeneities of tufa carbonates, inferring similar 2D models and capturing the full petrophysical variations of the studied tufas is not possible.Despite these limitations however, correlations between various analytical approaches are fairly strong.This study confirms the increase in porosity from vent and rimstone facies to slope and cascade facies, and is able to tie acoustic, petrophysical and EM F I G U R E 1 2 Comparison of measured gas porosity against porosity derived from the methods applied in this study: (A, black data points) EM porosity calculated from the CRIM model (Equation 10), (B, blue data points) porosity from V p (Raymer-Hunt-Gardner Equation 6), (C, red data points) porosity from V p (linear regression of tufa and travertine data in Figure 6, Equation 9).
properties to measurable rock fabric parameters.At least the two endmember facies (1 and 2 above) can be differentiated even in the absence of petrographic data.Additional research on and other tufa complexes, using more samples, will be needed to reinforce these conclusions.
The approach followed in this study thus allows ping properties from wells or samples, such as porosity or saturation, onto 3D subsurface geometries (Chan & Knight, 1999;Knight, 2001;Neal et al., 2008), which is particularly useful in situations where no direct measurements of subsurface properties are possible.This ability is critical in many groundwater, geotechnical, environmental and reservoir applications (Ghose & Slob, 2006;Greaves et al., 1996;Knight, 2001), and with inclusion of a GPR system in the Mars Perseverance mission, can also constrain hydrogeological conditions on other planets (Pellicer et al., 2014;Rossi et al., 2008).

| CONCLUSIONS
1.A combined petrophysical (gas porosity, permeability, acoustic velocity, GPR-derived EM porosity) and petrographic (rock fabrics, pore shape parameters) study was carried out on a Pleistocene tufa carbonate mound.Rocks form a stiff crystalline framework consistent with other examples of continental carbonates.Petrographically, they vary between two endmembers: (1) a peloidal-thrombolitic rock framework with small and geometrically more complex intercrystalline-intergranular pores, and (2) a network of upright sparry shrubs with large and geometrically simpler framework, vegetation mouldic and vuggy pores.2. Rock fabric and degree of cementation are the dominant controls on petrophysical and acoustic properties, which accordingly vary between two endmembers.The strongly cemented micritic peloidal and thrombolitic frameworks are more strongly cemented and have low porosity and permeability.The polygonal crystalline framework and cement bridges allow better transmission of acoustic waves.Sparry shrubs have a more open framework with larger pores, limited cementation and reduced crystal-to-crystal contacts, which cause lower acoustic velocities and higher porosity and permeability.3. Porosity determined from GPR measurements positively correlates with gas porosity, although it will be lower for a given sample.This is probably an artefact due to inconsistent water and gas saturation of samples during measurements.4. In highly heterogeneous materials such as tufa, EM and petrophysical approaches yield compatible data across a range of observational scales.The study demonstrates that radar models can be compared with direct measurements of physical subsurface properties. 5.The results of this study support the spatial facies distribution for the tufa mounds derived in previous studies from GPR and resistivity, by confirming an overall concentrically decrease in porosity from the centre to the edges of the mounds.The highly heterogeneous and fragile nature of tufa samples only allowed collection of a reduced dataset that should be expanded through further studies.

ACKNO WLE DGE MENTS
An international mobility grant Jose Castillejo to J.

F
Fence diagram showing the interpreted GPR profiles and mound architecture in La Cassola (for location see top right inset and Figure 1B).Photographs 1 and 2 illustrate typical outcrop morphologies.Modified from Pellicer et al. (2014).
Perimeter over Area (median PoA) against median pore area.Colours represent V p , whereas the symbols reflect the observed lithotypes.Two pool stromatolite samples are excluded from the graph.

F
Linear dependence of solid dielectric permittivity on CRIM model porosity.Data are plotted according to facies.

Facies Sample Plug subsample He porosity He permeability (mD) V
p (m/s) V s (m/s) Abbreviations: N/A, analysis failed; n.d., no laboratory analysis carried out; n.d.*, no bulk density calculated; n.d.#, analysis carried out, but results considered unreliable.
Comparison of electromagnetic wave velocities for different morphological elements, estimated in the laboratory under completely dry and saturated conditions (this study), and estimates from fieldbased GPR from Pellicer et al. (2014).
; Regnet F I G U R E 1 0 P.C. (CAS19/00311) by the Spanish Ministry of Education, Culture and Sports (MECD) is gratefully acknowledged.The authors wish to thank Matthew McClellan and Katie Bonkowski (Davie) for GPR measurement laboratory support, and Steve May, Lee Paul and Ernie Rutter (Manchester) for acoustic-petrophysical measurement support.Constructive comments by Cathy Hollis, Giovanna Della Porta, and Michał Gradziński on an earlier version, as well as by two anonymous journal reviewers have greatly improved the quality of the manuscript.