Stéphane Popinet

Adaptive modelling of long-distance wave propagation and fine-scale flooding during the Tohoku tsunami

Adaptive modelling of long-distance wave propagation

and fine-scale flooding during the Tohoku tsunami S.

Introduction

On 11th March 2011 an unexpected very large earthquake struck the northeastern shore of Honshu: the Tohoku region. While the earthquake itself caused remarkably little damage -a testament to the high quality of earthquake preparedness in Japan -the subsequent tsunami caused the loss of more than 15,000 lives and extensively damaged infrastructure, overwhelming the extended network of sea defences. This disaster dramatically underlined the vulnerability of coastal communities to hazards such as tsunamis, even in countries as well prepared as Japan. A positive note is that thanks to the extensive, long-term investment in monitoring infrastructure in Japan, this event is certainly the best ever documented large earthquake and tsunami in history. Extensive datasets exist documenting ground deformation, offshore and inshore tsunami waves, global atmospheric waves, extent of inundation and damage etc... This wealth of data will hopefully lead to improved understanding of these extreme events, and of our capability to mitigate their consequences on society. One of the scientific challenges is to be able to construct a consistent model of the sequence of events leading to the observed timeseries for very different but tightly coupled processes. In this context, the main goal of this article is to demonstrate how adaptive tsunami modelling can effectively be used to obtain a consistent model of both long-distance tsunami wave propagation and local flooding.

In adaptive models the spatial resolution is not fixed in time but varies according to the local complexity of the solution. This is different from both constant-resolution models (which use a resolution constant in space and time) and variable-resolution models (which use a resolution variable in space but constant in time). The increased flexibility of adaptive methods is particularly valuable for processes which involve a wide and time-varying range of scales. Many fundamental physical processes, such as turbulence, belong to this category and tsunami waves are a clear geophysical example of such scalings. As demonstrated in [START_REF] Popinet | [END_REF]] exploiting the "gaps" in the scale distribution of typical tsunami waves leads to a computational cost (i.e. number of grid points) which scales roughly like ∆ -1.4 , with ∆ the mesh size, rather than ∆ -2 for constant-or variable-resolution models. This reduction in the effective dimension (or "fractal dimension") of the problem leads to increasingly large gains in computational efficiency with spatial resolution.

While adaptive methods are now relatively common for industrial computational fluid dynamics applications [START_REF] Yiu | Quadtree grid generation: information handling, boundary fitting and CFD applications[END_REF], Popinet (2003), Popinet (2009)], they are comparatively less developed for geophysical fluid dynamics [START_REF] Liang | Simulation of dam-and dyke-break hydrodynamics on dynamically adaptive quadtree grids[END_REF], [START_REF] Popinet | A tree-based solver for adaptive ocean modelling[END_REF], [START_REF] Liang | Adaptive quadtree simulation of shallow flows with wet-dry fronts over complex topography[END_REF], [START_REF] Popinet | A quadtree-adaptive spectral wave model[END_REF]]. The quadtreeadaptive tsunami model used in this article, implemented within the Gerris flow solver framework is one of the few adaptive tsunami models published [START_REF] George | High-resolution methods and adaptive refinement for tsunami propagation and inundation[END_REF], [START_REF] Harig | Tsunami simulations on several scales[END_REF], [START_REF] Popinet | [END_REF]]. A summary of the model is given in the next section. Section 3 describes a classical validation test case and section 4 gives a detailed comparison between field data and model results for the Tohoku tsunami.

Numerical method

This section gives an overall summary of the numerical method. A detailed description can be found in [START_REF] Popinet | [END_REF]] and references therein.

Tsunamis are most commonly modelled using a long-wave approximation of the mass and momentum conservation equations for a fluid with a free-surface. In this approximation, the slopes of both the free-surface and the bathymetry are assumed to be vanishingly small. Although these assumptions are often violated in the case of tsunamis, particularly close to shore, the long-wave approximation still gives reasonable results in practice. Better approximations include the Boussinesq wave equations which tend to be significantly more expensive to solve [START_REF] Grilli | Source constraints and model simulation of the December 26, 2004, Indian Ocean Tsunami[END_REF]].

The long-wave approximation can be obtained either by formal vertical averaging of the Navier-Stokes equations with a free-surface or by considering mass and momentum conservation through vertical slices in the water column, as first done by [START_REF] De Saint-Venant | Théorie du mouvement non-permanent des eaux, avec application aux crues des rivières et à l'introduction des marées dans leur lit[END_REF]]. The resulting equations of motion are often written in differential form as the set of non-linear advection equations

∂ t h + ∂ x (h u) + ∂ y (h v) = 0 ∂ t (h u) + ∂ x (h u 2 + 1 2 g h 2 ) + ∂ y (h u v) = -h g ∂ x z ∂ t (h v) + ∂ x (h u v) + ∂ y (h v 2 + 1 2 g h 2 ) = -h g ∂ y z
with h the free surface elevation, u = (u, v) the horizontal velocity vector, g the acceleration of gravity and z the depth of the bathymetry. This corresponds to a hyperbolic system of conservation laws which, using Stokes' theorem, can be expressed in integral form as

∂ t Ω q dΩ = ∂Ω f (q) • n d∂Ω - Ω h g ∇z, (1)
where Ω is a given subset of space, ∂Ω its boundary and n the unit normal vector on this boundary. For conservation of mass and momentum in the shallow-water context, Ω is a subset of bidimensional space and q and f are written

q =   h h u h v   , f (q) =     h u h v h u 2 + 1 2 g h 2 h u v h u v h v 2 + 1 2 g h 2     (2)
Dissipative terms (viscosity and/or friction) are also often added but they do not change the fundamental structure of the equations. This system of conservation laws is hyperbolic and admits an entropy inequality, which guarantees that the total energy h u 2 /2 + g h 2 /2 behaves physically. Another important property is that the steady state

u = 0, h + z = constant, (3)
is a solution of this system. This simply reflects the expected property that a lake at rest remains so. Numerical methods for hyperbolic systems of conservation laws have been extensively developed in the context of compressible gas dynamics with applications in aeronautics in particular.

Many of these schemes are suitable for solving the Saint-Venant equations provided two important specificities are taken into account: the solver needs to properly account for vanishing fluid depths (equivalent to a vanishing density i.e. a rarefaction wave) and needs to verify exactly the lake-at-rest equilibrium solution (i.e. be a "well-balanced" scheme). The Saint-Venant solver used in this article is based on the numerical scheme analysed in detail by [START_REF] Audusse | A fast and stable well-balanced scheme with hydrostatic reconstruction for shallow water flows[END_REF]] which verifies both properties. The scheme uses a second-order, slopelimited Godunov discretisation. The Riemann problem necessary to obtain Godunov fluxes is solved using an approximate HLLC (Harten-Lax-van Leer Contact) solver which guarantees positivity of the water depth. The choice of slope limiter is important as it controls the amount of numerical dissipation required to stabilise the solution close to discontinuities (i.e. hydraulic jumps and contact discontinuities). Based on the study in [START_REF] Popinet | [END_REF]], the Sweby limiter seems to be a good compromise between stability of short waves and low dissipation of long waves. A MUSCL-type discretisation is used to generalise the method to second-order timestepping. The timestep is constrained by a classical CFL condition

∆t < 1 2 min ∆ |u| + g h √
with ∆ the mesh size. To summarise, the overall scheme is second-order accurate in space and time, preserves the positivity of the water depth and the lake-at-rest condition and is volumeand momentum-conserving. This scheme is implemented within the Gerris Flow Solver framework [Popinet (2003), Popinet (2003Popinet ( -2011))] which provides a number of additional capabilities including: quadtreebased adaptive discretisation, generalised orthogonal coordinates [START_REF] Popinet | [END_REF]] (used for longitude-latitude discretisation in this article) and parallelism [START_REF] Agbaglah | Parallel simulation of multiphase flows using octree adaptivity and the volume-of-fluid method[END_REF]]. An example of the quadtree structure is represented in Figure 1 together with its logical (tree) representation. This tree can be conveniently described as a "family tree" where each parent cell can have zero or four children cells. An important parameter is the level of a given cell in the tree. The root cell has level zero by convention and the level increases by one for each successive generation in the tree. The level of the cells in the tree are also given.

Quadtree discretisations are a good compromise between the simplicity of (inflexible) regular Cartesian meshes and the flexibility of (complex) unstructured meshes. New cells are easily added or removed from the hierarchy and fine-grained control of the spatial resolution is possible (in contrast to block-structured AMR for example [START_REF] George | High-resolution methods and adaptive refinement for tsunami propagation and inundation[END_REF]). Interpolation on newly-refined or coarsened cells is also designed to ensure conservation of quantities.

In the context of tsunami modelling, using adaptive meshes also requires a technique for efficient reconstruction of the detailed bathymetry. With spatial resolutions varying over several orders of magnitude (e.g. from 250 kilometres down to 250 metres for the results in this article) it is generally not possible or efficient to reconstruct and store the bathymetry at the finest resolution in main memory. A terrain database system has been developed specifically for this task within Gerris. It also relies on a tree-like hierarchy for fast access to very large bathymetric datasets. In constrast to the quadtree structure stored in main memory, the bathymetry database relies on a 2d-tree data structure [START_REF] Bentley | Data structures for range searching[END_REF]] mostly stored on-disk so that the size of the database is disk-and not memory-limited. With this data structure the cost of constructing the average depth of a discretisation cell scales like O( N √ ) + O( m √ ) where N is the total number of samples in the database and m is the number of samples contained in the cell [START_REF] Popinet | [END_REF]]. In practice the typical cost of adaptive bathymetry reconstruction is roughly 10% of the total cost.

Tsunami runup onto a plane beach

The implementation of the Saint-Venant solver within Gerris has been validated both with simple test cases and complex tsunami models [START_REF] Popinet | [END_REF], Popinet (2010a), [START_REF] Popinet ; Popinet | Gerris -Monai tsunami benchmark[END_REF], [START_REF] Popinet | [END_REF]]. For the sake of completeness, this section describes a further test case which has not been previously published. The test case was proposed at the Third International Workshop on Long-Wave Runup Models [Liu (2004)] and seeks to validate the Saint-Venant solution of a plane wave running up a sloping beach. The geometry is chosen to replicate the spatial scales of the vertical cross-section of a typical (large) tsunami wave. The problem is thus one-dimensional for the (depth-averaged) Saint-Venant equations. The domain is 60 kilometres long, the slope of the beach is 1/10. The initial wave profile is given numerically [IWLWRM ( 2004)] and has maximum and minimum amplitudes of + 3 and -8.8 metres respectively, located approximately 20 km and 8 km offshore respectively. The evolution with time of the wave profile close to the shore is illustrated in Figure 2. The (semi)-analytical solution is obtained using the initial-value-problem (IVP) technique introduced by Carrier et al [START_REF] Carrier | Tsunami run-up and draw-down on a plane beach[END_REF]] and is given numerically on the web site [START_REF] Iwlwrm | Third international workshop on long-wave runup models: Tsunami runup onto a plane beach[END_REF]]. The initial wave minimum is located closer to the shore than the maximum so that the first wave is negative with a maximum dryout of the beach at around t = 175 seconds, the positive wave then runs up the beach and reaches a maximum elevation of approximately 15 metres at t = 220 seconds. For this test case only static mesh refinement is used. The spatial resolution is increased closer inshore according to

l(x) = l min + (l max -l min ) x/L
where l is the level of refinement, x is the distance to the shoreline, L = 50 km, l min = 7 and l max is varied between 11 and 14 to study the convergence of the solution with spatial resolution. It is clear from Figure 2 that the numerical solution at the maximum resolution is a good approximation of the analytical solution. A quantitative assessment of the quality of the numerical solution is given in Figure 3. The average (L1-norm) and maximum errors at t = 220 seconds are computed as

e 1 = i ∆ i |y i -y(x i )| i ∆ i e max = max i |y i -y(x i )|
where (x i , y i ) is the discrete numerical solution for the wave profile, y(x i ) is the (interpolated) analytical solution at position x i and ∆ i is the (variable) grid size. The average norm shows first-order convergence and the maximum error slightly less than first-order convergence. This is not unexpected given that errors are dominated by the contact discontinuities at the wet-dry transition, where the Godunov scheme (or any other shock-capturing scheme) reduces to firstorder accuracy. From a tsunami modelling perspective, this test case demonstrates that, even for very simple coastline geometries, spatial resolutions of less than a few tens of metres are necessary to accurately capture wave runup and inundation on the shoreline. All the data necessary to reproduce this test case are given on the Gerris web site [Popinet (2011a)].

Results for the Tohoku tsunami

The adaptive model was used to simulate the 11th March 2011 Tohoku tsunami. The simulation domain was chosen so as to validate both long distance tsunami propagation and high-resolution local inundation. The Saint-Venant equations are solved on a longitude-latitude grid centered on the epicentre of the earthquake (142.597 • E, 38.486 • N) and spanning 73 • in both latitude and longitude (Figure 4). The topography is obtained from two datasets: the GEBCO_08 30 arcsecond dataset [IOC (2011a)] and the SRTM 3 arc-second dataset [NASA ( 2011)] (for the Japanese islands only). Both datasets are stored using the 2d-tree (KDT) hierarchical terrain database system of Gerris with on-disk storage sizes of:

• GEBCO_08: global dataset, 933 million points, 21 GBytes, KDT index: 65 MBytes

• SRTM, 3 arc-second: Japan only, 58 million points, 1.4 GBytes, KDT index: 4.1 MBytes These on-disk datasets are accessed during the simulation to reconstruct the bathymetry as required as the spatial resolution evolves. The size of the KDT index gives an order of magnitude estimate of the amount of dynamic memory (RAM) required to access each dataset. A critical part of any tsunami model is the choice of an appropriate initial surface elevation (the "source model"). In this study, I have chosen not to use a source model based on inversion of DART buoys and tidal gauges timeseries but rather based directly on seismic inversions. The rationale for this choice is that this allows independent validation of the source model using the DART buoys and tide gauges data, and also that, in a realtime forecasting scenario, using only seismic information may allow a quicker estimate of the source (because seismic waves travel faster than gravity waves).

Several deformation models for the Tohoku earthquake have been made available by the seismic research community, and some of them were released only a few hours after the event. A selection of such models is available on the USGS Tohoku web site [ USGS (2011)]. Different methods can be used for the seismic inversion (combined with different weighting for the various seismic datasets) and the scatter in the resulting predicted deformations can be large. For this study I have used one of the deformation models obtained by Shao, Li and Ji, UCSB, more precisely their preferred "Preliminary model III" released on March 14th 2011 [START_REF] Shao | Preliminary Result of the Mar 11, 2011 Mw 9.1 Honshu Earthquake[END_REF]]. This model predicts vertical deformations which are significantly larger than any of the other models mentioned on the USGS Tohoku web site and was the only source model I tested which resulted in satisfying agreement with tsunami wave measurements. The initial vertical wave ele-vation was assumed to be equal to the total vertical bottom displacement predicted by the source model. This displacement was obtained by superposition of the 190 Okada subfaults displacements used in the seismic inversion. The fault model also includes the relative rupturing times of each subfault, but I chose to assume that the total rupturing time (3 minutes) was small enough to be negligible. The resulting initial vertical displacement is illustrated in Figure 5. For the adaptive model, the choice of spatial resolution is very flexible and can be based on several simultaneous criteria, depending on which specific information one would like to extract from the model. In this study I chose to illustrate how adaptivity can help obtain simultaneous solutions for two aspects of tsunami modelling which are generally considered as conflicting in term of resolution requirements: long-distance propagation and detailed inundation/damage predictions. In this study, both regimes are captured using a single criterion but different maximum spatial resolutions: the mesh size ∆ is adapted so that

|∇h| ∆ < ǫ (4)
where ǫ is a control parameter set to 2.5 cm. ǫ can be interpreted as the maximum truncation error in the free surface elevation of a spatially first-order discretisation scheme. This corresponds to the maximum error expected near sharp discontinuities (i.e. shocks). In smooth regions, the scheme is second-order and errors will be significantly smaller. There are primarily two regions where one can expect larger discretisation errors: wave fronts and "breaking" waves in shallow water i.e. the two regimes of long-distance wave propagation and inundation we are interested in. I have shown in a previous study that adequate propagation of long-wave energy over large distances with the non-linear, shock-capturing numerical scheme requires spatial resolutions of the order of one arc-minute in order to minimise numerical wave dissipation [START_REF] Popinet | [END_REF]]. Following this, the maximum resolution for deep water wave propagation is set to 12 levels of quadtree refinement, giving a maximum deep water resolution of close to 1 arc-minute (which is also comparable to the spatial resolution of the deep water GEBCO bathymetry). Detailed inundation modelling, on the other hand often requires higher spatial resolution, in order to capture fine-scale topographic effects as well as wave runup (see the previous section). In this study the maximum resolution of "inundable areas" (any part of the domain above sea level) is set to 15 levels of refinement (or 8 arc-seconds, or ≈ 250 metres) which is comparable to the spatial resolution of the SRTM topographic dataset. Note that this choice of adaptive refinement does not rely on any a priori assumptions on the location of inundated areas and/or propagation of the tsunami wave and could thus be applied to any other event. The variable resolution is also used to implement "sponge layers" on the open boundaries of the domain.

Within bands approximately 3.65 degrees wide along all four boundaries of the domain, the spatial resolution is coarsened to only 5 levels (approximately 2.3 • resolution). The resulting numerical dissipation dampens waves before they exit the domain, thus minimising spurious wave reflections at boundaries.

With this particular choice of spatial resolutions, the global timestep is constrained by the CFL condition applied to gravity waves propagating in deep water (not by the highly-resolved waves propagating in inundated areas) and is of the order of 2 seconds. The initial evolution of the wave and mesh is illustrated in Figures 6 and7. Figure 7 details the evolution of inundation and the associated adaptive mesh refinement in the Miyagi prefecture (Sendai) area which was highly impacted by the tsunami. The detail shown covers an area of roughly 220 × 180 km. Figure 8 gives the corresponding evolution in the total number of grid points over 10 hours of modelled physical time. The initial number of grid points is approximately 50,000 which corresponds to the number of elements necessary to resolve the source of Figure 5 with the accuracy specified by (4). After a rapid initial increase, the number of grid points reaches a max-imum of approximately one million at t = 5 hours, the time at which waves start exiting the computational domain. This can be compared with 2 24 , the total number of grid points which would be necessary if the whole domain was resolved at a constant resolution of ≈ 1 arc-minute (the resolution necessary to model deep-water wave propagation correctly), as well as an upper bound of 2 30 ≈ one billion, if the whole domain was resolved at the maximum inundation resolution of 8 arc-seconds. A somewhat more realistic upper bound can be estimated by assuming that a static "nested mesh" approach could be used to resolve only low-lying (altitude less than e.g. 50 metres) emerged land areas at the maximum resolution (8 arc-seconds) and the rest of the domain at a resolution of one arc-minute. This still gives about 50 million grid points, so that a fair estimate of the gain in number of grid points obtained with the adaptive method for this particular simulation is a factor of order 15 to 50. Note also that this factor is scale-dependent, so that as pointed out in [START_REF] Popinet | [END_REF]] it will increase with increasing maximum spatial resolution. The gain in runtime is more difficult to estimate as it depends on the details of the respective implementations and optimisations of different methods. Previous comparisons show that the current quadtree implementation in Gerris is roughly five times slower than a comparable numerical scheme implemented on a regular Cartesian grid so that the gain in runtime would be a factor from 3 to 10. Note however that this also assumes that the speed of the regular Cartesian grid scheme would not be affected by the large memory requirements, which is unlikely in practice. The results presented here were obtained on a single-processor standard desktop PC with an average runtime of roughly two hours for one hour of physical time (this ratio is much lower in the initial phases of the simulation where few grid points are used). This runtime could be decreased by using the parallel capabilities of Gerris (including dynamic parallel load-balancing) [START_REF] Agbaglah | Parallel simulation of multiphase flows using octree adaptivity and the volume-of-fluid method[END_REF]]. A general overview of the maximum wave elevation reached over 10 hours is given in Figures 9, 10 and 11. Maximum wave elevations in excess of 40 metres are predicted along the coast of the Oshika peninsula . These extreme values extend to the north of Ofunato, roughly between latitudes 38 • N and 39.5 • N. Large values are otherwise predicted along the Japanese coast from 35 • N to 42 • N (Figure 10). Offshore, the maximum amplitude is reached along an east-south-east direction with strong modulations due to the many small submerged and emerged topographic features in the Pacific (Figure 9). A first validation of the model is obtained by comparing the modelled wave elevation timeseries with observed sea level at DART buoys locations. Figure 12 illustrates these timeseries for all the DART buoys deployed in the simulation domain (see Figure 9) which were active during the event and which measured deviations larger than 10 cm. The DART timeseries were obtained directly from the raw data available in real-time on the DART website [NOAA (2011)] and were low-pass filtered via Fourier decomposition to remove the tidal components. The model correctly predicts both the arrival times and maximum wave elevations at all locations (keeping in mind that DART observations were not used to constrain the source model). The initial wave front is steeper than observed for most stations however. This was also observed in a previous study [START_REF] Popinet | [END_REF]] and may be due to excessive anti-diffusion of the limiter in the shock-capturing scheme. It is also clear that the model reproduces many of the high-frequency, secondary wave reflections. For example, the clear modelled and measured signal around t = 4 hours for DART 21419 is most likely due to the reflection of the primary wave off the Emperor seamount chain (the underwater range extending south-east from DART 21416 in Figure 4). While the DART buoys data validate the long-distance, deep-water propagation, tidal gauges can give information on the inshore wave propagation. In this study, I chose to use only data available within the Global Sea Level Observing System (GLOSS) [IOC (2011b)]. The GLOSS database system provides a unified way to access hundreds of sea level recording systems, in near real-time, which greatly simplifies the implementation of an automated tsunami forecasting system. Unfortunately, many existing tide gauges have not been integrated within GLOSS yet. For the area I consider only six sites are available for the Japanese islands. Another limitation of the GLOSS database is that the location of tide gauges is only recorded with an accuracy of one arc-minute. This can easily lead to tide gauges being incorrectly located in "dry" areas of the numerical model. For this study the locations of gauges have been adjusted to their most probable locations based on bathymetry and satellite images.

Figure 13 and 14 give a comparison between the numerical model and measurements for all the GLOSS stations in the simulation domain of Figure 4 which were operating at the time of the event and which measured deviations from tides larger than 10 cm. Predictions of tsunami timeseries at tidal gauge locations is notoriously more difficult than for open ocean buoys, as it is well known that measured tidal timeseries are often dominated by the resonant modes (seiches) of the small-scale bays and harbours where the gauge is located [START_REF] Takahashi | [END_REF]Aida (1963), Miller (1972)]. Keeping this in mind, it is clear that the arrival time and amplitude of the first wave are very well predicted for all locations of Figure 13. This is particularly true of the stations at Hanasaki, Naha, Saipan and Pago bay and less so for Ofunato, Omaezaki, Ishigakijima and Wake island. Interestingly, the Ofunato tide gauge measured a positive first wave before ceasing transmission at t ≈ 30 minutes. The model predicts a negative first wave which is clearly related to the negative initial amplitude distribution of the source model (Figure 5). The characteristics of the following waves are very site-specific but can be roughly classified according to their dominant frequency. Some sites such as Hanasaki or Ishigakijima have dominant waves at relatively long periods ( ≈ 1 hour) which are characteristic of large-scale resonances from the coast to the edge of the continental shelf. These sites are located on open coastlines and are thus well-represented by the model and associated bathymetry. Other sites such as Naha and Omaezaki are located within harbours protected by breakwaters which are not represented at the maximum model resolution of 250 m (and in the topographic data used). As demonstrated by [Abe (2011)], for sites such as Omaezaki the detailed geometry of the harbour is crucial in determining the amplification of short-period components. This could explain the under-estimation of wave amplitudes for these sites. Finally, the sites at Wake island and Pago bay are characterised by very high frequencies (periods smaller than 15 minutes). Both the characteristic frequencies and amplitudes of the measured waves are well captured by the model. Both sites are located on small, isolated islands surrounded by deep ocean. The origin of the observed and modelled high-frequency waves is readily apparent on Figure 6 (bottom row). They are due to multiple scattering of the primary wave off the many seamounts and atolls of the north-west Pacific. These high-frequency signatures are also present (at lower amplitudes) in the modelled and measured DART buoys timeseries (Figure 12). Prediction of the extent of tsunami inundation as a function of detailed topography is an important aspect of tsunami modelling, particularly when used for urban planning and risk assessment. Extensive survey work after the Tohoku tsunami done in very difficult conditions by the Tohoku Earthquake Tsunami Joint Survey Group [TTJT (2011b), TTJT (2011a), [START_REF] Mori | Survey of 2011 Tohoku earthquake tsunami inundation and run-up[END_REF]] as yielded an impressive amount of data on this event: more than 5000 individual GPS records of wave height have been collected along the entire eastern Japanese coastline. In addition to this field data, various satellites have also recorded different aspects of the tsunami, including detailed maps of the extent of flooding one day after the event. Figure 15 (left) displays a detail of the topography contours (red), survey data points (blue) and extent of flooding (green) estimated using analysis of high-resolution Synthetic Aperture Radar (SAR) imaging [DLR (2011)] for the Miyagi prefecture area.

In the southern half of the map, the flooding extent is clearly controlled by the relatively steep topography, and less so for the northern half. Figure 15 (right) shows a comparison between the satellite estimation of the flooding extent on 12th March with both the modelled maximum extent of the tsunami (green) and the modelled inundated area (depth larger than 10 cm) at t = 10 hours (black). Both sets of curves are a satisfactory approximation of the satellite data, with the maximum extent curve (green) lying beyond the satellite data as expected, and the 10 cm flooding curve (black) generally underestimating slightly the extent of flooding compared with the satellite estimate. These curves are available as Google Earth KML files for the entire simulation domain of Figure 4 [ Popinet (2011b)] and can thus be compared visually with satellite images of the damaged areas. This also confirms the good model predictions of the extent of tsunami damage. Finally Figure 16 illustrates a statistical comparison of the surveyed and modelled wave elevation for all locations in the survey dataset. The survey data covers an area extending from 32 • N to 44 • N. The comparison in Figure 16 assumes that the wave elevation is mostly correlated with latitude which is broadly confirmed by the reasonably narrow range between the lower and upper quartiles (shaded areas). The model predictions are in broad agreement with the survey data, both for the qualitative wave distribution along the coastline and the elevations reached. Very good agreement is obtained between 36 and 37.5 • N as well as north of 42 • N (on the coast of Hokkaido where the Hanasaki tide gauge is located, Figure 13). The extreme wave elevations (peak values in excess of 40 metres) observed between 38.5 and 40.5 • N are also captured by the model with the correct amplitudes, but the modelled extremes occur about one degree south of the observed values: this corresponds to the Ofunato and Miyako regions respectively (see Figure 11). Both regions are characterised by very steep and indented coastlines which, together with proximity to the source, explains the high wave elevations reached there. The sources of error for both datasets (surveyed and modelled) are of course different: the model will suffer from a representation bias due for example to the inaccurate representation of steep terrain (even at the maximum resolution of 250 metres); while the surveyed data is likely to suffer from a sampling bias due to the difficulty of surveying every location and/or due to particular choices in the sampling procedure (e.g. focusing on tsunami damage to infrastructure etc...). Such a sampling bias could explain the discrepancy between modelled and surveyed data for the Miyagi area (corresponding to the inundation maps of Figure 15). The surveyed data for this area seems to be biased low (with many samples even reporting negative wave elevations) while the modelled elevation distribution seem to be more consistent with the observed extent of inundation.

Model or survey biases alone cannot however explain the large latitudinal shift in maximum wave elevation between the Ofunato and Miyako regions. The modelled distribution is consistent with the source model: the Oshika peninsula and Ofunato where the maximum wave elevations are predicted lie on a line perpendicular to the fault direction and going through the maximum of initial wave amplitude (Figure 5). In contrast the Miyako area lies at the edge of the area of maximum deformation. A possible explanation for the discrepancy is thus that the source model is locally inaccurate and that a significant fraction of the rupture energy travelled in a more northerly direction toward the Miyako coastline. Assuming that the total energy released should be unchanged, this would also divert some of the energy away from the Ofunato area. The short record from the Ofunato tide gauge (top graph of Figure 13) also suggests that the source model may be locally inaccurate: the observed first wave is positive whereas the modelled first wave is negative (with an amplitude of -5 metres) which is consistent with the large negative amplitude area of the source model west of the epicentre. To try to clarify this issue, the model output is compared in Figures 17 and18 with wave height records from a further six wave buoys and two bottom pressure gauges located a few tens of kilometres offshore the Sanriku coastline (locations indicated on Figure 5). An obvious discrepancy is that all the modelled timeseries show a large negative first (or second) wave which is absent in the field records. As noted before for the Ofunato tide gauge, this is consistent with the large negative amplitude in the source model (south of TM1 in Figure 5). For most sites, the peak amplitude of the largest positive wave (and subsequent waves) is correctly captured by the model with the notable exceptions of buoys 804, 803 and 801. The peak wave amplitude is under-estimated by a factor of two for buoy 804 and over-estimated by a similar factor at buoys 803 and 801. This is consistent with the bias in predicted wave elevations of Figure 16: buoys 803 and 801 are located offshore the Ofunato region (positive bias) whereas buoy 804 is located offshore the Miyako region (negative bias). This confirms that the biases in predicted wave elevations on the coastline are likely due to local errors in the source model rather than an inaccurate representation of the dynamics of flooding on this steep terrain.

Note also that several other source models proposed for this event suffer from a similar north-south bias. Using for example the source model proposed by Fuji et al [START_REF] Fuji | Tsunami source of the 2011 off the Pacific coast of Tohoku, Japan earthquake[END_REF]], based on inversion of tide and pressure gauges and wave buoys (including all the sites of Figures 17 and18) still leads to a large under-prediction of the wave amplitude and runup in the Miyako region (but does reduce the amplitude in the Ofunato region). Petukhin et al also report a very similar bias for their source model based purely on seismic inversion (see Figure 3 of [START_REF] Petukhin | Simulation of tsunami and longperiod ground motions during the M9.0 2011 Tohoku-oki earthquake[END_REF]) and suggest that this may be due to an inaccurate representation of the bathymetry for this area. Based on the results presented here, I do no think that this is a sufficient explanation since both the northern and southern parts of the Sanriku coast are described with similar datasets and do not differ markedly in term of complexity. The discrepancy is more likely due to the interaction of the main wave -which looks correctly described given the agreement with far-field measurements -with the initial wave field between the epicentre and the coast.

Conclusion

The quadtree-adaptive Saint-Venant solver implemented within Gerris was shown to provide robust and accurate solutions for the full range of wave phenomena associated with the Tohoku tsunami. The adaptive spatial resolution ranged from approximately 250 metres to 250 km for a domain spanning 73 • of latitude and longitude. The maximum spatial resolution proved adequate to resolve inundation processes in most areas while adaptivity allowed to track deep-water wave fronts at a typical resolution of 1 arc-minute ( ≈ 1.8 km) which was sufficient to ensure good long-distance wave energy conservation. A conservative estimate of the gain in computational efficiency obtained with the adaptive model is a factor of ten compared to a discretisation using a spatially-variable but constant-in-time resolution. This estimate increases to several orders of magnitude when compared to simple, spatially-constant (e.g. Cartesian) discretisations. Furthermore, as pointed out in [START_REF] Popinet | [END_REF]] this gain is scale-dependent and will increase with maximum spatial resolution. The criteria used for adaptive mesh refinement are simple and do not rely on a priori assumptions about the solution. They are thus applicable to any tsunami event.

As in most tsunami models, the quality of the initial vertical displacement field is crucial to the accurate prediction of the subsequent wave evolution. The source model of Shao, Li and Ji was shown to lead to remarkably accurate predictions of initial wave amplitudes and arrival times as well as subsequent multiple reflections from bathymetry for all the DART buoys located within the simulation domain. Shao, Li and Ji only used information on seismic waves for their source inversion, so that comparisons between the modelled and observed surface waves are a true independent validation of both source and tsunami model consistency. This study demonstrates that accurate and consistent source models for tsunami wave propagation can be constructed using seismic waves only. Coupled with the high level of automation readily possible with the adaptive model, this opens the door to realtime forecasting of tsunami propagation and inundation, where the main time constraint is related to acquisition and processing of the fast (a few km/s) seismic waves signals rather than the slower (a few hundreds of m/s) surface gravity waves signals.

The results for deep-water waves are confirmed by comparisons with in-shore GLOSS stations timeseries. All stations show satisfactory agreement when allowing for the well-known dependencies of in-shore signals on (unresolved) topographic details such as breakwaters. Model results also compare well with satellite estimates of the extent of flooding as well as with extensive point survey data along the eastern Japanese coastline. However there is a significant discrepancy between surveyed and modelled wave elevations on the Sanriku coastline, where extreme wave elevations (up to 40 metres) were observed. These extremes values are reproduced in the model but further south than observed. Analysis of further inshore wave buoys and pressure gauge records indicate that this discrepancy is most probably due to an incorrect incoming wave rather than inaccurate modelling of inundation on the complex and steep shoreline. A similar discrepancy has been noted by other authors, using different source models [START_REF] Fuji | Tsunami source of the 2011 off the Pacific coast of Tohoku, Japan earthquake[END_REF]), Petukhin et al. (2011)]. This may indicate that our understanding of the local seafloor deformations during the Tohoku earthquake is incomplete. A potential mechanism which could explain the discrepancy is the occurence of local submarine landslides (i.e. plastic rather than elastic deformations) triggered by the earthquake.

Figure 1 .

 1 Figure 1. An example of quadtree discretisation (left) together with its logical representation (right).The level of the cells in the tree are also given.

Figure 2 .

 2 Figure2. Topography, analytical (lines) and numerical (symbols) wave profiles for the times indicated in the legend. The numerical results are for l max = 14 i.e. a maximum spatial resolution of ≈ 3.6 metres. The numerical solution is represented using only one in every six discretisation points for clarity.

Figure 3 .

 3 Figure 3. Convergence of the average and maximum errors with spatial resolution. The errors are for the predicted wave profile at t = 220 sec.

Figure 4 .

 4 Figure 4. Extent of the simulation domain and bathymetry (altitude in metres). DART buoys (black labels) and GLOSS stations (red labels).

Figure 5 .

 5 Figure 5. Initial vertical displacement (metres) derived from the subfaults model of Shao, Li and Ji, UCSB (Preliminary model III). The location of wave buoys (801-807) and bottom pressure gauges records (TM1 and TM2 located between TM1 and 802) used in this study are also indicated.

Figure 6 .

 6 Figure 6. Wave evolution over most of the modelled domain. Rows from top to bottom, t = 1, 2, 4 hours. Left column: evolution of wave elevation, dark blue < -1 metre, dark red > 2 metres. Right column: corresponding evolution of spatial resolution, dark blue 5 levels, ≈ 2.3 • , yellow 12 levels, ≈ 1 arc-minute, dark red 15 levels, ≈ 8 arc-seconds, ≈ 250 metres.

Figure 7 .

 7 Figure 7. Detail of the wave evolution in the Miyagi prefecture area. The black line in all figures is the coastline. Rows from top to bottom, t = 1, 2, 4 hours. Left column: evolution of wave elevation (colorbar in metres). Right column: corresponding evolution of spatial resolution, yellow 12 levels ≈ 1 arc-minute, dark red 15 levels ≈ 8 arc-seconds ≈ 250 metres.

Figure 8 .

 8 Figure 8. Evolution with physical time of the total number of grid points.

Figure 9 .

 9 Figure 9. Maximum wave elevation reached over 10 hours. The (logarithmic) colorbar gives the maximum wave elevation in metres. The low values close to the boundaries of the domain are due to the mesh coarsening used to limit spurious wave reflections.

Figure 10 .Figure 11 .

 1011 Figure 10. Detail of the maximum wave elevation reached over 10 hours. The (logarithmic) colorbar gives the maximum wave elevation in metres. The model predict maximum wave elevations of more than 40 metres on the coast of the Oshika peninsula (south of Ofunato).

Figure 12 .

 12 Figure 12. Timeseries of observed and modelled wave height at DART buoys locations, ordered by increasing arrival time from top to bottom. Horizontal axis is time in hours.

Figure 13 .

 13 Figure 13. Measured and modelled timeseries at GLOSS stations locations. Horizontal axis is time in hours. Distances from the epicentre are indicated in the legends.

Figure 14 .

 14 Figure 14. Measured and modelled timeseries at GLOSS stations locations (continued). Horizontal axis is time in hours. Distances from the epicentre are indicated in the legends.

Figure 15 .

 15 Figure 15. Maps of inundation extent for the Miyagi prefecture area. Left: 0 and 10 metres contours (red), SAR estimate of inundation extent on 12 March 2011 (green), field survey data points (blue). Right: SAR inundation estimate (green), maximum modelled inundation extent (red) and modelled inundated areas at t = 10 hours (black).

Figure 16 .

 16 Figure 16. Wave elevation reached as a function of latitude. For each dataset (modelled and surveyed) the curve is the median elevation in each 0.2 • -wide bin, and the shaded areas delimit the lower and upper quartiles of the wave elevation distribution. The vertical lines indicate the different regions displayed in Figure 11.

Figure 17 .

 17 Figure 17. Comparison between recorded and modelled wave amplitudes for the inshore locations indicated in Figure 5. The stations are ordered by latitude from top to bottom. The horizontal axis is time in minutes.

Figure 18 .

 18 Figure 18. Comparison between recorded and modelled wave amplitudes for the inshore locations indicated in Figure 5 (continued). The stations are ordered by latitude from top to bottom. The horizontal axis is time in minutes.

Numerical method

Section 4

Section 4

Section 4

Section 4


M A Herrada

J M Lopez-Herrera

Alfonso M Gañán-Calvo

E J Vega

J M Montanero

S Popinet

J M López-Herrera

Numerical simulation of electrospray in the cone-jet mode

Keywords: numbers: 47.61.Jd, 47.65.-d, 47.55.-t, 47.15.-x Electrohydrodynamic, EHD, Electrospray, VOF

ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

I. INTRODUCTION

The use of electrostatic fields to spray liquids is a procedure that has been used since the XVIIth century [START_REF] Gilbert | De Magnete[END_REF]. In the most typical configuration, the liquid is injected at a small flow rate through a metallic needle where an electric potential of the order of kilovolts is imposed with respect to a grounded electrode situated downstream. Although many electrospray modes have been identified [START_REF] Jaworek | [END_REF], the one which has attracted most attention is the cone-jet mode, mainly due to its success as a continuous source of ions [3]. In this mode, the meniscus adopts a stable quasi-conical shape from whose apex a stationary thin jet is ejected. The jet eventually breaks up into droplets due to the capillary forces, which gives rise to a charged mist.

All theories describing the physics behind the cone-jet mode distinguish three regions [4,5]: (i) the upstream meniscus region attached to the needle, where ohmic conduction is the dominant mechanism of charge transport and the voltage drop is negligible; (ii) the downstream jet region, where the relevant charge transport mechanism is surface convection; and (iii) an intermediate transition region (the cone neck), where the charge transport mechanism evolves from conduction-dominant to convection-dominant. It is generally accepted that both the electric current circulating through the system and the jet's diameter are essentially governed by this neck region. Two theoretical approaches have been proposed to model this critical region: one considers electrical relaxation phenomena in this region [6], and the other assumes that the electric charges are relaxed over the entire interface, including the cone neck region [4,7,8]. They propose different scaling laws for both the jet's size and issued electric current. Experimental data for the size of the resulting droplets indicate that the jet's diameter scales with the governing parameters as indicated by the later theory [9]. However, there have as yet been no experimental data on the neck clearly in favour of one or the other of these two approaches. There have been several numerical studies of electrospray in the literature. The first numerical simulation was due to Hartman et al. [10] who used a one-dimensional model to reach steady cone-jet configurations by successive iterations. Their study focused on describing the meniscus and the neck regions. Yan et al. [11] refined that model by solving the axisymmetric Navier-Stokes equations and incorporating grid adaptation into their calculations. An additional constraint imposed in these two studies, even though it is not common operational practice, is that the jet is forced to hit a flat grounded electrode located downstream. Higuera [12] examined numerically the neck region in order to analyze the electric current circulating through the system. The presence of the needle and the breakup region were obviated by using as up-and down-stream boundary conditions far-field asymptotic behaviors for the conical meniscus [START_REF] Taylor | Proc. Roy. Soc. Lond. A[END_REF][START_REF] Barrero | [END_REF] and the jet [4], respectively. In all these studies, the flow was assumed to be strictly steady, thus neglecting the growth of perturbations on the free surface and the subsequent jet breakup. Collins et al. [15] proposed a more sophisticated numerical scheme to study the entire electrohydrodynamic (EHD) tip streaming flow, including the formation of the cusp, the emitted jet, and the resulting droplets. Despite the disparity of the length scales present in this phenomenon, it was accurately described by this numerical approach. Reznik et al. [16] studied the temporal evolution of an electrified compound core-shell droplet towards the emission of a coaxial jet. All the aforementioned numerical studies have in common that, in the bulk, the electrical relaxation time was short as compared to any hydrodynamic time present in the process, and thus the liquid bulk was free of electrical charges, which were only allowed to exist at the interface. Consequently, the surface charge density σ e was regarded as an unknown governed by the surface balance equation [6] d(2πF

u s σ e ) ds = 2πF KE ni , ( 1
)
where s is the distance along the interface, F and u s are the interface radial position and velocity, respectively, K is the liquid electrical conductivity, and E ni is the normal component of the electric field on the inner side of the interface.

In the present work, a robust and computationally efficient numerical method is developed to simulate steady EHD atomization processes. As in previous approaches [10][11][12]15],

the electric charges are confined to the interface. The results are compared with those calculated with Gerris, an open source free software library developed by Popinet [17,18] and extended by López-Herrera et al. [19] to study EHD phenomena. Experiments are conducted to partially validate the numerical results. The cone-jet transition is examined very close to the stability limit to assess the validity of the hypothesis on the charge behavior.

The rest of the paper is organized as follows. The problem is formulated in Sec. II.

Sections III and IV give the descriptions of the numerical schemes and the experimental method, respectively. Finally, the results are presented and discussed in Sec. V.

II. FORMULATION OF THE PROBLEM

The geometrical and electrical configurations considered in this work are represented schematically in Fig. 1. A metallic needle of inner radius R i and outer radius R o is held at a constant voltage V . The needle is brought face to face up close to a planar grounded electrode located at a distance H ′ . A liquid of density ρ, viscosity µ, electrical conductivity K, and permittivity ε is injected through the needle at a constant flow rate Q. The flow is fully developed inside the capillary, so that there is a parabolic Hagen-Poiseuille velocity profile upstream at a distance L n from the needle's exit. The ambient medium is a perfect dielectric gas of permittivity ε o equal to that of the vacuum. The gas-liquid surface tension is γ. Because of the small values of the Bond number and gas/liquid density and viscosity ratios, the effects of gravity and the gas dynamics can be neglected.

The equations are made dimensionless with the inner radius R i , the liquid density ρ, the surface tension γ, and the applied voltage V , which yields the characteristic time and velocity scales t c = (ρR 3 i /γ) 1/2 and v c = R i /t c , respectively. The dimensionless axisymmetric incompressible Navier-Stokes equations for the liquid velocity v(r, z; t) and pressure p(r, z; t) fields are (ru) r + rw z = 0, (2)

u t + uu r + wu z = -p r + C[u rr + (u/r) r + u zz ] + χf r e , (3)
w t + uw r + ww z = -p z + C[w rr + w r /r + w zz ] + χf z e , ( 4
)
where t is the time, r/z is the radial/axial coordinate, u/w is the radial/axial velocity

component, χ = ε o V 2 /(R i γ)
is the electric Bond number, and C = µ(ργR i ) -1/2 is the Ohnesorge number. The subscripts t, r, and z here and henceforth denote the partial derivatives with respect to the corresponding variables. The terms f r e / f z e stand for the radial/axial volumetric electric force,

  f r e f z e   = ρ e   E r E z   - (E r + E z ) 2 2   β r β z   , ( 5
)
where ρ e is the volumetric charge density, E r /E z is the radial/axial electric field, and β = ε/ε o is the relative permittivity.

In addition to the Navier-Stokes equations for the liquid flow, the electric potentials ϕ i and ϕ o in both the inner (liquid) and outer (gas) domains must be calculated. They obey the Poisson and the Laplace equations

ϕ i zz + ϕ i rr + ϕ i r /r = ρ e /β, (6)
ϕ o zz + ϕ o rr + ϕ o r /r = 0, (7) respectively.
Taking into account the kinematic compatibility and equilibrium of tangential and normal stresses at the interface r = F (z), one gets the following equations:

F z w -u = 0, (8)
C (1 -F 2 z )(w r + u z ) + 2F z (u r -w z ) (1 + F 2 z ) 1/2 = χσ e E t , (9)
p + F F zz -1 -F 2 z F (1 + F 2 z ) 3/2 - 2C[u r -F z (w r + u z ) + F 2 z w z ] 1 + F 2 z = χ 2 [ (E no ) 2 -β(E ni ) 2 + (β -1)(E t ) 2 ] , ( 10
)
where σ e is the surface charge density, E t is the tangential component of the electric field, and E ni and E no are the normal components of the electrical field on the inner and outer sides of the interface, respectively. The right-hand sides of Eqs. ( 9) and ( 10) are the Maxwell stresses resulting from the accumulation of free electric charges at the interface and the jump of permittivity across this surface.

The electrical field and the surface charge density can be calculated from the inner ϕ i and outer ϕ o electric potentials, and the free surface shape F :

E ni = F z ϕ i z -ϕ i r √ 1 + F 2 z , ( 11
)
E no = F z ϕ o z -ϕ o r √ 1 + F 2 z , ( 12
)
E t = -F z ϕ o r + ϕ o z √ 1 + F 2 z = -F z ϕ i r + ϕ i z √ 1 + F 2 z , ( 13
)
σ e = E no -βE ni , ( 14
)
Note that the continuity of the electric potential across the interface, ϕ i = ϕ o , has been taken into account in Eq. ( 13). The interfacial equations are completed by imposing the surface charge conservation at the free surface r = F (z) [Eq. ( 1)] [20]:

F z u + w 1 + F 2 z σ ez - σ e 1 + F 2 z [ u r + F 2 z w z -F z (u z + w r ) ] = αE ni , ( 15
)
where

α = β [ρR 3 i K 2 /(γε 2 )]
1/2 is the dimensionless electrical conductivity.

We shall assume that all characteristic times of hydrodynamic nature are long as compared to the electrical relaxation time (Saville's leaky dielectric model Saville [21]). Therefore, the volumetric charge density ρ e is zero, and the liquid conductivity K is constant over the entire liquid domain. Also, the liquid electrical permittivity is assumed to be homogeneous, and thus the volumetric dielectrophoretic force proportional to the permittivity gradient is zero. In this case, the electrical effects on the liquid dynamics are exclusively associated with the Maxwell stresses at the interface.

As mentioned above, a Hagen-Poiseuille velocity profile is assumed at the entrance of the liquid domain z = 0:

u = 0, w = 2v e (1 -r 2 ), (16)
where

v e = Q/(πR 2 i v c
) is the dimensionless mean velocity. At the needle wall, the electric potential is fixed and the liquid does not slip, i.e.,

ϕ i = ϕ e = 1 and u = w = 0. ( 17
)
We impose the standard regularity conditions

ϕ i r = u = w r = 0 (18)
on the symmetry axis, and the outflow conditions

u z = w z = F z = (σ e ) z = 0 ( 19
)
at the right-hand end z e = (H + L n )/R i of the computational domain (the red rectangle in Fig. 1).

For the configuration considered in this work, Gañán-Calvo et al. [22] found the following analytical solution for the far-field electric potential:

ϕ 1 (r ′ , z ′ ) = -K v log(4H ′ /R i ) log { [r ′2 + (1 -z ′ ) 2 ] 1/2 + (1 -z ′ ) [r ′2 + (1 + z ′ ) 2 ] 1/2 + (1 + z ′ ) } , ( 20
)
where r ′ and z ′ are cylindrical coordinates with origin at the intersection between the symmetry axis and the grounded planar electrode (see Fig. 1), and K v is a dimensionless constant which depends on the ratio H ′ /R i . We take advantage of this solution to limit the calculations to the meniscus and the jet, leaving aside the charged spray. Specifically, we assume that the electric potential matches the analytical solution ϕ 1 at the boundary r = r e (r e ≡ R e /R i ). A logarithmic drop of voltage

ϕ 2 = 1 -[1 -ϕ 1 (r e , z ′ e )] log r/ log r e , z ′ e ≡ (H ′ + L n )/R i , ( 21
)
is applied at the boundary z = 0 and 1 < r < r e . Finally, the condition

ϕ z = (ϕ 1 ) z ( 22
)
is considered for the right-hand end z = z e of both the liquid and gas computational domains.

III. NUMERICAL PROCEDURES

We solved the problem formulated in the previous section with two numerical schemes: a tracking interface technique proposed in this work, and an adaptation of the Volume-of-Fluid (VoF) method developed by Popinet [17,18] and extended by López-Herrera et al. [19] to study EHD phenomena.

A. The proposed method

Since we search for steady solutions, we drop any time dependence in the Navier-Stokes equations ( 2)-( 4). As mentioned in Sec. II, we assume that the characteristic electric time is much shorter than any hydrodynamic time, and hence all the electrical charges are confined to the free surface (ρ e = 0). In this case, the electric conductivity K can be regarded as constant. Because the liquid electrical permittivity is also assumed to be homogeneous, the electric forces vanish in the bulk domain, and enter into the problem through the Maxwell stresses in the free surface boundary conditions.

A tracking interface technique was implemented to capture the steady liquid-gas interface position, r = F (z). We solved the nonlinear system of equations formed by the bulk equations ( 2)-( 7), the interface equations ( 8)-( 15), and the rest of the boundary conditions ( 16)- (17). In this way, one obtains the spatial dependence of the meniscus-jet shape F , the velocity and pressures fields (u, w, p) in the liquid domain, the surface charge density σ e ,

and the inner and outer potential fields, ϕ i and ϕ o . The current intensity I is the sum of the contributions due to the bulk conduction (I b ) and surface convection (I s ), calculated at any axial position z along the cone-jet as

I b (z) = 2π ∫ F (z) 0 αE iz (r, z)rdr, I s (z) = 2πσ e (z)v s (z), (23)
where v s is the dimensionless liquid velocity at the interface. All these quantities are determined as functions of the five dimensionless control parameters C, v e , β, χ, and α.

A boundary-fitted coordinate system [23] was used to calculate the solution. Both the liquid and the gas domains were mapped onto the fixed rectangular domain (0 ≤ η (l,g) ≤ 1, 0 ≤ ξ ≤ z e ) by means of the coordinate transformations (η l = r/F, ξ = z) and (η g = (r -F )/(r e -F ), ξ = z) for the liquid and gas domains, respectively. The liquid and gas domains were discretized using n η l and n η g Chebychev spectral collocation points [24] along the η l and η g axes, respectively, while the axial direction ξ was discretized using n ξ uniformly distributed points. We used fourth-order central finite differences to get the axial derivatives except close to the ξ-boundaries, where second-order upwind finite differences were applied.

All the simulations presented here were done with n η l = 11, n η g = 20, and n ξ = 160.

The equations in the discretized domain yielded a system of (4n η l + n η g ) × n ξ + 2n ξ nonlinear equations which was solved iteratively using a trust-region method [START_REF] Nocedal | Numerical Optimization[END_REF] implemented in the Matlab subroutine Fsolve. The initial guess for the iterations was obtained as follows. For given values of C, β, χ, and α, we first set a sufficiently high liquid velocity v e . Because a quasi-cylindrical jet is expected in this case, the equations were solved with u = 0, w = 2v e (1 -r 2 ), F = 1, and ϕ i = ϕ e = 1 as the initial guess. Once the algorithm converged to a solution, the control parameter space was explored using a continuation technique. The algorithm showed itself to be robust and computationally efficient. It converged to the sought solution for all the cases considered in less than 1 hour running on an Intel Xeon(R) X5670@2.93Ghz processor.

B. The VoF method

We also performed simulations using the open source fluid simulation engine Gerris [START_REF] Popinet | The Gerris flow solver[END_REF].

Gerris is a VoF method originally conceived to deal with the time-dependent incompressible

Euler equations [17]. It has since been improved and extended considerably to suit problems arising in such diverse fields as multiphase [18], geostrophic [START_REF] Popinet | [END_REF], granular [28], and EHD [19] flows. The VoF method treats immiscible fluids as a single one with spatially varying properties through the interface. Therefore, the balance equation ( 1) for the surface charge density can not be explicitly considered. To circumvent this difficulty, the VoF method solves the evolution equation (A.8) (see Appendix) for the volumetric charge density ρ e with K =const. It simulates a transient regime which eventually reaches the steady conejet state of electrospray (provided that this state is stable for the selected values of the governing parameters). We have verified that the steady states reached in the simulations were consistent with the leaky-dielectric model, i.e., all the charges were confined in a very thin layer next to the interface. In this regard, one can assert that the numerical method converged to the steady solution of the leaky-dielectric model.

We set as initial condition a planar interface located inside the capillary tube. As time proceeds, the meniscus protrudes and a jet emerges by the action of electric forces. The simulation is continued until a stable meniscus-jet is obtained. The grid is adapted at every time-step according to the charge distribution and the interface curvature. The minimum and maximum cell sizes at the last stage of the simulation were h min = 3.91 × 10 -3 and h max = 0.03125, respectively. More details of the method are given elsewhere [19]. The algorithm converged to the sought solution in about 72 hours running on an Intel Xeon(R) X5670@2.93Ghz processor.

IV. EXPERIMENTAL METHOD

A. Experimental setup

The experimental setup is shown schematically in Fig. 2. The liquid was injected at a constant flow rate Q by a stepper motor through a capillary (A) of inner (outer) radius

R i = 100 (R o = 110
) µm located in front of a metallic plate at a distance H ′ = 1 mm. The plate had an orifice (B) of 350 µm in diameter located in front of the capillary. The plate covered the upper face of a metallic cubic cell (C). We used a high-precision orientation system (D) and a translation stage (E) to ensure the correct alignment of these elements, and to set the capillary-to-orifice distance H ′ . An electric potential V was applied to the end of the feeding capillary through a DC high voltage power supply (bertan 205B-10R) (F), and the cubic cell was used as ground electrode. A liquid meniscus was formed in the open air and stretched by the action of the electric field. A microjet tapered from the meniscus tip and moved vertically towards the plate orifice. A prescribed negative gauge pressure (about 10 mbar) was applied in the cubic cell by using a suction pump to produce an air stream coflowing with the jet. Both the liquid jet and the coaxial air stream crossed the plate orifice, which prevented liquid accumulation on the metallic plate. In this way, all the electric charges were collected in the cubic cell. Because the full scale of typical electric current measurements was in the nanoampere range, special care was taken for electrical shielding and grounding. The electric current I transported by the liquid jet was measured using a picoammeter (keithley model 6485) (G) connected to the cell. It must be noted that the mechanical effects associated with the gas stream were negligible (the air speed was smaller than about 10% of that of the liquid at the orifice, and the density and pressure drops in the gas stream were smaller than about 2%). In addition, air ionization effects were less probable than in quiescent air because charges in the gas that could cause ionization cascades were flushed by the gas stream.

Digital images of 1280 × 960 pixels were acquired using a CCD camera (AVT Stringay To check that the fluid configuration was axisymmetric, we also acquired images of the liquid meniscus using an auxiliary CCD camera (not shown in Fig. 2) with an optical axis perpendicular to that of the main camera. All these elements were mounted on an optical table with a pneumatic anti-vibration isolation system (M) to damp the vibrations coming from the building.

Experiments were conducted with 1-octanol (Panreac 99% ps) (ρ = 827 kg/m 3 , µ = 0.0081 kg/ms, γ = 0.0266 N/m, K = 9 × 10 -7 S/m, β = 10). The surface tension was determined using the Theoretical Image Fitting Analysis (TIFA) method [29]. The electrical conductivity was measured by applying a voltage difference between the ends of a cylindrical borosilicate capillary full of the working liquid, and then measuring the resulting electric current. The viscosity, density, and electrical permittivity were taken from the manufacturer's specifications. By way of illustration, Fig. 3 shows the image of a jet produced by the experimental setup described above with Q = 0.5 ml/h and V = 1.6 kV.

B. Image processing technique

We shall here briefly describe the main aspects of the image processing technique used to precisely determine the free surface position. Details of the procedure are given elsewhere [30,31]. The contours of the free surface in the image were detected using a two-stage procedure. In the first stage, a set of pixels probably corresponding to the contour being sought was extracted using Otsu's method [32]. The accuracy of Otsu's method is limited to the pixel size. In the second stage, the accuracy of the detected contours was improved to the sub-pixel level by analyzing the gray intensity profile along the direction normal to the contour. Fitting the sigmoid (Boltzmann) function [33] to the gray intensity values allowed us to obtain a continuous function in the transient region of the gray profile. The contour point was found by applying the local thresholding criterion. The number of resulting contour points was approximately equal to the number of pixels of the image along the vertical direction. Once the left x l and right x r contours had been detected, they were rotated to their vertical position (the rotation angle was less than 2 • in all the cases analyzed), and the symmetry axis was found. Small non-axisymmetric perturbations of the free-surface were ignored by considering the axisymmetric contour F = (x r -x l )/2. This quantity was compared with the numerical results.

V. RESULTS AND DISCUSSION

Numerical simulations of the problem formulated in Sec. II were conducted with the tracking interface technique proposed in this work (Sec. III A), and a VoF method [17,18] (Sec. III B). The following set of geometrical parameters was considered in all the simulations: r e = 6.67, L n /R i = 1.77, H/R i = 9.4, and H ′ /R i = 40 (for which K v = 0.56).

In addition, the needle width was set to zero (R i = R o ) in the proposed method, while the value R o /R i = 1.22 were selected in the VoF technique to improve the algorithm convergence.

Figure 4 shows a comparison between the two numerical methods for the set of parameters C = 0.18, β = 10, α = 15.3 , χ = 10, and v e = 0.0182, for which a stable meniscus-jet configuration was reached. The two methods predict almost the same free surface shape, and show the existence of a recirculation cell inside the liquid meniscus. This recirculation pattern is caused by the action of the shear Maxwell stresses at the free surface, which drag along that surface an amount of liquid that can not be exhausted by the thin emitted jet.

This flow structure is characteristic of low-viscosity liquids electrosprayed at small flow rates [34], and has also been observed in coflowing systems for the low Q-range [35,36].

We also checked our tracking numerical method against experiment for two different cases. Figure 5 shows the results for the free surface position F (z) of a 1-octanol meniscusjet produced by applying a voltage V = 1600 V while injecting liquid at Q = 0.5 and 1 ml/h. These flow rates are close to the minimum value Q min ∼ γε o (ρK) -1 = 1.14 ml/h for which the cone-jet mode is stable [4,8,37]. Indeed, we were not able to produce a stable 1-octanol jet for Q < 0.5 ml/h. The corresponding dimensionless parameters are C = 0.182, β = 10, α = 15.3, χ = 10, and v e = 0.00913 and 0.0182 for Q = 0.5 and 1 ml/h, respectively. There is remarkable agreement between the experimental and numerical contours. The average deviation between them was less than 0.6 µm in both cases. The simulations captured perfectly the reduction of the jet radius when the flow rate was reduced by half. It is worth emphasizing that the agreement between the numerical and the experimental results for these limiting configurations shows the validity of the charge relaxation hypothesis [4,7,8,37] over the entire parameter space. FIG. 5. Free surface position F (z) of a 1-octanol meniscus-jet produced with V = 1600 V, and Q = 0.5 ml/h (circles and solid line) and Q = 1 ml/h (triangles and dashed line). The symbols and lines correspond to the experiments and proposed numerical method, respectively. Here, F and z stand for the corresponding dimensional quantities, and the z axis is the upward vertical axis with origin at the bottom of the experimental image.

Optical diffraction renders the present fluid configuration almost opaque. Therefore, it is quite difficult to get information on the flow pattern inside the meniscus by standard experimental techniques, such as micro-PIV or dye-injection. On the contrary, numerical simulations allow not only the visualization of the recirculation cell taking place for small Q, but also the quantification of its size as Q decreases. Figure 6 illustrates this capability of the proposed numerical method. One of the most relevant quantities in electrospray is the electric current I transported by the jet. A proof of the consistency of our simulations is the global conservation of electric charge, which implies that the total electric current I = I b + I s should remain constant along the meniscus-jet. Figure 7 shows the axial distribution of the current intensities associated with the bulk conduction and the surface convection, as well as the total value, for the experimental case with Q = 0.5 ml/h. The values were made dimensionless with the characteristic intensity I 0 = (ε 0 γ 2 /ρ) 1/2 [4,8]. As has been noted in theoretical works [4,8,12], bulk conduction is the dominant mechanism of charge transport in the meniscus, while surface convection is dominant in the jet region. The neck can be seen as the region where the transition between these two mechanisms takes place. The dimensionless intensity I/I 0 measured in the experiments was 1.4, reasonably close to the numerical prediction 1.2 and the scaling law I/I 0 = 2.6(Q/Q 0 ) 1/2 = 1.7 [8].

To summarize, a robust and efficient simulation method has been proposed to simulate steady EHD phenomena. Calculations were greatly simplified in this method by assuming that the electric charges are accumulated within an infinitely thin layer at the interface, and then electric forces result from the Maxwell stresses at that surface exclusively. The validity of the discretization scheme has been shown by comparing the results with those obtained with a much more sophisticated VoF method. We have also shown that the charge relaxation simplification is correct in the electrospray cone-jet mode. The free surface position of the meniscus-jet was measured very precisely in experiments with 1-octanol. The comparison between the numerical predictions for the liquid contour and the experimental results confirmed the validity of the proposed method. The numerical method also provided Here, n (k) , r (k) , z (k) , and ω (k) are the concentration, production rate by chemical reaction, valence, and mobility of the k th species, respectively, while k B , e, and T are the Boltzmann constant, the electron charge, and fluid temperature, respectively.

The charge density ρ e can be calculated as

ρ e = ∑ k ez (k) n (k) .
(A.6)

A conservation equation for this variable can be readily derived by adding the equations The production rates r (k) depend on the concentrations of the species and the chemical kinetics of the reactions occurring in the system. Note that these reactions can be extremely complex and, in some occasions, unknown. In fact, the rates r (k) depend on each other because the production/destruction of some species implies the destruction/production of others. Therefore, care in the modelling of the chemical reactions must be taken to prevent overdetermination. Nevertheless, three scenarios are simple enough to be easily modeled.

The first is the unipolar injection, where a single species generated at an upstream electrode migrates through the fluid with zero production rate. The second one is the evolution of strong electrolytes. In this case, neutral species are fully dissociated and chemical reactions are absent (r (k) = 0). The last case is the evolution of a binary z-z electrolyte analyzed by Saville [21]. Three species are diluted in the fluid: the neutral, the ionic, and the counterionic. In the forward reaction, the neutral species dissociate producing the same amount of ions and counterions, while in the recombination reaction the ionic and counter-ionic species disintegrate at the same rate to produce the neutral one. The source/sink term ∑ k ez (k) r (k) in Eq. (A.7) vanishes in the three cases described above. In the leaky-dielectric model [21], all characteristic times of hydrodynamic nature in the bulk are much longer than the electrical relaxation time. Therefore, the volumetric charge density ρ e is assumed to be zero, and electrical charges are only allowed to exist at the interface. Therefore, Eq. (A.8) reduces to ∇ • (KE) = 0 . (A.9)

It must be noted that Saville [21] keeps the bulk electric forces in Eq. (A.4) although, strictly speaking, they should be removed. Our VoF method solves Eqs. (A.1)-(A.4) and (A.9) with K =const.

  FIG. 1. (Color online) Sketch of the problem's formulation. The rectangle (red online) denotes the limits of the computational domain.

FFIG. 2 .

 2 FIG. 2. Schematic diagram of the experimental setup.

FIG. 3 . 1 -

 31 FIG. 3. 1-octanol jet produced by the experimental setup described in the text. The left (right) image was acquired with the auxiliary (working) CCD camera.

FIG. 4 .

 4 FIG. 4. Free surface position and streamlines calculated with the tracking interface technique [solid line and (a)] and VoF method [symbols and (b)] for C = 0.18, β = 10 , α = 15.3 , χ = 10, and v e = 0.02.

FIG. 6 .

 6 FIG. 6. Streamlines plotted on the corresponding experimental images for a 1-octanol meniscus-jet produced with V = 1600 V, and Q = 0.5 ml/h (a) and Q = 1 ml/h (b).

(A. 5 )z

 5 for the different species:(k) ω (k) ∇ 2 n (k) + ∑ k ez (k) r (k) . (A.7)In Eq. (A.7), migration due to thermal diffusion (the second term of the right-side of the equation) can be neglected as compared to electromigration, provided that ℓeE o /(k B T ) ≫ 1 (ℓ and E o are the characteristic length and electric field, respectively). In electrospray,ℓ ∼ (γε 2 o /ρK 2 ) 1/3 and E o ∼ (γ/ε o ℓ) 1/2 [4], which yields ℓeE o /(k B T ) ≃ 8 × 10 3 for 1-octanol. Therefore, thermal diffusion is negligible for the configurations considered in the present work.

  Taking into account these considerations, Eq. (A.7) reduces to the expressionDρ e Dt = ∇ • (-KE) , K = ∑ k ω (k) e 2 z (k)2 n (k) . (A.8)

ACKNOWLEDGMENTS

Partial support from the Ministry of Science and Education, Junta de Extremadura, and Junta de Andalucía (Spain) through Grants Nos. DPI2010-21103, GR10047, and P08-TEP-04128, respectively, is gratefully acknowledged.

correct values of the emitted electric current even for configurations very close to the cone-jet stability limit.

Appendix: Equations of the VoF method

In this appendix, we derive the equations solved by the VoF method from the Taylor-Melcher leaky-dielectric model [21] . The starting point is the Maxwell electrostatic equations

the Navier-Stokes equations

and the equations for the conservation of neutral and charged species .5)