A novel approach to model the hydrodynamic response of the surficial level of the Ionian multilayered aquifer during episodic rainfall events

The paper presents a modeling framework to analyze the effect of episodic rainfall supply on groundwater dynamics in the Ionian coastal plain multilayered aquifer. The focus is essentially on the short-term behavior of the shallowest layer, with a specific analysis on episodic rainfall events. In the studied aquifer, groundwater level responds sensitively to rainfall events, highlighting the presence of preferential recharge zones. The hydraulic head peak is a function of groundwater level antecedent to the rainfall event. A kinematic dispersion wave model was used to model infiltration processes via preferential pathways. A one-dimensional and non-linear particle based numerical model was developed. Particles with constant water volume travel according to celerity and hydraulic dispersion to simulate the infiltration rate wave through the vadose zone. The flow rate that reaches the water table represents the input function to determine groundwater level fluctuations along groundwater flow direction and according to the lithological features of the surficial levels of the multilayered aquifer, its storage capacity changes passing from unconfined to confined conditions. The model was validated with high time resolution time series of precipitation and groundwater level coming from Terra Montonata meteo-climatic and groundwater level monitoring station. The developed model represents a basis for evaluating and predicting groundwater discharge of the shallowest layers of the Ionian coastal multilayered aquifer under natural conditions including episodic rainfall.


Introduction
The understanding of hydrological behavior of an aquifer is fundamental to achieve a sustainable management of the groundwater resource. An increased knowledge of the groundwater dynamics and water balance is important for the prediction, evaluation, and control of the groundwater evolution from both qualitative and quantitative point of view (Brunner and Kinzelbach, 2008;Cherubini and Pastore, 2011;Lee et al., 2012). In shallow aquifers, where the water table quickly responds to water inputs, recharge events may be isolated and associated with individual precipitation events (Meinzer, 1923;Healy and Cook, 2002).
The term preferential flow usually refers to three physically distinct processes: macropore flow, finger flow, or funnel flow (Kung, 1993;Ogawa et al., 1999). Finger flow or fingering, results from air entrapment or water repellence causing an instability in the wetting front and its breakup into fingers (Glass et al., 1988, Wang et al., 1998, Hardie et al., 2011. Hill & Parange (1972) found that in a layered soil whose upper layer is finer and less conductive than the coarser layer beneath, the wetting front becomes unstable and breaks into narrow wetting columns, or "fingers," that move downward at a velocity equal to the saturated conductivity divided by the saturated volumetric water content of the bottom layer. Periodical intermittent flow events of persistence occurrence can occur periodically.
Threads along the preferential flow channel would snap, drain, and then reform again (Su et al.,1999).
The kinematic diffusion wave theory is a useful approach to model the vertical movement of the moisture under uniform and/or preferential flow (Germann and Beven, 1985;Germann, 1990;Rasmussen, 2013;Rossmann et al., 2014). Some authors derived a kinematic diffusion form of Richards equation, mathematically equivalent to the advection diffusion equation (Rasmussen et al., 2000;Bücker-Gittel et al. 2002;Cata and Mohrlok, 2006;Zehe and Jackish 2016;Pastore et al., 2017). Nimmo (2010) introduced an additional term to the advection diffusion form of Richards equation in order to consider preferential flow. Alternative approaches based on Newtons law found a power law relationship between the infiltration rate and mobile water content within the draining porosity (Germann, 1990;Chen and Wagenet, 1992;Germann and Karlen, 2016). The combination of this functional relationship and the continuity equation gives rise to the kinematic wave model (Germann and Beven, 1985). The latter is strictly convective and leads to an overestimation of the infiltration rate. A dispersion effect that attenuates the kinematic waves is observed (Di Pietro and Lafolie, 1991;Di Pietro and Germann, 2015). It is associated to several factors such as the contribution of mesopores in which capillary forces may be significant, as well as the presence of intricate pore paths that produce water flow dispersion within the draining porosity. The introduction of the dispersion component gives rise to the kinematic dispersion wave model (Di Pietro et al., 2003;Moradzadeh et al. 2020). Note that kinematic diffusion model and kinematic dispersion wave model are formally equivalent. Both present convective and diffusion/dispersion terms. The difference lies in the constitutive functions used to represent convective and diffusion/dispersion functions.
The implementation of preferential flow into the numerical model as well as the determination of model parameters remain the subject of discussion (Jarvis, 2016). Modeling of the infiltration processes through the vadose zone and the quantification of the recharge hydrograph of the aquifer are strictly connected with the basic characteristics of unsaturated and saturated zones (Šimůnek et al., 2003;Diamantopoulos et al., 2015). The assessment of groundwater recharge requires a numerical model able to represent all key process controlling the infiltration process and with a small computational cost.
In this paper an improved modeling framework and the relative numerical approach has been developed to represent infiltration processes, soil moisture dynamics, aquifer recharge and groundwater level hydrographs during episodic rainfall event. The infiltration processes are modelled like a one -dimensional flow along vertical direction using the kinematic dispersion wave model.
Power law constitutive relationships between the infiltration rate and the mobile water content have Preprints (www.preprints.org) | NOT PEER-REVIEWED | Posted: 2 August 2020 doi:10.20944/preprints202008.0039.v1 been used. Furthermore, dispersive attenuation of kinematic processes has been taken into account.
Generally, the kinematic dispersion wave equation is solved numerically by means of Eulerian methods (Di Pietro et al., 2003;Moradzadeh et al., 2020) that typically encounter numerical problems such as numerical dispersion and they suffer of computation burden when dealing with strong variation of boundary conditions. During intense rainfall events, the infiltration rate, the water content as well as the water table fluctuation change rapidly giving rise to numerical instabilities in Eulerian numerical model. In order to overcome these difficulties a nonlinear particlebased model has been developed to numerically solve the kinematic dispersion wave equation. Particles move according to convection and dispersion terms which are function of the particle density. The infiltration rate hydrograph that reaches the water table is determined and it has been used to simulate groundwater level fluctuation of the surficial aquifer according to its storage coefficient. The length of the vadose zone changes on the basis of the water table depth.
The developed method is applied to the case study of the surficial level of the Ionian coastal plain multi-layered aquifer (Italy). Previous studies carried out in the same aquifer, on the basis of monthly average precipitation and groundwater levels, show that recharge takes place mainly from the upward aquifer, whereas direct recharge is interdicted by the presence of surficial silty and clay layer (Polemio et al, 2003;Fidelibus et al., 2004;Doglioni et al., 2011). The analysis of precipitation and groundwater records at higher resolution scale puts in evidence a quick response to intense rainfall events that are not coherent with lateral or upward recharge supply. This implies the presence of direct recharge.
Then the quick groundwater level fluctuation response may be favored by the presence of several reclamation channels crossing the coastal plain which cut the surficial impervious layer creating a preferential flow path towards the aquifer. Moreover, according to the lithological features of the surficial aquifer, unconfinedconfined flow conversion may occur giving rise to a change in storage.
Hydrogeological parameters have been estimated according to the comparison between modeled results and observed data.
The specific goals of this paper are:

Geological, lithological and hydrogeological setting
The study site is located in a wide and flat area called Ionian (Metaponto) coastal plain with highest elevation of 15 m, facing in the Ionian Sea in the southernmost part of Basilicata region (Southern Italy) within the interfluvial area between River Cavone and River Basento, which represent two of the five main rivers (Agri, Cavone, Basento and Bradano) that run parallel to each other before reaching the shoreline (Fig. 1a). The study site lies in the southern part of the Bradanic Trough (Boenzi et al., 1976) (Fig. 1b). The Bradanic Trough is filled up by a thick Pliocene to Pleistocene sequence mainly of Sub-Appennines Clays (Cherubini et al., 1989;Cherubini and Lupo, 2002;Galeandro et al., 2017), passing upwards to coarse-grained coastal and continental deposits known as "regressive coastal deposits of the Bradanic Through" giving rise to the well-known staircase of the marineterraced deposits of the Metaponto. Magri (1967, 1969), detected seven orders of terraces. Successively, other authors identified more terraces -up to eleven-and that each marineterraced deposit is not simply related to a single transgressiveregressive cycle, but it presents an internal cycle of higher frequency (Boenzi et al., 1976;Brückner, 1980). Guerricchio and Melidoro (1986), Bentivenga et al. (2004) and Doglioni and Simeone (2020), have shown that this lower part of the trough is crossed by tectonic structure and distensive faults. The lower level terraced deposits correspond to a part of marine -terraced deposit of the seventh order (Vezzani, 1967), which were recently named by Cilumbriello et al. (2010) as "Sabbie e Conglomerati di Policoro" unit and in part "Sabbie e Conglomerati di Masseria Ricotta" unit.
Alluvial, transitional, coastal and marine deposits characterize the plain. Several authors (Polemio et al., 2003;Pescatore et al., 2009;Cilumbriello et al., 2010)  The upper unit presents a thickness of 30 m deepening locally in correspondence of the paleovalley and is characterized by offshore transition silty-clay, then to deltaic silty-sands and fluvial sands.
According to Polemio et al. (2003) and Cilumbriello at al. (2010), the lithological features of the study area gave rise to the presence of two main aquifers. The former contains the aquifers of marine terraces and their continuity is interrupted by the presence of the river valleys. Groundwater flows generally in unconfined conditions, in the coarse grained deposits of the terraces. The direct recharge is ensured by the limited and nonhomogeneous extension of low hydraulic conductivity levels placed mostly at the aquifer bed. The latter is an aquifer hosted in the coastal plain deposits, where groundwater flows in permeable level. Figure 2a shows the schematic cross geological section of interfluve area. Figure 2b shows the strata sequence detected in correspondence of the Terra Montonata monitoring station.
The multilayered aquifer lies on clays and silty clays of the Subappenine Clays (Argille Subappennine) formation ( Fig. 2a). Close to the coast and along the strike direction the lower part of the aquifer system is partitioned, and groundwater flows in the paleovalleys within sandy silty and sandy gravelly layers which alternate with low permeable silty layers. Moving landward, sandy silty and sandy gravelly layers become thicker giving rise to a continuous permeable porous body. Differently they become thinner moving seaward up to termination within shelftransition clays and silt clays. Then, inside the paleovalley and close to the coast the multilayered aquifer shows artesian features. Above the paleovalleys, groundwater splits into several permeable layers and the shallowest of which can be characterized by free piezometric oscillations. As showed in Figure  Well test analyses conducted in the whole Metaponto coastal plain by Polemio et al. (2003), show hydraulic conductivity values in the range 3.47×10 -6 -5.69×10 -3 ms -1 with an average value of 2.28×10 -4 . The hydraulic conductivity increases from the coast to the inland while moving from Bradano to Sinni catchments.
The increase of agricultural activities during the last decades encouraged the overexploitation of groundwater from coastal aquifer system. Fidelibus (2004) estimated a value of ten m 3 y -1 of groundwater pumped from the wells drilled by the Consortium for Reclamation of Bradano and Metaponto. Anyway, private and illegal wells are also present in the area, extracting groundwater from the coastal aquifer especially during the dry season. The occurrence of many years of drought together with the overexploitation of the groundwater resource led to qualitative and quantitative deterioration of this water resource.

Hydrodynamic observations
The present work refers to time series data collecting by Protezione Civile Basilicata -Centro  A time lag of 3-4 months exists between the monthly average rainfall and groundwater level peaks.
It shows that a large part of groundwater supply is due to lateral/upward recharge probably through the high permeability marine -terraced deposits of Conglomerati di Masseria Ricotta unit more lithologically and hydraulically connected with the surficial aquifer.
Anyway, the analyses of precipitation and groundwater level at daily time scale evidence the presence of an impulsive response to intense rainfall. Water table responds immediately to intense rainfall events. This implies there are quite short paths of recharge along the vadose zone that contribute to groundwater supply. This behaviour is evident from autumn to spring period where the more significant rainfall events in terms of intensity, magnitude and amount of precipitation are present. Figure 4 and Figure 5 show four daily time series during the wet season (1 th October -30 th April) in which the quick response of groundwater level due to intense rainfall is evident.
These quick responses are characterized by two aspects. The first aspect concerns the fact that similar rainfall events in terms of cumulative precipitation give rise to different groundwater level fluctuations in terms of peak values. In other words, the surficial aquifer presents a different hydraulic behavior during the rainfall event depending on the initial value of groundwater level. This behavior seems to be interrelated to the depth of the bed of the silty and clay unit (gray bar on Figure 4 and The second aspect regards the fact that groundwater responds very sensitively. The time lag between significant rainfall event and the relative groundwater fluctuation is about 2-3 days. This behavior is not coherent with a groundwater supply only due the high permeability marine-terraced deposits.
Direct recharge is also present.
However direct recharge is not the main recharge mechanism of the surficial aquifer. This is demonstrated through the analysis of the master recession curves of the significant episodic rainfall event. Groundwater level falling after rising due to significant rainfall events decreases exponentially according to the following recursion expression (Heppner and Nimmo, 2005;Crosbie et al., 2012): is the groundwater level that would attain if no preferential flow recharge occur,

Conceptual model of aquifer during episodic rainfall events
The key observations described in the previous sections permit to implement the following conceptual model of the hydrogeological processes that involve the surficial aquifer in the study area (Fig. 6).
Long term lateral recharge combines with direct recharge that takes place through the preferential flow mechanism in vadose zone. Preferential flow occurs under various conditions at different spatial and temporal scales. Moreover, hydromechanic processes at larger scale are the consequence of the ones at pore scale. In the study area, two class of macropores are observed at the pore scale: biopores created by earthworm and shrinkage cracks formed by drying of swelling phenomena. Earthworms presents a diameter of 1-3 mm and can extend up to 6 m in vertical length. Shrinkage cracks give rise to a rapid movement of depth of around 1 m (Hardie et al., 2011). As consequence, water accumulated backfilling the void space. At darcian scale preferential flow mechanisms can be due to unstable flow.
Under certain conditions, the wetting front moving downwards breaks into fingers. They represent preferential pathways that facilitate recharge flow. Compression of air below the wetting front causes instability (Wang et al.,1998) and the wetting front breakup into fingers.
At areal scale the numerous reclamation channels in the study area represent local surface depressions where surficial flow concentrates and is transmitted to the underlying shallow aquifer.
Preferential flow recharge mechanism is modelled as onedimensional flow in z direction through the vadose zone. Given that the interfluvial area is mainly covered by arable crops with a topsoil more or less of 0.30 m and the investigate period regards only the wet season, evapotranspiration processes have not been taken into account.
When the precipitation p(t) [LT -1 ] exceeds the infiltration capacity of the topsoil qc [LT -1 ], overland surficial flow is generated causing the infiltration rate at top soil q(0, t) [LT -1 ] along the preferential flow paths. Moreover, the infiltration rate at topsoil is limited above according to the following expression: Where qcrit [LT -1 ] represents the critical value of the infiltration rate, the precipitation in excess which does not feed the preferential flow paths (Nimmo, 2010).
The infiltration rate travels along the vadose zone reaching the water table giving rise to the infiltration rate hydrograph at water table q(L, t) . L [L] is the water table depth from topsoil.
The total rate of the ground water level fluctuation is the sum of the recession and accretion rate (Nimmo, 2010): Where H [L] is the height of the groundwater level above the base level

Kinematic dispersion wave model
Kinematic dispersion wave model is used to determine the infiltration rate hydrograph at water table.
Fluid flow through the one -dimensional variably unsaturated medium occurs through preferential pathways, no exchange exists between impervious stratum and the preferential pathway. Assuming water density as constant and neglecting the inertial terms in the momentum equation, according to Di Petro et al. (2003) the conservation laws can be written as: Where w D [L 2 T -1 ] is the hydraulic dispersion and w c [LT -1 ] is the celerity.
The infiltration rate wave propagates through convection and dispersion processes governed by the changes of the mobile volumetric water content along the preferential pathways. Celerity is defined as the derivative of the infiltration rate respect to the mobile volumetric water content: Whereas hydraulic dispersion depends on the celerity and the water dispersivity coefficient: Neglecting the second term on the right side of equation 5, the celerity can be written as: The infiltration processes depend on three parameter a, b and α. The following initial and boundary conditions have been imposed: Equation 6 has been solved numerically using nonlinear random walk method. The Numerical scheme is described in the next section. The program was written in Matlab environment.

Numerical scheme
Randomwalk method uses particle tracking to solve kinematic process whereas dispersion processes are simulated by adding random displacement to each particle in addition to advective displacement. It is known that the spacetime distribution of particle can be represented by the Fokker -Planck equation which is not identical to Equation (6). In analogy to solute transport problem, the celerity term cw is replaced by: For each time step t the specific volume of water which enter in the subsoil at z = 0 is . The accumulated specific water volume W [L] will be updated as t t t in W W w + =+. The volumetric water content at z = 0, according to the Equation (5) will be equal to: Then N particles having each a specific water volume ( ) 0, / tN  are released at z = 0. The depth zi of each i-th particle at time t+t is updates as: Where Z is a normally distributed random number, cw,i' and Dw,i represent the celerity and the hydraulic diffusivity associated to i-th particle. They are function of the infiltration rate which change along the depth.
The onedimensional domain along the depth between the topsoil and the water table of length L is discretized with n cells. The infiltration rate and the volumetric water content is assumed constant in space within each cell.
For each cell j (j = 1,…,n) the specific volumetric water content is determined as: Finally, the infiltration rate that reaches the water table q(t+t, L) is determined as: The groundwater level fluctuation is determined using the recursion form of equation 2: Finally, the domain length L is updated:

Preliminary analysis
Groundwater level time series at daily scale show a quick response of about of 2-3 days at significant rainfall event. On the other hand, average monthly time series highlight a long-term recharge behavior with a peak delay between average monthly groundwater levels and precipitation of 3-4 months. This long-term recharge behavior is due to lateral recharge processes, in which marine-terraced deposits feed the surficial aquifer. In order to validate this thesis, an estimation of groundwater wave propagation from the marine-terraced deposits to the monitoring station can be made. The Groundwater flow in surficial aquifer can be represented as a one-dimensional flow: Where T [LT -2 ] is the transmissivity of aquifer and x is the horizontal direction. Marine -terraces deposits present a minimum distance from the monitoring station equal to x0 = 1200 m. The coastal plain aquifer is characterized by an average hydraulic conductivity value of K = 2.2×10 -4 ms -1 with an average thickness of 20 m (Fig 2). Then the transmissivity will be equal to T = 4.4 ×10 -3 m 2 s -1 .
Assuming that the coastal aquifer is subject an instantaneous unit rise pulse at x = 0, groundwater level wave will be propagated according to the following analytical solution for diffusion: The average value of storage coefficient S can be determined as the ratio between the cumulative infiltration rate at topsoil and the cumulative accretion rate at the end of the rainfall event: In order to determine the cumulative infiltration, a value for the infiltration capacity qc and the critical rainfall rate qcrit must be assumed. According to the permeability of the surficial thin clayey matrix, qc has been assumed equal to 0.5 mmh -1 . qcrit is estimated equal to 6 mmh -1 , representing the value that allows to have a value for the average storage coefficients S (determined for each rainfall event) closest to 0.2.
The calculated values assumed by S put in evidence a fundamental aspect of the episodic recharge behavior in the study area (Fig. 7). S reduces systematically when the groundwater level exceeds the bottom of the silty clay unit (5.0 m AMSL) reaching a value of 0.07 for the event #2, 0.08 for event #4 and 0.06 for the event #6. For the event characterized by a groundwater peak below the bottom of the impervious layer the storage parameter reaches a value in the range 0.11 -0.16. Figure 8 shows the comparison between the cumulative precipitation and the accretion rate for the event #1 and event #2.
Successively, for each event the lag time between cumulative infiltration rate at water table and cumulative accretion rate as the difference between the respective residence time is estimated: The calculated lag times confirm the presence of direct recharge mechanism in the study area through preferential pathways. They vary in the range between 2.19 -9.11 d, whereas the respective value of the average celerity w c obtained as the ratio between the average depth of the water table and the time lag assume a values in the range between 491.77 -2267.12 mmd -1 .

Comparison of model prediction with observed data
The kinematic dispersion model has been used to predict groundwater level fluctuation starting from the hourly precipitation time series. Simulation results have been compared with the observed groundwater level time series. The kinematic dispersion equation has been solved using the developed particle -based model in order to determine the infiltration rate hydrograph at water table q(L, t).
Then groundwater level fluctuation H(t) has been obtained. Finally, the simulated groundwater level is determined as the sum of H(t) and the value of groundwater level in absence of direct groundwater recharge Zw0. Note that due to the presence of long-term lateral groundwater recharge mechanism Zw0 varies in time as showed in Figure 4 and Figure 5, whereas the recession time τ is constant equal to 18 d.
In order to control the efficiency and performance of the simulation, the time step t and the grid cell size z for the evaluation of the particle density have been chosen so as to satisfy a small Courant The storage parameter reduces as water table depth decreases according to the estimated step function (Fig. 9).
The fitting results are showed in Figure 10  A comparative analysis between the infiltration rates at the topsoil and water table was made, in order to investigate the effects of the characteristics of rainfall events on the infiltration processes along the preferential pathways. The cumulative infiltration rate curve indicates the amount of water that crossed a certain surface (topsoil or water table) at a specific time. In order words they indicate the time required to reach a determined amount of water that crossed a surface of vadose zone. Then, given the topsoil and the water table surface, it is possible to determine their respective times needed to be crossed by a determined amount of water. Water table surface will have a time lag respect to topsoil surface ( Fig. 12 and Fig. 13).
This time lag represents the travel time of water needed to reach the water table from topsoil. Then, the average velocity of the wetting front (average celerity) can be determined as the ratio between the depth of the water table from the topsoil and the time lag. The average celerity curve is almost specular to time lag curve. In correspondence of a rainfall event, the average celerity increases rapidly according to the rainfall intensity. When a rainfall event passes, the celerity decreases in a potential way reaching a minimum value of 0.25 md -1 . The maximum values reached by celerity are strictly dependent on rainfall intensity, resulting equal to 500 mmd -1 for lower intensity rainfall event and 1800 md -1 for higher intensity rainfall event.
The numerical simulation results are consistent with the findings of the preliminary analysis. Both disclose that the storage coefficient decreases as water table increases according to the presence of the surficial silty clay level. Moreover, the dependence of the celerity on the rainfall event has been highlighted, finding consistent values to the preliminary analysis .

Discussion
The present study investigates the direct recharge of the surficial level of the Ionian coastal multi layered aquifer by means of an improved modeling framework. The analysis of time series at monthly scale (precipitation and groundwater level) follows the expected patterns in which the marineterraced deposits feed the surficial aquifer through a long-term lateral recharge mechanism.
Nevertheless, the analysis at daily scale or less shows a behavior not explainable by this recharge mechanism alone. Long term lateral recharge combines with direct recharge through the vadose zone giving rise to a quick response of the groundwater level.
The study area is susceptible to preferential flow due to different physical mechanisms involving the infiltration processes in the vadose zone at different spatial and temporal scales. The kinematic dispersion model captures the impact of preferential flow mechanism at field scale of site. The model response is governed by three parameters. According to the kinematic theory, the preferential diffusion index a should vary between 2 and 3. According to these two limits the conductance terms b varies in the range between 3.6×10 3 mmh -1 and 3.6×10 4 mmh -1 .
The parameter a and b govern the celerity representing the velocity of the wetting front. As a decreases, b should increase in order to maintain the same order of magnitude of the celerity. A value of water dispersivity αw equal to 200 mm is found. This parameter attenuates the infiltration wave, leading both to an infiltration rate hydrograph at water table and rise limb of groundwater level more distributed in time.
The values found for the conductance term and water dispersivity are consistent with those found by Di Pietro et al. (2003) and Moradzadeh et al. (2020) at laboratory scale. In the present study the maximum infiltration rate is equal to 6 mmh -1 resulting lower than those used by the authors presenting a minimum value of 30.35 mmh -1 . As a result, the conductance term is lower in our study according to the expected theory value. Anyway, the water dispersivity is higher. This is due to several factors linked to the scale dependence of dispersion phenomena. In analogous manner to solute dispersion theory, water dispersion results scale dependent. As the depth of vadose zone increases, the probability that the wetting front moving downwards breakup in more and more fingers increases. On the other hands, the capillary contribution on the wetting front movement can increase attenuating the infiltration rate propagation. Moreover, as the depth increases the conductance of the preferential flow pathways can reduce. As a result, the infiltration rate hydrograph at water table is more attenuated.
The parameter qcrit indicate how much rain, occurring at high intensity, becomes recharge. The critical rate is 6 mmh -1 limiting mainly the significant rainfall event occurred between 2003-2004 and 2004-2005. According to source responsive theory (Nimmo, 2010), qcrit is equal to the product between a constant parameter that assumes a value of 2.1×10 -6 m 2 h -1 at temperature of water equal to 20 °C, and the smallest value of the facial area density Mmin (L -1 ) which characterizes the preferential flow path.
Then in the present study a value of Mmin equal to 2857 m -1 is found. It represents the minimum fraction of the total specific surface area of the vadose zone on which preferential flow takes place. This parameter measures the susceptibility of a field site to preferential flow. Cuthbert et al. (2013) found a value of Mmin between 250 and 750 m -1 at field site in Shropshire (UK). Nimmo (2010) found a value of Mmin equal to 4000 m -1 for Masser site in Pennsylvania.
Local surface depressions due to the presence of several reclamation channels in the study area play an important role on preferential flow mechanism at areal scale (Hendrickx and Flury, 2001). During rainfall events water accumulates in the reclamation channels increasing the opportunity of preferential flow.
The velocity of the wetting front is strictly correlated to the rainfall intensity. Its maximum value results in the range between 1500 and 2000 mmd -1 consistently with the value reported for preferential flow (German and Hansel, 2006;Hincapié and Germann, 2009a;Kung et al., 2006;Hardie et al., 2011). The comparison between the observed groundwater level and simulated ones highlighted a delay on the rising period of the simulated groundwater level with respect to the observed ones that systematically occurs at the first rainfall event of each investigated time series. This can be ascribed to the fact that the value of mobile water content is assumed equal to zero at the beginning of simulation, underestimating the speed of the infiltration wave. Anyway, since the time series begin just after the dry season, the imposed initial condition should not be very different from the real case.
Another explanation can be attributable to the fact that the velocity of the wetting front is greater when the soil is dry in contrast with the theory that higher antecedent soil moisture condition hydraulically actives the preferential pathways (Granowsky et al., 1994;Kung et al., 2000;Jaynes et al., 2001;Jarvis, 2007). However, several authors support the theory that when the soil is dry, preferential flow is more evident (Shipitalo and Edwards, 1996;Lin et al., 1998;Greco et al., 2002;Hardie et al., 2011). Under lower antecedent soil moisture condition shrinkage cracks play an important role in rapid and deep water movement through dry soil. Water backfills the shrinkage cracks increasing the opportunity of preferential flow via preferential flow pathways due to wetting instability due to the compression of air below the accumulated water into the shrinkage cracks. Under higher antecedent soil moisture preferential flow is more dominated by the stable preferential pathways. In this circumstance, lateral flow from macropore to soil matrix reduces, as a consequence the size and number of the preferential flow increases.
When the groundwater level affects less permeable strata, a lower value of the storage coefficient is found due to the presence of aquitards represented by the surficial silty and clay unit that locally confine the shallow aquifer. Changes in storage coefficient have great effect in the amplitude of the groundwater level fluctuation. The conducted analysis disclosed the presence of a transition zone at depth from topsoil equal to 5 m that corresponding to bottom of the surficial silty clay unit detected at the site. Storage coefficient decreases passing from 0.3 to 0.08 evidencing a transitional behavior where confined and unconfined condition coexist. Then the surficial aquifer passes from unconfined condition to weak confined condition and vice versa.

Summary and conclusions
The study presents an improved modeling framework based on the kinematic dispersion theory to simulate water flow through preferential paths and predict groundwater level fluctuation in a surficial aquifer, that occasionally can flow in a confined condition. The prediction accuracy evaluation and comparison indicated that the kinematic dispersion model and its numerical solution by the proposed particle-based model is able to capture the preferential flow recharge mechanism in the study area, showing a good agreement with the observed groundwater level time series.
The local stratigraphic sequence is characterized by a surficial thin silty and clay level covering the loamy sandy levels hosting the aquifer. Therefore, the direct recharge appears quite difficult for the presence of a low permeability level covering the high permeability levels hosting the groundwater.
The recharge could be attributed partly to the coarsegrained deposits associated to Conglomerati di Masseria Ricotta unit outcropping about 1200 m upstream that are hydraulically connected with the coastal aquifer system. Anyway, the response of groundwater level is almost immediate, due to a sort of direct recharge. The study area is subject to infiltration through preferential flow paths due to different physical processes at spatial and temporal scale. At areal scale, the presence of several reclamation channels widespread in the plain cutting the thin layer of low permeability deposits outcropping all over the plain, leading to local topographic depressions that encourage the opportunity of the preferential flow.
Moreover, the conversion from unconfined to weak confined groundwater flow highlighted. When the groundwater level rise above the bottom of the silty and clay unit, weakly confined condition occurs. Storage coefficient reduces from 0.3 to 0.08. These estimated values are coherent with a critical rainfall rate qcrit of 6 mmh -1 .
Recharge via preferential flow is strictly correlated to the precipitation characteristics in terms of duration, magnitude and intensity. In the study area, high intense rainfall events favor direct recharge of surficial aquifer via preferential flow dynamics, but at the same time reduces the travel time of the mobile water in vadose zone increasing the risk of surficial groundwater contamination as a consequence of preferential flow. The comparison between model prediction and observed data indicates that the susceptibility of preferential flow in the study area results relatively high. The results evidences the influence of moisture antecedent conditions on preferential flow mechanisms. Antecedent dry conditions seem to favor preferential flow via shrinking cracks. Anyway, more detailed analysis must be conducted.
Future development of this study will regard the implementation of precise weighting lysimeters in the study area to observe and analyze the preferential flow processes at several depths in vadose zone.
Furthermore, transient well test analysis is recommended in order to deepen the mixed unconfinedweakly confined groundwater behavior.

RMSE = 57.84 m
Preprints (www.preprints.org) | NOT PEER-REVIEWED | Posted: 2 August 2020 doi:10.20944/preprints202008.0039.v1 Figure 12. a) Daily precipitation (mm), cumulative infiltration rate at topsoil (mm) (black curve) and cumulative infiltration rate at water table (mm) (gray curve). b) variation of the time lag (grey curve) (d) and celerity (md -1 ) (black curve) as function of the cumulative infiltration rate. The two graphs permit to highlight the relationship between the rainfall characteristics in terms of amount and intensity.
Preprints (www.preprints.org) | NOT PEER-REVIEWED | Posted: 2 August 2020 doi:10.20944/preprints202008.0039.v1 Figure 13. a) Daily precipitation (mm), cumulative infiltration rate at topsoil (mm) (black curve) and cumulative infiltration rate at water table (mm) (gray curve). b) variation of the time lag (grey curve) (d) and celerity (md -1 ) (black curve) as function of the cumulative infiltration rate a water table. The two graphs permit to highlight the relationship between the rainfall characteristics in terms of amount and intensity.