H Manzanarez

J P Mericq

P Guenoun

D Bouyer
email: denis.bouyer@umontpellier.fr

Modeling the interplay between Solvent Evaporation and Phase Separation Dynamics during Membrane Preparation

Keywords: Phase-field simulation, Cahn-Hilliard model, mobility, membranes, evaporation process

Introduction

Phase separation processes are widely used in industry for manufacturing various types of products, from usual metals to polymer solutions. Preparation of polymeric membrane is one of these applications of great interest [START_REF] Van De Witte | Phase separation processes in polymer solutions in relation to membrane formation[END_REF][START_REF] Ulbricht | Advanced functional polymer membranes[END_REF]. In the last four decades, research has been dedicated to the polymeric membrane formation mechanisms in order to better control the final membrane morphology [START_REF] Kesting | Synthetic polymeric membranes: a structural perspective[END_REF][START_REF] Mulder | Basic Principles of Membrane Technology[END_REF]. Starting from a homogeneous polymeric solution composed of a polymer dissolved in a good solvent, a thermodynamic demixing process is induced by a temperature change (Temperature Induced Phase Separation or TIPS process) [START_REF] Caneba | Polymer membrane formation through the thermal-inversion process. 1. experimental study of membrane structure formation[END_REF][START_REF] Lloyd | Microporous membrane formation via thermally induced phase separation. i. solid-liquid phase separation[END_REF][START_REF] Lloyd | Microporous membrane formation via thermally-induced phase separation. ii. liquid-liquid phase separation[END_REF][START_REF] Shang | Preparation and membrane performance of poly(ethyleneco-vinyl alcohol) hollow fiber membrane via thermally induced phase separation[END_REF] or by the intrusion of a non-solvent of the polymer (Non-solvent Phase Separation or NIPS process) [START_REF] Reuvers | Formation of membranes by means of immersion precipitation: Part i. a model to describe mass transfer during immersion precipitation[END_REF][START_REF] Kim | Formation of polyurethane membranes by immersion precipitation. ii. morphology formation[END_REF]. Starting from a ternary system composed of a polymer, a good solvent and a small amount of a non-solvent of the polymer, a faster evaporation rate of the solvent comparing to that of non-solvent can also induce the phase inversion (Dry-casting process) [START_REF] Shojaie Saeed | Dense polymer film and membrane formation via the dry-cast process part i. model development[END_REF]. During the demixing process, two phases will be created: a polymer-rich phase mainly composed of polymer and a polymerlean phase mainly composed of solvent (and/or non-solvent depending on the process). The polymer-rich phase form the membrane matrix after extraction of the polymer-lean phase which will form membrane pores.

One of most important challenges in membrane manufacturing concerns the control of the final morphology that will strongly affect the membrane performances towards the targeted applications. For instance, asymmetric structures characterized by a pore structure that gradually changes from very large pores to very fine pores at the membrane surface [START_REF] Ulbricht | Advanced functional polymer membranes[END_REF], will be targeted for pressure driven membrane processes (filtration in water treatment applications for example). The upper selective layer, responsible for membrane selectivity, should be as thin as possible, while the pore size strongly increases beneath this selective layer to maximize the filtration flux through the membrane. On the contrary, symmetric membranes with uniform structures through the entire membrane thickness could be interesting for applications such as dialysis and electrodialysis, but also microfiltration [START_REF] Ladewig | Fundamentals of Membrane Processes[END_REF].

Controlling the whole membrane structure is therefore the key point in membrane preparation, but it still remains a goal hard to achieve since the membrane formation mechanisms are quite complex and particularly difficult to simulate, and hence to predict. Phase separation can be described using the equations of Cahn and Hilliard [START_REF] Cahn | Free energy of a nonuniform system. i. interfacial free energy[END_REF] for polymeric systems, where the free energy of mixing of the polymeric system is derived from Flory Huggins theory [START_REF] Flory | Principles of Polymer Chemistry[END_REF] and the mobility term has to be described using a specific equation.

In a recent paper, [START_REF] Manzanarez | Modeling phase inversion using cahn-hilliard equations -influence of the mobility on the pattern formation dynamics[END_REF] [START_REF] Manzanarez | Modeling phase inversion using cahn-hilliard equations -influence of the mobility on the pattern formation dynamics[END_REF] investigated the influence of this mobility term on the phase separation dynamics for a closed binary polymeric system [START_REF] Barton | Dynamics of spinodal decomposition in polymer solutions near a glass transition[END_REF][START_REF] Mino | Three-dimensional phase-phield simulations of membrane porous structure formation by thermally induced phase separation in polymer solutions[END_REF].

However, additional features have to be described to simulate the membrane formation since the phase separation is coupled with transfer phenomena occurring at membrane interfaces. Indeed, mass exchanges often occurs between the membrane and the external environment simultaneously with phase separation. For instance, solvent extraction and non-solvent intake occur during NIPS process [START_REF] Zhou | Phase field simulations of early stage structure formation during immersion precipitation of polymeric membranes in 2d and 3d[END_REF][START_REF] Tree | A 34 multi-fluid model for microstructure formation in polymer membranes[END_REF], while solvent and non-solvent evaporation will be involved during dry-casting process. Recently, we also exhibited how the phase inversion performed by LCST-TIPS process for water soluble polymer systems was coupled with solvent evaporation [START_REF] Bouyer | Modeling the membrane formation of novel pva membranes for predicting the composition path and their final morphology[END_REF]. Focusing on the modeling of TIPS process, a wide literature exists and various types of models have been developed during the last 30 years (phase field methods, dissipative particle dynamics méthods, Coarse grain simulation, Monte-Carlo simulation. . . ). Caneba and Song [START_REF] Caneba | Polymer membrane formation through the thermal-inversion process. 1. experimental study of membrane structure formation[END_REF] were one of the first to develop a 1D phase field model to simulate the TIPS process. They used Cahn-Hilliard equations for spinodal decomposition, Flory-Huggins for thermodynamics and Vrentas models for the description of the mobility terms. Later, Barton and Mc Hugh [START_REF] Barton | Dynamics of spinodal decomposition in polymer solutions near a glass transition[END_REF][START_REF] Barton | Modeling the dynamics of membrane structure formation in quenched polymer solutions[END_REF] added a temperature gradient due to heat transport to simulate the droplets growth during demixing process, in 1D geometry yet.

Using phase field methods, the impact of a temperature gradient was also investigated by Lee and coworkers [START_REF] Lee | Morphology development and characterization of the phase-separated structure resulting from the thermal-induced phase separation phenomenon in polymer solutions under a temperature gradient[END_REF] in 2D geometry and by Chan [START_REF] Chan | Effect of concentration gradient on the thermal-induced phase separation phenomenon in polymer solutions[END_REF] in 1D geometry to better understand the formation of anisotropic morphologies by TIPS process. Later the Cahn-Hilliard equations were solved in 3D geometry for modeling TIPS process [START_REF] Mino | Three-dimensional phase-phield simulations of membrane porous structure formation by thermally induced phase separation in polymer solutions[END_REF]. Using different modeling method, He et al. [START_REF] He | Dissipative particle dynamics simulation on the membrane formation of polymer-diluent system via thermally induced phase separation[END_REF], Tang and coworkers simulated the TIPS process by Dissipative Particle Dynamics simulation (DPD) (2013,2015,2016) [START_REF] Tang | Three-dimensional analysis of membrane formation via thermally induced phase separation by dissipative particle dynamics simulation[END_REF][START_REF] Tang | Investigation on the membrane formation process of polymer-diluent system via thermally induced phase separation accompanied with mass transfer across the interface: Dissipative particle dynamics simulation and its experimental verification[END_REF][START_REF] Tang | Multiscale simulation on the membrane formation pro-cess via thermally induced phase separation accompanied with heat transfer[END_REF]. Even more recently, Tang and coworkers [START_REF] Tang | Formation of polyethersulfone membranes via nonsolvent induced phase separation process from dissipative particle dynamics simulations[END_REF] used DPD and Coarse Grain methods to simulate the coupling between phase separation and mass transfer when the UCST-TIPS process is conducted by immerging a hot polymer solution into a cold water bath. In the latter case, mass exchanges are expected to be very rapid since they occur in liquid phase. However, the coupling between mass exchanges by solvent evaporation and phase separation was less investigated.

Mino et al. [START_REF] Mino | Three-dimensional phase-phield simulations of membrane porous structure formation by thermally induced phase separation in polymer solutions[END_REF] only considered in their simulations an initial concentration gradient that could be due to an initial solvent evaporation but their simulations did not involve the direct coupling between both phenomena. However, for LCST-TIPS process the coupling between solvent evaporation and phase separation dynamics is crucial since both dynamics are slowed and concomitant [START_REF] Bouyer | Modeling the membrane formation of novel pva membranes for predicting the composition path and their final morphology[END_REF][START_REF] Hanafia | Fabrication of novel porous membrane from biobased water-soluble polymer (hydroxypropylcellulose)[END_REF]. Whatever the process aforementioned, it is of prime importance to elucidate how the mass exchanges affect the phase separation dynamics, and hence the final membrane morphology. Surprisingly, to the best of our knowledge, few theoretical studies have considered this coupling [START_REF] Tree | Mass-transfer driven spinodal decomposition in a ternary polymer solution[END_REF].

This paper focuses on the coupling between the phase separation induced by TIPS process and mass transfer phenomena for a simplified binary polymer/solvent system. More specifically, the solvent evaporation occurring at the upper membrane interface will be simulated and its interplay with the phase separation dynamics will be investigated.

In the first part of the paper, the coupling between phase separation and solvent evaporation will be simulated and discussed in a horizontal 2D plan, and then in the second part, the simulations will be performed in a 2D vertical cross-section of the membrane in order to investigate the possible formation of concentration gradients.

Theory

The diffusion equation that describes the phase separation is a modified time-dependent Ginzburg-Landau theory for a conserved order parameter.

∂φ (r, t) ∂t = ∇ • Λ (φ, r) ∇ δF δφ (r) + σ (φ, r) (1)
where φ (r, t) is the volume fraction of the polymer, Λ (φ, r) the mobility term and F is the free energy functional of the system can thus be expressed as [START_REF] Cahn | Free energy of a nonuniform system. i. interfacial free energy[END_REF] :

F [φ (r, t)] = dr f (φ) + C 2 |∇φ| 2 (2)
where f (φ) is the free energy of mixing per lattice site for polymer solutions described by the Flory-Huggins theory [START_REF] Flory | Principles of Polymer Chemistry[END_REF]:

f (φ) = k B T v 0 φ N lnφ + (1 -φ) ln (1 -φ) + χφ (1 -φ) (3)
were the degree of polymerization N = 150, v 0 is the volume of the monomer and χ is the interaction parameter. This parameter is supposed to be here an inverse function of temperature. All quenches studied here are made at a constant χ and are consequently isothermal. C|∇φ| 2 /2 is the gradient energy contribution which describes the cost of an interface between the two phases resulting from the phase separation. The gradient parameter C follow here the Debye approximation [START_REF] Debye | Angular dissymmetry of the critical opalescence in liquid mixtures[END_REF]:

C (φ) ≈ k B T v 0 χ 3 R 2 g (4)
R g is the radius of gyration of the polymer (R 2 g ≈ a 2 N/6 were a is the Kuhn length and N the degree of polymerization). The term C can be related to the polymer chain length and the surface tension. Thus, from a purely numerical point of view, this parameter imposes a space discretization (typically a mesh must be lower than C) which makes it possible to determine the final size of the simulated field (L = 250 C). Finally, the source term σ (φ, r) is added to the continuity equation in order to add evaporative transfers at the interface between the system and the outside.

In a recent paper, we investigated the influence of the mobility term on the phase separation dynamics [START_REF] Manzanarez | Modeling phase inversion using cahn-hilliard equations -influence of the mobility on the pattern formation dynamics[END_REF] using Fourier transforms and Minkowski descriptors. Fast, Slow and Vrentas model were compared for various conditions of initial quenching with a polymeric system described by the Flory-Huggins theory. For binary systems, those models were shown to give somewhat different results in terms of growth law. Experimental data found in literature however shown to mostly validate the Fast model [START_REF] Kramer | Interdiffusion and marker movements in concentrated polymer-polymer diffusion couples[END_REF][START_REF] Sillescu | Relation of interdiffusion and self-diffusion in polymer mixtures[END_REF]. This model was consequently chosen here and writes as:

Λ (φ, r) = φ (1 -φ) [φD s + (1 -φ) N D p ] (5)
In this study, constant values are taken to estimate the tracer diffusion of the solvent D s and the polymer D p = D s /N .

To describe the solvent evaporation during phase inversion process, a solvent flux was added at the upper interface, from the polymeric system to the external environment:

J evap (φ, r upper ) = k∆g (φ) (6)
where ∆g (φ) = a i s (T i ) -a ∞ s (T ∞ ) is the difference in activities between the upper polymer solution at the interface and the solvent in the external environment [START_REF] De | Solvent evaporation of spin cast films: "crust" effects[END_REF][START_REF] Dušková-Smrčková | Solvent activity changes and phase separation during crosslinking of coating films[END_REF][START_REF] Ovejero | Thermodynamic modeling and simulation of styrene-butadiene rubbers (sbr) solvent equilibrium staged processes[END_REF][START_REF] Bouyer | Validation of mass-transfer model for vips process using in situ measurements performed by near-infrared spectroscopy[END_REF].The solvent activity at the upper interface was calculated using the mathematical description of the free energy of mixing, herein the Flory-Huggins theory:

a s = exp 1 k B T ∂∆G m ∂n s (7)
In equation ( 6), k represents the mass transfer coefficient that mainly depends on the hydrodynamic conditions in the air above the interface and can be calculated using semi-empirical correlations. The evaporation regime are described by the convective Biot number:

Bi = k D 0 l(t = 0) (8)
For the simulations of this work, a simplified approach was used and four different values of k were chosen for testing different regime of evaporation: Bi = 0.01, 0.1 and 0.5. The solvent activity in the external environment was assumed to be null.

It must be noted that the simplified model of evaporation used here is based on the assumption that diffusion of species is the limiting factor of evaporation [START_REF] Tsay | Mass transfer dynamics of the evaporation step in membrane formation by phase inversion[END_REF]. This model can be complemented by taking into account the gelification of the polymer in the evaporating crust [START_REF] De | Solvent evaporation of spin cast films: "crust" effects[END_REF][START_REF] Ozawa | Diffusion process during drying to cause the skin formation in polymer solutions[END_REF][START_REF] Rabani | A phase-field model for the evaporation of thin film mixtures[END_REF],

a phenomenon which is ignored here. Moreover, another approach, alternative to diffusion-limited models, was recently proposed for dealing with non-equilibrium situations [START_REF] Cummings | Modeling solvent evaporation during thin film formation in phase separat-39 ing polymer mixtures[END_REF].

Methods

The non-dimensionnal system of partial differential equations was numerically solved in two dimensions using finite element software: COMSOL Multiphysics 5.4 with the scaling parameters

L 0 = l x = 1[µm], t 0 = L 2 0 /D 0
and the Biot number (equation 8).

A structured moving mesh was used with 64x64 nodes and it was refined until no change in the numerical result was obtained. A variable time step was used to improve the numerical resolution. For the results in the YZ plane, the resolution of mass transfers process induces a displacement of the upper domain boundary that smoothly displaces the mesh nodes at the surface and inside the bulk of the domain.

The images resulting from phase separation-evaporation simulation were analyzed using Fast Fourier Transform (FFT) performed by Image J (NIST), as well as by calculating the Minkowski descriptors for estimating the topological indicators such as the volume fraction, the connectivity, and the Euler characteristics. More details on the procedure to obtain such data, can be found in our previous paper [START_REF] Manzanarez | Modeling phase inversion using cahn-hilliard equations -influence of the mobility on the pattern formation dynamics[END_REF]. For all pictures, the FFT returned a reciprocal space image exhibiting a ring, confirming the existence of a characteristic size of phase separation L m (t). For determining this size with precision, these rings were radially averaged to provide a I(q) curve-also called structure factor-where I is proportional to the square of the Fourier transform of concentration correlations and q is the wavevector of reciprocal space.

These structure factors were compared at different distances (see below) of the evaporating surface for assessing the phase separation homogeneity.

The thresholding method transforms the grayscale image to black/white images by sweeping all thresholds. Then, three topological descriptors are extracted from those images: F is the ratio of the area occupied by the white pixels divided by the total pixel number, U is the interface density obtained by counting the pixel number of black/white neighbors and Ec is the Euler characteristic which can be assimilated in this work to the difference between the number of black non-connected domains and the number of white nonconnected domains. All images were normalized in size before analysis to ensure a clear comparison from one image to another

Results

Two sets of results are presented herein: in a first part, simulations were carried out in a 2D XY plane (horizontal plane) as reported in Figure 1. In a second part, simulations in a 2D vertical YZ plane are presented to simulate the phase separation in the membrane cross-section. Starting from this phase diagram (Figure 2) and this quenching point, in a closed system the phase separation would lead to the formation of a continuous polymer-lean phase and a disperse polymer-rich phase. At equilibrium, each phase would be expected to tend to its equilibrium concentration: φ a = 1.4.10 -4 for the polymer-lean phase and φ b = 0.386 for the polymerrich phase (cf. Figure 2). In absence of evaporation, the global polymer volume concentration φ is constant, equal to 0.08 and the volume fractions of each phase keep constant values, equal to φ rich = 0.21 and φ lean = 0.79, respectively, as calculated by the lever rule. Now, when considering the coupling between phase separation and solvent evaporation, the global polymer concentration is expected to gradually increase in the system in such a way that the system follows a composition path that will ultimately reach the right of the phase diagram, following the dotted line at constant χ = 0.7 (Figure 2).

For the first set of simulations in the 2D XY plane, the nondimensional mass transfer coefficient k was fixed at 0.1. In Figures 3 a-j In a previous paper [START_REF] Manzanarez | Modeling phase inversion using cahn-hilliard equations -influence of the mobility on the pattern formation dynamics[END_REF], we analyzed the patterns using both Fourier transform and Minkowski descriptors. The latter method was used in this work to analyze more deeply the influence of the coupling between solvent evaporation and phase inversion dynamics: the patterns were binarized using a chosen threshold, and then the binarized images were analyzed with three Minkowski descriptors: volume fraction, connectivity, Euler characteristics.

The use of Minkowski descriptors requires performing a prior binarization of the patterns. The choice of the threshold is not trivial since it could significantly affects the curves interpretation. In our previous paper, the binarization threshold was chosen equal to the initial polymer concentration in such a way that during the phase separation, the regions characterized by higher polymer concentration than the initial polymer concentration were represented in white color, while the regions with lower concentrations than the initial polymer concentration were represented in black color. In this way, it was easy to catch the formation of the polymer-lean and polymer-rich phases as soon as the demixing process started [START_REF] Manzanarez | Modeling phase inversion using cahn-hilliard equations -influence of the mobility on the pattern formation dynamics[END_REF]. Although the problem is different in this work since the solvent evaporation induces a continuous increase of the polymer concentration, the same thresholding procedure was F is shown to tend towards F ∞ = 0.21 with regards to the level rule. After t = 50, at the end of the simulation, Figure 4.a exhibits that the system is almost at equilibrium in terms of volume fraction of each phase (F is close to 0.21). On the contrary, when the evaporation is coupled to the phase separation, the system is continually in non-equilibrium state and the dynamics of 1 Actually, as shown by volume fraction temporal evolution, the average domains concentration is very quickly different from the initial value : this latter choice is then still suitable for defining separating domains. However, as the interface is rather steep between phases, the error is believed to be negligible since 0.08 is far from actual concentrations of evaporating phases phase separation is different. Due to continuous solvent loss, the descriptor F gradually increases until reaching unity, corresponding to pure polymer, around t = 20. Note that no slope change was observed once the composition path passes through the binodal curve, corresponding to a polymer concentration close to 0.386. At this time (around t = 20), the system is composed of a continuous phase characterized by a high concentration in polymer with a covered area fraction of 1.

As reported in a previous paper [START_REF] Manzanarez | Modeling phase inversion using cahn-hilliard equations -influence of the mobility on the pattern formation dynamics[END_REF], U is the interface density (boundary When coupling solvent evaporation and phase separation (solid curves), the mean polymer concentration continuously increases, in agreement with the displacement of the composition path in the phase diagram toward higher polymer concentrations. For the lean phase, the curve is slightly above the curve corresponding to the closed system until around t = 12, which suggests that the relaxation dynamics due to the phase separation is too slow to compensate the continuous solvent loss. In other words, the system is continuously forced to be in non-equilibrium state because of solvent evaporation.

length
The same conclusion can be drawn considering the polymer-rich phase:

dl dt = J evap (9)
J evap is determined by the flux of solvent, expressed by the equation [START_REF] Lloyd | Microporous membrane formation via thermally-induced phase separation. ii. liquid-liquid phase separation[END_REF]. Simulations were performed assuming three different quenching points at three initial polymer concentrations (φ init = 0.08, 0.14 and 0.20, respectively), and four nondimensional values of the mass transfer coefficient k that correspond at the Biot number (Bi = 0.01, Bi = 0.1 and Bi = 0.5).

At early times, a surface directed phase separation is evidenced due to breaking of symmetry caused by evaporation. This is evidenced as a surface composition wave limited in extension to a few wavelengths (see top of figures 7 at t=1 where horizontal domains are evidenced) and whose wavelength value is close to the bulk phase separation characteristic length.

These surface directed phase separation patterns have been evidenced

both experimentally due to preferential wetting constraints [START_REF] Guenoun | Dynamics of wetting and phase separation[END_REF] or theoretically because of thermal gradients on one side [START_REF] Ball | Spinodal decomposition and pattern formation near surfaces[END_REF] or, more recently, because of solvent replacement in a ternary solution [START_REF] Tree | Mass-transfer driven spinodal decomposition in a ternary polymer solution[END_REF]. The evaporating surface can be viewed as a wetting constraint breaking isotropy, contrary to the neutral sides and bottom surfaces where isotropy is maintained by periodic boundary conditions and Dirichlet condition respectively.

At later times, when the surface directed phase separation wave has disappeared, a dense layer is formed in the vicinity of the upper interface.

In Figure 8 are reported the patterns obtained for the three initial quenching points, without evaporation and with solvent evaporation for the interme- diate evaporation rate (Bi = 0.1) at late times (t >= 6). With evaporation, rich phase when φ init ≈ 0.20. In all cases the dense phase (or skin) on top is clearly inhomogeneous along z and forms a quasi-planar interface with the phase-separating region below for φ init = 0.08 whereas the skin is continuously linked with the phase-separating region in the other cases.

A qualitative visual observation of the snapshots suggests that the formation of the upper dense polymer layer weakly affects the dynamics of phase separation beneath the dense layer: the patterns are very similar with and without the solvent evaporation. The simulations launched with φ init = 0.14 that show interconnected structures, exhibit slight differences near the dense layer region. At φ init = 0.08 and φ init = 0.20 the differences are even weaker with and without solvent evaporation (see below Figure 12 for a quantitative analysis).

To refine the previous observations, we reported in Figure 9 the patterns obtained at t = 1 of the twelve simulations (3 concentrations, 4 evaporation rates). The first, second and third lines correspond to quenching point at initial concentrations (φ init = 0.08), φ init = 0.14 and φ init = 0.20, respectively. The evaporation rate taken into account for those simulations were Bi = 0.01, Bi = 0.1 and Bi = 0.5 a, for columns 1, 2 and 3 respectively.

Not surprisingly, increasing the mass transfer coefficient leads to faster evaporation, to a faster decrease of domain height. A dense layer is formed which seems thicker when increasing not only the initial polymer concentration, but also the mass transfer coefficient Bi.

The dense layer is fairly easy to define when the continuous phase is the polymer-lean phase (case for φ init = 0.08) since an interface is clearly formed. For the other cases (φ init = 0.08 and φ init = 0.20), the dense layer thickness should be carefully evaluated (Figure 9). We decided to define the lower boundary of the dense layer as the Z-coordinate where the polymer concentration reaches the equilibrium concentration (0.386): actually a concentration larger than 0.386 is very quickly reached as soon as the thick layer is detectable. The dense layer thickness was reported in Figure 10.b for the aforementioned initial conditions of simulation (φ init = 0.08, 0.14 and 0.20

and Bi = 0.1 at t=2).

As visible in Figure 9, the thickness of the dense layer increases more rapidly at higher initial polymer concentration (φ init = 0.20). This is due to the fact that at higher initial polymer concentration, the binodal of the dense phase is reached earlier compared to lower initial polymer concentration when solvent evaporation occurs. Since the lower interface of the dense layer was not perfectly flat in the simulated patterns due to the phase separation and the presence of droplets or interconnected structures, its thickness was estimated using an average along the Z-axis. The difficulty of identifying the dense layer is exemplified by the wavy curve shape in Figure 10.b, especially when droplets of polymer-lean phase are dispersed in a continuous phase of polymer-rich phase (φ init = 0.20).

Below the dense (skin) layer, a visual observation of the polymer-rich and polymer-lean phases indicates that no gradient exist in polymer concentration in both phases. In other words, the patterns presented above suggest that the relaxation dynamics below the skin layer are sufficiently fast to maintain the polymer-rich and polymer-lean phases close to the equilibrium values despite the solvent evaporation at the upper interface.

In order to better evidence this absence of gradients, we report in Figure • to focus on the polymer-lean continuous phase (images d), e) and f) in Figure 10), a threshold was fixed at φ s2 = φ a + 0.05 ≈ 0.05 in such a way that when the concentration exceeded 0.05, its value was fixed to a white color. On the contrary, for concentrations between 0 and 0.05, the color scale is respected and exhibits the absence of color gradient, i.e. of concentration gradient in the z-direction in the polymer-lean phase.

The Z-average concentrations over the continuous phase were calculated along a vertical Z line at each Y-coordinate (Figure 11). the plots confirms that below the dense layer and at t = 2, the concentration in the polymer- lean phase for φ init = 0 is very close to 1.10 -4 and the concentration in the polymer-rich phase for φ init = 0.37 is also very close to 0.386. No concentration gradient was observed from the bottom to the dense layer interface, whatever the case.

To complete the quantitative image analysis, 2D Fourier transform (FFT)

were calculated for the images to check to what extent the evaporation did affect the phase separation dynamics. For the three pictures at the latest times of separation and the three different initial concentrations, the structure factors were calculated in two rectangular windows. Each window is rectangular of width 446 pixels and height 152 pixels and one is located close to the interface and the other down close to the lower border (see Figure 12 for details). Structure factors are shown to be similar close and far to the interface, proving homogeneity. In particular, for each pair of curves at equal time and concentration, a peak at small wavevectors is evidenced de-spite the large typical distance between domains which pushes the peak to the y-axis. This peak is however distinguishable and representative of the distance between domains, which is similar at the top and bottom.

To summarize the previous simulation results obtained in the YZ place, we exhibited that the continuous solvent loss at the top surface due to evaporation induces the formation of a skin layer (note here that no change of dynamics such as gelation or glass transition was assumed to take place in this layer [START_REF] Ovejero | Thermodynamic modeling and simulation of styrene-butadiene rubbers (sbr) solvent equilibrium staged processes[END_REF]), which suggests that mass transfer localized at the upper interface is faster than the potential inflow of solvent from deeper layers by molecular diffusion. The system relaxes to minimize its free energy in such a way that an equilibrium is reached between the lower part of the dense layer and the adjacent bottom separating layer composed of polymer-lean phase or polymer-rich phase. The thickness of the gradient zone was shown to increase during time whatever the Biot number, i.e. the relative evaporation rate driven by the air flow conditions. In this way, beneath this dense layer, the phase separation in the bulk solution was shown not to be affected by the mass transfer occurring at the upper system interface.

Conclusion

Herein, we developed a model that coupled the demixing process and the solvent evaporation during the membrane formation by TIPS process.

Simulations have been performed in a 2D geometry in XY plane (membrane surface) and YZ plane (cross-section). The simulations in the X-Y plane clearly predicts the existence of an evaporation regime where an initially minority phase rich in polymer will be turned into a majority phase (percolation inversion). This is confirmed in Y-Z simulations where a "polymer skin layer" appears after some time on top of an initially dilute phase. Surprisingly and interestingly enough, our results predict that the skin layer is a gradient zone in concentration between a value that tends to one at the top surface and the equilibrium value of the polymer-rich phase which is maintained throughout the evaporation process. This result demonstrates that beneath the skin layer, the phase separation was not affected by the solvent loss at the top surface and stays homogeneous through the entire bulk volume, which was not expected a priori. In this context, the simulation results presented in this work allow a better understanding of the interplay between the solvent evaporation and demixing process, especially for predicting the skin layer thickness, which depend on the evaporation rate. The formation of this skin layer (whose porosity is often facilitated by porogen additives in industry) is of crucial importance for the membrane preparation since it controls the membrane selectivity. The thickness of this layer also plays an important role in the membrane permeability since the main resistance to mass transfer is localized in it. Furthermore, the model provides an insight on the interplay between the solvent evaporation and the demixing process deeper in membrane. Indeed, this work demonstrated that the dynamics of phase separation below the skin layer was not affected by the solvent evaporation, meaning that the pore size within the membrane bulk is not affected by the solvent evaporation. This suggests that the global membrane porosity (the void ratio) would not be affected by the solvent evaporation. On the theoretical side, our results are foreseen to be extended to a tridimensional geometry for coupling X-Y and Y-Z processes.

Figure 1 :Figure 2 :

 12 Figure 1: Schematic representation of the evaporation process in the XY plane

Figure 3 :

 3 Figure 3: Time evolution of the patterns obtained in closed system (a-e) and with a coupling between phase separation and solvent evaporation (f-j). φ is the polymer concentration.

Figure 4 : 1 .

 41 Figure 4: Time evolution of the Minkowski descriptors in closed system (dash line) and for the coupling between phase inversion and solvent evaporation (blue line) for k = 0.1. a. represents the variation of the covered area F occupied by the rich phase with time, b. the interface density U and c. the Euler characteristic E c .

Figure 4 .

 4 Figure 4.a represents the variation of the covered area occupied by the rich phase with time. Without evaporation the volume fraction of each phase is not expected to change during demixing process, so the Minkowski descriptor

  in 2D) and Ec is the Euler characteristic, useful for analyzing the connectivity of domains. Without evaporation, Figure 4.b exhibits that the boundary length U continuously decreases. The small droplets formed at initial stage are expected to grow and they coalesce with other droplets or disappear due to ripening effect, thus decreasing the total interface length. In presence of solvent evaporation, the curve of the interface density U shows the same trend as the curve without evaporation during the first time steps, and then the interface density U is shown to decrease steeper until zero due to the disappearance of solvent droplets. Without solvent evaporation, the Euler characteristic E c en Figure 4.c was shown to sharply decrease at the beginning of the phase separation because of the creation of numerous dispersed droplets, and then E c slightly increases due to the reduction of the droplets number. The Euler characteristics E c was shown to keep negative values in a closed system, indicating the absence of percolation inversion: the polymer-rich phase is always the dispersed phase into a continuous polymer-lean phase. On the contrary, the

Figure 4 .Figure 5 :

 45 Figure 4.c exhibits that a percolation inversion is detected using this descrip-

ZFigure 6 :

 6 Figure 6: Schematic representation of the evaporation process in YZ plane

Figure 7 :

 7 Figure 7: Early times phase separation patterns under evaporation of solvent by the top most surface for Bi = 0.1. A surface composition wave is established more rapidly than the bulk PS. This wave eventually transforms in the top skin layer which is characteristic of late times evolution.

Figure 8 :

 8 Figure 8: Patterns obtained at increasing time steps for different values of φ init with Bi = 0.1.

Figure 9 :

 9 Figure 9: Patterns obtained at different initial polymer concentrations φ init and different values of the mass transfer coefficients Bi for t = 1.

Figure 10 :

 10 Figure 10: Bulk polymer concentration following the color code in the concentrated phase (a), b) ,c)) and the lean phase (d),e), f)) respectively for different quenching concentrations: left φ init = 0.08, middle φ init = 0.14 and right φ init = 0.20 at t=2s. Thresholds are chosen to render the lean (resp. the concentrated phase) white to reveal the other one (see text)

Figure 11 :

 11 Figure 11: Z-average concentrations over the continuous phase for the three different concentrations after the thresholding procedure of Figure 10.

Figure 12 :

 12 Figure 12: (left) Superposition of structure factors at two different distances from the dense evaporating layer and for the three concentrations at the latest time studied. Oscillations are non-physical and merely due to the pixellisation inhomogeneity whereas the dip at a wavevector of ca. 251 is due to the sudden drop in analyzed pixels number above the maximum inscribed circle radius. (right) Position of the chosen analysis rectangles.

the formation of a dense layer and the decrease of the total height were clearly 380 observed, whatever the initial polymer concentration φ init .

381

When the initial polymer concentration (φ init ) is lower than 0.14, droplets 382 of polymer-rich phase are dispersed into a continuous polymer-lean phase.

383

Besides, symmetrical interconnected phases are observed when φ init ≈ 0.14 384 and droplets of polymer-lean phase are dispersed into a continuous polymer-


M.-O Bristeau
email: marie-odile.bristeau@inria.fr

B Di Martino
email: dimartinob@univ-corse.fr

A Mangeney
email: anne.mangeney@gmail.com

J Sainte-Marie
email: jacques.sainte-marie@inria.fr

F Souille
email: fabien.souille@edf.fr

Some quasi-analytical solutions for propagative waves in free surface Euler equations

This note describes some quasi-analytical solutions for wave propagation in free surface Euler equations and linearized Euler equations. The obtained solutions vary from a sinusoidal form to a form with singularities. They allow a numerical validation of the free-surface Euler codes.

Résumé

Solutions quasi-analytiques d'ondes propagatives dans les équations d'Euler à surface libre Cette note décrit des solutions quasi-analytiques correspondant à la propagation d'ondes dans les équations d'Euler et d'Euler linéarisées à surface libre. Les solutions obtenues varient d'une forme sinusoïdale à une forme présentant des singularités. Elles permettent de valider numériquement les codes de simulation des équations d'Euler à surface libre.

Introduction

The water wave problem described by the Euler equations with a free surface has been widely studied in the literature, see e.g. [START_REF] Craik | The origins of water wave theory[END_REF][START_REF] Houghton | Nonlinear shallow fluid flow over an isolated ridge[END_REF][START_REF] Lannes | The water waves problem (book, 321 p)[END_REF][START_REF] Longuet-Higgins | A Theory of the Origin of Microseisms[END_REF][START_REF] Stoker | The formation of breakers and bores the theory of nonlinear wave propagation in shallow water and open channels[END_REF]. This paper proposes some quasi-analytical solution of these equations that allow, for example, to validate the efficiency of the numerical tools. These analytical solutions exhibit singularities of the free surface when the wave amplitude increases. As far as the authors know, it is the first analytical solutions having such a behavior corresponding to an existence result given by W. Strauss, see [START_REF] Constantin | Exact steady periodic water waves with vorticity[END_REF][START_REF] Strauss | Bound on the Slope of Steady Water Waves with Favorable Vorticity[END_REF] and references therein. Some other explicit solutions have been presented in the literature, for example by Boulanger et al. [START_REF] Boulanger | Analytical solutions for the free surface hydrostatic Euler equations[END_REF] and Daboussy et al. [START_REF] Daboussy | On explicite solutions of the free-surface Euler equations in the presence of gravity[END_REF] for the steady state. Following the methodology of Constantin and Strauss [START_REF] Constantin | Exact steady periodic water waves with vorticity[END_REF], Kalimeris [START_REF] Kalimeris | Analytical approximation and numerical simulations for periodic travelling water waves[END_REF] proposes an asymptotic expansion of the Euler system reducing the problem resolution to a cascade of ODEs. On the one hand, the result of Kalimeris is not reduced to flows without vorticity, on the other hand the proposed solutions -also exhibiting singularities of the free surface -are not analytical because obtained through an iterative numerical process.

Solutions presented in Section 2 are solutions of the Euler linearized system up to a negligible term. The proposed solutions are irrotationnal and are compared in Section 3 to the solutions of Airy and third order Stokes waves. In Section 4 this result is extended to the nonlinear Euler system through an additional pressure term on the free surface. Same type of quasi-analytical solutions are proposed in Section 5 for the stationnary waves.

We consider the Euler system and the linearized Euler system over a flat bottom for x ∈ R and 0 ≤ z ≤ h(t, x) given respectively by (1)-( 3) and ( 4)- [START_REF] Daboussy | On explicite solutions of the free-surface Euler equations in the presence of gravity[END_REF], where u(t, x, z), w(t, x, z) are the two components of the velocity in the (x, z) domain, h(t, x) is the water depth, p(t, x, z) is the pressure and ρ 0 is the density assumed to be constant:

∂u ∂x + ∂w ∂z = 0, (1)
∂u ∂t + u ∂u ∂x + w ∂u ∂z + 1 ρ 0 ∂p ∂x = 0, (2) ∂w ∂t + u ∂w ∂x + w ∂w ∂z + 1 ρ 0 ∂p ∂z = -g, (3)
Euler system

∂u ∂x + ∂w ∂z = 0, (4)
∂u ∂t + u 0 ∂u ∂x + 1 ρ 0 ∂p ∂x = 0, (5)
∂w ∂t + u 0 ∂w ∂x + 1 ρ 0 ∂p ∂z = -g, (6)

Linearized Euler system

These systems are completed by initial conditions (u(0, x, z) = u 0 (x, z), w(0, x, z) = w 0 (x, z), p(0, x, z) = p 0 (x, z)), a dynamic boundary condition at the free surface

p s = p(t, x, h(t, x)) = p a (t, x), (7)
a kinematic boundary condition at the free surface and a non-penetration condition at the bottom

∂h ∂t + u s ∂h ∂x -w s = 0, w b = 0. ( 8
)
where the subscript s (resp. b) denotes the considered quantity at the free surface (resp. at the bottom).

Remark 1 For the sake of simplicity we have used the same notations for the solution of the Euler and linearized Euler system but it is clear that they correspond to different solutions.

The linearized Euler system (4)-( 6) is obtained by assuming that the velocity components u and w are such that u = u 0 + O(ε), w = O(ε) with ε 1 and u 0 = cst. Around the solution (u 0 , 0, p a + ρ 0 g(h -z)), the solution (u, w, p) of ( 4)-( 6) yields a remainder term of order O(ε 2 ) in ( 1)- [START_REF] Constantin | Exact steady periodic water waves with vorticity[END_REF]. It is important to notice that in most cases the linearized Euler system does not admit any energy balance. However, simple computations show that when the quantity

e s,b = |U s | 2 2 ∂η ∂t + |U s | 2 2 u 0 ∂η ∂x ,
where U = (u, w)

and |f | 2 = |(f 1 , f 2 )| 2 = f 2 1 + f 2 2
, can be written under the conservative form e s,b = ∂α s,b ∂x , with α s,b = α(h, u 0 , U s ), the linearized Euler system (4)-( 6) completed with ( 7) and ( 8) admits an energy balance of the form

∂ ∂t η z b (E + p a ) dz + ∂ ∂x η z b u 0 E + u(p + gz) dz + α s,b = h ∂p a ∂t , (9)
with

E defined by E = u 2 +w 2 2 + gz.
The water depth h(t, x) does not appear directly in systems ( 1)-( 3) and ( 4)-( 6), it can be obtained by integrating equation ( 1) from z = 0 to z = h(t, x)

∂h ∂t + ∂ ∂x h(t,x) 0 u(t, x, ξ)dξ = 0. ( 10
)
We consider the linearized Euler system (4)-( 6) completed by the boundary conditions ( 7) and [START_REF] Houghton | Nonlinear shallow fluid flow over an isolated ridge[END_REF]. A large part of the results are based on the properties of the LambertW functions (LW). The LambertW functions is the main branch of a set of functions corresponding to the inverse relation of the function f (z) = ze z where z is any complex number [START_REF] Corless | On the LambertW function[END_REF]. LW(x) is the unique real solution of LW(x)e LW(x) = x.

For x ≥ -1/e and for x = 0 and x = -1/e we have LW (x) = LW(x)

x(1 + LW(x)) . Moreover, LW (x) ∼ x→0 x.
Then the following proposition holds.

Proposition 2.1 Let (b, u 0 , h 0 , k, a, c, ω) ∈ R 2 × R 5 + such that |a| < 1/e and h 0 k/c > LW(-|a|) and f : x → - c k LW a cos x c + b , (12)
where LW represents the LambertW function. Let p a (t) be any given function.

The functions h, u, w and p defined by

h(t, x) = h 0 + f (kx -(ω + ku 0 )t), (13)
u(t, x, z) = u 0 - ωa k e k c (z-h 0 ) cos kx -(ω + ku 0 )t c + b , (14)
w(t, x, z) = - c k ∂u ∂x , (15)
p(t, x, z) = p a (t) + ρ 0 g(h 0 -z) - ρ 0 gca k e k c (z-h 0 ) cos kx -(ω + ku 0 )t c + b , (16)
are quasi-analytical solutions of the linearized Euler system (4)-( 6) completed by the boundary conditions [START_REF] Dingemans | Water wave propagation over uneven bottoms[END_REF] and (8) if and only if the following relation holds

ω k = gc k . ( 17
)
More precisely, equations (4)-( 6), dynamic pressure condition [START_REF] Dingemans | Water wave propagation over uneven bottoms[END_REF] and the kinematic free surface condition in equation (8) are verified exactly while. For Eqs (8) and (10), we get:

w b = O e -h 0 k c and ∂h ∂t + ∂ ∂x h(t,x) 0 u(t, x, ξ)dξ = O e -h 0 k c . ( 18
)

Sketch of proof

To make the proof more readable, we set u 0 = 0, b = 0 and we take Y = a cos ((kx -ωt)/c).

From (4) and using the definition of u given in ( 14) we obtain easily the expression of w. From ( 14)-( 16), by simple derivation of u and p, equation ( 5) gives

- ω 2 a k + ga e k c (z-h 0 ) sin(Y ) = 0.
This equation is verified for all Y only if the relation (17) holds. This relation is similar to the dispersion relation in the Airy theory when kh 0 1. Expression of p is obtained by integrating (6) from z to h. It is easy to verify a posteriori that ( 6) is verified by taking the derivative of w and p appearing in this equation. From (16), we observe that

p(t, x, h(t, x)) = p a (t) -gf (kx -ωt)) - gc k e k c f (kx-ωt)) Y, = p a (t) + gc k LW(Y ) - gc k e -LW(Y ) Y = p a (t).
The main difficulty is to verify the surface evolution equation [START_REF] Lannes | The water waves problem (book, 321 p)[END_REF]. From ( 13) we obtain

∂h ∂t = ωa k sin kx -ωt c LW(Y ) Y (1 + LW(Y )) . ( 19
)
With the expression of u given in (14) we deduce using [START_REF] Longuet-Higgins | A Theory of the Origin of Microseisms[END_REF] that

h(t,x) 0 u(t, x, ξ)dξ = - ω k Y h(t,x) 0 e k c (ξ-h 0 ) dξ, = - ωc k 2 Y e -k c (h(t,x)-h 0 ) -e k c h 0 , = ωc k 2 Y e -k c h 0 -LW(Y ) .
Then, using the same expression of the derivative of the LW function, we have

∂ ∂x h(t,x) 0 u(t, x, ξ)dξ = - ωa k sin kx -ωt c + b e -k c h 0 + LW(Y ) Y (1 + LW(Y ))
.

(20) We deduce that equation [START_REF] Lannes | The water waves problem (book, 321 p)[END_REF] is verified up to a term ωa k sin kx-ωt c + b e -k c h 0 . In the same way, vertical velocity at the bottom is given by w b = -ωa c e -k c h 0 . The solutions proposed in prop. 2.1 for the linearized Euler system are not exactly analytical solutions in the sense that additional terms in O(e -h 0 k c ) appear. But when h 0 k 1, considering e.g. h 0 = 100 m, k = 0.2 m -1 and c = 1 gives e -h 0 k c ≈ 10 -9 1. These solutions satisfy the energy balance (9) with

α s,b = - gω 4k c 2 LW a cos kx-(ω+ku 0 )t c + b 2 k 2 cos kx-(ω+ku 0 )t c + b 2 .

Comparison with the Airy and Stokes wave theories

The form of the free surface of our quasi-analytical solution depends of the parameter a. When a is small, surface elevation is close to a sinusoidal function since we have seen that LW (x) ∼ x→0

x. When |a| is near to 1/e, the top of the wave is narrower than its bottom. This property is presented in Figure 1 with parameters (c = 1, k = 1, t = 0 and b = 0). For |a| = 1/e, the function f is not differentiable in (2m + 1)π, ∀m ∈ Z. The non-differentiable character of the solution was studied from a theoretical point of view by Strauss [START_REF] Constantin | Exact steady periodic water waves with vorticity[END_REF]. Here we give an explicit expression of this solution.

In the literature, some analytical solutions of free surface simplified models have been proposed [START_REF] Craik | The origins of water wave theory[END_REF]. The most known are the Airy wave and the thirdorder Stokes wave, and we propose here some numerical comparisons.

By setting θ = (k 0 x -ω 0 t + b 0 ), in the Airy wave theory, surface elevation is given by h(t, x) = h 0 + a 0 cos(θ), (21) and in the third-order Stokes wave on deep water, surface elevation is given by [START_REF] Dingemans | Water wave propagation over uneven bottoms[END_REF]

h(t, x) = h 0 +a 0 1 - 1 16 (k 0 a 0 ) 2 cos(θ) + 1 2 k 0 a 0 cos(2θ) + 3 8 (k 0 a 0 ) 2 cos(3θ) .
(22) To make a comparison with water depth given by equation ( 13), we have to set a = a 0 k 0 (using again LW (x) ∼ x→0 x), k = ck 0 , b = b 0 and ω = cω 0 . We have plotted in Figure 2 a comparison between the free surface obtained by Airy theory, third order Stokes theory and equation [START_REF] Strauss | Bound on the Slope of Steady Water Waves with Favorable Vorticity[END_REF] for the given set of parameters (a = 1/e, k = 1, c = 1, t = 0, b = 0). With this value of a we observe numerically the maximum of differences between the three waves. Of course, for small value of a, solution of ( 13) can be very close to the two other waves.

Fig. 1. Free surface given by the function f : x → -LW(a cos(x)) for three values of the parameter a : 0.10 (dot line), 0.30 (dash-dot line), 1/e (dash line), and comparison with function (-1/e) cos(x) (solid line). Fig. 2. Free surface given by Airy theory (solid line), third order Stokes theory (dash-dot line) and equation (13) (dot line) for the given set of parameters :

(a = 1/e, k = 1, c = 1, t = 0, b = 0).

Propagating waves for the Euler system

We consider the Euler system (1)-( 3) completed by the boundary conditions [START_REF] Dingemans | Water wave propagation over uneven bottoms[END_REF] and [START_REF] Houghton | Nonlinear shallow fluid flow over an isolated ridge[END_REF]. The previous solution can be extended by considering a small pressure term at the free surface and the following proposition holds. Proposition 4.1 Under the same conditions of the proposition 2.1 the functions h, u, w and p defined by

h(t, x) = h 0 + f (kx -(ω + ku 0 )t), u(t, x, z) = u 0 - ωa c e k c (z-h 0 ) cos kx -(ω + ku 0 )t c + b , w(t, x, z) = - c k ∂u ∂x , p(t, x, z) = g(h 0 -z) + ga 2 c 2k 1 -e 2k c (z-h 0 ) - gca k e k c (z-h 0 ) cos kx -(ω + ku 0 )t c + b ,
are quasi-analytical solutions of the Euler system (1)-(3) completed by the boundary conditions [START_REF] Dingemans | Water wave propagation over uneven bottoms[END_REF] and (8) iff the relation (17) holds. At the free surface, the pressure p is such that p s = p a (t, x) with

p a (t, x) = ga 2 c 2k 1 -e 2k c (h(t,x)-h 0 ) . (23)
More precisely, equations (4)-( 6), dynamic pressure condition [START_REF] Dingemans | Water wave propagation over uneven bottoms[END_REF] and the kinematic free surface condition in equation ( 8) are verified exactly. For Eqs. ( 8) and (10), we get

w b = O e -h 0 k c and ∂h ∂t + ∂ ∂x h(t,x) 0 u(t, x, ξ)dξ = O e -h 0 k c . (24)
Pressure at the surface is equal to zero when h(t, x) = h 0 . If h(t, x) -h 0 is small, this pressure is of the order a 3 . The proof is similar to the one given in prop. 2.1 but is more tedious due to the nonlinearity.

Remark 2 The solutions proposed in props. 2.1 and 4.1 are irrotational. In [START_REF] Amick | Bounds for water waves[END_REF], C. J. Amick proved that for any irrotational wave, the angle of inclination of the free surface with respect to the horizontal must be less than 31.15˚. In our analytical solutions, the angle of inclination is less or equal to 45˚, but we have an additional source term (23) that can justify this inclination.

Standing waves

Now we consider the situation of standing waves that occur when two progressive waves of same amplitude travel in opposite direction. The results depicted in this paragraph are based on the following remark: for small values of the parameter a, one has the Taylor expansion

LW a cos kx -ωt c = a cos kx -ωt c -a 2 cos 2 kx -ωt c + O(a 3 ).
Proposition 5.1 Under the same conditions as in proposition 2.1 the functions h, u, w and p defined by

h(t, x) = h 0 + f (kx -ωt) + f (kx + ωt), (25)
u(t, x, z) = ω c e k c (z-h(t,x)) f (kx -ωt) -f (kx + ωt) , (26)
w(t, x, z) = - c k ∂u ∂x , (27)
p(t, x, z) = p a (t) + g(h -z) + h(t,x) z ∂w ∂t dz, (28)
are quasi-analytical solutions of the linearized Euler system (4)-( 6) completed by the boundary conditions (7) and (8) iff the relation (17) holds. At the free surface, the pressure p is such that p s = p a (t, x) with

p a (t, x) = ga 2 c 2k 1 -e 2k c (h(t,x)-h 0 ) . (29)
More precisely, equations (4) and (6), dynamic pressure condition (7) and the kinematic free surface condition in equation (8) are verified exactly. For Eqs. (5), ( 8) and (10), we get

∂u ∂t + u 0 ∂u ∂x + 1 ρ 0 ∂p ∂x = O a 3 , (30)
w b = O ae -h 0 k c and ∂h ∂t + ∂ ∂x h(t,x) 0 u(t, x, ξ)dξ = O ae -h 0 k c . (31)
Since a Taylor expansion of Eq. (25

) gives h(t, x) = h 0 -2ac k cos kx c cos ωt c + O a 2
k , the proposed solution corresponds, up to terms in O a 2 k , to a standing wave.

The proofs rely on simple but tedious computations similar to those performed in the proof of prop. k , the proposed solution corresponds, up to terms in O a k , to a standing wave.

The proofs rely on simple but tedious computations similar to those performed in the proof of prop. 2.1.

2 . 1 . 5 . 2

 2152 Proposition Under the same conditions as in proposition 4.1 the functions h, u, w and p defined byh(t, x) = h 0 + f (kx -ωt) + f (kx + ωt),(32)u(t, x, z) = ω c e k c (z-h(t,x)) f (kx -ωt) -f (kx + ωt) , , x, z) = g(h -z) + ga 2 c 2k 1 -e 2k c (z-h 0 ) + g(h -z) + h(t,x) z ∂w ∂t dz,(35)are quasi-analytical solutions of the Euler system (1)-(3) completed by the boundary conditions (7) and (8) iff the relation (17) holds. At the free surface, the pressure p is such that p s = p a (t, x) with p a (t, x) = ga 2 More precisely, equations (1) and (3), dynamic pressure condition (7) and the kinematic free surface condition in equation (8) are verified exactly. For Eqs. (1), (8) and (10), we get x, ξ)dξ = O a 2 e -h 0 k c . (38) Since a Taylor expansion of Eq. (25) gives h(t, x) = h 0 -2ac k cos kx c cos ωt c + O a 2

Acknowledgements

The work presented in this paper was supported in part by the ANR Mimosa (ANR-14-CE01-0012) and ERC contract ERC-CG-2013-PE10-617472 SLIDE-QUAKES. The authors would like to thank Professor E. Stutzmann for fruitful discussions that have helped to realise this work.