Analytical Solution for the Pressure Oscillations Caused by Trains of Solitary Waves within Confined Coastal Aquifers

An exact analytical solution is proposed for the pressure oscillations within deep coastal aquifers under the action of tidal level time-variations attributable to train of solitary waves originating off-shore. The purpose of the study is to relate the characteristics of the response of the system to amplitude, steepness, and asymmetry of the soliciting waves, in order to assess its vulnerability to events like violent seaquakes and consequent tsunamis. The time needed by the forcing perturbations, approximated by consecutive triangular impulses, to attain their maximum is assumed to be always smaller than the aquifer diffusive time, in order to evaluate the consequences of the sudden raise of water level along the shoreline, typical of those quite extreme phenomena.


Introduction
Common waves propagating on the sea surface are a direct consequence of the wind.When it starts to blow, any even weak irregularity in the flow pattern will produce corresponding local variations of the atmospheric pressure on the water table, triggering off the formation of series of ripples amplified and pushed away by the wind itself.That phenomenon of gradual growth and associated coalescence can lead, off-shore and under the action of particularly intense fluxes, to waves of considerable height, which in extreme cases can measure up to 8 meters in the Mediterranean sea and over 15 meters in the oceans.Due to their very slow attenuation, mainly attributable to the effect of the friction, those sea waves keep on propagating when the forcing mechanism stops, reaching zones that can be even very far from the origin.Totally different is the genesis of the solitary waves, also known as "tsunami" (a Japanese term meaning "port wave"), which are induced by impulsive phenomena like submarine earthquakes (e.g., [1,2], see also Figure 1) or the fall in the sea of big rock and ice masses.
The tsunami waves travel on the sea surface according to peculiar laws, with almost no decay, over distances of thousands of kilometres.Their rate of propagation (Lagrange velocity) is given by v = g y, (1) where y indicates the local sea depth and g is the specific gravity.For instance, a sea depth of 5 km would correspond to a tsunami velocity of about 800 km/h (comparable to the regime speed of an aircraft).Characterized by a limited amplitude when they propagate in open sea (about 1 m), tsunamis' waves become gigantic water walls (10-30 m) while breaking on the coast, with predictable ruinous effects, which sometimes can be even catastrophic.Two eloquent examples of the tsunami destructive power are represented by the Lisbona seaquake in 1755 and the event determined by the 2004 Sumatra Island earthquake, which released an enormous quantity of energy (about 3 × 10 18 joule) and caused the death of more than 200,000 people in different countries along the coasts of the Indian ocean.The sudden and disastrous increment of the "port-wave" amplitude near the shoreline can be easily understood based on (1): where the sea depth diminishes, waves velocity diminishes too, but due to the weak dissipation and the almost total conservation of the energy carried out by the perturbations, their amplitude must compensate the deceleration.The purpose of present study consists in the analysis of the effects induced by trains of solitary waves-like those produced by a tsunami-in terms of pressure oscillations within a deep confined coastal aquifer, connected to the sea along one of its free boundaries, in order to assess the impact of the phenomenon on the groundwater dynamics according to the degree of steepness and asymmetry of the forcing perturbation.Elementary solutions of the diffusion equation for the piezometric head in porous formations were proposed by Carslaw and Jaeger [3]; applicative contributions, essentially focusing on the drainage of unconfined aquifers, came from Hooghoudt [4,5], Edelman [6], Charni [7], Glover [8][9][10], Dumm [11,12], Von Schilfgarde [13,14], Brooks [15], and Hammad [16].An exhaustive review of the results is discussed in Bear [17].A variety of solutions for confined and unconfined aquifers subject to harmonic oscillations of the sea level were provided by Williams et al. [18], by resorting to mathematical, hydraulic, and electrical analogical models.The analysis of the propagation of periodic water-table waves in unconfined aquifers through an extension of the firstorder Dupuit-Forchheimer approach is presented by Nielsen et al. [19].A fully nonlinear, diffusive, and weakly dispersive equation is derived by Liu and Wen [20] to describe the propagation of gravity waves in shallow porous media through asymptotic expansions of the target unknowns in terms of formation depth normalized with respect to the horizontal length-scale, as well as in terms of free-surface fluctuations normalized with respect to the average water level.Further examples of analysis of aquifers response to tidal fluctuations can be found in Serfes [21], Liu [22], Jiao and Tang [23].Recently, Wang et al. [24] have analyzed the water well resonance induced by preearthquake signals in a confined aquifer resorting to an analytical linearized approach.In the present paper, an exact analytical solution is proposed for the pressure waves propagating within a homogeneous confined coastal aquifer as a consequence of a train of generally asymmetric, triangular-like sea waves approximating the sequence of impulses to which can be reduced a train of solitons (e.g., El et al. [25]), conventionally assuming an indefinite flow domain along the shoreline.

Mathematical Formulation
The purpose of present work is represented by the derivation of a solution in terms of time-dependent piezometric head h = z + p/γ (with z indicating the altitude, p the pressure, and γ the specific weight) within a confined, homogeneous, uniform thickness coastal aquifer delimited by an inland constant-level basin which is hydraulically tributary of the sea (Figure 2), during the breaking of a train of solitary waves.The porous formation, characterized by conductivity K and specific storage S, is considered unbounded along the transverse direction and large L in the longitudinal one.Note that the assumption of uniform thickness makes sense when dealing with big aquifers, whose horizontal dimensions dominate its depth.On the other hand, the homogeneity of the porous matrix can be justified considering that the coastal formations are usually located below wide alluvial plains or estuaries characterized by fine and quite homogeneous sediments.The effect of the tidal oscillation at the shoreline (x = L), is approximated by a series of triangular impulses (Figure 3); indeed, the solitary wave, experimentally studied for the first time in the nineteenth century, is a very peculiar case of wave, fully elevated above the undisturbed water level, that can cover very large distances in open sea without appreciable deformations and decay.The analysis will be carried out in the x-h plane, with reference to the generic cross-section.Governing diffusion equation (e.g., [17]) is then expressed by (2) associated with the following initial and boundary conditions: where F(t) indicates the external forcing, h 0 − h 1 the steady difference of level between the inland basin and the undisturbed sea surface, and the origin of the reference frame is placed in correspondence of the latter.The problem will be handled in dimensionless form assuming suitable space and time scales, λ s and λ t .A reasonable and almost imposed choice for λ s is represented by the porous formation length.About λ t , in order to simplify the resulting dimensionless governing equation, one can assume or with D representing the porous medium diffusivity.Notice that the selected time-scale turns out to be coincident with the characteristic aquifer dissipative time.With reference to the only first impulse or "triangle", the three subsystems to be solved are (for the sake of simplicity, the dimensionless variables will be indicated with the same symbols used for the dimensional ones).
(A) For 0 ≤ t ≤ θ, with θ = T max /λ t : with G indicating the temporal asymmetry index of the tidal wave at x = 1 (see Figure 3, where G is conventionally assumed >1), h max the peak, and T max the time needed by the perturbation to attain it.Note that for G = 1 raising and decay period have the same length, whereas G / = 1 may account for the specific generation mechanism of the solitons and the orientation of the fault plane or for the effect of the forthcoming wave breaking.Furthermore, the detrending of the h B and h C initial condition has been performed in order to avoid the superposition with the effect of the permanent upstream boundary constraint.
Resorting to the Fourier representation of the space functions for confined domains, the following piecewise solution is found in the unified time reference frame: where Due to the linearity of the governing equation, the extension of the basic solution represented by ( 9) and ( 10) to the whole series of perturbations is straightforward: where the subscript 0 indicates the above solution, N indicates the total number of solitons, h maxs = R s h max , and G s , θ s , and R s indicate the coefficient of asymmetry, the dimensionless raise length, and the reduction coefficient of the single wave following the leading component, respectively; plus, with τ i expressing the time-lag between consecutive impulses.

Discussion
Figures 4-9 show the solution at fixed times and positions through the aquifer for N = 2 and R 1 = 1, and three different combinations of parameters (θ 0,1 = 0.01; G 0,1 = 1; h max 0,1 = 10), (θ 0,1 = 0.01; G 0,1 = 3; h max 0,1 = 5) and (θ 0,1 = 0.1; G 0,1 = 1; h max 0,1 = 1), respectively.The steady difference of inland basin-undisturbed sea surface level h 0 − h 1 has been set at 0.1.The time-lag τ 1 between the two waves of the pair has been assumed equal to 2/5 of their length, while the total simulation time is given by twice the length plus the timelag.Practically, the two solitons are identical in shape with no appreciable decay (that could be the case for the first impulses determined by particularly intense earthquakes, with the only cause of attenuation, related to the porous matrix capacity of dissipation, affecting their propagation through the aquifer).The three combinations of parameters characterizing the forcing perturbation have been selected in order to keep the input waves area constant and, therefore, to represent sequences of tidal time oscillations characterized by the same budget of energy.The only differences accounted for are represented by the sea depth near the beach, making the disturbances more or less steep (indeed, as discussed above, shallower waters correspond to narrower and steeper perturbations) and by the mechanism of wave generation, which can affect its degree of asymmetry.As one can see, everything else being the same, a steeper solicitation reduces the portion of the aquifer subject to the increment of pressure (Figures 4, 5, and 6) and the part of it involved into the drainage process toward the sea, which goes from the local maximum of the single curve to the shoreline (x = 1).That kind of behavior is definitely in contrast with that exhibited by nondissipative systems governed by hyperbolic partial differential equations.As an example, consider the somehow similar case of the elastic oscillations within the terminal part of a hydroelectric power plant, where the pressure perturbations propagating through the forced conduit, kept at a constant piezometric head in correspondence of its upstream extreme section, are determined by the gradual plugging of a valve.In the classic approach to the problem, the fluid is assumed to be nonviscous and the two coupled governing equations are where a indicates the rate of propagation of the perturbation and U the section averaged conduit velocity.For a fast closure, which is completed in a time smaller than that needed by the elementary disturbance originating at the valve to cover there and back the conduit, the portion of it interested by the maximum overpressure is given by where L c is the conduit length and t represents the duration of the closing manoeuvre.Thus, faster the blocking of the flux, larger the part of the system that will undergo the strongest solicitation.The response of the aquifer under the conditions examined in the present study may be explained through the virtual increase of its diffusivity (the ratio D = K/S) as inverse function of θ, for θ < 1 (or T max < L 2 /D).That conclusion is strengthened by the analysis of Figures 7, 8,  and 9, showing the piezometric oscillations in time at x = 0.8, 0.5, and 0.2.It is clearly seen that steeper input waves enhance the porous medium tail effect as a consequence of the relative reduction of the storativity and, therefore, of the deceleration of the drainage phase.Moreover, shorter raising periods are correlated with more pronounced shifts between the peaks at different locations (for θ = 0.01 and G = 1, at x = 0.5 and 0.2, there are no peaks at all, at least for the given simulation time).The asymmetry of the sea wave does not seem to affect the degree of symmetry of its porous medium counterpart; on the contrary, a longer decay phase (G > 1) induces more peaked, and relatively less tailed piezometric head oscillations (see Figure 7 as compared to Figure 8).Finally, a steeper input wave induces a relatively smaller porous medium first peak at a fixed location in time: for θ = 0.01 and h max = 10, at x = 0.8 that first peak is about 20% of h max , whereas for θ = 0.1 and h max = 5 the percentage goes up to about 55% (see Figure 7 as compared to Figure 9).All the above mentioned characteristics have well-defined consequences in terms of global response of the system to a train of solitons; indeed, looking at Figures 7, 8, and 9 one can easily infer the hysteretic effect of the periodicity (or pseudoperiodicity, if τ s does not keep constant), manifesting itself through the progressive decay of the relative maxima and the presumable large-time saturation (constant increase without oscillations).That means that the reduction of the porous medium storage capacity during each cycle is not automatically recovered at the end of it and before the next, with a typical cumulative effect which is more pronounced for steeper forcing shock-waves.
When analyzed in terms of environmental impact, the proposed solution shows that the coasts surrounded by shallower waters, and virtually associated with relatively steeper solitary waves are characterized by a more elevated risk of overpressure going beyond the resistance limits of the aquifer confining layers in the near-shore zone; nevertheless, the presence of deeper coastal contours, associated with apparently more innocuous and flatter surface waves, determines a considerable risk in terms of strong overpressures as well as in terms of salt/pollutants intrusion even for the inland zones, far from the shoreline.

Conclusions
The paper has dealt with the analytical treatment of the pressure oscillations taking place within deep homogeneous coastal aquifers and caused by the breaking of a series of solitary waves.The aim of the study was represented by the assessment of the impact that catastrophic natural events like seaquakes and consequent tsunamis produce on the groundwater dynamics and, therefore, on the subsurface stability as well as on the water quality preservation.The simplifying hypotheses underlying the mathematical  formulation concern the shape of the forcing time-function, assimilated to a series of triangular impulses characterized by different dimensions and different degrees of asymmetry, in order to account for the progressive attenuation of the peaks of the train of solitons, as well as for the role played by the specific genesis of the waves and the orientation of the deep fault plane.The resulting plots obtained for the first two components of the series in absence of appreciable decay, and drawn at fixed locations in time and at fixed times in space, show that: (i) Steeper the soliciting perturbation (and, therefore, shallower the near-shore contour), smaller and more peripheral the portion of the aquifer interested by conspicuous overpressures and involved in the process of sea-to-aquifer flow, with associated salt/pollutants intrusion.
(ii) Steeper the soliciting perturbation, slower the drainage phase, and more pronounced the hysteretic effect due to the apparent reduction of storativity.
Those characteristics are in definite contrast with what happens in nondissipative systems governed by hyperbolic equations (see, e.g., the case of the forced conduit of an hydroelectric power plant subject to the more or less fast closure of the downstream valve, traditionally handled assuming the fluid as nonviscous), in which the length of the domain affected by the maximum disturbance is inversely proportional to the duration of the forcing perturbation.
Finally, everything else being the same (same budget of energy carried out by the train of solitons), the degree of temporal asymmetry of the tidal oscillations does not seem to influence their porous medium counterpart, which exhibits a quite regular shape and a relatively smaller tail in case of asymmetrical impulses.
When analyzed in terms of environmental impact, and due to the inverse proportionality between steepness of tidal waves and local sea depth, the proposed solution shows that the coasts surrounded by shallower waters are characterized by a more elevated risk of overpressure confined to the nearshore zone, while the presence of deeper coastal contours determines a considerable risk in terms of strong overpressures as well as in terms of salt/pollutants intrusion even for the inland zones.Once the characteristics of the expected, most catastrophic events in a given place were statistically known, (9) to ( 12) could be used in order to predict the peaks reached by the overpressures soliciting the confining layers of the involved aquifer, as well as to determine the extent of it potentially free from the sea contamination and, therefore, safely usable for potable/agricultural destinations without need for remediation interventions.

Notation D:
Aquifer diffusivity F(t): Forcing sea wave g: Specificgravity G: Forcing wave asymmetry index K: Aquifer hydraulic conductivity h: Piezometric head h 0 − h 1 : Difference of level between inland basin and undisturbed sea surface h max : Peak of the forcing wave L: Aquifer length p: Pressure R: Coefficient of attenuation for consecutive impulses S: Aquifer specific storage t: T ime T max : Time at which the forcing peak is reached v: Solitons rate of propagation x: Longitudinal coordinate y: Local sea depth z: Aquifer points altitude γ: Specificweight λ s : Lengthscale λ t : Timescale θ: Dimensionless T max .

Figure 2 :
Figure 2: Sketch of the aquifer system (the vertical scale has been exaggerated).