J Wolanin

L Michel

D Tabacchioni

J M Zanotti

J Peters

I Imaz

B Coasne

M Plazanet
email: marie.plazanet@univ-grenoble-alpes.fr

C Picard
email: cyril.picard@univ-grenoble-alpes.fr

Heterogeneous Microscopic Dynamics of Intruded Water in a Superhydrophobic Nanoconfinement: Neutron Scattering and Molecular Modeling

With their strong confining porosity and versatile surface 2 chemistry, Zeolitic Imidazolate Frameworks -including the 3 prototypical ZIF-8 -display exceptional properties for vari-4 ous applications. In particular, the forced intrusion of water 5 at high pressure (∼ 25 MPa) into ZIF-8 nanopores is of in-6 terest for energy storage. Such a system reveals also ideal to 7 study experimentally water dynamics and thermodynamics 8 in ultra-hydrophobic confinement. Here, we report on neu-9 tron scattering experiments to probe the molecular dynam-10 ics of water within ZIF-8 nanopores under high pressure up 11 to 38 MPa. In addition to an overall confinement-induced 12 slowing down, we provide evidence for strong dynamical 13 heterogeneities with different underlying molecular dynam-14 ics. Using complementary molecular simulations, these het-15 erogeneities are found to correspond to different microscopic 16 mechanisms inherent to vicinal molecules located in strongly 17 adsorbing sites (ligands) and other molecules nanoconfined 18 in the cavity center. These findings unveil a complex micro-19 scopic dynamics, which results from the combination of sur-20 face residence times and exchanges between the cavity sur-21 face and center.

Introduction 23

Interfacial and confined water has become a scientific field 24 in itself at the crossroad of physics, chemistry and biology 25 with significant implication for applied science such as in 26 geophysics, mechanical/fluid engineering, etc. [START_REF] Lynch | Water in nanopores 485 and biological channels: A molecular simulation per-486 spective[END_REF][START_REF] Kavokine | Fluids at the Nanoscale: From Continuum to Subcontinuum Transport[END_REF] Such wa-27 ter, which displays marked differences with its bulk counter-28 part, is of importance for situations where physicochemical 29 interactions with surfaces are at play: protein folding, bio-30 logical nanochannels, nanofiltration/phase separation, catal-31 ysis, etc. [START_REF] Lynch | Water in nanopores 485 and biological channels: A molecular simulation per-486 spective[END_REF][START_REF] Gravelle | Optimizing water permeability through the hourglass shape of aquaporins[END_REF][START_REF] Jungwirth | Biological Water or Rather Water in Biology?[END_REF] Despite the microscopic origin and nanomet-32 ric range of these physicochemical interactions, the impact 33 of interfaces on the water structure and dynamics -with a 34 pronounced dependence on confinement size, geometry and 35 chemistry -manifests itself at the macroscopic scale as illus-36 trated by considering water flow properties. [START_REF] Knight | Water properties under nano-scale confinement[END_REF][START_REF] Sendner | Interfacial water at hydrophobic and hydrophilic surfaces: Slip, viscosity, and diffusion[END_REF][START_REF] Kolesnikov | Anomalously Soft Dynamics of Water in a Nanotube: A Revelation of Nanoscale Confinement[END_REF][START_REF] Chen | Role of Hydrogen Bonding in Hysteresis Observed in Sorption-Induced Swelling of Soft Nanoporous Polymers[END_REF] Hydrophilic, i.e. wetting, surfaces are typically characterized by a no slip 38 boundary condition while hydrophobic, i.e. non-wetting, 39 surfaces are known to exhibit slippage. [START_REF] Bocquet | Nanofluidics, from bulk to interfaces[END_REF][START_REF] Siboulet | Hydrophobic Transition in Porous Amorphous Silica[END_REF] Such slip bound-40 ary conditions are broadly considered as a mean to overcome 41 the permeability-selectivity trade-off in membrane science such as for medical applications, water desalination, blue en-ergy harvesting [START_REF] Zhang | Surface slip on rotating graphene membrane enables the temporal selectivity that breaks the permeability-selectivity tradeoff[END_REF] (with the goal to compete with high performance biological selective nanochannels [START_REF] Gravelle | Optimizing water permeability through the hourglass shape of aquaporins[END_REF] ).

Despite their significant impact at the macroscopic level, the molecular specificities of water near solid surface remain only partly elucidated. At the microscopic level, a large body of experimental and numerical studies is available on the dynamics of water confined within wetting, i.e. hydrophilic, nanoporous media. In this case, the strong molecular interactions with the hydrophilic surfaces combined with severe in-pore confinement cause a decrease in water translational diffusion coefficient -i.e. a slowing down of water dynamics which gets more pronounced with tighter confinement and/or stronger surface interactions. [START_REF] Bellissent-Funel | Status of experiments probing the dynamics of water in confinement[END_REF][START_REF] Siboulet | Water self-diffusion at the surface of silica glasses: effect of hydrophilic to hydrophobic transition[END_REF][START_REF] Romero-Vargas Castrill Ón | Evolution from surfaceinfluenced to bulk-like dynamics in nanoscopically confined water[END_REF][START_REF] Osti | Characteristic features of water dynamics in restricted geometries investigated with quasi-elastic neutron scattering[END_REF][START_REF] Baum | Dynamical and Structural Properties of Water in Silica Nanoconfinement: Impact of Pore Size, Ion Nature, and Electrolyte Concentration[END_REF] In partially wetting materials, available data on moderately hydrophobic materials such as carbon nanotubes suggest a more complex scenario with strongly enhanced water transport as probed using permeability measurements. [START_REF] Hummer | Water conduction through the hydrophobic channel of a carbon nanotube[END_REF] While such faster transport can be rationalized by elucidating the very specific interactions between water and the host surface (weakened hydrogen-bond network, curvature-induced friction decrease, etc.), [START_REF] Falk | Molecular origin of fast water transport in carbon nanotube membranes: superlubricity versus curvature dependent friction[END_REF] it does not necessarily imply an increase in the diffusion coefficient and associated faster selfdynamics. As a matter of fact, even if the solid surface impact is less pronounced with hydrophobic surfaces, a decreased self-diffusion coefficient with respect to that of bulk water has been also reported in the literature by different groups. [START_REF] Siboulet | Water self-diffusion at the surface of silica glasses: effect of hydrophilic to hydrophobic transition[END_REF][START_REF] Jensen | The hydrophobic effect: Molecular dynamics simulations of water confined between extended hydrophobic and hydrophilic surfaces[END_REF] Nevertheless, in all these studies dedicated to hydrophilic and so-called hydrophobic materials, the affinity of water for the solid matrices is always strong enough to promote a spontaneous water condensation within the pores at a bulk vapor pressure lower than the saturation pressure.

Beyond hydrophilic and moderately hydrophobic materials, the microscopic dynamics in superhydrophobic nanoporous solids -such as all-silica zeolites and some Metal Organic Frameworks for which a pressure much larger than the water saturation pressure is needed to force water into the porosity -remains to be elucidated. Despite a well-established non-wetting thermodynamic behavior, [START_REF] Fraux | Forced intrusion of water and aqueous solutions in microporous materials: from fundamental thermodynamics to energy storage devices[END_REF] the molecular mechanisms ruling water transport in these materials involve complex phenomena which are often overlooked in such superhydrophobic environment (residence steps, in-pore relocation, etc.). In particular, the elemental dynamical mechanisms remain almost unexplored using microscopic experiments (e.g. neutron scattering, nuclear magnetic resonance). In this context, Zeolitic Imidazolate Frameworks (ZIF) are attractive candidates owing to their tunable chemistry, variety of topology and nanoscale pore size. [START_REF] Chen | Zeolitic imidazolate framework materials: recent progress in synthesis and applications[END_REF][START_REF] Ortiz | What makes zeolitic imidazolate frameworks hydrophobic or hydrophilic? The impact of geometry and functionalization on water adsorption[END_REF][START_REF] Grosu | Stability of zeolitic imidazolate frameworks: effect of forced water intrusion and framework flexibility dynamics[END_REF][START_REF] Michelin-Jamois | Giant osmotic pressure in the forced wetting of hydrophobic nanopores[END_REF] ZIF-8 is made of tetrahedral zinc ions linked by imidazolate organic groups (see insert in Figure 1). These two different moieties form a cubic structure of spherical cages, 11.6 Å in diameter, interconnected by 3.4 Å narrow gates which fill with water at a pressure of the order of 25 MPa at room temperature. The ZIF-8 structure, which is superhydrophobic, presents local heterogeneities as regularly encountered in biological confinement. The large microporous volume of ZIF-8 (∼ 0.4 cm 3 /g) [START_REF] Troyano | Colloidal metal organic framework particles : the pioneering case of ZIF-8[END_REF] is experimentally of interest to reach -after intrusion -a significant content of confined water with respect to bulk water. While water within ZIF-8 has been largely explored from a structural viewpoint, [START_REF] Ortiz | What makes zeolitic imidazolate frameworks hydrophobic or hydrophilic? The impact of geometry and functionalization on water adsorption[END_REF][START_REF] Calero | Underlying adsorption mechanisms of water in hydrophobic and hydrophilic zeolite imidazolate frameworks: ZIF-71 and ZIF-90[END_REF][START_REF] Nalaparaju | Molecular understanding for the adsorption of water and alcohols in hydrophilic and hydrophobic zeolitic metal-organic frameworks[END_REF][START_REF] Fraux | Dynamics, and Thermodynamics of Intruded Electrolytes in ZIF-8[END_REF][29] its dynamics has received only very little attention. On the one hand, Han et al. [START_REF] Han | Assessing the Im-597 pact of Point Defects on Molecular Diffusion in ZIF-8 598 Using Molecular Simulations[END_REF] considered the impact of defects on the dynamics of gas water molecules in ZIF-8, on the other hand Fraux et al. [START_REF] Fraux | Dynamics, and Thermodynamics of Intruded Electrolytes in ZIF-8[END_REF] reported a ten-fold decrease in liquid water self diffusivity. Despite these pioneering works, we note that the details of the molecular mechanisms leading to such a dynamics slowing down remain to be identified.

Here, we report on a coupled experimental and molecular simulation investigation of the microscopic dynamics of intruded water in superhydrophobic confinement. To the best of our knowledge, we performed the first direct measurement of the diffusion coefficient of liquid water in such superhydrophobic environment by means of quasi-elastic neutron scattering (QENS). Our experiments are combined with a dual Grand Canonical Monte Carlo/molecular dynamics simulation strategy which sheds light on the heterogeneous dynamics of water confined in such hydrophobic environment. While we confirm the average slowing down of the dynamics for confined water, careful analysis of our experimental data in the light of molecular simulations unveils a complex microscopic dynamics related to two different dynamical mechanisms: reduced dynamics for water located close to strongly adsorbing sites (ligands) and bulk-like water in the cavity center.

Material and Methods

The molecular simulation approach first involves Grand for neutrons, the signal arising from hydrogen motions of water dominates the scattered intensity, therefore providing 152 a measurement of the self-dynamics of water in the sam-153 ple. Descriptions of the experimental environment together 154 with calibration and adequate background measurements 155 are given in the Supporting Information (S2). In brief, an 156 optimized aluminium high pressure cell was designed 32 to 157 prepare a sub-millimetric uniform and dense layer of ZIF-8 158 powder (BASF, Basolite R Z1200), with similar volume frac-159 tions in intra-grain pores and inter-grain bulk voids. Mea-160 surements were performed under hydrostatic pressures of 161 4, 20, 28 and 38 MPa which are below and above the water 162 intrusion pressure of 25 MPa (see Figure 2a) at ambient tem-163 perature. [START_REF] Michelin-Jamois | Giant osmotic pressure in the forced wetting of hydrophobic nanopores[END_REF] The two lower pressures correspond to a matrix 164 empty of liquid whereas the two higher pressures corre-165 spond to a matrix filled with liquid of interest for the study 166 of water dynamics in this superhydrophobic framework. 167 The reproducibility of the measurements was assessed by 168 performing several short measurements during successive 169 intrusion cycle from low to high pressure.

170

Results and Discussion

171

We first investigate the water dynamics, using molecular 172 simulation, from the mean square displacements (MSD) for 173 bulk and confined water with respect to time (Figure 1). In 174 the long time limit t > 1 ns, the short-wave vector displace-175 ments (large distances) are dominated by translational mo-176 tion where both the O and H atoms diffuse over the same 177 distance on a given time scale. While the Fickian regime is 178 reached on a short time scale for bulk water (t ∼ 100 ps), we 179 found that much longer simulations are required to investi-180 gate this normal diffusion regime for confined water (typi-181 cally, t ∼ 10 ns). The diffusivity of water confined in ZIF-182 8 at 300 K is found to be smaller than its bulk counterpart; 183 typically, D T ∼ 3.4 × 10 -9 m 2 /s for bulk water while D T 184 ∼ 2.6 × 10 -10 m 2 /s for confined water. These two values 185 obtained with the TIP4P water model being 30% larger than 186 the ones obtained with the SPC/E model (extended simple 187 point charge model) [START_REF] Fraux | Dynamics, and Thermodynamics of Intruded Electrolytes in ZIF-8[END_REF] . However, importantly, the ten-fold 188 decrease in the water self diffusivity between bulk and con-189 fined values is independent of the force field or its scaling (as detailed in the Supplementary Information).

In the short time limit (t < 10-100 ps), for both bulk and confined water, the O and H mean square displacements depart from each other due to rotational effects. Interestingly, while the mean square displacements for bulk water is timescale independent in the sense that it exhibits a single regime once rotational effects are no longer noticeable, the data for confined water indicate a heterogeneous dynamics with different regimes. In more detail, between the time range dominated by rotation/translation coupling (t < 100 ps) and the normal Fickian regime (t ∼ 1-10 ns), a subdiffusive regime where ∆r 2 ∼ t α with α < 1 is observed.

At this stage it is useful to comfort these initial numerical results with experimental values of the diffusion coefficient D T . The confined water dynamics is extracted from the dynamical structure factor S(q, ω). [START_REF] Bée | Quasielastic neutron scattering[END_REF] The experimental S(ω) spectra (summed over the whole q-range) at the dif- nal is therefore a direct measurement of the water dynamics 226 in such superhydrophobic confinement. The limited ZIF-8 227 compressibility, [START_REF] Ortiz | Anisotropic elastic properties of flexible metal-organic 615 frameworks: How soft are soft porous crystals?[END_REF] which leads to minor deformations of the 228 matrix on the considered pressure range -essentially related 229 to a slight dilation of the order of 0.5% occurring within the 230 narrow pressure window associated to the filling process [START_REF] Tortora | Giant Negative Compressibility by Liquid 620 Intrusion into Superhydrophobic Flexible Nanoporous 621 Frameworks[END_REF] -231 is expected to not alter the liquid water dynamics within the 232 pores after intrusion. The superimposition of signals mea-233 sured at 28 MPa and 38 MPa confirms that once the matrix 234 is filled, the pressure increase does not alter the water dy-235 namics, neither directly through a potential change of water 236 structure nor indirectly because of the flexibility of the ma-237 trix. The dynamical structure factor, S(q, ω), in the inves-238 tigated energy range, can be approximated by the product 239 of the Debye -Waller factor, accounting for the vibrational 240 motions, with the decoupled diffusive motions, such as rota-241 tional, S R (q, ω), and translational, S T (q, ω), ones:

33 242 S(q, ω) = [S T (q, ω) + S R (q, ω)] exp - q 2 u 2 3 ( 1
)
with u 2 the mean-square vibrational amplitude of the hy-243 drogen atoms. In a first approach, both translational and 244 rotational motions were described using a Lorentzian func-245 tion of unconstrained half-width at half-maximum (hwhm), 246 Γ. Γ T and Γ R refer to the translational and rotational mo-247 tions, respectively. In the Fickian regime, the variation of Γ T 248 is linear with q 2 as Γ T = D T Fick q 2 where D T Fick is the transla-249 tional self-diffusion coefficient. For the localized rotational 250 motion, the width is q-independent: . The solid lines are attributed to the results obtained with the analysis 2 (Eq. 2) using the jump diffusion model. We note the presence of two solid blue lines associated to the heterogeneous behavior of the intruded highpressure water obtained with the analysis 2.

Γ R = 2D R with D R
ter dynamics. Moreover, we assume that inter-grain water is extracted parameters are given in Table 1.

265

The translational dynamics for bulk water and inter-grain 266 water obey to Fick's law for q 2 < 3 Å-

281

In contrast to inter-grain water, the translational dynamics 282 for intra-grain water follows Fick's law up to q 2 ∼ 2 Å-2 only.

283

In the low q Fickian regime, the confinement effect leads to a 284 decrease of both translational and rotational diffusion coef- ulation, to a linear Fickian regime that extends much beyond q 2 ∼ 2 Å-2 . As detailed in the Supporting Information (S3), the deviation from Fick's law beyond this q value, marked by the onset of a plateau for Γ T , provides a signature of the non homogeneous behavior of intra-grain water, due to the influence of the discrete atomic structure surrounding water (also highlighted in Figure 1 with the subdiffusive regime at short times). To further analyze this departure from an effective average behavior, a more rigorous model is needed to cover and accurately describe both the Fickian and non-Fickian dynamical regimes.

Table 1: Numerical values of the extracted fitting parameters for the analysis 1. D T Fick was calculated with the Fick's law on the appropriate q-range. We performed therefore a deeper analysis based on a formal description of water translation through the use of the jump diffusion model 33 coupled to isotropic rotational diffusion in a sphere of radius R (equal to the OH distance in a water molecule ∼ 1 Å). Each population will then be described by the new fitting model:

S(q, ω) = exp - q 2 u 2 3 S(q, ω) T⊗R (2)
with

S(q, ω) T⊗R = 1 π Γ T ω 2 + Γ 2 T ⊗ A 0 (q).δ(ω) + ∞ ∑ i=1 A i (q). 1 π Γ i ω 2 + Γ 2 i ( 3
)
where Γ T is defined within the jump-diffusion model:

Γ T = D T Jump q 2 1 + τD T Jump q 2 (4)
with τ the residence time, A i (q) = (2i + 1)j 2 i (qR) where j i is the spherical Bessel function of order i and

Γ i = i(i + 1)D R .
As expected, a single water population only is not sufficient to accurately fit our experimental data on the entire (q, w) range. On the other hand, a good fit is obtained using two water populations (ratio of 70/30 for fast and slow population, respectively). Each of these two populations are described with the model given in Eqs. 2 and 3 (fitting curves are presented in the Supporting Information (S2)). The resulting fitting parameters are shown in Table 2. The fast population is similar to bulk water while the slow population has a translational diffusion coefficient one order of magnitude lower as inferred from the MSD obtained using molecular simulation (Figure 1). As detailed in the Supporting Information (S3), at low q, the average dynamics of the two populations is similar to the one presented in Analysis 1 with an effective translational diffusion coefficient ∼ 30% than that of bulk water. This moderate deviation at low q is related to the small proportion of the slow water population with respect to the fast one. At high q, however, the slow water popula- As shown in Figure 4a, the self (incoherent) intermediate 348 scattering function, f (q, t) (inverse time Fourier transform of S(q, ω)) extracted from molecular dynamics simulations for bulk water conforms to the classical behavior expected for Fickian diffusion. The superimposition of the curves plotted with respect to q 2 t reveals that the rotational contribution is restricted to a very short time scale. This is further confirmed by the large value of D R ∼ 1 ps -1 obtained from the fit of f -inset of Figure 4a -using the inverse time Fourier transform of the formal model Eq. 2 (see Supporting Information S4). The dominant translational contribution of the form exp(-D T Jump q 2 t) is characterized by a diffusion coefficient D T Jump which agrees with the one inferred from MSD calculations. In contrast, as previously observed from our experimental data and molecular simulation MSD, f (q, t) for confined water display a non-Fickian behavior which can be approximated again with the contribution of a fast and a slow population (ratio 60/40) -each of them modeled with the inverse time Fourier transform of Eq. 2 (see inset of Figure 4a). Nevertheless, D T Jump for the fast population remains four times smaller than that for bulk water; this is a Figure 4: (a) Self (incoherent) intermediate scattering function f (q, t) as a function of q 2 t for different scattering vector q = | q | for bulk water (open circles) and for water confined in ZIF-8 (full dots). The color code is as follows: q = 0.7 Å-1 (red), q = 1.0 Å-1 (green), q = 1.5 Å-1 (blue), and q = 2.0 Å-1 (grey). Inset: f (q, t) according to t for bulk (open circle) and confined water (full dots). Dashed lines, bulk water: fit with a single water population using the inverse Fourier transform of Eq. 2 with D T Jump = 4.2 × 10 -9 m 2 /s, τ = 0.9 ps and D R = 0. At the microscopic level, the strong variations of water density as obtained from GCMC simulation (see Figure 4b and Figure 4c) reveal and confirm a complex heterogeneous water structure within the ZIF-8 porosity. These data, which are consistent with previous molecular simulation in the field, [START_REF] Ortiz | What makes zeolitic imidazolate frameworks hydrophobic or hydrophilic? The impact of geometry and functionalization on water adsorption[END_REF] show that water is adsorbed on specific sites, close

to the ligands at a distance about 0.4 -0. the propagators for the different water populations confined in ZIF-8 (note the use of a log scale for improved visualization purpose). In more detail, we measured the propagator for molecules initially located at time t = 0 in the strongly adsorbing sites close to the ZIF-8 ligands (see density contour plot shown in Figure 4b) and that for molecules initially located in the cavity center not directly in contact with ZIF-8 ligand/surface. Such data confirm that the fast diffusion mode (i.e. subpeak in the propagator located at distance r ∼ 3.5 Å) corresponds to molecules located within the cavity center while the slow diffusion model corresponds to molecules initially located in the strongly adsorbing sites. Interestingly, for both contributions, some molecules undergo displacements corresponding to the other mode. This is due to the fact that, during the characteristic time t = 10 ps considered here, some molecules leave the cavity center to migrate towards the strongly adsorbing sites and reciprocally.

Additionally to the analysis of the translational motions, we extracted the rotational relaxation time distribution for 434 water in the different regions of the ZIF-8 cavities to clarify 435 more specifically the impact of the surface heterogeneities 39 436 (see Figure 5b). The timescale of rotational motions enables 437 indeed a more local analysis since molecules, over such a 438 short timescale, do not significantly exchange between dif-439 ferent sites. Owing to their intrinsically slow dynamics, the 440 rotational motion of the molecules located in the strongly ad-441 sorbing sites is slower than for molecules in the pore center 442 again with a typical rotational relaxation time a few times 443 larger than for water in the pore center. Such a slow rota-444 tional motion only applies to molecules in the strong adsorb-445 ing sites. For molecules at the ZIF-8 cavity surface, away 446 from adsorbing sites, we indeed confirm that the relaxation By means of a coupled experimental and molecular simula-

  Canonical Monte Carlo simulations to determine the number of adsorbed water molecules as a function of pressure. Then starting from configurations obtained at a pressure where the ZIF-8 porosity is filled with water, we performed molecular dynamics simulations to compare bulk and confined water dynamics. A supercell consisting of a 2 × 2 × 2 ZIF-8 unit cell (the cubic lattice parameter is approximately 1.7 nm) was built to conduct these calculations. Periodic boundary conditions were used along the 3 space directions to avoid finite size effects. Following previous works on water in ZIF, the TIP4P ( Transferable Intermolecular Potential-4 Point) allatom model is used to describe the water molecule while both a nominal and rescaled force field proposed by Amrouche et al. 31 is used to model the water/ZIF-8 molecular interactions. All details on the molecular simulation strategy can be found in the Supporting Information (S5). QENS experiments were performed at Institut Laue-Langevin (ILL, Grenoble) on the time-of-flight spectrometer IN6-Sharp with an incident wavelength of 5.1 Å -leading to an elastic energy resolution of 70 µeV (dynamics faster than ∼ 10 ps) over the momentum transfer q-range of 0.4-2 Å-1 . Because of the large hydrogen incoherent cross section

Figure 1 :

 1 Figure 1: Mean square displacements ∆r 2 as obtained using molecular dynamics for the O and H atoms of water in bulk water (dashed lines) and confined in ZIF-8 (solid lines). The black dashed line is a guide to the eye to indicate the Fickian regime where ∆r 2 ∼ t. The insert shows the ZIF-8 structure.

  ferent investigated hydrostatic pressures are shown in Figure 2b. The contribution from the ZIF-dry matrix, which is purely elastic, confirms the absence of QENS signal that could arise from imidazole ring rotation on the investigated time scale. Moreover, complementary measurements have been performed on the spectrometer IN13 (ILL), covering a shorter (∼ 100 ps) timescale than IN6-SHARP. The results show that the rotational motion of the 2-methyl-imidazole rings is not influenced by the intrusion process (as it occurs on a timescale that is included in the elastic line of the present experiment). Details are presented in the Supporting Information (S1). We also observe that spectra below P int =25 MPa (4 and 20 MPa) are very similar as well as those above P int (28 and 38 MPa). Below 25 MPa, only inter-grain water is present outside the ZIF-8 grains as the nanopores are empty. Thus, the observed intensity increase between 20 and 28 MPa (Figure 2c) provides evidence for the presence of water within the nanopores above P int . The obtained high-pressure sig-

Figure 2 :

 2 Figure 2: (a) Experimental water intrusion/extrusion cycle at room temperature in a ZIF-8 material. V(P) indicates the intruded and extruded volume V of water per g of ZIF material as a function of the exerted pressure P. (b) S(ω) summed over q [from 0.4 to 2 Å-1 ] of the quasi-elastic spectra S(q, ω) obtained for dry ZIF-8 (black), low-pressure water in contact with ZIF-8 (red) and high-pressure water intruded in ZIF-8 (blue). Inset: schematics of these three situations. (c) Zoom on the elastic peak showing the overall increase of the intensity upon water addition.

Figure 3 :

 3 Figure3: Half-width at half-maximum Γ T of the translational contribution versus q 2 for bulk water (black), low-pressure water in contact with ZIF-8 (red) and high-pressure water intruded in ZIF-8 (blue). Error bars are smaller than the symbol size. The symbols refer to the results obtained with the analysis 1 (Eq. 1). The dashed segments indicate the Fickian regime where Γ T ∼ q 2 . The solid lines are attributed to the results obtained with the analysis 2 (Eq. 2) using the jump diffusion model. We note the presence of two solid blue lines associated to the heterogeneous behavior of the intruded highpressure water obtained with the analysis 2.

258 36 261

 36 identical at low and high pressures, as justified by the very 259 low-pressure dependence of the bulk translational diffusion 260 coefficient in the investigated pressure range (< 40 MPa). More information on this analysis and curve fitting quality 262 are given in the Supporting Information (S2). The obtained 263 q 2 dependence of Γ T is shown in Figure 3 (symbols) while the 264

285

  ficients with respect to their bulk counterpart. The slowing 286 down measured experimentally is only of the order of 30% 287 with respect to the bulk values. Such a limited departure 288 from the bulk dynamics should lead, for a single water pop-289

Analysis 1 :

 1 eq. 1 D T Fick ×10 -5 [cm 2 /s] D R [ps -1 ] Bulk water 2.01 ± 0.03 2.06 ± 0.01 ZIF-8: P<25 MPa 1.92 ± 0.04 2.09 ± 0.02 ZIF-8: P>25 MPa 1.22 ± 0.02 1.75 ± 0.02

333S4)

  for the two populations, the translational term expresses 334 as the sum of two time-decaying exponentials. At long time 335 this translational term should be controlled by the slower ex-336 ponential decay only associated to the slow population. This 337 means that at long time a Fickian regime with a diffusion co-338 efficient of the order of 1.7 × 10 -10 m 2 /s is expected. This 339 analysis of experimental data confirms, at long time, the ten-340 fold decrease observed from numerical MSD and reveals, at 341 short time, that the heterogeneous dynamics of confined wa-342 ter in ZIF-8 can be interpreted as the interplay between a slow 343 and a main bulk-like water population. 344 Attempting the same two-population decomposition us-345 ing numerical data, a similar trend -even if to a less ex-346 tent -is also observed when focusing on short time scales.

  347

  Figure4: (a) Self (incoherent) intermediate scattering function f (q, t) as a function of q 2 t for different scattering vector q = | q | for bulk water (open circles) and for water confined in ZIF-8 (full dots). The color code is as follows: q = 0.7 Å-1 (red), q = 1.0 Å-1 (green), q = 1.5 Å-1 (blue), and q = 2.0 Å-1 (grey). Inset: f (q, t) according to t for bulk (open circle) and confined water (full dots). Dashed lines, bulk water: fit with a single water population using the inverse Fourier transform of Eq. 2 with D T Jump = 4.2 × 10 -9 m 2 /s, τ = 0.9 ps and D R = 0.8 ps -1 . Full lines, confined water: fit with a fast and a slow population (ratio 60/40) with D T Jump = 1.1 × 10 -9 /3.4 × 10 -10 m 2 /s, τ = 0.25/20 ps and D R = 2/0.1 ps -1 (b) Contour plot showing the density of water in a ZIF-8 material as obtained using Grand Canonical Monte Carlo simulations. The density increases with color as follows: blue, green, yellow, orange and red. (c) Simulated density profiles ρ(r) for the O and H atoms of water intruded in ZIF-8. r is the distance from the center of the ZIF cages shown in the contour plot in (b). The dashed line indicates the bulk density of water at room temperature.)

  5 nm from the pore center [see marked density spots in Figure 4b], and in less adsorbing sites located closer to the pore center [see less marked density spots at a distance ∼ 0.2 nm]. Although these adsorbing sites are hydrophobic enough to not trap water at ambient pressure, their relative affinity for the intruded water contribute to the marked slowing down its dynamics observed at long time.For times much longer than the residence time, exchanges between the different adsorbing sites contribute to the complex heterogeneous dynamics and hence to the effective behavior accessible experimentally. Such a contribution from this molecular intermittence can be investigated by means of the diffusion propagator G(r, t) (inverse space and time Fourier transform of S(q, ω)) which corresponds to the probability that a molecule gets located at a position r at a time t while being at r = 0 at a time t = 0. At short times (t < 10 ps), G(r, t) shown in Figure5afor both O and H atoms of confined water is monomodal and reflects local displacements (r ∼ 0.1 nm) which are similar to those found for bulk water on the same time scale (although of a slightly smaller amplitude in confinement). For larger times (t > 10 ps), the propagator becomes bimodal for O and H atoms, therefore indicating two characteristic translational modes coexisting on the same time scale. We show in the insert of Figure5a

Figure 5 :

 5 Figure 5: (a) Diffusion propagator G(r, t) as obtained using molecular dynamics at room temperature for the O (red lines) and H (blue lines) atoms of water confined in ZIF-8. The dashed line is a guide to the eye to indicate that the data are estimated at increasing times t = 1, 5, 10 and 20 ps. As discussed in the main text, the insert shows the different contributions to the G(r, t) data for the O atoms at t = 20 ps: all molecules in the ZIF-8 cage (red), molecules in the cage center (blue) and molecules in strongly adsorbing sites (green). (b) Rotational autocorrelation time θ(t)θ(0) showing the change in a water molecule time at a time t = 0 and a time t later. The dashed and solid lines are for bulk water and water confined at room temperature in ZIF-8. The insert shows the same data for confined water but for different subpopulations: molecules in the cage center (blue), molecules at the cage surface (purple), and molecules in strongly adsorbing sites (green).

451

  tion study, we provide microscopic insights into the dynam-452 ics of liquid water intruded at high pressure within the ZIF-8 453 superhydrophobic nanoporous material. For the first time, 454 the direct experimental measurement of water in superhy-455 drophobic cavities enables a comparison with a molecular 456 simulation investigation. Qualitatively, a slowing down of 457 water is evidenced in both approaches. From an analysis 458 of experimental data based on the contribution of two water 459 populations, we show that at short distances, the main pop-460 ulation presents a bulk-like behavior while the slow popula-461 tion has a translational diffusion coefficient ten times smaller 462 than the bulk coefficient. This arbitrary subdivision into two 463 populations is sustained by the microscopic details obtained 464 from our molecular simulation approach. A fast population 465 is indeed identified in the cavity center while a slower pop-466 ulation is found close to the different adsorbing sites at the 467 cage surface (especially at the vicinity of imidazole ligands). 468 Even if water has a long residence time in these adsorbing 469 sites, exchanges with the pore center are also observed. Such 470 intermittence is found to contribute to the overall heteroge-471 neous water dynamics and could explain the discrepancy 472 between experimental and numerical values of the transla-473 tional diffusion coefficient of the fast population. We believe 474 that our findings represent an important step forward in the 475 understanding of complex dynamical behavior in superhy-476 drophobic systems which are relevant to many industrial ap-477 plications.

Table 2 : Numerical values of the extracted fitting parameters of the confined water in ZIF-8: P>25 MPa for the Analysis 2

 2

		Analysis 2: eq. 2	D T Jump ×10 -5 [cm 2 /s]	τ [ps]	D R [ps -1 ]
		ZIF-8: P<25 MPa	2.09 ± 0.05	0.55 ± 0.01 1.83 ± 0.01
		ZIF-8: P>25 MPa: Population 1	1.87 ± 0.02	0.25 ± 0.05 2.18 ± 0.04
		ZIF-8: P>25 MPa: Population 2	0.167 ± 0.007	2.55 ± 1.6	0.17 ± 0.01
		tion, despite its limited proportion, becomes a leading term
	330	which contributes to the plateau observed for Γ T in Figure 3
	331	in the framework of analysis 1. Furthermore, moving now
	332	to the temporal domain using the inverse time Fourier trans-
		form of the formal model Eq. 2 (see Supporting Information

time distribution is close to the one found for molecules in 448 the pore center.

Conclusion Supporting Information Available

479

The following files are available free of charge.

480

IN13 data; IN6 Experimental details and fitting curves;

Complementary information on the analysis 1 and 2 in the 482 experimental part; Computational details Acknowledgement The authors thank the Institut Laue-Langevin (Grenoble, France) for the allocation of neutron beam time. The neutron data set is available at DOI: http://10.5291/ILL-DATA.EASY-537. This work was supported by the French Research Agency (ANR LyStEn 15-CE06-0006, ANR TWIST 17-CE08-0003). Some of the computations in this paper were performed using the Froggy platform of the GRICAD infrastructure (https://gricad.univgrenoble-alpes.fr), which is supported by the Rhne-Alpes region (GRANT CPER07-13 CIRA) and the Equip@Meso project (reference ANR-10-EQPX-29-01) of the programme Investissements dAvenir supervised by the French Research Agency. We thank Claude Payre and Jér ôme Giraud for high pressure equipment and the design of high pressure cell. We are grateful to David Farruseng who initiated fruitful collaboration for this work.

TOC Graphic


Denys Dutykh

Frédéric Dias

Energy of tsunami waves generated by bottom motion

Keywords: water waves, shallow-water equations, tsunami energy

In the vast literature on tsunami research, few articles have been devoted to energy issues. A theoretical investigation on the energy of waves generated by bottom motion is performed here. We start with the full incompressible Euler equations in the presence of a free surface and derive both dispersive and non-dispersive shallowwater equations with an energy equation. It is shown that dispersive effects only appear at higher order in the energy budget. Then we solve the Cauchy-Poisson problem of tsunami generation for the linearized water wave equations. Exchanges between potential and kinetic energies are clearly revealed.

Introduction

Oceanic waves can be devastating as shown by recent events. Whilst some areas are more vulnerable than others, the recent history shows that catastrophic waves can hit even where they are not expected. The tsunami waves generated by the huge undersea earthquake in Indonesia on 26 December 2004 caused devastation across most of the coasts of the Bay of Bengal. The tsunami waves generated by the massive submarine landslide in Papua-New Guinea on 17 July 1998 as well as the 17 July 2006 Java tsunami and the 2 April 2007 Solomon Islands tsunami also caused devastation, but on a smaller scale. Unfortunately, such cataclysmic tsunamis are likely to be generated again by earthquakes, massive landslides or volcano eruptions [START_REF] Synolakis | Tsunami science before and beyond Boxing Day 2004[END_REF].

Information on tsunami energy can be obtained by applying the normal mode representation of tsunami waves, as introduced by [START_REF] Ward | Relationship of tsunami generation and an earthquake source[END_REF]. For example, [START_REF] Okal | Normal mode energetics for far-field tsunamis generated by dislocations and landslides[END_REF] considers the total energy released into tsunami waves. He obtains expressions for the energy of tsunamis (see his expressions (31) for a tsunami generated by an earthquake and (36) for a tsunami generated by a landslide). In the case of a landslide, he computes the ratio between tsunami energy and total change in energy due to the slide. In the present paper, we use the incompressible fluid dynamics equations. Tsunamis have traditionally been considered as non-dispersive long waves. However various types of data (bottom pressure records [START_REF] Ritsema | Inversion of deep-ocean tsunami records for 1987 to 1988 Gulf of Alaska earthquake parameters[END_REF]; satellite data [START_REF] Kulikov | Dispersion of the Sumatra tsunami waves in the Indian Ocean detected by satellite altimetry[END_REF]; hydrophone records [START_REF] Okal | Quantification of hydrophone records of the 2004 Sumatra tsunami[END_REF]) indicate that tsunamis are made up of a very long dispersive wave train, especially when they have enough time to propagate. These waves travel across the ocean surface in all directions away from the generation region. Recent numerical computations using dispersive wave models such as the Boussinesq equations show as much as 20% reduction of tsunami amplitude in certain locations due to dispersion [START_REF] Dao | Tsunami propagation modelling -a sensivity study[END_REF]. But one has to be careful with the interpretation of satellite data: as indicated by [START_REF] Kânoglu | Initial value problem solution of nonlinear shallow water-wave equations[END_REF], the mid-ocean steepness of the 2004 Sumatra tsunami measured from satellite altimeter data was less than 10 -5 . Nonlinear dispersive theory is necessary only when examining steep gravity waves, which is not the case in deep water.

The wavelength of tsunamis and, consequently, their period depend essentially on the source mechanism. If the tsunami is generated by a large and shallow earthquake, its initial wavelength and period will be greater. On the other hand, if the tsunami is caused by a landslide (which happens less commonly but can be devastating as well), both its initial wavelength and period will be shorter as indicated for example by [START_REF] Kulikov | The landslide tsunami of November 3, 1994, Skagway Harbor, Alaska[END_REF]. From these empirical considerations one can conclude that dispersive effects are a priori more important for submarine landslide and slump scenarios than for tsunamigenic earthquakes.

Once a tsunami has been generated, its energy is distributed throughout the water column. Clearly, the more water is displaced, the more energetic is the tsunami (compare for example the December 2004 and March 2005 Sumatra tsunamis). Due to the large scale of this natural phenomenon and limited power of computers, tsunami wave modellers have to adopt simplified models which reduce a fully three-dimensional (3D) problem to a two-dimensional (2D) one. This approach is natural, since in the case of very long waves the water column moves as a whole. Consequently the flow is almost 2D. Among these models one can mention the nonlinear shallow water equations (SWE), Boussinesq type models, the Green-Naghdi and Serre equations. There is a wide variety of models, depending on whether or not the effects of run-up/run-down, bottom friction, turbulence, Coriolis effects, tidal effects, etc, are included.

Today scientists can easily predict when a tsunami will arrive at various places by knowing source characteristics and bathymetry data along the paths to those places. Unfortunately one does not know as much about the energy propagation of such waves. Obviously tsunami amplitude is enhanced over the major oceanic ridges. [START_REF] Titov | The global reach of the 26 December 2004 Sumatra tsunami[END_REF] clearly describe the waveguide type effect from midocean ridges that has funnelled the 2004 megatsunami away from the tip of Africa. As emphasized by [START_REF] Kowalik | The tsunami of 26 December, 2004: Numerical modeling and energy considerations[END_REF], travel-time computation based on the first arrival time may lead to errors in the prediction of tsunami arrival time as higher energy waves propagate slower along ridges. At the beginning, the energy is essentially potential, although it depends on the generation mechanism. Then it redistributes itself into half kinetic and half potential energies. Finally, it converts its potential component into kinetic energy. How do these conversions take place? The purpose of this study is to shed some light on this topic and to see if the importance of dispersion in tsunamis can be studied by looking at the energy rather than at wave profiles.

Previous researchers have considered the topic of tsunami wave energy. [START_REF] Dotsenko | The effect of frequency dispersion on plane waves generated as a result of bed motions[END_REF], [START_REF] Kajiura | Tsunami source, energy and the directivity of wave radiation[END_REF], [START_REF] Velichko | Amplitude-energy characteristics of tsunami waves for various types of seismic sources generating them[END_REF] studied the energy exchange between the solid bottom and the overlying water associated with the bottom deformation. There were recent attempts to obtain equations for tsunami energy propagation. We can mention here the work of [START_REF] Tinti | Energy of water waves induced by submarine landslides[END_REF] Article submitted to Royal Society devoted to idealized theoretical cases and the work of [START_REF] Kowalik | The tsunami of 26 December, 2004: Numerical modeling and energy considerations[END_REF] using the energy flux point of view to study the changes in the 2004 tsunami signal as it travelled from Indonesia to the Pacific Ocean. We believe that these models can be improved, given the present state of the art in wave modelling.

A point of interest is that some of the equations used for wave modelling have an infinite number of conserved quantities. There has been some confusion in the literature on which quantities can be called energy. Indeed there is here an interesting question. In incompressible fluid mechanics, the internal energy equation is decoupled from the equation of continuity and from the fundamental law of dynamics. It is used only when one is interested in computing the temperature field once the velocity distribution is known. In addition to the internal energy equation, one can write a total energy (internal energy + kinetic energy) equation, or a total enthalpy equation. The confusing part is that for perfect fluids one usually defines the total energy differently: it is the sum of internal energy, kinetic energy, and potential energies associated to body forces such as gravitational forces and to the pressure field. If in addition the fluid is incompressible, then the internal energy remains constant. In the classical textbooks on water waves [START_REF] Johnson | A modern introduction to the mathematical theory of water waves[END_REF][START_REF] Stoker | Water waves, the mathematical theory with applications[END_REF], one usually introduces the energy E as the sum of kinetic and potential energies and then looks for a partial differential equation giving the time derivative dE/dt (incidentally the meaning of d/dt is not always clearly defined). In any case, when one uses a depth-integrated model such as the nonlinear SWE, one can compute the energy a posteriori (the potential energy is based on the free-surface elevation and the kinetic energy on the horizontal velocity). But one can also apply the nonlinear shallow water assumptions to the full energy equation to begin with. Then one obtains a nonlinear shallow water approximation of the energy equation. Are these two approaches equivalent? We show that the answer is no.

First we present the energy equation in three different forms: full water wave equations, dispersive SWE and non-dispersive SWE. Surprisingly, the energy equation is the same for dispersive and non-dispersive SWE at leading order. Then we present some numerical computations over a flat bottom. It allows us to concentrate on the generation process and the energy transfer through the moving seabed. We refer to Duytkh et al. [2008] for simulations of some real world events including energy pumping and its transformation over uneven bathymetry, using the SWE with energy. Finally we present the energy equation in a fourth form: the linearized dispersive water wave equations. We solve the Cauchy-Poisson problem of tsunami generation. Exchanges between potential and kinetic energies are clearly revealed.

Derivation of the energy equation

Consider the 3D fluid domain shown in Figure 1. It is bounded above by the free surface z * = η * (x * , y * , t * ) and below by the solid boundary with prescribed motion z * = -h * (x * , y * , t * ). A Cartesian coordinate system with the z * -axis pointing vertically upwards and the x * Oy * -plane coinciding with the still-water level is chosen.

The fluid is assumed to be inviscid. Its motion is governed by the 3D Euler equations, written here in their incompressible form (see for example [START_REF] Gisler | Tsunami simulations[END_REF] for the compressible counterpart, after replacing ρ by p in the pressure term of his Eq. 3):

∇ • u * = 0, (2.1) ∂u * ∂t * + ∇ • u * ⊗ u * + p * ρ * Id = g, (2.2) ∂e * ∂t * + ∇ • e * + p * ρ * u * = 0, (2.3)
where ρ * is the fluid density, u * = (u * , v * , w * ) the velocity vector, e * the sum of kinetic energy density e * K = 1 2 |u * | 2 and potential energy density e * P = gz * , p * the pressure and g the acceleration due to gravity. In the present study g = (0, 0, -g). For incompressible flows, the energy equation (2.3) is redundant. Indeed it can be obtained from Eq. (2.2). However we keep it since it is not equivalent to derive shallow water equations with or without the energy equation.

Equations (2.1)-(2.3) have to be completed by the kinematic and dynamic boundary conditions. Since surface tension effects are not important for long waves, the dynamic boundary condition on the free surface reads

p * = p * s , at z * = η * . (2.4)
Later we will replace the surface pressure p * s by 0 but we keep it arbitrary for now. The kinematic boundary conditions on the free surface and at the seabed are, respectively,

w * = ∂η * ∂t * + u * ∂η * ∂x * + v * ∂η * ∂y * , z * = η * , (2.5) w * = - ∂h * ∂t * -u * ∂h * ∂x * -v * ∂h * ∂y * , z * = -h * . (2.6)
Below we denote the horizontal gradient by ∇ ⊥ and the horizontal velocity by u * ⊥ . After a few manipulations and integration across the water column from bottom to Article submitted to Royal Society top, one can write the following global energy equation:

∂E * ∂t * + ∇ ⊥ • Φ * + P * = 0, (2.7)
where E * is the sum of kinetic and potential energies in the flow, per unit horizontal area, Φ * the horizontal energy flux vector, and P * the net energy input due to the pressure forces doing work on the upper and lower boundaries of the fluid. They are given by the following expressions:

E * = η * -h * 1 2 ρ * |u * | 2 + ρ * gz * dz * , (2.8) Φ * = η * -h * u * ⊥ 1 2 ρ * |u * | 2 + p * + ρ * gz * dz * ,
(2.9)

P * = p * s η * t * + p * b h * t * , (2.10)
where p * s is the pressure exerted on the free surface and p * b the bottom pressure. In the case of a stationary bottom boundary and of a free surface on which the pressure vanishes, then as expected the net energy input P * is identically zero. Energy can be brought to the system by a moving bottom or by a pressure disturbance on the free surface. From now on, we take p * s = 0.

(a) Dimensionless equations

The problem of tsunami propagation possesses two characteristic length scales: the average water depth h 0 for the vertical dimension and a typical wavelength l for the horizontal dimensions. It is classical to introduce the following dimensionless variables. The scaling for the independent variables is

x = x * l , y = y * l , z = z * h 0 , t = √ gh 0 l t * .
In order to introduce the dimensionless dependent variables we need one more parameter, the typical wave amplitude a:

u = h 0 a √ gh 0 u * , v = h 0 a √ gh 0 v * , w = h 0 l h 0 a √ gh 0 w * , η = η * a , h = h * h 0 , π = p * + ρ * gz * ρ * ga , p = p * ρ * gh 0 , e = e * gh 0 .
The hydrostatic pressure -ρ * gz * has been incorporated into π.

The following dimensionless parameters, which are assumed to be small, are introduced:

ε := a/h 0 , µ := h 0 /l.
The parameter ε represents the relative importance of nonlinear terms and µ measures dispersive effects. Note that

e = 1 2 ε 2 (u 2 + v 2 ) + 1 2 ε 2 µ 2 w 2 + z.
Article submitted to Royal Society

The Euler equations of motion (2.1)-( 2.3) become in dimensionless form

µ 2 (u x + v y ) + w z = 0, (2.11) µ 2 u t + εµ 2 (u 2 ) x + εµ 2 (uv) y + ε(uw) z + µ 2 π x = 0,
(2.12)

µ 2 v t + εµ 2 (uv) x + εµ 2 (v 2 ) y + ε(vw) z + µ 2 π y = 0, (2.13) µ 2 w t + εµ 2 (uw) x + εµ 2 (vw) y + ε(w 2 ) z + µ 2 π z = 0, (2.14)
µ 2 e t + εµ 2 ((e + p)u) x + εµ 2 ((e + p)v) y + ε ((e + p)w) z = 0.

(2.15)

In dimensionless form the boundary conditions (2.4)-(2.6) become

π = η, z = εη, (2.16) w = µ 2 η t + εµ 2 uη x + εµ 2 vη y , z = εη, (2.17) εw = -µ 2 h t -εµ 2 uh x -εµ 2 vh y , z = -h. (2.18)
In the case of a static bottom h = h(x, y), the order of magnitude of the vertical velocity w at the bottom is O(µ 2 ). With a moving bathymetry the behaviour is different:

w| z=-h = - µ 2 ε h t + O(µ 2 ).
(b) Integration over the depth

Next we reduce the above 3D problem (2.11)-(2.18) into a 2D one by integrating all the equations over the water column. We begin with the continuity equation (2.11) which we integrate with respect to z from -h to εη. Taking into account the boundary conditions (2.17) and (2.18), one obtains

η t + ∂ ∂x εη -h u dz + ∂ ∂y εη -h v dz = - 1 ε h t .
(2.19)

A source term appears in equation (2.19) due to the moving bathymetry. For tsunami generation, the function h * can be represented as follows:

h * (x * , y * , t * ) = h * 0 (x * , y * ) -ζ * (x * , y * , t * ),
where h * 0 (x * , y * ) is the static sea bed profile and ζ * (x * , y * , t * ) the bottom displacement due for example to coseismic displacements or landslides. In nondimensional form this representation takes the form:

h(x, y, t) = h 0 (x, y) -εζ(x, y, t),
since the bottom variation must be of the same order of magnitude as the typical wave amplitude a: ζ * (x * , y * , t * ) = aζ(x, y, t). Differentiating with respect to time yields 1 ε

∂h ∂t = - ∂ζ ∂t = O(1).
Below we will replace 1 ε h t by -ζ t .

Article submitted to Royal Society

Integrating the vertical momentum conservation equation (2.14) yields

π| z=-h = η + ∂ ∂t εη -h w dz + ε ∂ ∂x εη -h uw dz + ε ∂ ∂y εη -h vw dz.
(2.20)

A more general expression for the pressure can be obtained if we integrate equation (2.14) from z to εη and use the boundary conditions (2.16) and (2.17):

π = η + ∂ ∂t εη z w dz + ε ∂ ∂x εη z uw dz + ε ∂ ∂y εη z vw dz - ε µ 2 w 2 .
(2.21)

The vertical velocity w is obtained by integrating the continuity equation (2.11) from -h to z and applying the seabed kinematic condition (2.18):

w = µ 2 ζ t -µ 2   ∂ ∂x z -h u dz + ∂ ∂y z -h v dz   . (2.22)
Finally, the integration of the equation of energy conservation (2.15) yields

∂ ∂t εη -h e dz + ε ∂ ∂x εη -h (e + p)u dz + ε ∂ ∂y εη -h (e + p)v dz -ε p| z=-h ζ t = 0, (2.23)
which is nothing else than the dimensionless counterpart of Eq. (2.7).

All the equations derived above are exact and no assumption has been made about the orders of magnitude of ε and µ.

(c) The nonlinear shallow-water equations with energy equation

In Appendix A, we briefly summarize the derivation of various systems of shallow-water wave equations. The non-dispersive SWE are obtained by taking the limit µ → 0. Here we provide the dispersive and non-dispersive SWE in their conservative forms [START_REF] Eskilsson | Discontinuous Galerkin spectral/hp element modelling of dispersive shallow water systems[END_REF] with dimensions, based on the depth-averaged horizontal velocity ū * . The total water depth h * + η * is denoted by H * (x * , y * , t * ). The definition of E * is given by Eq. (2.8).

(i) SWE with dispersion and energy equation

∂H * ∂t * + ∇ • (H * ū * ) = 0, (2.24) ∂(H * ū * ) ∂t * + ∇ • H * ū * ⊗ ū * + 1 2 gH * 2 + h * 3 6 ∇ ∇ • H * ū * h * - h * 2 2 ∇ (∇ • (H * ū * )) t * = gH * ∇h * , (2.25) ∂E * ∂t * + ∇ • ū * E * + 1 2 ρ * gH * 2 = ρ * gH * ∂ζ * ∂t * . (2.

26)

Article submitted to Royal Society (ii) SWE without dispersion and with energy equation

∂H * ∂t * + ∇ • (H * ū * ) = 0, (2.27) ∂(H * ū * ) ∂t * + ∇ • H * ū * ⊗ ū * + 1 2 gH * 2 = gH * ∇h * ,
(2.28)

∂E * ∂t * + ∇ • ū * E * + 1 2 ρ * gH * 2 = ρ * gH * ∂ζ * ∂t * .
(2.29)

Between (i) and (ii), only the equation for the evolution of the horizontal velocity (2.28) differs. In particular, it is interesting to note (this is the first main result of the present paper) that the energy equations are the same in the dispersive and non-dispersive cases. Differences appear only at higher order (terms of order O(εµ 2 )). The physical meaning is that the vertical velocity as well as the nonuniform structure of the horizontal velocity appear only at the next order. An additional remark can be made about the hyperbolic structure of the system (2.27)-(2.29). We restrict our observations to one space dimension. Let us assume that a shock wave propagates from left to right at velocity s > 0. The states before and after the discontinuity are denoted by (H * l,r , ū * l,r , E * l,r ) correspondingly. From entropy considerations one can conclude that for an admissible state H * l > H * r . For a general system of conservation laws

v t + ∂ x f (v) = 0, the Rankine-Hugoniot relations have the simple form f (v r ) -f (v l ) = s(v r -v l ).
Some simple algebraic calculations yield the following relations:

s = H * l ū * l -H * r ū * r H * l -H * r , ū * r = ū * l ± 1 2 g(H * 2 l -H * 2 r ) 1 H * r - 1 H * l , E * r H * r = E * l H * l - 1 2 ρ * g 1 H * r - 1 H * l H * 2 r ū * r -H * 2 l ū * l ū * r -ū * l .
(2.30)

These formulas relate left and right states connected in the (H * , ū * , E * ) space by a shock wave. The first two relations are well known. What is new to our knowledge is the formula (2.30) which gives an insight into the energy states in a shock wave. In practice, they can be used for the theoretical analysis of bores and as a validation test for nonlinear SWE codes (with energy equation).

Simulations of energy

Next we illustrate the main features of energy evolution in tsunami generation. The importance of dispersive effects strongly depends on the extent of the source area (the smaller the source the stronger the dispersive effects) and the ocean depth in the source area [START_REF] Kervella | Comparison between three-dimensional linear and nonlinear tsunami generation models[END_REF]. We restrict our study to the nondispersive SWE (2.27)-(2.29). We solve these equations numerically with a finite volume method [Duytkh et al., 2008].

While the common practice in modeling tsunami generation consists in translating the initial sea bottom deformation to the water surface, thus neglecting all dynamical effects, we prefer to include some dynamics in the process in an effort to be closer to what happens in reality [START_REF] Dutykh | Linear theory of wave generation by a moving bottom[END_REF]. We construct the bottom motion by multiplying Okada's static solution ζ * OK (x * , y * ) † by a function of time [START_REF] Hammack | A note on tsunamis: their generation and propagation in an ocean of uniform depth[END_REF]:

h * (x * , y * , t * ) = h * 0 (x * , y * ) -(1 -e -α * t * )ζ * OK (x * , y * ).
The parameter α * is related to the characteristic time t * 0 under consideration. We chose

1 -e -α * t * 0 = 2 3 ⇔ α * = log 3 t * 0 .
The various parameters used in the computations are given in Table 1.

Other time laws are possible and we refer to [START_REF] Dutykh | Water waves generated by a moving bottom[END_REF] for more details. In the present numerical computations we chose h * 0 (x * , y * ) = h 0 = const. This choice is not only made for the sake of simplicity. Another reason is that Okada's solution is derived within the assumption of an elastic half-space which does not take into account the bathymetry. In order to be coherent with this solution we assume the bottom to be flat before deformation.

With our definition of potential energy, the total energy is nonzero both at time t * = 0 and at time t * → ∞. For practical applications it is important to isolate the energy available to the tsunami wave. One possibility is to define the wave energy as follows:

E * wave = E * + 1 2 ρ * gh * 2 . (3.1)
Clearly E * wave = 0 both at time t * = 0 and at time t * → ∞ when the wave has left the generation region. From Eq. (2.29), one finds that the energy equation satisfied by E * wave is

∂E * wave ∂t * + ∇ • ū * E * wave + 1 2 ρ * gH * 2 - 1 2 ρ * gh * 2 = ρ * g(H * -h * ) ∂ζ * ∂t * . (3.
2) † Okada's solution is a steady analytical solution for the seafloor displacement following an underwater earthquake, based on dislocation theory in an elastic half-space [START_REF] Okada | Surface deformation due to shear and tensile faults in a half-space[END_REF]. [START_REF] Ben-Menahem | Amplitude patterns of tsunami waves from submarine earthquakes[END_REF], tsunami energy radiates primarily at right angles to a rupturing fault (see also [START_REF] Kajiura | Tsunami source, energy and the directivity of wave radiation[END_REF]). The distribution of potential energy makes sense when one compares the potential energy plots with the free-surface plots. The total energy spreads in a more uniform manner across the area affected by the waves.

Article submitted to Royal Society

It is also of interest to see how the wave energy E * wave , integrated over the whole fluid domain, varies during the generation process. This evolution is shown in figure 10. Three curves are plotted in figure 10(a). The solid curve simply is the energy (2.8) integrated over the water surface. The middle curve will be explained in Section 4 (linearized theory). The bottom curve is obtained as follows. Imagine that one is solving the SWE without the energy equation. In order to look at the of numerical method we use. In figure 10(b), one clearly sees the exchange between potential and kinetic energies until equipartition is reached.

Energy in the framework of the dispersive linearized equations

In the case of tsunamis generated by earthquakes, nonlinear effects are not important during the process of generation and propagation. This is why it is valid to use the linearized water-wave equations. [START_REF] Dutykh | Linear theory of wave generation by a moving bottom[END_REF] and others showed that taking an instantaneous seabed deformation is not equivalent to instantaneously transferring the seabed deformation to the ocean surface, except in the framework of the linearized shallow water equations (very long waves). The difference comes from the vertical velocities and dispersion. In this case we must go back to the initial set of equations (2.1)-(2.3). Since the motion starts from the state of rest, it can be considered as irrotational (potential flow) and one can introduce the velocity potential u * = ∇φ * .

We perform the linearization of the equations (2.1) and (2.2), and of the boundary conditions (2.4)-(2.6). It is equivalent to taking the limit of the equations as ε → 0. For the sake of convenience, we switch back to the physical variables. The linearized problem in dimensional variables reads [START_REF] Dutykh | Water waves generated by a moving bottom[END_REF]

∆ ⊥ φ * + ∂ 2 φ * /∂z * 2 = 0, (x * , y * , z * ) ∈ R 2 × [-h * , 0], (4.1) ∂φ * ∂z * = ∂η * ∂t * , z * = 0, kinematic condition, (4.2) ∂φ * ∂t * + gη * = 0, z * = 0, dynamic condition. (4.3)
Article submitted to Royal Society

Within linear theory the forces that cause perturbations are so weak that the boundary condition at the bottom (2.6) is also simplified:

∂φ * ∂z * = ∂ζ * ∂t * , z * = -h * .
(4.4)

The bottom motion appears in the right-hand side of Eq. (4.4). The Laplace equation (4.1) together with the boundary conditions (4.2), (4.3) and (4.4) determine the boundary-value problem for the velocity potential φ * within the linear theory. In order to solve the equations for a prescribed bottom motion, one can use Fourier and Laplace transforms (this approach is followed here) or Green's functions. Three scenarios are considered for the bottom motion [START_REF] Dutykh | Linear theory of wave generation by a moving bottom[END_REF]: the passive generation in which the deformation of the sea bottom is simply translated to the free surface (one is then solving an initial value problem) and two dynamical processes

ζ * (x * , y * , t * ) = T (t * )ζ * OK (x * , y * )
, where ζ * OK (x * , y * ) is given by Okada's solution. The two choices for T are the instantaneous deformation with

T (t * ) = T i (t * ) = H(t * )
, where H(t * ) denotes the Heaviside step function, and the exponential law used above in Section 3:

T (t * ) = T e (t * ) = 0, t * < 0, 1 -e -α * t * , t * ≥ 0, with α * > 0. Let ζ * OK be the Fourier transform of ζ * OK , ω 2 = g|k * | tanh(|k * |h *
) the dispersion relation and x * = (x * , y * ). We provide the general integral solution for the free surface elevation in the three cases, and the velocity potential in the passive case (the expressions for the other two active cases are a bit cumbersome) [START_REF] Dutykh | Linear theory of wave generation by a moving bottom[END_REF],

η * o (x * , t * ) = 1 (2π) 2 R 2 ζ * OK e ik * •x * cos(ωt * ) dk * , (passive) (4.5) η * i (x * , t * ) = 1 (2π) 2 R 2 ζ * OK e ik
φ * o (x * , t * ) = 1 4π 2 R 2 -g ω ζ * OK e ik * •x * (cosh |k * |z * + tanh |k * |h * sinh |k * |z * ) dk * .
(4.8) Then one can easily compute the kinetic and potential energies:

E * K = 1 2 ρ * R 2 η * -h * |∇φ * | 2 dx * dz * , E * P = 1 2 ρ * g R 2 η * 2 dx * . (4.9)
Results are shown in figures 11 and 12. Even though there are differences during the first few seconds, the three mecanisms lead to the same almost exact equipartition between kinetic and potential energy once the dipolar waves start to propagate. The simplest estimate proposed for the energy of tsunamis generated by a dislocation source is that given by [START_REF] Okal | A theoretical comparison of tsunamis from dislocations and landslides[END_REF]. They compute the increase in potential energy of the ocean by displacing a volume of water S × δh * from the bottom to the surface of the ocean. This also represents the work of the pressure forces displacing the ocean bottom. Then they explain that the center of mass of the displaced water, initially at height δh * /2 above the ocean floor, is transferred to the ocean surface, so that the change in potential energy is not as much. The difference between the two is the energy available to the tsunami wave:

E * wave dx * = 1 2 ρ * gS(δh * ) 2 . (4.10)
Incidentally, this expression does not depend on the sign of δh * and is also valid for a sudden subsidence of a section of the ocean floor. It can be extended to a more realistic sea floor deformation, such as the one used in this paper, ζ * OK (x * , y * ): (4.11) This quantity corresponds to the square in figure 11(b).

E * wave dx * = 1 2 ρ * g (ζ * OK ) 2 dx * .

Concluding remarks

In this article we provided a formal derivation of the total energy equation in the framework of the nonlinear SWE, both for dispersive and non-dispersive waves. We also made an attempt to better understand the energy transfer from a moving bottom to the water above. The importance of this topic is clear, given the serious hazard that tsunamis represent for coastal regions.

Tsunami energy can be studied at several levels. A simple formula was given by [START_REF] Okal | A theoretical comparison of tsunamis from dislocations and landslides[END_REF]. In the present paper, we extended it to more spatially realistic sea floor deformations. But this formula does not involve any dynamics, which can play an important role in tsunami generation. A somewhat counterintuitive consequence of this simple estimate is that the energy carried by the full tsunami wave is practically independent of depth [START_REF] Okal | Source discriminants for near-field tsunamis[END_REF]. One would think that the energy involved in lifting a 4 km column of water above a rupture is different from lifting only 100 m of water (whether instantaneously or in a few seconds). But no difference can be found in wave height. The assumption of incompressibility may in fact no longer be valid, especially for a very deep ocean.

The emphasis of this paper has been on tsunami generation. At the runup or inundation stage, energy is also quite important. Potential energy is being transferred into kinetic energy and the study of these exchanges is left for future work. Our approach can also be used to analyze the structure of the wave field in caustics more accurately than with ray theory [START_REF] Berry | Focused tsunami waves[END_REF]. The present paper can be considered as a first step towards a better understanding of tsunami wave energy, in order to provide scales for tsunami magnitudes for example. More profound mathematical and physical analysis is needed.

Another extension is the study of other mechanisms. For example [START_REF] Okal | Normal mode energetics for far-field tsunamis generated by dislocations and landslides[END_REF] showed that the combination of a lesser absolute level of excitation and a more

Article submitted to Royal Society

We see that the horizontal velocities vary quadratically with respect to z over the depth at leading order O(µ 2 ).

Following [START_REF] Ursell | The long-wave paradox in the theory of gravity waves[END_REF], we introduce a number which measures the relative importance of nonlinear and dispersive effects in long waves:

S := ε/µ 2 .
In order to simplify the computations we now assume that S = O(1). Dispersive terms can be neglected if S ≫ 1.

An expression for the pressure can be obtained by substituting (A 5) and (A 3) into (2.21), integrating and retaining leading order terms:

π = η + µ 2 2 z 2 ∇ • u bt + µ 2 z ∇ • (hu bt ) -ζ tt + O(µ 4 ). (A 6)
The equation for the free-surface evolution is derived by substituting (A 5) into the depth-integrated continuity equation (2.19) and integrating:

η t + ∇ • (h + εη)u b -µ 2 ∇ • h 2 2 ∇ ∇ • (hu b ) -ζ t - h 3 3 ∇ ∇ • u b -ζ t = 0. (A 7)
The equation for the evolution of the horizontal velocity is obtained by substituting (A 5) and (A 6) into (2.12) and (2.13):

u bt + ε(u b • ∇)u b + ∇η + µ 2 h 2 2 ∇(∇ • u b ) -h∇ ∇ • (hu b ) -ζ t t = 0.
Using the irrotationality condition this equation can be rewritten as The higher-order terms are of order O(εµ 2 ). Another possibility for the choice of variables is to introduce the depth averaged velocity. The corresponding standard Boussinesq-type equations were obtained by Peregrine [1967] in the case of a fixed seabed. We extend the results to a moving bathymetry and add the energy equation. The main advantage of the depth averaged velocity consists in the fact that the continuity equation (or equivalently the equation for the free-surface elevation) is very simple and exact in this variable.

u bt + ε 2 ∇ |u b | 2 + ∇η + µ 2
Let us rewrite all equations in terms of the depth-averaged velocity defined by ū = 1 h + εη Since the energy equation is redundant for incompressible flows, the linear dispersion relation is unaffected by the inclusion of the energy equation. As is wellknown, it can be improved by defining the horizontal velocity at an arbitrary level.

Figure 1 .

 1 Figure 1. Sketch of the 3D fluid domain for wave generation by a moving bottom.

Figure 2 .

 2 Figure 2. Tsunami generation leading to a dipolar wave form; t * = 4 s. The same scales are used in Figures2-9: for the free surface, red denotes a rise of sea level; for the total energy, the scale goes from 0 (blue) to 5 × 10 4 (red); for the potential energy, the scale goes from 0 (blue) to 10 4 (red).

Figure 3

 3 Figure 3. t * = 6 s

FigureFigure 6 Figure 9

 69 Figure 4. t * = 10 s

Figure 11 .Figure 12 .

 1112 Figure 11. Time evolution of kinetic and potential energies (4.9) for three mechanisms of tsunami generation: passive generation, instantaneous bottom motion and exponential bottom motion. The square in plot (b) indicates the estimate (4.11).

  u b ) -h∇ ∇ • (hu b ) -

  The depth-integrated continuity equation (2.19) yields immediatelyη t + ∇ • (h + εη)ū -ζ t = 0. (A 11)In order to derive equations for the horizontal velocity and the energy we need a relation between u b and ū. The desired relation is deduced directly from the definition (A 10) by substituting (A 5) in it:ū = u b -µ 2 2 h∇ ∇ • (hu b ) -ζ t + µ 2 3 h 2 ∇(∇ • u b ) + O(ε 2 + εµ 2 + µ 4 ).Inverting the last equation yieldsu b = ū + µ 2 h 2 ∇ ∇ • (hu b ) -ζ t -h 2 3 ∇(∇ • u b ) + O(ε 2 + εµ 2 + µ 4 ). (A 12)Substituting the relation (A 12) into equation (A 8) gives the standard Boussinesq equations for a moving bottom:ūt + ε 2 ∇ |ū| 2 + ∇η + µ 2 h 2 6 ∇(∇ • ū) -h 2 ∇(∇ • (hū) -ζ t ) t = 0. (A 13)The energy equation is obtained by substituting the relation (A 12) into equation (

Table 1 .

 1 Values of physical parameters used for the energy density computations.

	Article submitted to Royal Society

Table 2 .

 2 Values of physical parameters used for the Cauchy-Poisson analysis of tsunami generation.

	Article submitted to Royal Society

Article submitted to Royal Society

pronounced shift of the spectral energy towards higher frequencies characterized by strong dispersion makes landslide sources significantly deficient far-field tsunami generators, as compared to classical dislocations.

The authors would like to thank Jean-Michel Ghidaglia who suggested the idea of considering Boussinesq equations with energy. We are also grateful to him for numerous fruitful discussions on numerical and fluid mechanics topics. The second author acknowledges the support from the EU project TRANSFER (Tsunami Risk ANd Strategies For the European Region) of the sixth Framework Programme under contract no. 037058.

Appendix A. Derivation of dispersive shallow-water equations with variable bathymetry

We derive the Boussinesq equations following the method used, for example, by [START_REF] Nwogu | Alternative form of Boussinesq equations for nearshore wave propagation[END_REF], [START_REF] Villeneuve | Nonlinear, dispersive, shallow-water waves developed by a moving bed[END_REF], Yoon and Liu [1989]. A moving bathymetry was considered in the framework of the Boussinesq equations by Villeneuve and Savage [1993], but they dealt with a 2D problem leading to a 1D system of equations.

We need to know the depth dependence of the horizontal velocity u ⊥ in order to reduce the problem to a 2D one. We expand u ⊥ in a Taylor series in the vertical coordinate z about the seabed z = -h:

From now on, the horizontal velocity at the bottom is denoted by

If the flow is assumed to be irrotational, there is an additional relation which closes the system:

Substituting (2.22) into (A 2) clearly shows that (u ⊥ ) z (x, y, -h, t) = O(µ 2 ). Substituting (A 1) into (2.22) and integrating then yields

The vertical velocity varies linearly with respect to z over the depth at leading order O(µ 2 ). The horizontal velocities can be found by integrating the irrotationality condition (A 2) from -h to z: