Interaction between Perched Epikarst Aquifer and Unsaturated Soil Cover in the Initiation of Shallow Landslides in Pyroclastic Soils

Abstract

A physically based mathematical model of the slope of Cervinara (southern Italy), which is characterized by a shallow pyroclastic soil cover laying upon a limestone fractured bedrock, has been developed. Previous and current ongoing monitoring suggested that leakage through the soil–bedrock interface occurred, with leaking water temporarily stored in a perched aquifer located in the upper part of the fractured limestone (epikarst). This aquifer supplied several springs, and recharge to the deeper groundwater circulation occurred. Hence, in the proposed model, the unsaturated water flow taking place within the soil cover is coupled with the saturated water flow in the perched aquifer. The application of the model to the simulation of the slope hydrologic behavior over a period of 11 years, between 2006–2017, provides realistic results in terms of soil storage, epikarst storage, spring discharge, and groundwater recharge. The different response times of soil and epikarst aquifer to precipitation input allow distinguishing the hydrological predisposing causes of potential landsliding (i.e., a few months of persistent rainfall that is capable of filling the epikarst aquifer) from the triggers, which are represented by single intense rainfall events. The application of the model offers a key of interpretation of the hydrological processes leading to the landslide that occurred on 16 December 1999.

1. Introduction

Rainfall-induced landslides in shallow granular soil covers are among the most hardly predictable landslide types. In fact, the initiation of such landslides strongly depends on the characteristics of the triggering rainfall event (e.g., duration, mean intensity, peak intensity), as well as on the moisture conditions of the potentially unstable soil layer at the onset of the event e.g., Crozier et al. [1,2,3,4]. The triggering of shallow landslides requires an increased storage of water within the soil cover, leading to pore pressure increase [5]. Hence, the correct assessment of the water balance of the soil cover is mandatory, considering infiltration as well as drainage mechanisms that are both strongly and non-linearly affected by the initial soil moisture state e.g., Tromp-van Meerveld et al., Nieber and Sidle, Damiano et al. [6,7,8]. Particularly, the hydraulic behavior of the boundaries of the soil cover, through which it exchanges water with the surrounding hydrological systems, needs to be understood [5]. In fact, it has also often been recognized in shallow soil covers that soil moisture follows a long-term trend e.g., Cascini et al., Comegna et al. [9,10], and so the initiation of landslides has been related by many authors to rainfall accumulated over various time intervals e.g., Glade et al., Chleborad et al. [11,12], which are in any case much longer than the triggering event duration. This confirms the role played by large-scale (in space and time) hydrological processes as the predisposing causes of landslides, with the triggering rainfall being only the last push for the slope to fail [13].

In Campania (southern Italy), steep slopes covered with a few meters of pyroclastic loose granular deposits, which are mainly ashes (sand to loamy sand) with some layers of pumices (sandy gravel) laying upon fractured limestone bedrock, are sometimes subjected to fast destructive flowslides e.g., Picarelli et al. [14]. Although it is commonly recognized that slope failure occurs when rainfall infiltration leads to the vanishing of the contribution to shear strength offered by soil suction in unsaturated conditions e.g., Lu and Likos, Greco and Gargano [15,16], the process leading to the establishment of such a condition is still debated. Many studies indicate that a prominent role is played by the hydraulic behavior of the soil–bedrock interface. In some cases, the interface is considered as poorly pervious, owing to the presence of a layer of fine-textured altered soil (silty sand) with low hydraulic conductivity [17,18]. In other studies, the interface is modeled as a capillary barrier due to the coarse dimension of limestone fractures, under the hypothesis that they contain air at atmospheric pressure [19]. Conversely, other authors consider the soil–bedrock interface pervious, and focus on the formation of perched aquifers within the fractured limestone, which affect the water potential of the unsaturated soil cover [20], and may give rise to exfiltration, thus favoring slope failure [21].

The fractured limestone rocks of Southern Apennines are indeed characterized by the presence of large karst aquifers, for which regional mean annual groundwater recharge up to 500 mm has been estimated, ranging between 30–200 mm in the part of Campania where the limestone is covered with pyroclastic deposits [22]. The exchange of water between the soil cover and the karst aquifer occurs through the epikarst, which is the upper part of the karst system, and is weathered and densely fractured, more permeable, and porous than the underlying carbonate rock [23]. Often, a perched aquifer develops in the epikarst, from which percolation through vertical fractures and the conduits ensure the recharge to the deep groundwater system [24,25]. The highly conductive and porous pyroclastic covers of Campania favor diffuse rainwater infiltration and the formation of a thick (up to 10–15 m) epikarst [26], storing water that is drained by temporary springs [27]. The fractures of the epikarst are usually filled with the overlying soil, are penetrated by tree roots, and contribute to evapotranspiration e.g., Bakalowicz, Jukic and Denic-Jukic [23,28], so that the hydraulic connection between the perched aquifer and the unsaturated soil cover is possible e.g., Fu et al. [29].

In this paper, based on the evidence from field monitoring carried out at the slope of Cervinara (about 40 km east of Naples), a simple physically based mathematical model is proposed that describes the interaction between the unsaturated pyroclastic soil cover and the underlying perched aquifer temporarily storing infiltrating water. Coupled modeling of unsaturated and saturated flows has been indeed shown to provide a better description of water exchanges between the vadose zone and shallow groundwater systems e.g., Xu et al., Szymkiewicz et al. [30,31] than traditional uncoupled approaches with static boundary conditions at the base of the unsaturated zone. The model is applied to the simulation of the hydrological response of the slope to precipitations for a 11-year period, showing how the water level in the perched aquifer may affect the water potential of the soil cover, allowing to reproduce the typical seasonal soil moisture variability observed in the field. Finally, the simulation of the response in 1999, when a catastrophic flowslide occurred on the night between 15–16 December, offers a possible key of interpretation of the process leading to the triggering of the landslide.

2. Materials and Methods

This study refers to the northeast slope of Mount Cornito, near the town of Cervinara, in Campania (southern Italy), which is located in the mountainous area about 40 km northeast of Naples (Figure 1). The slope was involved in a series of shallow landslides that occurred in an area of about 4 km² during the night between 15–16 December 1999. The investigated area is part of the southern Apennines zone. The basal rock is Mesozoic–Cenozoic limestone, and the pyroclastic cover is mainly constituted by ashes (sand to loamy sand), with layers of pumiceous lapilli (sand with gravel). Near the interface with the bedrock, a layer of weathered ashes (silty sands) is present [18]. The slope ground surface has a regular inclination between 36° and 40° (Figure 1c). The deposit has the typical features found in the area, with the pumice layers only a few decimeters thick, while ashes reach a few meters. However, apart the foot of the slope, where accumulated colluvium reaches the depth of several meters of remolded material, the total thickness of the primary deposits rarely exceeds 2.5 m along the slope.

Table 1. Main physical characteristics of the pyroclastic soils of the cover of Cervinara (n, total porosity; γ_dry, dry unit weight; K_sat, saturated hydraulic conductivity; ф’, effective friction angle; c’, effective cohesion).

Soil	Clay %	Silt %	Sand %	Gravel %	n (%)	γdry (kN/m3)	Ksat (m/s)	ф’ (°)	c’ (kPa)
coarse pumices	0	0–4	48–60	38–56	50–55	1.2–1.3	5 × 10−6–1 × 10−5	-	-
ashes	2–5	4–14	75–85	0–8	70–75	0.7–0.8	1 × 10−6–6 × 10−6	38	0
fine pumices	0–4	2–15	55–60	25–40	50–55	1.2–1.3	-	-	-
weathered ashes	2–10	10–40	40–80	2–8	60–65	0.9–1.0	8 × 10−7–1 × 10−6	31	11

summarizes the main characteristics of the soils constituting the cover.

The climate is Mediterranean, with a mean annual precipitation of about 1600 mm, and total potential evapotranspiration, calculated with the empirical formula of Thornthwaite [32], between 700 mm (at the elevation of 750 m a.s.l.) and 800 mm (360 m a.s.l.).

Table 2. Mean potential evapotranspiration ($PET$ (L)), calculated with the empirical formula of Thornthwaite, and precipitation depth (L), based on meteorological data measured at Cervinara between 2001 and 2017.

Period	$\mathit{P}\mathit{E}\mathit{T}\left(\mathbf{mm}\right)\mathbf{at}320\mathbf{m}\mathbf{a}.\mathbf{s}.\mathbf{l}.$	$\mathit{P}\mathit{E}\mathit{T}\left(\mathbf{mm}\right)\mathbf{at}760\mathbf{m}\mathbf{a}.\mathbf{s}.\mathbf{l}.$	Mean Precipitation Depth (mm)
January	14	10	226
February	16	12	168
March	31	26	201
April	51	45	127
May	85	76	87
June	121	107	71
July	145	127	54
August	135	118	43
September	98	86	94
October	64	57	124
November	30	26	214
December	17	14	193
Year	808	703	1602

reports the mean monthly values of potential evapotranspiration and precipitation measured between 2001–2017. The vegetation cover is mostly constituted by deciduous chestnut woods (Castanea sativa), which are cultivated for fruit production. In late spring and summer, a dense understory of brushes and ferns (Pteridium aquilinum) flourish. There are many perennial and intermittent springs all around the slopes of the area, which are mainly located near the foot of the slopes, but in some cases also along the slopes, at mid-altitude (Figure 2). The springs supply a network of small streams and creeks, all reaching the nearby river Isclero.

2.1. The Hydrometeorological Monitoring Station

Hydrological monitoring activities have been carried out at the slope of Cervinara since 2001. Initially, this was done by means of several tensiometers installed at various locations and depths throughout the slope, which were read manually every couple of weeks, and a tipping bucket rain gauge with sensitivity of 0.2 mm [18]. Afterwards, between November 2009 and June 2012, slope monitoring went on by means of an automatic station equipped with the tipping bucket rain gauge, tensiometers, and Time Domain Reflectometry (TDR) probes. Soil suction and water content data were collected every two hours [10].

Since November 2017, new monitoring activities started. A hydrometeorological station was installed at the elevation of 585 m a.s.l., near the scarp of the landslide that occurred in December 1999 (Figure 2). It is equipped with a tipping bucket rain gauge with a sensitivity of 0.2 mm, a thermohygrometer, a pyranometer, an anemometer, and a thermocouple for soil temperature measurement. Six tensiometers (depths between −0.3 m and −3.2 m below the ground surface) and seven TDR probes (depths between −0.5 m and −2.5 m), which were connected to a Campbell Scientific TDR100 device, were also installed in an area of about 2 m², in order to instrument a vertical alignment of the soil cover. All of the data were collected hourly. Aiming at collecting information that is useful for the assessment of the water balance of the area, the water level in two streams close to the monitoring station is also monitored (Figure 2). One of the monitored stream sections is located at the foot of the slope (Section 1, elevation 535 m a.s.l.), and is supplied by several semi-perennial springs, so that some small discharge is sometimes found also in summer. The other monitored stream (Section 2, elevation 570 m a.s.l.), in contrast, receives water only during the wet season.

2.2. The Mathematical Model of the Slope

As will be shown in the next section, experimental data from field monitoring suggest that the hydrological regime of the unsaturated soil cover of the slope of Cervinara is affected by the storage of water within the underlying upper part of the limestone fractured bedrock (the epikarst). Hence, a simplified physically based mathematical model of the two interacting hydrological systems (soil cover and epikarst) has been developed. The regular geometry of the slope allows considering a two-dimensional 2D flow field with a constant inclination of the ground surface and constant thickness of the soil cover. The base of the epikarst is assumed to be an inclined plane, too (Figure 3).

The water fluxes through the soil cover, which are supposed as homogeneous for simplicity [33], are modeled with the 2D Richards equation written in the space coordinates $s$ (positive upslope) and $n$ (positive upward), respectively parallel and orthogonal to the slope:

$$\frac{d\theta }{d\psi }\frac{\partial \psi }{\partial t}=\frac{\partial }{\partial s}\left[k\left(\psi \right)\left(\sin \alpha +\frac{\partial \psi }{\partial s}\right)\right]+\frac{\partial }{\partial n}\left[k\left(\psi \right)\left(\cos \alpha +\frac{\partial \psi }{\partial n}\right)\right]-q_r$$

(1)

In Equation (1), $\psi$ (L) is soil water potential; $\theta$ (−) is soil volumetric water content; $\alpha$ (rad) is the slope inclination angle; $\theta \left(\psi \right)$ represents soil water retention curve; and $k\left(\psi \right)$ (L/T) represents the hydraulic conductivity function of the soil. The hydraulic characteristic curves of the soil, $\theta \left(\psi \right)$ and $k\left(\psi \right)$, are modeled with the van Genuchten–Mualem expressions [34]. The sink term $q_r$ (1/T) represents the root water uptake, which evenly distributes the evapotranspiration flux $E$ (L/T) throughout the soil cover according to the local value of water potential [35]:

$$\begin{array}{ccc}{q}_r=\lambda \left(\psi \right)E& \begin{array}{c}\begin{array}{c}\psi >{\psi }_1\\ {\psi }_1\ge \psi >{\psi }_2\end{array}\\ \begin{array}{c}{\psi }_2\ge \psi >{\psi }_3\\ \psi \le {\psi }_3\end{array}\end{array}& \begin{array}{c}\begin{array}{c}\lambda =\psi /{\psi }_1\\ \lambda =1\end{array}\\ \begin{array}{c}\lambda =1-\left(\psi -{\psi }_3\right)/\left({\psi }_2-{\psi }_3\right)\\ \lambda =0\end{array}\end{array}\end{array}$$

(2)

In Equation (2), the parameters ${\psi }_1=-0.5$ m, ${\psi }_2=-15.0$ m, and ${\psi }_3=-150$ m (i.e., the permanent wilting point of the plants) have been assigned according to values suggested in the literature [20,35,36]. The mean potential evapotranspiration flux $E$ of each month is obtained by dividing the values of $PET$ (L) of Table 1 by the soil cover thickness $H_c$ (L) and the duration of the month $T_m$, $E=PET/\left(T_m{H}_c\right)$.

At the ground surface ($n=H_c$), the infiltration capacity of the soil is calculated according to the current soil water potential value and gradient near the surface. Then, the rainfall intensity excess over the infiltration capacity is turned into overland runoff flow, and the difference is set as infiltration into the soil.

The water level $h$ (L) in the perched aquifer developing in the epikarst has been modeled by assuming the validity of the 1D Darcy equation:

$$n_e\frac{\partial h}{dt}=\frac{\partial }{ds}\left[K_{e}h\left(\sin \beta +\frac{\partial h}{ds}\right)\right]+q_b-q_s$$

(3)

In Equation (3), $n_e$ (−) is the effective porosity of the epikarst; $K_e$ (L/T) is the (saturated) hydraulic conductivity of the epikarst; and $\beta$ (rad) represents the inclination angle of the epikarst bed. The source and sink terms $q_b$ (L/T) and $q_s$ (L/T) represent, respectively, the flow crossing the soil–bedrock interface, and the outflow from the epikarst aquifer partly supplying the spring discharge and the deep groundwater recharge. For $q_s$, a simple linear dependence on water level is assumed, while $q_b$ is calculated as the flow rate crossing the interface at the base of the soil cover:

$$\begin{array}{c}{q}_s=K_e\frac{h}{L}\\ q_b=-k\left({\psi }_b\right)\left(\cos \alpha +{\frac{\partial \psi }{\partial n}\rfloor }_{n=n_b}\right)+dn\left\{{\frac{\partial \theta }{\partial t}\rfloor }_{n=n_b}+\frac{\partial }{\partial s}\left[k\left({\psi }_b\right)\left(\sin \alpha +{\frac{\partial \psi }{\partial s}|}_{n=n_b}\right)\right]\right\}\end{array}$$

(4)

In Equation (4), $L$ (L) represents the total length of the slope (Figure 3), as it has been assumed as a characteristic length that is suitable for estimating the mean gradient of the outflow from the epikarst; $n_b$ (L) represents the coordinate of the soil–bedrock interface in the adopted reference system; and ${\psi }_b$ (L) is the water potential at the base of the soil cover ($n=n_b$).

For the soil cover, a no-flow boundary condition has been assumed upslope, and a seepage face boundary conditions has been assumed downslope. For the epikarst aquifer, a free flow condition $\partial h/\partial s=0$ has been assumed upslope, and $h=0$, i.e., outflow at atmospheric pressure, has been assumed at the foot of the slope.

Besides the water exchange represented by $q_b$, the two flow fields have been coupled also by introducing the following equation at the soil–bedrock interface, neglecting the water flow through the upper unsaturated part of the epikarst and directly linking the water potential at the base of the soil cover with the water level in the epikarst aquifer:

$${\psi }_b=h-H_e$$

(5)

In Equation (5), $H_e$ (L) represents the epikarst thickness, $H_e$=$H_0+s\cos \alpha \left(\sin \alpha -\sin \beta \right)$, with $H_0$ being the epikarst thickness at the foot of the slope.

The system of differential equations of the model has been numerically solved with a fully implicit finite differences technique. At every time step, an iterative procedure has been implemented, so as to get the coupled solution of Equations (1) and (3). An upwind formulation has been adopted for the first-order spatial derivatives, while centered schemes have been used for the second-order derivatives. The non-linear algebraic system has been iteratively solved with an alternate direction implicit technique. To improve the stability of the numerical solver, water fluxes and hydraulic conductivity have been evaluated on a staggered grid [37].

3. Results and Discussion

In this section, first, some results from field monitoring are reported, which are at the basis of the development of the mathematical model described in Section 2.2, and which give some insight about the choice of the values of the parameters of the model. Then, the results of the application of the model to the simulation of the hydrological response of the slope to actual precipitation records are presented. The obtained results show that the coupled modeling of the unsaturated flow in the soil cover and the perched aquifer developing in the epikarst allows for a realistic simulation of the water potential regime in the soil cover, which is the result of the superposition of long-term seasonal fluctuations and short variations directly related to the sequence of rainfall events and dry intervals.

3.1. Slope Water Balance

The hydrometeorological monitoring activities carried out at the slope of Cervinara enable approximately assessing the water balance of the hydrological system consisting of the slope cover, the epikarst, and the deep groundwater system. More specifically, the measurement of rainfall and soil water content at various depths, the estimation of the evapotranspiration based on the measurements of meteorological variables, and the estimation of the water drained by the springs based on the measurements of water level in the streams, have been used to calculate the various terms of the water balance equation between 1 September 2017 and 6 June 2018:

$$V_s+V_e+V_g=\underset{0}{\overset{t}{{\displaystyle \int }}}\left(R-E-q_s\right)dt$$

(6)

Equation (6) is written in terms of specific storage, i.e., the terms $V_s$, $V_e$ and $V_g$ (L) on the left-hand side represent, respectively, the water volume stored in a column of unit area of soil, epikarst, and groundwater. The variation of specific storage in time is related, on the right-hand side, to the cumulated difference between rainfall intensity $R$ (L/T), evapotranspiration $E$ (L/T), and specific spring discharge $q_s$ (L/T). The overland runoff has been neglected in Equation (6), as it occurs only sporadically, during the highest rainfall intensity peaks, and does not significantly affect the long-term water balance.

The measurements of soil water content $\theta$, provided by the TDR probes installed at the depths of −0.30 m, −0.90 m, −1.40 m, and −2.20 m, allow estimating soil specific storage $V_s$ in the entire soil cover, for which a total thickness of 2.60 m was assumed:

$$V_s=\sum {\theta }_i\mathsf{\Delta }{z}_i$$

(7)

In Equation (7), ${\theta }_i$ and $\mathsf{\Delta }{z}_i$ represent the water content recorded by the i-th TDR probe and the soil profile height attributed to that probe, depending on the probe depth, so that $\sum \mathsf{\Delta }{z}_i=H_c$. The data of only four probes have been used, as the other installed probes experienced several out of order intervals during the considered period.

The water level measurements carried out in the two stream sections plotted in Figure 4, show that only at the end of November 2017 did water start flowing in the downslope stream (Section 1), and that after this date, there was always some discharge in the stream. In Section 2, instead, the establishment of flow was more intermittent. To estimate the stream discharge $Q_s$ from the measured water level $H$ (L), a stage–discharge relationship is needed. To establish such a relationship, the cross-sectional area of the stream Section 1 was measured on 22 December 2017, a day with low water level, and on 23 March 2018, when the highest ever observed water level was attained in both streams. On the same days, the flow mean velocity was estimated by timing small wood pieces floating within the stream. The obtained $\left(H,Q_s\right)$ couples were fitted with the following monomial stage–discharge relationship, in which the exponent assigned was equal to 2, so as to be close to the value of the relationship holding for uniform motion in a shallow trapezoidal section:

$$\frac{Q_s}{Q_{max}}={\left(\frac{H}{H_{max}}\right)}^2$$

(8)

In Equation (8), $H_{max}$ and $Q_{max}$ represent, respectively, the maximum measured water level and the corresponding maximum estimated flow discharge. The same relationship has been assumed valid also for the other monitored stream, in which the discharge was estimated only once, on 23 March 2018.

Table 3. Measured water levels, estimated flow discharges in the two monitored stream sections, and estimated specific spring discharge from the whole drained area. Bold values represent discharge estimated by measuring the cross-sectional area and estimating the mean flow velocity. The other discharge values are obtained by means of Equation (8).

Measurement Date	Section 1		Section 2		${\mathit{q}}_{\mathit{s}}(\mathbf{mm}/\mathbf{h})$
	$\mathit{H}\left(\mathbf{cm}\right)$	${\mathit{Q}}_{\mathit{s}}(\mathbf{L}/\mathbf{s})$	$\mathit{H}\left(\mathbf{cm}\right)$	${\mathit{Q}}_{\mathit{s}}(\mathbf{L}/\mathbf{s})$
7 December 2017	0.0	0.0	0.0	0.0	0.0
12 December 2017	2.0	0.2	0.0	0.0	0.0
22 December 2017	11.0	13.8	0.0	0.0	1.7
11 January 2018	14.0	25.3	0.0	0.0	3.1
18 January 2018	28.0	143.1	7.0	7.4	18.7
5 February 2018	18.0	47.4	0.0	0.0	5.9
8 March 2018	35.0	250.0	10.0	18.1	33.4
23 March 2018	35.0	250.0	15.0	50.0	37.3
20 April 2018	15.0	30.1	0.0	0.0	3.7
11 May 2018	14.0	25.3	0.0	0.0	3.1
6 June 2018	12.5	19.1	0.0	0.0	2.4

summarizes the measured water levels and the corresponding estimated discharges, in the two monitored sections. The obtained stream discharge values should be considered just as a rough estimate, which surely would change if the measurements were carried out with a more sophisticated technique. However, the obtained values are suitable to get a realistic assessment of the water balance of the system. The sum of the two measured discharges was divided by the entire underground watershed extension, which was estimated to be around 2.5 km² from the hydrogeological map of the area, in order to obtain the specific discharge of the springs draining the epikarst aquifer.

In Figure 4, the estimated trend of the specific storage of Equation (6) is reported. More specifically, as no information about the leakage from the epikarst toward the deep groundwater was available, the terms $V_e$ and $V_g$ could not be separated, and so only their sum $V_e+V_g$ could be estimated from the other measured/estimated terms of the balance equation. However, the spring discharge estimated from the stream water level was small, and it would further decrease during summer, when it usually becomes negligible. Hence, it can be argued that the epikarst storage was already very small at the end of the monitoring period, so that the estimated value of around 300 mm is mostly representative of deep groundwater recharge. Such a value is in reasonable agreement with the mean annual regional groundwater recharge estimated in the area [22].

Although the TDR probes were out of order for most of the rainy season, it is clear how the soil storage peaks shortly after the most intense rain storms (e.g., the value of more than 530 mm reached on 7 December), while the variations of storage in the fractured limestone are much smoother. Conversely, the drainage of water from the soil (reduction of soil storage) results were slower than the decrease of limestone storage.

3.2. Estimation of the Parameters of the Slope Model

The parameters of the hydraulic characteristic curves of the van Genuchten–Mualem model had been calibrated in previous studies by back-analysis of field measurements of soil water content and suction [33]. Although a single homogeneous soil layer is assumed, the parameters proved to be able to effectively reproduce the observed trends. The adopted parameters are reported in

Table 4. Parameters of the coupled mathematical model of unsaturated soil cover and epikarst aquifer.

Soil cover	Cover thickness, $H_c$ (m)	2.0
	Slope length (m)	200
	Saturated water content, ${\theta }_{sat}$ (−)	0.75
	Residual water content, ${\theta }_{res}$ (−)	0.01
	Air entry point, ${\alpha }_{vG}$ (m−1)	6.0
	Van Genuchten retention curve shape parameter, $n_{vG}$ (−)	1.1
	Saturated hydraulic conductivity, $k_{sat}$ (m/s)	3.0 × 10−5
Epikarst	Epikarst thickness, $H_e$ (m)	14
	Effective porosity, $n_e$ (−)	0.015
	Epikarst hydraulic conductivity, $K_e$ (m/s)	1.1 × 10−6

.

For the parameters of the epikarst aquifer, namely the effective porosity $n_e$, the epikarst thickness $H_e$, and the hydraulic conductivity $K_e$, the values reported in the literature results were widely spread, owing to the great differences exhibited by karstic systems in terms of density and dimension of fractures, which are affected by climate as well as by rock and soil cover mineralogy. The epikarst thickness is very much related to the characteristics of the soil cover and the yearly precipitation amount [23]. In Mediterranean climate, an epikarst thickness of 10 m was suggested [38]. In a karst aquifer covered by the pyroclastic soil of Campania (very close to the site of Cervinara), epikarst thickness up to 15 m was reported, depending on the soil cover thickness (the thicker the soil cover, the thicker the epikarst) and on vegetation cover (thicker epikarst under woodlands, thinner under grasslands) [26]. For simplicity, a constant thickness $H_e=14 \mathrm{m}$ was assumed. For $K_e$, values between 10⁻⁷ m/s and 10⁻³ m/s have been suggested e.g., Freeze and Cherry [39], and such a wide interval was confirmed by more recent studies e.g., White, Worthington, Allocca et al., Binet et al. [40,41,42,43]. For $n_e$, the reported values range from less than 1% [44] up to 7% [29,45]. In this study, the values of $n_e$ and $K_e$ have been chosen so that in an 11-year long simulation with real rainfall records, the perched aquifer resulted in dynamic equilibrium (i.e., it returned nearly to the same water level at the end of every hydrological year, with seasonal fluctuations consistent with the typical observed water potential regime of the overlying soil cover). The results of these simulations, which are given in the next subsection, led to the assumption of $n_e=1.5\%$ and $K_e=1.1\times {10}^{-6}$ m/s.

3.3. Long-Term Simulations

The coupled model has been applied to the simulation of the response of the slope to precipitations from 1 September 2006 to 31 August 2017. For this simulation, the rainfall height series with a sampling interval of 10 min, recorded by the rain gauge of Cervinara of the meteorological alert network of the Civil Protection Agency of Campania, has been used.

Table 5. Annual precipitations recorded by the rain gauge of Cervinara of the meteorological alert network of the Civil Protection Agency of Campania. Every year begins on 1 September and ends on 31 August.

Year	Annual Rainfall Height (mm)
2006–2007	1235
2007–2008	1431
2008–2009	2119
2009–2010	1739
2010–2011	1928
2011–2012	1122
2012–2013	2129
2013–2014	2000
2014–2015	1427
2015–2016	1524
2016–2017	1090
Mean	1613

summarizes the yearly precipitation depths of the series. The initial conditions on 1 September 2006 have been set by running a two-year long simulation with the precipitation record of 2013–2014 rescaled to have a total yearly precipitation that is equal to the yearly mean value (1613 mm).

Figure 5 reports the predicted trends of soil storage, epikarst aquifer water level, and yearly cumulated epikarst outflow during the entire simulation period. The monthly precipitation histogram is also plotted. It looks clear how the aquifer fluctuations are affected by the seasonal precipitation amount, especially in spring, when the evapotranspiration subtracts a significant water amount from the soil cover (Table 2). The trends predicted in all of the years resemble those observed in 2017–2018, confirming that the model is capable of reliably reproducing the actual response of the slope. The predicted fluctuations of soil storage correspond to the mean volumetric water content in the cover between 0.3 and 0.5. The fluctuations of epikarst storage correspond to soil water potential at the base of the cover, oscillating between −12.0 m to −14.0 m in summer, and −1.0 m to −3.0 m in winter. Both the modeled seasonal trends of soil water potential and water content result in good agreement with field observations available for similar slopes e.g.,Napolitano et al., Cascini et al., Comegna et al. [4,9,10].

It is interesting to note the different response time to precipitations of epikarst storage, which shows a smooth behavior, compared to soil storage, the graph of which is characterized by abrupt peaks following rainy periods. The effects of potentially triggering events (heavy rain storms) are clearly distinguishable from those of the hydrological predisposing cause, which are represented by the seasonal soil moisture fluctuations related to the water level in the aquifer. It becomes clear why in some cases (e.g., 2008–2009 and 2010–2011), although the occurred long-lasting precipitations (in both years, monthly precipitation higher than 600 mm was recorded) led to unusual filling of the epikarst aquifer, this predisposing condition did not result in any landslide triggering. In fact, during the periods of maximum water level in the aquifer, none of the occurred rain storms was large enough to cause a significant increase of soil storage, i.e., the infiltration rate was nearly compensated by the leakage towards the aquifer. The maximum soil storage values—slightly higher than 1000 mm—were attained in January 2009, November 2010, and February 2014, and correspond to a mean water content in the profile of about 0.5, which is still far from saturation.

3.4. Simulation of Year 1999

The importance of the predisposing conditions induced in the soil by the high water level developing in the underlying epikarst aquifer is confirmed by the results of the simulation of the year 1999. The total yearly precipitation depth was 1807 mm, i.e., not very high, although larger than the average. However, a landslide was triggered along the slope in the early morning of 16 December, after a rain storm of about 325 mm in 50 h. The initial conditions at the beginning of the year have been assigned, in terms of soil water potential and aquifer water level, by taking the average of the values, obtained on 1 January at 0:00, of all 11 years considered for the long-term simulations (i.e., from 2007 to 2017).

Figure 6 reports the predicted time history of the water level in the aquifer, together with the daily hyetograph, for the whole year 1999. It looks clear that during the first half of the year, the predicted aquifer water level remained always at least 5 m below the soil–bedrock interface, so that dry conditions at the base of the soil column were present during the entire period. Nonetheless, after the typical summer decrease, which led between 10–13 October to ${\psi }_b<-10 \mathrm{m}$ at the depth of −2.0 m, during the following period of frequent precipitations (412 mm in less than two months, between 15 October and 13 December), the aquifer level rose nearly 5 m. In the same period, the mean soil water content grew only from 0.39 to 0.46. Hence, the process of water level build-up in the perched aquifer, developing in about two months, represents the hydrological cause that created the predisposing conditions for the slope failure. The triggering event consisted of nearly 325 mm of rainfall fallen in 48 h, between 6 a.m. on 14 December and the same hour on 16 December. Such a rain storm effectively induced slope failure thanks to the existence of the predisposing conditions.

In fact, in an unsaturated cohesionless shallow cover characterized by an inclination angle close to the friction angle as the considered case, the triggering of a landslide requires the attainment of a nearly saturated condition. In fact, under the infinite slope hypothesis, which can be assumed valid for a large part of the considered slope, the safety factor assumes the following expression:

$$F_S=\frac{{\tau }_{lim}}{{\tau }_{\alpha }}=\frac{c^{\prime }+\left({\sigma }_{\alpha }-u_a\right)\tan {\varphi }^{\prime }+\left(u_a-u_w\right)\chi \tan {\varphi }^{\prime }}{\gamma \left(H_c-n\right)\sin \alpha \cos \alpha }$$

(9)

In Equation (9), ${\sigma }_{\alpha }$ and ${\tau }_{\alpha }$ represent, respectively, normal and shear stress components acting along a plane parallel to the slope; $u_w$ and $u_a$ are pore water and air pressure (the latter assumed as equal to atmospheric pressure); $\chi$ is the Bishop coefficient, which was introduced to calculate the contribution of soil suction to the shear strength of an unsaturated soil, and is related to the degree of saturation of the soil e.g., Greco and Gargano, Lu et al. [16,46]; $\gamma$ is the mean specific weight of the soil layer between the ground surface and the depth $H_c-n$. As it is $\left({\sigma }_{\alpha }-u_a\right)=\gamma \left(H_c-n\right){\left(\cos \alpha \right)}^2$, in a soil with $c^{\prime }=0$, similar to the ashes of Cervinara, the vanishing of suction implies that for $\alpha \cong {\varphi }^{\prime }$, it results $F_S\cong 1$, i.e., landslide-triggering conditions.

Given the high porosity of the considered deposit and the typical soil water storage during the wet season rarely exceeding a mean water content of 0.5, a large water volume has to be stored in the deposit to approach saturation, which is namely about 500 mm for the considered soil cover thickness of 2.0 m. As the triggering rainfall depth was “only” about 325 mm, this implies that saturation was reached only locally within the soil cover, and that the drainage from the soil cover during the triggering event was ineffective.

When the water level in the aquifer is high, the water potential gradient required for the establishment of an intense flow rate $q_b$ through the soil–bedrock interface may lead to soil saturation (suction vanishing), and consequently to the triggering of a landslide along a slope with an inclination angle close to the soil effective friction angle, as the considered slope.

Figure 7 reports the volumetric water content profiles throughout the soil cover that were simulated at various times between 14–18 December. It appears clear how with the given initial conditions, the intense infiltration led to the establishment of a high soil moisture gradient at the base of the cover, so that most of the infiltrating water remained stored in the upper 1.5 m until the early morning of 16 December. This is due to the dry state of the soil at the base of the cover before the triggering event (i.e., $\theta \cong 0.26$, corresponding to ${\psi }_b\cong -5 \mathrm{m}$) caused by the relatively low water level in the perched aquifer. With the low unsaturated hydraulic conductivity of such a dry soil layer, only under large water potential gradient water could penetrate the base of the cover and cross the soil–bedrock interface. Through the soil water retention curve, large water potential gradients correspond to large soil moisture gradients. The simulated soil moisture profiles show that the maximum soil moisture gradient at the base of the cover developed on 16 December, between midnight and 5 am (which was the actual triggering time), leading to water content as high as 0.68 between 1.0–1.3 m below the ground surface. As said, the attainment of this value is related not only to the characteristics of the triggering rainfall event, but also to the initial soil moisture profile and to the water content at the soil–bedrock interface, which in turn depend on the water level in the epikarst aquifer. Both these latter features are the result of the predisposing hydrological phase.

The effect of the infiltration features on the potential triggering of a landslide is shown in Figure 8, which reports the predicted evolution of the factor of safety at various depths, as well as the trend of the mean water content of the soil cover during the triggering event. The worst conditions for slope stability are predicted around midnight on 16 December, when $F_S\cong 1.3$ was attained at depths between 1.0–1.5 m below the soil surface (respectively, $n=1.0 \mathrm{m}$ and $n=0.5 \mathrm{m}$). Such depth is consistent with the reported depth of the failure surface, which was well above the interface with the bedrock [47]. Although the predicted minimum value of $F_S$ would not imply slope failure, it is small enough to be considered indicative of potential instability. In fact, heterogeneities play a big role in the triggering of landslides, which cannot be completely captured by the simplified modeling approach adopted here. The real slope presents variable inclination, soil cover thickness, and hydraulic and geotechnical soil properties. Furthermore, the considered cover, which is here schematized for simplicity as homogeneous, is actually characterized by soil layers with different properties (Table 1). The presence of layers with low hydraulic conductivity may exacerbate the mechanism suggested by the modeling results, enhancing the moisture gradient required for the infiltrating water to reach the soil–bedrock interface. This would lead to higher water content within the soil profile, and consequently to a smaller factor of safety. Indeed, it has been recently shown that in the pyroclastic covers of Campania, either coarse layers (i.e., pumices) in dry conditions, or fine layers (i.e., weathered ashes) in wet conditions may hinder the infiltration process [8]. Although the detailed description of the infiltration process falls beyond the aim of this simplified approach, the model results shed some light on the concurrence of long-term predisposing processes and short triggering mechanisms for landslide initiation.

4. Conclusions

The triggering of rainfall-induced landslides in shallow soil covers depends on the characteristics of the triggering rainfall, as well as on the soil moisture conditions at the onset of the event. In this paper, the seasonal dynamics of soil moisture typically observed in the slopes of Campania (southern Italy), which were characterized by shallow pyroclastic covers laying upon a fractured and karstified limestone bedrock, is linked to the fluctuations of the water level in a perched aquifer located in the upper part of the bedrock (epikarst). A simplified physically-based mathematical model of the coupled interaction between the unsaturated flows within the soil cover and the saturated flow in the aquifer is proposed. The parameters describing the hydraulic properties of the soil cover are chosen based on the results of previous laboratory and field investigations. The hydraulic parameters of the epikarst aquifer are assigned according to literature indications, and to the results of currently ongoing field monitoring.

The model has been applied to simulate the hydrological response to precipitations of the slope of Cervinara (about 40 km NE of Naples). The simulation results of the period 2006–2017 show that the model is capable of realistically reproducing the typical seasonal soil moisture trends observed in the area. The model clearly identifies the different response times to precipitations of the soil cover storage, which are rapidly affected by rainy days, and of the perched epikarst aquifer water level, which changes in response to precipitations at a monthly scale. This allows the distinction between the predisposing causes of landslides, i.e., the relatively slow processes of accumulation of water in the perched epikarst aquifer (and its effects on the hydraulic behavior of the soil–bedrock interface), and the triggers, which are represented by single intense rain storms that are too short for allowing effective drainage from the slope cover. The identified cause–trigger interplay offers a sound and convincing interpretation of the landslide that occurred on 16 December 1999.