Modeling Rain Isotopic Composition under Orographic Control: A Landscape Approach for Hydrogeological Applications

Abstract

Oxygen isotopic composition is useful for individuating recharge areas of groundwater bodies by the comparison with those of local rainfalls. While on a global scale general relationships, such as the isotopic vertical gradient or continentality effects, efficiently describe spatial variations of the isotopic signature, hydrogeological applications need spatial models that are more focused on the effects of local topographic structures and/or subsoil geology. This work presents a case study in northeastern Sicily (Italy) characterized by complex geological and orographic structures, in which isotopic composition of rainfalls is governed by orographic effects and the varying initial composition of humid air masses. We used a black box approach, comparing the average isotopic composition of rain collected from a network of eight samplers with their spatial descriptors (elevation, latitude and longitude). We obtained the best correlation with the simultaneous use of all these variables, applying their multiple linear correlation equation to transform the 1 × 1 km digital elevation model (DEM) of the study area into a digital isotopic model (DIM). The reliability of the DIM was confirmed by its good agreement with the oxygen isotopic composition contour map of the local groundwater.

1. Introduction

The isotopic signature of natural waters, expressed in terms of δ¹⁸O and δ²H, is an important tool for understanding the complex relationships between rainfall and groundwater [1], particularly with reference to the identification of the recharge areas of aquifers.

If contributions from external hydrogeological structures, and/or re-evaporation of infiltrated water, due to high temperature/low humidity climatic regimes are negligible, the average groundwater composition of an aquifer is equal to the long-term average of the rain isotopic composition of its recharge area. In making this assumption, the main problem is the estimation of the error in extrapolating and interpolating data acquired in from limited number of sampling points (rain gauges) to entire catchment areas.

To correctly address this issue, the main factors controlling the isotopic signature of precipitation have to be considered, such as:

• the initial composition of condensable atmospheric water vapor according to Rayleigh distillation processes [2];
• the dynamic and thermodynamic factors affecting rainfall events, such as internal cloud or cloud-base temperatures [3,4];
• the physical state of precipitation and rainfall rates [3]; and
• the condensation history of air masses, expressed in terms of distances from the sea and/or altitude gradients [3].

This study proposed a simplified method for describing the complex isotopic fractionation processes affecting rain under orographic control by using a “black box” model based on few territorial variables and able to match the precision level needed in the hydrogeological approach.

The target area was the northeastern sector of Sicily (Italy) (Figure 1), for which isotopic information for rain and groundwater are available as the result of a study carried out in the years 2004–2005 under an agreement between the INGV (Italian National Institute of Geophysics and Volcanology), Sezione di Palermo and the Water Authority of the Sicilian Regional Government [5]. The acquired data were further discussed in a previous work [6] and extended to the whole Sicily.

We re-elaborated these data, recalculating the vertical isotopic gradient and applying a more articulated spatial modeling algorithm to better describe the orographic drive on the isotopic composition of precipitations in the study area.

2. Theoretical Background

The isotopic signature of meteoric waters is expressed in terms of δ¹⁸O and δ²H, generally correlated by a strong linear relationship [3]. Thus, we can perform spatial modeling of the isotopic signature using only one of these two variables; this simplification was adopted in the development of the presented method and in particular all the following calculations will only concern δ¹⁸O.

The oxygen isotopic composition of precipitation is influenced by the kinetics and thermodynamics of the condensation processes, which are in turn driven by interaction with orography [3]. This interaction is particularly relevant in Sicily, where the isotopic composition of hydrometeors is ruled by the provenance of humid air masses, i.e., Atlantic weather systems from the northwest or African perturbances from the southern quadrants, and their interaction with the orographic structure [7].

Although on a global scale the vertical isotopic gradient and the continentality effect [3] efficiently describe spatial isotopic variations in natural precipitations, the complex processes described above make these variables ineffective for hydrogeological applications at detailed spatial scales. If condensation and precipitation of hydrometeors are controlled by orography, the direction of the vector representing the vertical gradient may coincide, totally or partially, with that of the “distance from the sea”, and the reverse of these vectors with the same or the opposite, leading to constructive or destructive spatial phase interferences between the different effects. In other words, spatial isotopic models should be based on more effective “geometric factors”, able to represent both horizontal (distance from the sea) and vertical (successions of valleys and mountain chains) variations, reciprocally interacting in all the possible combinations.

The simplest way to implement such a geometric factor is a search for multiple correlations between the average (weighed to rain heights) isotopic composition of precipitation and the triplets of spatial coordinates (east, north, elevation) topologically defining the points of a sampling network. The necessary condition for making this approach work is that the numerosity and geometry of the network must be properly designed to describe the orographic structure. Once a statistically significant correlation equation, linking isotopic composition to one or more spatial variables, is selected, isotopic models can be generated to transform a digital elevation model (DEM) into a digital isotopic model (DIM) [8]. This is the approach that we followed in this work.

3. Study Area Settings

The study area was located in the eastern sector of the Sicilian northern mountain chain, namely the Peloritani and Nebrodi mounts (Figure 1). Their physiography is characterized by mountains and hills formed from massive rocks, such as arenites and alpine metamorphic rocks [9]. In this area, the slopes are very steep and have accentuated V-shaped valleys and landslides are common, particularly along the coastal slopes [9].

The climate is typically Mediterranean to subtropical, characterized by warm-dry summers and cool-wet winters [9]. The mean annual temperature ranges between 19–20 °C in the coastal area, 14–16 °C inland and 8–10 °C on the highest mountains. The annual precipitation is usually lower than 700 mm, concentrated in the late autumn and early spring. The northern (Tyrrhenian Sea) and western coasts (Ionian Sea) are slightly wetter, with an annual precipitation that can reach 800–900 mm near the mountains [9].

From the geostructural point of view, the studied mountainous chain represents a complex south-southeast vergent fold and thrust belt build-up from the Neogene–Quaternary. It is the result of both postcollisional convergence between Africa and a complex European crust, and a coeval roll-back of the subduction hinge of the Ionian lithosphere [10,11,12]. Here we can distinguish:

• “European” Palaeozoic igneous and metamorphic rocks, which outcrop in the easternmost sector of the chain (Peloritani Mts.);
• “Tethyan” terrains (Upper Jurassic–Oligocene deep-water carbonates, sandy mudstones and Upper Oligocene–Lower Miocene flysch deposits) of the Sicilide domain; and
• “African” tectonic units (Sicilian domain) [13].

From the hydrogeological point of view, this geotectonic assemblage determined the formation of main aquifers allocated in fractured rocks (carbonates, sandstones and secondarily metamorphic and igneous rocks). The overthrust of clayey, low-permeable terrains, interposed between fractured rocks, locally interrupts the general hydraulic continuity of these aquifers, generating springs at medium-high altitudes. The coastal plain hosts huge underground bodies, exploited by hundreds of wells, fed by a lateral recharge from the inland mountain chain and only to a minor extent by direct precipitation (unpublished data from INGV-Palermo database).

4. Materials and Methods

Isotopic compositions of meteoric and groundwater used in this work refer to data acquired in the years 2004–2005 on behalf of the Piano Tutela Acque della Regione Siciliana (Water Protection Plan of the Sicilian Regional Government) [5,6].

A network of nine rain gauges (Figure 1 and

Table 1. Coordinates, yearly amounts and measured and calculated oxygen isotopic compositions of rain sampled from the eight-point network illustrated in Figure 1. δ¹⁸O values refer to yearly averages weighted to monthly rainfall amounts, expressed as the deviation per mil (δ‰) from the reference VSMOW. Calculated values were obtained by applying the first equation shown in Table 2.

ID	East UTM33N	North UTM33N	Elevation (m a.s.l.)	Rainfall (mm)	δ18O Measured	δ18O Calculated
SAM	467,973	421,3944	18	762	−6.1	−6.1
CCL	494,846	4,224,085	236	970	−6.4	−6.5
SBA	480,777	4,207,739	737	1223	−7.6	−7.3
SSC	442,871	4,208,031	102	731	−6.2	−6.3
ALF	473,815	4,208,401	391	1083	−7.0	−6.7
CES	474,378	4,188,851	1112	788	−8.1	−8.1
NDS	511,434	4,207,778	620	1159	−6.6	−7.0
ALI	537,113	4,206,322	22	826	−5.5	−5.5
CPE	557,072	4,235,966	21	616	−5.9	−5.8

) were installed and operated in the period from May 2004 to May 2005 in the studied area, collecting samples on an approximately monthly basis. The rain gauges were bulk collectors, filled with approximately 250 cm³ of pure Vaseline oil to avoid evaporation. The sites were chosen to cover the entire area from the coast to the maximum elevations. We excluded the sites snow is a significant part of the winter precipitations from our analyses because rain collectors were not designed for sampling snow, which avoids isotopic fractionation due to evaporation. Groundwater were sampled from a network of 121 sites, composed of springs, wells and drainage galleries. All the samples (rain and groundwater) were collected in double-capped PET bottles and kept refrigerated at 4 °C prior to analysis.

Water samples were analyzed for their oxygen isotopic composition in the laboratory facilities of Istituto Nazionale di Geofisica e Vulcanologia, Sezione di Palermo (Italy), using an Analytical Precision AP 2003 mass spectrometer. The isotope ratio was expressed as the deviation per mil (δ‰) from the reference VSMOW, with an uncertainty of ±0.1‰.

The DEM used for the spatial model had a dimension of 1000 × 1000 m and was extracted from a 100 × 100 m grid covering Italy [14]. Multiple linear regression analysis was carried out using the online tool MLRA [15]. Contour maps were drafted by implemented the kriging algorithm in Golden Software Surfer, release 16.

5. Results and Discussion

The bases for all calculations were the average yearly isotopic compositions of rain, weighed to monthly rainfall amounts, sampled in the network of eight points (Figure 1) and reported in Table 1.

The simplest way of identifying the possible relationship between oxygen isotopic composition and territorial variables is plotting δ¹⁸O versus the altitudes of the sampling points (Figure 2). The vertical isotopic gradient from our data was −0.002 δ‰ m⁻¹, intermediate between those found by Liotta et al. [6] (−0.0013 δ‰ m⁻¹) and by Longinelli and Selmo [16] (−0.0013 δ‰ m⁻¹).

The lowest gradient in [6], who used the same data as this work, was determined by computations made using all the rain gauges, including those where snow was a significant part of the precipitation and which we excluded due to the uncontrolled evaporation of snow after collection (see Materials and Methods Section). The highest gradient found by [16] was related to the diverse sampling network, located in the Mt. Etna area (circa 3300 m a.s.l. high) under very different orographic conditions. The other feature evidenced in Figure 2 is an additional horizontal drift of the isotopic composition, highlighted by the sampling sites located close to sea level (SAM, CPE and ALI). Moving eastward and southward (Figure 1), δ¹⁸O values drift toward less negative values (at constant elevations), suggesting a different balance between the amount of rain precipitated from Atlantic and African humid air masses (see discussion in Chapter 2).

The possible existence of an orographic control extended to the 3D space, as suggested by Figure 2, led to testing not only the vertical isotopic gradient, but also the possible relationship between δ¹⁸O and the relative positions of the sampling sites in the horizontal plane. We calculated the parameters of the multiple linear regression equation of the type:

δ¹⁸O = ax + by + cz + d,

(1)

introducing as independent variables easting UTM, northing UTM and elevations of the sampling sites. For the calculations, the results of which are reported in

Table 2. Parameters of the multiple linear correlation equations of the type δ¹⁸O = ax + by + cz + d, calculated using as independent variables the elevation (m a.s.l.), false easting and northing (km) as the differences between the easting and northing UTM values of the sites and the origin of the local x,y plane set at UTM coordinate point E = 273,000, N = 4,055,000, corresponding to the minima of longitude and latitude UTM for Sicily.

Independent Variables	a	b	c	d	r2
Elevation(x), E(y), N(z)	−0.002266	0.006524	−0.01566	−4.8	0.93
Elevation(x), E(y)	−0.001925	0.004339		−6.9	0.91
Elevation(x), N(y)	−0.00207	−0.003402		−5.3	0.88
Elevation(x)	−0.001989			−5.9	0.88
N(x), E(y)					0.39

, we used false eastings and northings (in km) as the differences between the easting and northing UTM values of the sites and the origin of the local x,y plane set at the UTM coordinate point E = 273,000, N = 4,055,000, corresponding to the minima of longitude and latitude UTM for Sicily.

As shown in the table, the best correlation was obtained when the triplets of coordinates were simultaneously used, with a determination coefficient of r² = 0.93. Looking at the determination coefficients of the other solutions, vertical gradient was the main factor driving the isotopic composition of rain, followed by longitude UTM and, to a minor extent, latitude UTM. No correlation (r² = 0.39) was found using the simple position on the horizontal plane.

Once the first equation of Table 2 was selected as the transformation function from the triplets of geographical coordinates of a given point to its oxygen isotopic composition, we transformed the digital elevation model of the studied area into a digital isotopic model (DIM); we adopted a grid of 1 × 1 km to avoid introducing isotopic variations at a high spatial frequency, deprived of any physical consistency. The related contour map is presented in Figure 3.

As shown in the map, the oxygen isotopic contribution fit well with the orographic structure of the area, with minima under −9 δ‰ in its westernmost and highest sector (Mt. Soro, about 1800 m high), and relative maxima (−5.5 δ‰) in the southeastern sector. The positive isotopic drift was due to the combined effect of decreasing maximum altitudes of the mountainous chain and the major weight of African air masses.

A similar map has been presented in a previous study of this area [6], also based on a DEM, but using only elevations for retrieving isotopic data by applying the simple vertical isotopic gradient. The main differences between our new elaboration and the previous map consist of:

• the values of δ¹⁸O of our map which are progressively lower than in the previous one as elevation increases, due to the steeper vertical gradient we calculated (excluding sampling points affected by snow evaporation); and
• the absence in the previous map of the SE positive isotopic drift revealed by our map, because the previous elaboration did not consider horizontal variations.

The validation of the proposed DIM was obtained from the oxygen isotopic composition contour map of groundwater, presented in Figure 4. A comparison of the maps reveals a good accordance between the two isotopic distributions, with δ¹⁸O of groundwater slightly shifted toward more negative values compared to those of related rainfalls. This difference is logical, from the hydrogeological point of view, because the recharge areas of aquifers develop uphill from them, thus at major elevations where rainfalls are more negative.

Finally, a comparison between the DIM and a classical contour map draft with the triangulated irregular network (TIN) method, i.e., with simple linear interpolation between measured values, further strengthened the proposed methodology. Figure 5a shows the contour map based on TIN interpolation, while Figure 5b shows the difference between the DIM and TIN; these maps adopt metric coordinates (UTM) instead of the angular ones of the previous plots due to the needs of performing metric calculations. As expected, the TIN map (Figure 5a) was much more simplified than the one based on the DIM (Figure 3), especially in the easternmost area where the absence of sampling points at high elevations caused the overestimation of δ¹⁸O values (Figure 5b). Conversely, the missing estimation of the east and south drift of the isotopic composition, due to the different prevalences of Atlantic and African humid air masses, led to the underestimation of δ¹⁸O values in the southernmost area (Figure 5b).

6. Conclusions

The proposed method was based on a geostatistical approach (MLRA) through which a DEM can be transformed into an equivalent DIM of natural precipitations. It was successfully applied to northeastern Sicily (Italy), where the δ¹⁸O of rain is controlled by both orography and different initial compositions of precipitable water vapor, which can be simultaneously described by a generic 3D point matrix. The results obtained in this specific case study can be extended to other areas, characterized by similar conditions, by just changing (and statistically testing) the transferring function from the easting, northing and elevation coordinates into the triplet easting, northing and δ¹⁸O. The transferring function, in the present study, was identified as a multiple linear regression equation, but it is important to note that its form can also be different (exponential, polynomial, etc.).