Louis-Alexandre Couston

Daniel Lecoanet

Benjamin Favier

Michael Le Bars

The energy flux spectrum of internal waves generated by turbulent convection

. Specifically, we demonstrate that the wave energy flux spectrum scales like k 4 ⊥ f -13/2 for weakly-damped waves (with k ⊥ and f the waves' horizontal wavenumbers and frequencies), and that the total wave energy flux decays with z, the distance from the convective region, like z -13/8 .

Introduction

Internal gravity waves can be found in planetary atmospheres and oceans, radiative envelopes of stars, and planetary fluid cores. Numerous studies have shown that internal gravity waves can drive strong large-scale flows [START_REF] Grisouard | Forcing of oceanic mean flows by dissipating internal tides[END_REF][START_REF] Bordes | Experimental observation of a strong mean flow induced by internal gravity waves[END_REF][START_REF] Van Den Bremer | The wave-induced flow of internal gravity wavepackets with arbitrary aspect ratio[END_REF][START_REF] Couston | Order out of chaos: Slowly reversing mean flows emerge from turbulently generated internal waves[END_REF]) and enhance turbulence and mixing [START_REF] Staquet | Internal gravity waves: From instabilities to turbulence[END_REF][START_REF] Kunze | Internal-wave-driven mixing: Global geography and budgets[END_REF][START_REF] Thorpe | Models of energy loss from internal waves breaking in the ocean[END_REF]. Through their effects on mean flows and properties of the fluid in which they propagate, internal waves can have a significant impact on the long-term evolution of global systems, such as Earth's climate [START_REF] Alexander | Recent developments in gravity-wave effects in climate models and the global distribution of gravity-wave momentum flux from observations and models[END_REF] or the differential rotation of stars [START_REF] Rogers | Internal Gravity Waves Modulate the Apparent Misalignment of Exoplanets Around Hot Stars[END_REF].

The energy flux of internal waves is required to accurately predict the influence of waves on global dynamics, but is difficult to estimate, especially when the generating source is turbulent. Nevertheless, turbulent sources account for a significant fraction of the internal wave energy in the atmosphere, stars and planetary interiors, so they can hardly be neglected. Turbulent motions in the bottom Ekman layer [START_REF] Taylor | Internal gravity waves generated by a turbulent bottom ekman layer[END_REF] or the surface mixed layer of the ocean [START_REF] Munroe | Internal wave energy radiated from a turbulent mixed layer[END_REF], turbulent buoyant plumes in the atmosphere [START_REF] Ansong | Internal gravity waves generated by convective plumes[END_REF], fingering convection [START_REF] Garaud | Double-diffusive convection at low prandtl number[END_REF]) and penetrative convection in stars (Pinçon et al. 2016) are examples of turbulent sources that can generate internal gravity waves.

Past studies of turbulence-driven internal waves have focused on extracting scaling laws for the total wave flux with respect to problem parameters. Information on the energy flux spectrum, i.e. the distribution of energy for different wavenumbers k and frequencies f of internal waves, is often beyond reach or difficult to interpret because of experimental or numerical constraints. For instance, the simulations of [START_REF] Rogers | Internal gravity waves in massive stars: Angular momentum transport[END_REF] were restricted to only 2D because of this difficulty. However, knowledge of the energy flux spectrum is essential to predicting wave effects, since both wave propagation, damping and interaction with the background state depend strongly on k and f . Several theoretical approaches have been put forward to predict the wave spectrum due to turbulent motions. Under the linearized rapid distortion theory (e.g. [START_REF] Carruthers | Velocity fluctuations near an interface between a turbulent region and a stably stratified layer[END_REF], the wave amplitude is proportional to the convective velocity. We instead assume the wave amplitude is linearly proportional to the Reynolds stress (i.e., the square of the convective velocity), which can accurately predict waves in direct numerical simulations (DNS) [START_REF] Lecoanet | Numerical simulations of internal wave generation by convection in water[END_REF]. Two theoretical descriptions of the Reynolds stresses in turbulent convection are the "eddy" and "plume" theories. In the "eddy" theory, the convection is decomposed into eddies which follow Kolmogorov statistics and eddy stresses are used as generating sources of internal waves in the stable region [START_REF] Goldreich | Wave generation by turbulent convection[END_REF][START_REF] Lecoanet | Internal gravity wave excitation by turbulent convection[END_REF]). In the "plume" theory, the convection is seen as an ensemble of uncorrelated plumes impinging at the base of the stable layer (Pinçon et al. 2016). The "plume" theory predicts the wave energy flux spectrum to be Gaussian in k and f ; the "eddy" theory predicts power laws.

Here we validate for the first time the prediction of the wave flux spectrum based on the "eddy" theory in the limit of strong stable stratification [START_REF] Lecoanet | Internal gravity wave excitation by turbulent convection[END_REF] using 3D state-of-the-art DNS of a convective-stably-stratified fluid model. The DNS and analyses are challenging as they require a fully turbulent flow interacting with a strongly stratified layer and high spatial/temporal resolution to construct the spectra.

DNS model

The fluid model we use is similar to the one used in [START_REF] Couston | Dynamics of mixed convective-stably-stratified fluids[END_REF], though here it is 3D (figure 1). We solve the Navier-Stokes equations under the Boussinesq approximation, i.e. the fluid density is ρ = ρ 0 + δρ with ρ 0 the (constant) reference density and δρ ρ 0 are small variations due to temperature fluctuations. We consider constant kinematic viscosity ν, acceleration due to gravity g and thermal diffusivity κ. The governing equations for the velocity u = (u, v, w) and temperature T in a Cartesian (x, y, z) frame of reference (ẑ is the upward-pointing unit vector of the z-axis) read

∂ t u + (u • ∇)u = -∇p + ν∇ 2 u -gδρ/ρ 0 ẑ -u/τ,
(2.1a)

∂ t T + (u • ∇)T = κ∇ 2 T, (2.1b) ∇ • u = 0, (2.1c)
with p the pressure and -u/τ is a damping term used to prevent wave reflection from the top boundary. The density anomaly δρ is related to T through the equation of state

δρ/ρ 0 = -α(T )(T -T i ) = -α s (T -T i )H(T -T i ) + α s S(T -T i )H(T i -T ), (2.2)
with H the Heaviside function; α s > 0 is the thermal expansion coefficient for T > T i , with T i the inversion temperature at which the density anomaly δρ/ρ 0 = 0 is maximum, and S > 0 is the stiffness parameter. We assume periodicity in the horizontal directions and impose free-slip and fixed temperatures conditions on the top (T = T t ) and bottom boundaries (T = T b ), with T b > T i > T t . Then, because equation (2.2) is nonmonotonic, the lower part of the fluid is convectively unstable, the upper part is stably stratified and S indicates how sharp or "stiff" the transition is at the interface [START_REF] Couston | Dynamics of mixed convective-stably-stratified fluids[END_REF].

We use the characteristic length scale H, the thermal diffusion time H 2 /κ, and the temperature difference ∆T = T b -T i in order to define the non-dimensional variables:

(x, ỹ, z) = (x, y, z) H , t = tκ H 2 , ũ = uH κ , T = T -T i ∆T , p = pH 2 κ 2 , ρ = δρ ρ 0 α s ∆T . (2.3) Substituting (2.3) in (2.1)-(2.
2) and dropping tildes we obtain the dimensionless equations

∂ t u + (u • ∇)u = -∇p + P r∇ 2 u -P rRaρẑ -u/τ, (2.4a) ∂ t T + (u • ∇)T = ∇ 2 T, (2.4b) ∇ • u = 0, (2.4c) ρ = -T H(T ) + ST H(-T ),
(2.4d )

with P r = ν/κ the Prandtl number and Ra = α s g∆T H 3 /(νκ) the Rayleigh number. We note L x , L y , L z the dimensionless lengths of the domain in x, y, z directions and T top < 0 the dimensionless temperature at z = L z (T = 1 at z = 0 and the inversion temperature is T = 0 in dimensionless space). The dimensionless z-dependent damping coefficient used to prevent wave reflection from the top boundary is τ -1 (z) = 0.5f N * {tanh[(z -L z + 0.15)/0.05] + 1} with f N * the target buoyancy frequency (cf. next section). We solve (2.4) using the open-source pseudo-spectral code DEDALUS [START_REF] Burns | Dedalus: A Flexible Pseudo-Spectral Framework for Solving PDEs[END_REF]. We are interested in DNS results at statistical steady-state and in the limit of fully turbulent convection (large Ra) and strong stable stratification (large S). For ease of presentation, we want the convection to extend from z = 0 to z = h ≈ 1 at thermal equilibrium, i.e. such that H is approximately the size of the convection zone in dimensional variables and Ra is the effective Rayleigh number of the convection. h is an output of our simulations, but it can be controlled by iteratively adjusting T top (all other parameters fixed). At leading order, h ≈ 1 when T top is chosen such that the measured heat flux Q = wT -T z in the convection (overbar denotes the horizontal average) equals the expected diffusive heat flux of the stable region ∼ -T top /(L z -h) with h = 1. Because Q in the convection zone changes with T top , finding T top is carried out iteratively. We fix all problem parameters except T top , make an initial guess for T top , and run a simulation for a few convective turnover times, initialized with a temperature of T = (1 -z)H(1 -z) + T top (z -1)H(z -1) plus low-amplitude noise. We then update our guess for T top using T top = -Q(L z -h) with h = 1, where Q is the average heat flux at z = 0, and we reiterate until T top doesn't change by more than about 10%. The final T top gives an estimate for the buoyancy frequency f N * = -P rRaST top /(L z -1)/(2π) used in the expression for the damping term τ in (2.4a). The actual buoyancy frequency f N = N/(2π) is computed using N = √ P rRaSQ with Q averaged from z = 0 to z = 1 and over the simulation time ∆t. In all the simulations, we set L x = L y = L z = L with L = 2 or L = 3. Therefore, the ratio between the horizontal and vertical extent of the convective zone is larger than 2, which ensures that confinement effects associated with the horizontal periodic boundary conditions are not dominating the dynamics. We set P r = 1 for simplicity. Input and output parameters of the eight simulations included in this paper are reported in table 1.

The variables of interest include the kinetic energy spectrum and the energy flux spectrum. The kinetic energy per unit surface area and unit (thermal) time reads

K(z) = |u| 2 2 dxdydt L 2 ∆t ≈ fi,k ⊥j | û| 2 2δf δk ⊥ δf δk ⊥ ≡ fi,k ⊥j dKδf δk ⊥ ≈ ∂ 2 K ∂f ∂k ⊥ df dk ⊥ , (2.5)
with hat denoting the Fourier transform in (x, y, t), ∆t the period over which results are integrated in time, k ⊥ and f are the horizontal wavenumber and frequency, δk ⊥ and δf are the unit increments in spectral space, and dK is the discrete version of the kinetic energy spectrum ∂ 2 K ∂f ∂k ⊥ . Similarly, we define the z-dependent vertical energy flux as

F (z) = wp dxdydt L 2 ∆t ≈ fi,k ⊥j R( ŵ p * ) δf δk ⊥ δf δk ⊥ ≡ fi,k ⊥j dF δf δk ⊥ ≈ ∂ 2 F ∂f ∂k ⊥ df dk ⊥ , (2.6)
where * denotes the complex conjugate, dF is the discrete energy flux spectrum and R denotes the real part. Note that throughout this paper, we assume isotropic motions in x, y directions and express variables in spectral space in terms of the average horizontal wavenumber

k ⊥ = (k 2 x + k 2 y ) 1/2 (with k x , k y the wavenumbers in x, y directions). Thus, | û| 2 ≡ | û(k ⊥ )| 2 is the azimuthal average of | û(k x , k y )| 2 and R( ŵ p * ) is the azimuthal average of R[ ŵ(k x , k y )p * (k x , k y )]
based on the modulus of vector (k x , k y ), i.e. averaging over all orientations. We apply a Hann windowing function to all variables before taking the time Fourier transform and we discard the modes f = 0 and k ⊥ = 0 from the spectra.

Theoretical prediction

We now describe a heuristic model which predicts the wave energy flux spectrum at the interface, F i (f, k ⊥ ), between the convective and stably-stratified layers. This is based on the calculation in [START_REF] Lecoanet | Internal gravity wave excitation by turbulent convection[END_REF] (specifically, the "discontinuous" N 2 case), which itself is based on the earlier calculation in [START_REF] Goldreich | Wave generation by turbulent convection[END_REF]. We first assume that turbulence driven by the convection can be decomposed into "eddies," each of which has a characteristic length h and velocity u h . The largest such "eddy" is the large-scale circulation of the convection, which has length c = 2π/k c and velocity u c . We imagine this large scale circulation represents the injection scale of a Kolmogorov turbulent cascade and to recover the k -5/3 energy spectrum (k is the three-dimensional wavenumber), we require that the typical velocity of eddies with size h scales like u h /u c = (h/ c ) 1/3 . We assume each eddy is coherent for an eddy turnover time τ h = h/u h , which implies that the power spectrum in the convection zone is peaked along the curve τ h f c = (h/ c ) 2/3 . This ignores the fact that small eddies can be advected over their size h on timescales shorter than τ h by large-scale motions. This set of assumptions is very simplistic, but allows a theoretical derivation of the wave energy flux spectrum.

In this hypothesis, the waves are generated by Reynolds stresses within the convection zone. We calculate the waves generated by a single eddy with size h and corresponding turnover time τ h , and then sum up the wave energy flux over all eddies. An eddy with size h and turnover time τ h can only efficiently excite a wave with frequency f and wavenumber

k ⊥ if both f τ -1 h , c (k ⊥ /2π) c /h = (τ h f c ) -3/2 . (3.1)
This is because we assume the spectral representation of an eddy's velocity field is peaked at wavenumbers around 2π/h and decreases exponentially at higher wavenumbers due to analyticity (and similarly for the temporal spectrum with timescale τ h ). Waves with f > τ -1 h or k ⊥ > 2π/h then have exponentially small coupling to eddies of size h. To connect the properties of the convection to the waves, one can write an equation for the wave velocity as a function of turbulent stresses as [START_REF] Lecoanet | Internal gravity wave excitation by turbulent convection[END_REF],

∇ 2 ∂ 2 t w + N 2 (z)(∂ 2 x + ∂ 2 y )w = S, (3.2)
where S is related to the Reynolds stresses associated with convective eddies. Note that we have neglected the effects of diffusion in the excitation process. We also neglect thermal stresses within the convection zone (as in [START_REF] Lecoanet | Numerical simulations of internal wave generation by convection in water[END_REF]. [START_REF] Lecoanet | Internal gravity wave excitation by turbulent convection[END_REF] derives the internal waves' Green's function from this equation assuming N 2 (z) is constant in the convective and stably-stratified region, and discontinuous at the interface. One can then estimate the convolution of this Green's function with the Reynolds stress of an eddy to find that the wave energy flux from eddies of size h is

dF i d log k ⊥ d log f d log h = Au 3 h f c f N f f c k ⊥ h 2π 4 = Au 3 c f c f N f f c h c k ⊥ h 2π 4 , (3.3)
with A a multiplicative constant, assumed universal, i.e. independent of all problem, wave or eddy parameters. To calculate the wave energy flux for a given frequency and wavenumber, we must sum up excitation from all eddies which satisfy equation (3.1).

Because of the strong dependence on the eddy size in equation (3.3), the wave energy flux is dominated by the largest possible eddies which can excite a wave. The wave energy flux depends on whether the wave frequency is higher or lower than the frequency of the largest generating eddy

f * (k ⊥ ) = f c [ c k ⊥ /(2π)] 2/3 . Waves with (k ⊥ , f ) satisfying f < f * (k ⊥ )
are not observed in the bulk of the stable region in DNS because they are strongly dissipated (see discussion after (3.5)), so we neglect their contribution to the energy flux in the theory. The expression for the flux of the waves we keep, i.e. with frequency

f > f c [ c k ⊥ /(2π)] 2/3 , is dF i d log k ⊥ d log f = Au 3 c f c f N f f c -13/2 k ⊥ c 2π 4 . (3.4)
We recall that equation (3.4) is a prediction for the energy flux spectrum of inviscid propagating wave modes at the convective-stably-stratified interface. In the DNS, however, the energy flux at the interface has contributions from both inviscid and viscid wave modes and overshooting plumes, such that it is not expected to match the prediction (3.4). In order to avoid this discrepancy, we compare the theory with DNS away from the interface, where overshooting plumes and strongly-viscous modes can be neglected.

In order to extend the expression for the flux (3.4) at the interface to an expression valid in all of the stable layer we derive a wave decay rate from the linear equation

(∂ t -∇ 2 ) 2 ∇ 2 w + N 2 ∂ 2 x + ∂ 2 y w = 0, (3.5)
valid for P r = 1. Note this differs from the inviscid equation (3.2) used to derive the wave flux at the interface. Substituting a plane wave solution w ∼ e i(kxx+kyy-ωt) e (ikz-γ)z with ω > 0 in (3.5) yields a cubic complex equation for (ik z -γ) 2 whose six roots can be written in closed form. We assume that motions in the stable layer are dominated by waves propagating upward, i.e. with group velocity C gz = ∂ω/∂k z > 0, and whose amplitude decrease with z. One of the three roots (k zj , γ j ) with C gzj > 0 has a much smaller decay rate γ j than the other two, such that it is expected to be the dominant contribution to upward-propagating and decaying waves in DNS. Unsurprisingly, this solution is the "full-dissipation" equivalent of the simpler WKB solution

k z = ∓k ⊥ N/ω, γ = ±k 3 ⊥ N 3 /ω 4 , (3.6)
with k z < 0 for upward-propagating waves, derived in the weak dissipation limit and valid when ω k 2 O(1) and ω 2 N 2 . We neglect the other two more-rapidly decaying upward-propagating wave solutions of (3.5). Note that whether they capture some of the turbulent energy is beyond the scope of this work, and neglecting them is consistent with prediction (3.4) since the theory only considers the WKB mode (3.6) in its inviscid form. Given that a broad range of frequencies and wavenumbers is excited in the DNS, it may be expected that (3.6) does not accurately predict the decay rate of the dominant waves. However, we find that the WKB solution and the "full-dissipation" solution for γ yield similar results for the theoretical spectra presented here, so for simplicity we will use (3.6) in the remainder of this paper (note that this suggests that the discrepancy between the "full" and the WKB decay rate is significant only when both decay rates are so large that it affects waves whose amplitudes are so small that they can be neglected). Note for future discussions that the decay rate is strongest for low-frequency high-wavenumber waves, as clearly shown by (3.6).

Finally, the z-dependent theoretical wave flux spectrum becomes

dF th d log k ⊥ d log f = dF i d log k ⊥ d log f e -2γ(z-zi) = Au 3 c 2πf c N f f c -13 2 k ⊥ c 2π 4 e -2γ(z-zi) , (3.7)
with γ given by (3.6) and z i the interface height. A useful approximate expression can be obtained for the total flux as a function of z by integrating (3.7) in the limit f N f c (z -z i )k 3 c /π and sufficiently far from the interface, i.e. z -

z i [γ(k c , f c )] -1 .
The derivation is tedious but straightforward, yielding (cf. Appendix A)

F th = f N fc kc( f fc ) 3 2 kc dF th df dk ⊥ dk ⊥ df ≈ Au 3 c f 15 2 c f 47 8 N 5Γ 5 8 39 8 c 28(8π 2 ) 13 8 (z -z i ) -13 8 ≡ F (z -z i ) - 13
8 .

(3.8)

The domain of validity of the asymptotic solution is expected to encompass all of the wave region in DNS. In the next section we show that the wave flux spectrum in DNS follows the trends predicted by the theory, i.e. equations (3.4) and (3.8).

Results and discussion

Statistics of turbulent convection

The theory described in section 3 assumes that the spatial correlations of the turbulent convection follow the standard Kolmogorov prediction for homogeneous and isotropic turbulence. Figure 2(a) shows that the frequency-integrated kinetic energy spectrum fi dKδf indeed has a slope relatively close to k -5/3 ⊥ for all simulations, thus validating the use of Kolmogorov hypothesis. The theory also requires prescription of the temporal correlations of the turbulence fluctuations. Here, the theory uses the eddy turnover time as the relevant time scale, which is known as the straining hypothesis. However, one could argue that the temporal correlations of small-scale turbulent features seen at a fixed point in space are affected by the advection by large-scale motions, and that the sweeping hypothesis should be considered instead [START_REF] Tennekes | Eulerian and lagrangian time microscales in isotropic turbulence[END_REF][START_REF] Chen | Sweeping decorrelation in isotropic turbulence[END_REF]. The wavenumber integrated kinetic energy spectrum in 2(b) shows a frequency scaling somewhat weaker than f 5/3 , which is the typical scaling for turbulent flows in Eulerian frame with sweeping [START_REF] Chevillard | Intermittency of velocity time increments in turbulence[END_REF][START_REF] Canet | Spatiotemporal velocityvelocity correlation function in fully developed turbulence[END_REF], in agreement with convection experiments [START_REF] Sano | Turbulence in helium-gas free convection[END_REF]. The straining mechanism would instead lead to a steeper f -2 scaling, which is recovered when computing temporal correlations along Lagrangian trajectories in homogeneous turbulence [START_REF] Chevillard | Intermittency of velocity time increments in turbulence[END_REF] and in Rayleigh-Bénard convection [START_REF] Liot | Simultaneous temperature and velocity lagrangian measurements in turbulent thermal convection[END_REF]. Figure 2(c) shows k ⊥ f dK for simulation C 100 9 in (k ⊥ , f ) space, thus providing a more detailed picture of the turbulence statistics of the most turbulent simulation (similar results are obtained for other simulations). Again, it can be seen that the line of maximum k ⊥ f dK is reasonably well approximated by a power fit of the form f ∼ k a ⊥ (solid black line), with a close to 1 for all simulations, suggesting strong sweeping effects. Isocontours at higher frequencies, however, have a milder slope f ∼ k b ⊥ with b < 1, possibly suggesting that the eddy dominant frequency does scale like ∼ k 2/3 ⊥ (solid white line) and that sweeping does not affect all eddies. We will show that there is a good agreement between DNS results and the theory based on the straining hypothesis despite the presence of a strong large-scale flow. It should be noted that the theory in §3 changes with the choice of temporal correlations, as is already the case for other wave generation mechanisms by turbulence (for sound waves, see e.g. [START_REF] Zhou | Sweeping and straining effects in sound generation by high reynolds number isotropic turbulence[END_REF][START_REF] Favier | On space and time correlations of isotropic and rotating turbulence[END_REF].

Integrated wave energy flux

In this section we compare the total energy flux F (z) obtained in DNS with the theoretical prediction F th (z) (cf. (3.8)), and we infer the free parameters of the theory A and f c from a best fit. The theoretical total flux F th has a steep scaling with the buoyancy frequency f N , the interface height z i and the turbulence parameters c , f c and u c . While f N and z i are relatively well-constrained, the turbulence parameters are challenging to estimate from DNS of turbulent convection, making F th potentially highly sensitive to the definitions used. A slight variation of f c by a factor of 2, due to a change of definition, can lead to a change of F th by a factor of ≈ 200 at all heights. Here, instead of estimating the values of the turbulence parameters using ad-hoc definitions, we use a best-fit approach. For simplicity, we assume that c = 1 in all cases and thus take u c = f c . This leaves f c as the only free parameter in the theory, besides A, which is assumed to be the same for all simulations. A posteriori, we check that the best-fit f c is not too different from the rms velocity u rms = √ 2K at z = 0.7 in the convection zone (cf. equation (2.5) and note that u rms is relatively depth-invariant in the convective bulk).

We obtain f c for each simulation case by requiring that the asymptotic prediction for the total flux (3.8) equals the total flux in DNS at some height z * in the wave region, i.e.

f 21 2 c = F (z * ) A f 47 8 N 28(8π 2 ) 13 8 5Γ (5/8) (z * -z i ) 13 8 . (4.1)
This definition of f c may appear degenerate since f c can take different values depending on A or z * . However, as discussed below (cf. also figure 3a), once A is chosen, f c does not depend on z * provided that z * is sufficiently far from the interface (i.e. z * > 1.1). With A ∼ O(10), we find for all simulations that f c is almost equal to u rms (to within a factor of 2), which can be considered as a reference value for the turnover time of the largestscale eddies [START_REF] Lecoanet | Numerical simulations of internal wave generation by convection in water[END_REF]. From now on we will set the multiplicative constant arbitrarily to A = 30. Other values could be chosen, but do not lead to significant changes as long as f c is relatively close to u rms . We compare DNS results for F (z) (symbols) with the theoretical prediction F th (z) (solid lines) in figure 3(a) for the eight simulation cases, using the parameter values listed in table 1 (i.e. with A = 30). The total wave flux decreases with z in all cases and, as expected, is larger for higher Ra (for fixed S) and lower S (for fixed Ra) [START_REF] Couston | Dynamics of mixed convective-stably-stratified fluids[END_REF]. With f c obtained from a best-fit at just one particular height, we find an excellent agreement between DNS and the theoretical prediction in all of the bulk of the wave region (z ∈ [1.1, 1.8]) for all simulations, which demonstrates that f c would be similar for any z * ∈ [1.1, 1.8]. The (z -1) -13/8 scaling predicted by the asymptotic solution is well verified in figure 3(b), which shows the normalized wave flux F/ F (cf. equation (3.8)). All DNS results collapse on the asymptotic solution (z -1) -13/8 shown by the dashed line for (z -1) ∈ [0.1, 0.8], demonstrating that the theoretical flux spectrum integrated over the full range of waves excited compares well with the total flux in DNS.

Wave energy flux spectrum

We now investigate in details the scalings of the wave flux spectrum k ⊥ f dF in DNS with k ⊥ and f , and discuss how they compare with the theoretical prediction ∼ k 4 f -13/2 . We first show in figure 4(a)-(c) the overall wave flux spectrum k ⊥ f dF in (k, f ) space at heights z = 1, 1.1, 1.5 for simulation C 100 9 , which is the most turbulent case. The wave toward higher frequencies (resp. lower wavenumbers) at higher z as a result of damping. We have checked that similar trends in frequency space and wavenumber space are obtained for other low wavenumbers (within the inertial subrange) and high frequencies (f u rms ), suggesting an overall good agreement between DNS and theory for C 100 9 . We now demonstrate that a good agreement is obtained between DNS and the theory not only for C 100 9 but for all simulations listed in table 1. In figures 6(a)-(c), we show the wave flux spectrum at z = 1.3 as a function of f for k ⊥ = 3.7 with Ra = 4 × 10 7 , 2×10 8 , 10 9 in (a), (b), (c). Clearly, the DNS results (solid lines) overlap well with the theoretical predictions (dashed lines) obtained after estimating the convective frequency f c from a best-fit approach (cf. §4.2 and listed values in table 1). A good agreement is similarly obtained for the wave flux spectrum at z = 1.3 as a function of k ⊥ and for f /f N = 0.4 between DNS and the theory (cf. figures 6(d)-(f)). We have checked that similar overlaps between DNS and theory are obtained for other fixed wavenumbers in the inertial subrange and frequencies f N > f f c .

Concluding remarks

We have shown that the scaling of the energy flux spectrum of internal waves generated by turbulent convection in DNS follows the k 4 ⊥ f -13/2 prediction in the high S 1 regime [START_REF] Lecoanet | Internal gravity wave excitation by turbulent convection[END_REF], and consistently, the total flux decreases as (z -1) -13/8 for all simulations presented. As discussed in the previous section, we can match the flux predicted by the theory to the DNS if we multiply the analytical prediction by a universal constant A and use a best-fit for f c . It is unsurprising that the theory needs be rescaled since numerical factors have been dropped off along the derivation, and it is equally unsurprising that the best-fit f c equals a few times the reference frequency based on the rms velocity u rms . It is worth pointing out, however, that there are some variations of f c /u rms across cases, which may mask some effects not included in the theory. Future work should explore the sweeping hypothesis, and include longer integration time with additional simulations or laboratory experiments in the high turbulence, high S regime.

It is remarkable that there is such a satisfactorily agreement between the DNS and an approximate-at-best theory based on second-order statistics of turbulent motions. The theory, as simple as it is, may prove useful in assessing the detectability of internal waves in geophysical and astrophysical fluids such as stars and in estimating wave-driven changes of global dynamics from the energy in turbulent sources. with Ei the exponential integral and the final expression is the leading-order approximation in the limits f N /f c 1, z * 1 and z * (f c /f N ) 4 1. In the same limits, I 2 ≈2 21 Γ(5/8)(z * ) -13/8 , such that adding the leading-order approximations for I 1 and I 2 yields equation (3.8). Equation (3.8) is valid for roughly all simulation parameters and (z -z i ) ∈ [0.01, 10].

Figure 1 . 8 S

 18 Figure 1. Snapshots of (a) w(y = 0), (b) Tz -T z at y = 0 (overbar denotes x average), (c) w(z = 0.7), (d) w(z = 1.3) for case C 400 8 . Variables in the wave region (z > 1) in (a), (b) have been multiplied by 10 4 , 10 3 , respectively. A movie is available as Supplementary Material.

Figure 2 .

 2 Figure 2. Integrated kinetic energy at z = 0.7 as a function of (a) k ⊥ and (b) f . Solid, dash-dash, and dash-dot lines correspond to Ra = 4 × 10 7 , 2 × 10 8 , 10 9 and different colors correspond to different stiffnesses. Black solid lines are power laws with slope -5/3. (c) Kinetic energy spectrum k ⊥ f dK at z = 0.7 for C 100 9 . The white star shows the wavenumber-frequency pair ( -1 c , fc) of the large-scale circulation that allows a good agreement of the DNS with the theory. The pair ( -1 c , fc) is not unreasonably far from the line f (k ⊥ ) of maximum energy, which is shown by the black dash-dash line (with solid line showing the linear fit). Dashed purple lines show fc, fN . The solid white line shows slope +2/3.

Figure 3 .

 3 Figure 3. (a) Total flux in the wave region in DNS (triangles, circles, stars show results for Ra = 4 × 10 7 , 2 × 10 8 , 10 9 ) along with the predictions F th (solid lines) using the best-fit values for A and fc reported in table 1. (b) The normalized flux F/ F collapses for all simulations on the prediction F th / F ≈ (z -1) -13/8 (dashed line) for (z -1) ∈ [0.1, 0.8] (cf. Eqn (3.8)).

Figure 4 .

 4 Figure 4. Wave flux spectrum k ⊥ f dF in DNS for case C 100 9 at (a) z = 1.0, (b) 1.1, (c) 1.5. The color scale is shifted to lower amplitudes from (a) to (c) because the wave flux decreases with height. Solid black lines show isocontours of the decay rate 2γ. Dashed purple lines show fc, fN . White bins indicate k ⊥ f dF values that are negative.

Figure 5 .

 5 Figure 5. Wave flux spectrum k ⊥ f dF for fixed k ⊥ /(2π) = 1.7, 8.2 [(a)-(b)] and fixed f /fN = 0.2, 0.4 [(c)-(d)] for case C 100 9 . Lighter colors show results at higher heights. Solid lines show DNS results and dashed lines show theoretical results using the convection parameters listed in table 1. The blue lines show the theoretical scalings.

Figure 6 .

 6 Figure 6. Wave flux spectrum k ⊥ f dF at z = 1.3 for fixed k ⊥ /(2π) = 3.7 [(a)-(c)] and fixed f /fN = 0.4 [(d)-(f)] with Ra = 4 × 10 7 , 2 × 10 8 , 10 9 increasing from (a) to (c) and (d) to (f). The different colors correspond to different stiffnesses S. Solid lines show DNS results and dashed lines show theoretical results using the convection parameters listed in table 1.

  upper incomplete gamma function. Integrating I 1 with respect to ξ yields

Table 1 .

 1 Simulation

	1	9.8	1.6	19

parameters: number of Fourier modes nx, ny in x, y directions, Chebyshev modes n zl /nzu in lower/upper regions of the compound z basis (stitched at z = 1.2), typical time steps δt, best-fit convective frequency fc for A = 30 (cf. §4.2), time fc∆t of data collection, convective velocity urms and Brunt-Väisälä frequency fN . P r = 1 for simplicity.

† louisalexandre.couston@gmail.com

z * I2

Acknowledgements. The authors acknowledge funding by the European Research Council under the European Union's Horizon 2020 research and innovation program through Grant No. 681835-FLUDYCO-ERC-2015-CoG. D. L. is supported by a PCTS fellowship and a Lyman Spitzer Jr fellowship. Computations were conducted with support by the HPC resources of GENCI-IDRIS (Grant No. A0020407543 and A0040407543) and by the NASA High End Computing (HEC) Program through the NASA Advanced Supercomputing (NAS) Division at Ames Research Center on Pleiades with allocations GID s1647 and s1439.

Appendix A

Here we present the key steps required to derive (3.8). Let us use the variables


Géraldine Pichot

Patrick Laug

Jocelyne Erhel

Romain Le Goc

Caroline Darcel

Philippe Davy

Jean-Raynald De Dreuzy

Flow simulations in geology-based Discrete Fracture Networks

Motivation

The underground is a reservoir of natural resources (water, oil and gas, heat,...) and a potential warehouse storage solution. Using these resources and storage facilities in a sustainable way requires a good understanding of the physical, chemical and biological processes happening there. Also, the geometry of the subsurface couples these processes together. Here, numerical models are very useful: they reduce the costs and risks of in situ experiments and allow long-term predictions.

Geometric model

At first, we need to build a geometric model of the subsurface. Observations and experiments (investigations of drill cores and drill cuttings, lineament characterization, measurements of reflection seismics, transient electromagnetic soundings, ... [START_REF]Svensk Kärnbränslehantering AB (SKB), Site description of Forsmark at completion of the site investigation phase. SDM-Site Forsmark[END_REF]) give data for soil characterization and fracture locations. Fractures cannot be neglected as they critically impact the response of subsurface systems through their complex organisation over a broad range of scales. However, the information collected via the measurements are local. To derive an image of the full geometry, numerical models come into play. It is quite challenging to process tremendous amounts of noisy data, scattered and at different scales. Nevertheless, hydrogeologists are able to derive models for the subsurface, and especially for fracture systems [START_REF] Maillot | Connectivity, permeability, and channeling in randomly distributed and kinematically defined discrete fracture network models[END_REF]. The one we rely on is based on a spatial organization of fractures satisfying two main regimes: a dilute regime for the smallest fractures, where they can grow independently of each other, and a dense regime for which the density distribution is controlled by the mechanical interactions between fractures [START_REF] Davy | A likely universal model of fracture scaling and its consequence for crustal hydromechanics[END_REF][START_REF] Davy | A model of fracture nucleation, growth and arrest, and consequences for fracture density and scaling[END_REF]. In these models, called UFM (likely Universal Fracture Model), the fracture size distribution matches the observations and large fractures inhibit the smaller ones, creating T-termination configurations. Models mix spatial dimensions as fracture networks consist in a large number of 2D fractures interconnected in the 3D space [START_REF] Martin | Modeling fractures and barriers as interfaces for flow in porous media[END_REF]. These models are implemented in a dedicated software, called UFMLab. UFM networks may contain millions of fractures.

Flow problem

Let us focus on a simple process, which is the flow of water within the fractures. The surrounding rock matrix is assumed impervious. The associated problem is governed by the Poiseuille's law and the mass conservation equation in each fracture combined with continuity conditions at the intersections between fractures. The transmissivity within the fractures is assumed heterogeneous and modeled by random fields whose parameters are given by the experiments.

Mesh generation

The next step is to build a mesh of the domain. To easily enforce the continuity conditions at the intersections, we request a mesh which takes the intersections into account, either with matching or non-matching grids at the intersections. However, even in the non-matching cases, the intricate configurations of the fracture systems make it difficult for standard planar meshers to generate a mesh of the geometry. The idea is then to discretize first the border and intersections to build a geometric model, valid for the planar mesher. It implies automatic corrections to build valid curve discretizations. These corrections are implemented in the software BLSURF FRAC1 . Then a planar mesher is called, for example, BL2D [START_REF] Borouchaki | Parametric surface meshing using a combined advancing-front -generalized-Delaunay approach[END_REF] or BAMG [START_REF] Hecht | BAMG: Bidimensional Anisotropic Mesh Generator[END_REF]. Up to now, the largest network successfully meshed with this technique contains 508,339 fractures which generate 1,031,231 intersections. It is meshed in 15 minutes using the combination of BLSURF FRAC and BL2D, with 4 threads on a PC Intel Core i7 4 cores CPU @ 2.90 GHz, 32GB RAM memory. The final mesh contains 8,112,299 of triangles and 1,630,682 of intersection edges. Notice that the computational time can be improved on a computer with more RAM memory as memory swapping was detected during this numerical experiment. Also, we identified some sequential parts of BLSURF FRAC that can be parallelized.

Numerical methods and implementation

The next step is to solve the flow problem efficiently. We propose to solve this problem using a Mixed Hybrid Finite Elements Method (MHFEM) [START_REF] Arnold | Mixed and nonconforming finite element methods: postprocessing, and error estimates[END_REF][START_REF] Brezzi | Mixed and Hybrid Finite Element Methods[END_REF][START_REF] Roberts | Mixed and hybrid methods, Handbook of numerical analysis[END_REF]. Possible non-matching grids at the intersections between fractures are handled by using suitable Mortar conditions [START_REF] Arbogast | On the Implementation of Mixed Methods as Nonconforming Methods for Second-Order Elliptic Problems[END_REF][START_REF] Arbogast | Mixed finite element methods on non-matching multiblock grids[END_REF][START_REF] Maryška | Numerical simulation of fracture flow with a mixed-hybrid FEM stochastic discrete fracture network model[END_REF][START_REF] Pichot | A mixed hybrid Mortar method for solving flow in discrete fracture networks[END_REF][START_REF] Pichot | A generalized mixed hybrid mortar method for solving flow in stochastic discrete fracture networks[END_REF]. In both, the matching and non-matching grids cases, we end up with a linear system whose matrix is symmetric positive definite. It is a sparse matrix with a L-shape structure. The linear system can be solved either with a direct solver (Cholesky factorization) or with an iterative solver (preconditioned conjugate gradient or multigrid). A way to improve the computational time of iterative solvers is to use domain decomposition techniques [START_REF] Poirriez | Study and implementation of a domain decomposition method for the modelisation of flow in 3D fracture networks[END_REF][START_REF] Pichot | A Mortar BDD method for solving flow in stochastic discrete fracture networks[END_REF]. We developed a Matlab code, called NEF-Flow, which implements the MHFEM either for matching or non-matching meshes, with sink/source terms and contrasts in transmissivities. Efficiency is obtained using Matlab vectorization. For example, we solved the flow problem in the largest network described above in 21 minutes on a PC Intel Core i7 4 cores CPU @ 2.90 GHz, 32GB RAM memory with a direct solver. Here again memory swapping was detected. Further numerical experiments are now performed with iterative solvers to reduce memory requirements. Also to keep an accurate solution with less triangles, we are investigating the use of a higher order hybridizable discontinuous Galerkin technique, the HHO method [START_REF] Di Pietro | A Hybrid High-Order locking-free method for linear elasticity on general meshes[END_REF], combined with a posteriori error estimates [START_REF] Ern | Polynomial-degree-robust a posteriori estimates in a unified setting for conforming, nonconforming, discontinuous Galerkin, and mixed discretizations[END_REF] in order to refine only the fractures that carry most of the flow.

http://pages.saclay.inria.fr/patrick.laug/logiciels/logiciels.html

Acknowledgments. The author thanks the organizers of the workshop for their kind invitation. The author thanks the MFO and all its staff for the extremely good working and living conditions provided at Oberwolfach.


Guillaume Miquelard

Alain Guinault

Cyrille Sollogoub

A Amor

N Okhay

G Miquelard-Garnier

M Gervais

Combined compatibilization and plasticization effect of low molecular weight poly(lactic acid) in poly(lactic acid)/ poly(3-hydroxybutyrate-co-3-hydroxyvalerate) blends

Keywords: polymer blends and alloys, biopolymers, toughening, compatibilization, barrier properties

    Combined compatibilization and plasticization effect of low molecular weight poly(Lactic acid) in poly(lactic acid)/poly(3-hydroxybutyrate-co-3-hydroxyvalerate) blends

Introduction

Biopolymers have received great attention both in industry and in academia due to an increased concern toward the environmental impact of plastic waste and the saving of limited fossil energy. Poly (lactide) (PLA), one of the most used biobased and biodegradable polymers, is synthesized from lactic acid obtained from renewable resources (corn, sugar, potato). Its properties make it a potential candidate to be used in various applications, such as packaging, medicine, agriculture, textile, and automotive industries [START_REF] Auras | An overview of polylactides as packaging materials[END_REF]. However, despite its good optical properties and high tensile strength, its inherent brittleness and high gas permeability significantly limit its applications, especially in food packaging [START_REF] Auras | An overview of polylactides as packaging materials[END_REF][START_REF] Garlotta | A literature review of poly(lactic acid)[END_REF][START_REF] Nampoothiri | An overview of the recent developments in polylactide (PLA) research[END_REF][START_REF] Drumright | Polylactic acid technology[END_REF]. Many strategies have been recently developed and proposed to improve these properties, such as the addition of modifiers, copolymerization or blending [START_REF] Anderson | Toughening polylactide[END_REF][START_REF] Liu | Research progress in toughening modification of poly(lactic acid)[END_REF][START_REF] Kfoury | Recent advances in high performance poly(lactide): From 'green' plasticization to super-tough materials via (reactive) compounding[END_REF][START_REF] Rasal | Poly(lactic acid) modifications[END_REF].

Blending PLA with another biobased and biodegradable polymer has attracted much attention since it offers simply the opportunity to widen the range of physical properties, without compromising environmental features of the final material. Among polyhydroxyalcanoates (PHAs), which are aliphatic polyesters obtained by microbial fermentation, poly(3hydroxybutyrate) (PHB) and its copolymers (poly(3hydroxybutyrate-co-3-hydroxyvalerate) (PHBV), poly (3-hydroxybutyrate-co-3-hydroxyoctanoate (PHBHO), poly(3-hydroxybutyrate-co-3-hydroxyhexanoate) (PHBHHx)…) have been largely used for PLA-based blends [START_REF] Avella | Review properties of blends and composites based on poly(3-hydroxy)butyrate (PHB) and poly(3-hydroxybutyrate-hydroxyvalerate) (PHBV) copolymers[END_REF]. Some of the resulting blends exhibit improved properties compared to neat PLA: higher toughness and increased elongation at break [START_REF] Noda | Polymer alloys of nodax copolymers and poly(lactic acid)[END_REF][START_REF] Ma | Toughening of poly (lactic acid) by poly (β-hydroxybutyrateco-β-hydroxyvalerate) with high β-hydroxyvalerate content[END_REF][START_REF] Bartczak | Tough blends of poly(lactide) and amorphous poly([R,S]-3-hydroxy butyrate) -Morphology and properties[END_REF][START_REF] Lim | Effect of composition ratio on the thermal and physical properties of semicrystalline PLA/PHB-HHx composites[END_REF][START_REF] Zhang | Blending polylactic acid with polyhydroxybutyrate: The effect on thermal, mechanical, and biodegradation properties[END_REF][START_REF] Ferreira | Films of PLLA/PHBV: Thermal, morphological, and mechanical characterization[END_REF][START_REF] Nanda | The effects of process engineering on the performance of PLA and PHBV blends[END_REF], improved gas barrier properties [START_REF] Boufarguine | PLA/PHBV films with improved mechanical and gas barrier properties[END_REF][START_REF] Zembouai | A study of morphological, thermal, rheological and barrier properties of poly(3-hydroxybutyrate-co-3-hydroxyvalerate)/polylactide blends prepared by melt mixing[END_REF], increased thermal stability [START_REF] Nanda | The effects of process engineering on the performance of PLA and PHBV blends[END_REF], enhanced biodegradability [START_REF] Zhang | Blending polylactic acid with polyhydroxybutyrate: The effect on thermal, mechanical, and biodegradation properties[END_REF].

Those improvements are partly due to the fact that, although immiscible, the two polymers have a certain degree of chemical affinity. Miscibility and compatibility of PLA/PHB(V) have been thoroughly studied and have been found to depend classically on the molecular weight of both components, the blending method, the blend composition and the crystallinity of both components. Blümm and Owen [START_REF] Blümm | Miscibility, crystallization and melting of poly(3-hydroxybutyrate)/poly(L-lactide) blends[END_REF] observed miscibility of low molecular weight poly(L-lactide) (PLLA) ( -M n = 1800 g/mol) and immiscibility of high molecular weight PLLA ( -M n = 160 000 g/mol) blended with PHB ( -M n = 222 000 g/mol) over the whole composition range. Similarly, Koyama and Doi [START_REF] Koyama | Miscibility of binary blends of poly[(R)-3-hydroxybutyric acid] and poly[(S)-lactic acid][END_REF] investigated the miscibility of PHB ( -M w = 650 000 g/mol) blended with PLA of various molecular weights and determined that the blend is immiscible above a critical value of molecular weight of PLA around 20 000 g/mol. More generally, all literature agrees on the immiscibility of the blend as soon as the molecular weight of both components is high. Zhang et al. [START_REF] Zhang | Miscibility, crystallization and morphology of poly(β-hydroxybutyrate)/poly (D,L-lactide) blends[END_REF] found that melt blended samples exhibited greater compatibility than those prepared by solvent casting at room temperature, possibly due to a transesterification reaction between PHB and PLA occurring at high temperature. This transesterification, or more generally the presence of specific interactions between the two polymers, has been evidenced by many other authors [START_REF] Zhang | Blending polylactic acid with polyhydroxybutyrate: The effect on thermal, mechanical, and biodegradation properties[END_REF][START_REF] Ferreira | Films of PLLA/PHBV: Thermal, morphological, and mechanical characterization[END_REF][START_REF] Gerard | Morphology and molten-state rheology of polylactide and polyhydroxyalkanoate blends[END_REF][START_REF] Gérard | Polylactide/poly(hydroxybutyrate-co-hydroxyvalerate) blends: Morphology and mechanical properties[END_REF][START_REF] Zhang | Crystallization behaviors of poly(3-hydroxybutyrate) and poly(L-lactic acid) in their immiscible and miscible blends[END_REF], especially for blends containing a low amount of PHB(V). For such PLA/PHBV blends (typically weight fractions of 90%/10%), Gerard and Budtova [START_REF] Gerard | Morphology and molten-state rheology of polylactide and polyhydroxyalkanoate blends[END_REF] observed a peculiar morphology, with very small PHBV nodules (about 400 nm) well dispersed in PLA matrix, with increased adhesion at the interface, presumably explaining the increase of ductility (elongation at break of 200% when tested immediately after sample preparation) measured for this composition by the authors. All those studies have highlighted the importance of the interface on the final properties of PLA/PHB(V) blends and some attempts have been undertaken to modify it, using different strategies. Yoon et al. [START_REF] Yoon | Effect of poly(ethylene glycol)-block-poly(Llactide) on the poly[(R)-3-hydroxybutyrate]/poly(L-lactide) blends[END_REF] compatibilized a mixture of PHB and PLA ( -M n = 56 000 g/mol) using synthesized Poly(ethylene glycol) -b -PLLA (PEG-b-PLLA) (PEG is known to be compatible with PHB) and poly(vinyl acetate) (PVAc, compatible with both PLA and PHB). A small improvement in toughness was obtained with the addition of 2 wt% of block copolymers or PVAc, but the overall mechanical properties were not improved much. One possible reason proposed by the authors could be the solubilization of the compatibilizer molecules in either PHB or PLA, or in both, instead of placing themselves at the interface between the PHB and PLA. Recently, an alternative approach was proposed by Yang et al. [START_REF] Yang | Transesterification induced mechanical properties enhancement of PLLA/PHBV bio-alloy[END_REF] by adding transesterification catalysts in order to induce chemical reaction between the two components. They found a significant effect on the mechanical properties by adding a small amount of zinc acetate (0.1 wt% in the blend). Finally, some authors have formulated PLA/PHB(V) blends in order to improve their mechanical properties and/or processability, adding different components as plasticizers, either commercial or synthesized ad-hoc: PEG [START_REF] Wang | Mechanical, thermal and degradation properties of poly (D,L-lactide)/poly(hydroxybutyrate-co-hydroxyvalerate)/poly(ethylene glycol) blend[END_REF][START_REF] Arrieta | Combined effect of poly(hydroxybutyrate) and plasticizers on polylactic acid properties for film intended for food packaging[END_REF], ATBC [START_REF] Arrieta | Combined effect of poly(hydroxybutyrate) and plasticizers on polylactic acid properties for film intended for food packaging[END_REF], limonene [START_REF] Arrieta | Ternary PLA-PHB-limonene blends intended for biodegradable food packaging applications[END_REF], polyester plasticizer (Lapol 108) [START_REF] Abdelwahab | Thermal, mechanical and morphological characterization of plasticized PLA-PHB blends[END_REF], lactic acid oligomer (oLA) [START_REF] Armentano | Processing and characterization of plasticized PLA/PHB blends for biodegradable multiphase systems[END_REF] and tributyrin [START_REF] D'amico | Fully bio-based and biodegradable polylactic acid/poly(3-hydroxybutirate) blends: Use of a common plasticizer as performance improvement strategy[END_REF]. It has been accepted that concentrations from 10 to 20 wt% and up to 30 wt% are required to provide a substantial plasticization effect, which strongly modifies the mechanical properties: if the ductility is usually enhanced, the Young's modulus and the tensile strength are decreased. Besides, the gas barrier properties of the plasticized blends are often deteriorated as a consequence of the incorporation of the plasticizer, which induces a gain of mobility of the polymer chains. Nevertheless, some authors noticed that, when the plasticizer has good miscibility with both of the components, it can have a beneficial effect on the compatibility between them. For instance, Armentano et al. [START_REF] Armentano | Processing and characterization of plasticized PLA/PHB blends for biodegradable multiphase systems[END_REF] observed a reduction of the size of PHB domains with the addition of oLA ( -M n = 957 g/mol) in PLA/PHB (85/15 wt%), indicating improved interfacial properties. The enhanced ductility can therefore be interpreted as the consequence of a combined effect of plasticization and compatibilization. This oLA was added from 15 to 30 wt%. In this study, we develop a different strategy to compatibilize PLA/PHBV blends, while preserving their tensile strength, gas barrier properties and bio-based character. It consists in introducing a much smaller amount (between 0.1 and 1 wt%) of hydrolyzed PLA having an oligomer-type molar mass ( -M n = 5 000 g/mol) to the PLA/PHBV blends. It is anticipated that the low molecular weight PLA (oLA) will improve the compatibility between PLA and PHBV without a significant drop in the mechanical properties. The resulting thermal, morphological, mechanical and gas barrier properties of the obtained biobased PLA/PHBV blends are investigated.

Experimental part 2.1. Materials

2003D PLA, an extrusion grade, was purchased from Natureworks, USA. This PLA is a Poly(D,Llactide) with a percentage of D-lactic acid units of 4.3% [START_REF] Li | Synergistic effects of PEG and MWCNTs on crystallization behavior of PLLA[END_REF]. Its melt temperature is close to 152 °C with a melt flow index of 6 g/10 min (210°C, 2.16 kg) as given by Natureworks. The molecular weight as determined by SEC is -M w ≈ 150 000 g/mol with a dispersity Ð of 2.5. SEC was performed at 40 °C using THF (GPR Rectapur, VWR, France) as eluant on a Waters apparatus (USA) equipped with three columns Styragel HR0.5, HR3 and HR4 and with a Waters 2414 refractive index detector at an elution rate of 1 mL/min. Polystyrenes were used as standards. ENMAT Y1000P Poly(3-hydroxybutyrate-co-3-hydroxyvalerate) PHBV was obtained from Tianan Biologic, China. Y1000P is an additivated PHBV containing 8 mol% of HV groups, and -M w ≈ 340 000 g/mol (Ð ≈ 2.5) [START_REF] Corre | Morphology and functional properties of commercial polyhydroxyalkanoates: A comprehensive and comparative study[END_REF], with a melt temperature close to 170 °C. The melt flow index is reported to be 2.4 g/10 min at 170°C (2.16 kg) [START_REF] Javadi | Microcellular poly(hydroxybutyrate-cohydroxyvalerate)-hyperbranched polymer-nanoclay nanocomposites[END_REF]. The chemical structure of the two polymers is shown Figure 1.

Synthesis of low molecular weight PLA

Low molecular weight PLA <20 000 g/mol is known to be miscible with PHB [START_REF] Blümm | Miscibility, crystallization and melting of poly(3-hydroxybutyrate)/poly(L-lactide) blends[END_REF][START_REF] Koyama | Miscibility of binary blends of poly[(R)-3-hydroxybutyric acid] and poly[(S)-lactic acid][END_REF]. The low molecular weight PLA (5000 g/mol) that will be added to PLA/ PHBV is obtained by the hydrolysis of 2003D PLA. PLA has indeed the advantage of being highly hydrolysable. This hydrolysis can be accelerated by many factors, such as temperature and pH of the reaction mixture [START_REF] Zhang | Morphological behaviour of poly(lactic acid) during hydrolytic degradation[END_REF]. The hydrolysis is realized in distilled water at neutral pH and at a temperature of 100°C. About 40 g of PLA are introduced into 4 L of water in a mount of reflux containing a round-bottom flask of 5 L. The use of a large amount of water will keep a pH close to 7 during degradation. Indeed the production of lactic acid during the degradation tends to acidify and catalyze the hydrolysis reaction which may cause a rapid degradation of the macromolecular chain or a change in the degradation mechanism.

The obtained samples are dried to a constant weight in a ventilated oven at 65 °C and their molecular weight is measured by size-exclusion chromatography (SEC) using the same protocol as previously described. SEC analyses of the PLA after 2 to 43 hours of hydrolysis have two molecular chains distributions (one around 5000 g/mol and one with a higher average molar mass). After 45 hours, a single distribution around 5000 g/mol is obtained. Therefore the oLA ( -M n ≈ 5000 g/mol) with Ð close to 1.25 is selected for this study.

Film samples preparations

PHBV and PLA were used as received, in pellet form, but dried with desiccant air before use (4 h at 90 °C for the PHBV, 4 h at 80 °C for the PLA). PLA, PHBV and oLA pellets were then dry-blended in proportions of (90 - x) wt%, 10 wt% and x wt% (with x = 0, 0.1, 0.5 or 1) respectively and then extruded in a single-screw extruder SCAMEX Rheoscam (France) (with a screw diameter of 20 mm and a length over diameter of 20, a screw speed of 60 rpm, and an operation temperature of 175-185°C). A flat die allowed obtaining films (thickness of 0.35 mm and width of 100 mm). The final composition of the blend will be given in PLA/ PHBV/oLA form in the following. PHBV and PLA additivated with oLA were also realized for some experiments using the following solvent cast procedure: oLA was introduced in PHBV (or PLA) dissolved in chloroform (GPR Rectapur, VWR, France) with different ratios (0, 1, 5 and 10 wt%) at 55 °C. The obtained solution was then transferred into aluminum pans which were placed in a vacuum oven at 65 °C to obtain films by solvent evaporation. Those blends will be referred as PHBV or PLA (solvent cast) and PHBV (or PLA)/oLA y wt% with y = 1, 5 and 10.

Determination of interfacial tension

A relevant parameter to discuss the compatibility of polymer blends is the interfacial tension between two polymers. To determine the interfacial tension, several methods can be used [START_REF] Xing | Critical experimental comparison between five techniques for the determination of interfacial tension in polymer blends: Model system of polystyrene/polyamide-6[END_REF], among which the relaxation drop method was chosen. A Linkam optical rheology system (UK) was used for shearing polymers under an optical microscope equipped with a camera. A small quantity of PHBV was placed between two PLA sheets and sheared until individual PHBV droplets can be isolated under the microscope at 185 °C. Then, a small shear rate was applied for a very short time on the sample: the droplet deforms into an ellipsoid. When shear is stopped, the droplet relaxes back to its initial spheroidal shape. The evolution of the long axis of the ellipsoid as a function of time follows an exponential decay with a characteristic relaxation time related to the viscosity ratio of the polymers and the interfacial tension following Equation ( 1):

(1) where τ d is the relaxation time due to the relaxation of the interface, η m is the viscosity of the matrix, γ 12 is the interfacial tension of the blend, r 0 is the mean radius of the drop, and K = η d /η m is the zero shear viscosity ratio of the droplet and matrix. The shear viscosity was deduced from dynamic shear measurements carried out at 185 °C on an Anton Paar MCR 502 rheometer using 25 mm diameter parallel plates (1 mm gap) under a shear strain of 1%. For assessing reproducibility, measurements were carried out on at least five different samples. In order to see the effect of oLA on the PLA/PHBV interfacial tension, the same measurement had been performed with PHBV additivated with oLA realized by solvent cast (see details above).

Samples characterization 2.5.1. Differential scanning calorimetry (DSC)

measurements Glass transition temperatures and melting enthalpies of the polymer samples were measured using a TA Instruments Q10 machine (USA), which was equipped with a DSC Refrigerated Cooling System to achieve low temperatures. Polymer samples were weighed (5 to 10 mg) and placed into non-hermetic aluminum pans individually. DSC analysis was done under nitrogen (50 mL/min) by first heating the sample to 200 °C at 10 °C•min -1 to erase its thermal history, then equilibrating the samples at -20 C, followed by ramping up the temperature to 200 °C at a rate of 10 °C per minute. The crystallinity ratio of the polymers c was calculated by Equation ( 2):

(2) where c [%] is the crystallinity ratio, ∆H m [J•g -1 ] the melting enthalpy, ΔH cc [J•g -1 ] the cold crystallization enthalpy, w the weight fraction of the polymer in the blend and ∆H m 0 the melting enthalpy of a 100% crystalline polymer, reported to be 93 J•g -1 for PLA [START_REF] Fischer | Investigation of the structure of solution grown crystals of lactide copolymers by means of chemical reactions[END_REF] and 109 J•g -1 for PHBV [START_REF] Scandola | Polymer blends of natural poly(3hydroxybutyrate-co-3-hydroxyvalerate) and a synthetic atactic poly(3-hydroxybutyrate). Characterization and biodegradation studies[END_REF].

Dynamic mechanical thermal analysis (DMTA) experiments

The thermo-mechanical behavior has been investigated using a dynamic thermomechanical analyzer DMTA from BOHLIN INSTRUMENTS (UK). A temperature scan from -40 up to 200 °C was performed at the rate of 2 °C•min -1 while a dynamic tensile test was performed at a frequency of 1 Hz with a deformation of 0.5%.

Tensile tests

Films mechanical properties (Young's modulus E, tensile strength σ y and elongation at break ɛ b ) were measured on an Instron 4507 (UK) with a 100 N load cell with no extensometer. At least five dogbone shaped samples, with 58 mm of length and 5 mm in width, were taken in the center of the film in the extrusion direction. The deformation rate was fixed at 5 mm•min -1 for all samples, force and deformation were measured and stresses and strains were calculated from these measurements. For each sample, average values of elongation at break and Young's modulus (calculated by taking the initial slope of the stress-strain curves) were collected and compared.

Atomic force microscopy (AFM) studies

The morphology and the size of the dispersed phase in the polymer blends were investigated via AFM. Samples were sectioned in the cross section of the film with an ultramicrotome 2088 UltrotomeV (LKB, Sweden) at a cutting speed of 1 mm/s. A VEECO, MultiMode Nanoscope-model V ( 2008) AFM (USA) working in Tapping mode, in air using silicon nitride tips was used in this study. Images were recorded using scan speeds on the order of 10 mm/s and the applied force was the minimum possible (~1 nN). The obtained images were processed by ImageJ software to determine the area, width and length of PHBV (or additivated PHBV) nodules within the PLA matrix. 5 images for each of two blends PLA/PHBV and PLA/PHBV/oLA 1 wt%, were analyzed, corresponding to a number of nodules higher than 100 (147 and 111 respectively).

r K K K 40 1 19 16 2 3 d m 12 0 $ x c h = + + + Q Q Q V V V X X w H H H 100 X c m m cc 0 $ D D D = -

Gas permeability

Helium permeability was measured by a specific homemade analyzer at 23 °C and 0% relative humidity, based on the ISO 15105-2:2003 method. Circular portions cut from the films (surface = 23.75 cm 2 ) were inserted between two hermetically sealed compartments drained using nitrogen. A helium constant flow (80 mL•min -1 ) was introduced in the downstream part of the cell and is measured in the upstream part, using an helium detector (mass spectrometer).

Measuring helium permeability rather than oxygen permeability presents several advantages: because helium molecules are smaller than oxygen, experimental time is reduced from 24 hours to 30 min. Moreover, helium is a neutral gas, which prevents from any possible interaction of the permeant gas within the polymer matrix and appears thus as a relevant probe of the blend morphology. Permeability was determined from the transmission rate. For assessing reproducibility of permeability, measurements were carried out on at least two different samples obtained from the same film.

Results and discussion

3.1. Characterization of PHBV and PLA additivated with oLA (obtained by solvent casting) Differential scanning calorimetry (DSC) analysis of different additivated PHBV and PLA formulations was conducted in order to determine the effect of the PLA short chains introduction on the thermal properties of PHBV and PLA. DSC thermal properties of the neat PHBV, PLA, oLA and of different PHBV/ oLA and PLA/oLA formulations are summarized in Table 1. Note that the ratio of oLA compared to PHBV in the formulations here is the same as in the final PLA/PHBV/oLA blends. Only the second heating has been considered for a better comparison with what happened after single screw extrusion. The T g value for neat PLA during the first and second heating scan is about 60 °C. However the thermogram of neat PHBV does not show any T g peak. Indeed DSC is not sensitive enough for the detection of this T g , given the high crystallinity of PHBV. From melting enthalpy, the degree of crystallinity of the different formulations can be determined. oLA has, as expected, a lower T g compared to PLA as well as a lower T cc and T m and a higher crystallinity ratio. Addition of oLA to PLA from 1 to 5 wt% does not influence the T g of PLA nor the degree of crystallinity. Therefore the addition of oLA by itself in PLA does not induce any modification of the behavior of PLA with the amount added to the system. oLA is not expected to have any influence in PLA/PHBV/oLA blends on the PLA part as the amount used will be lower or equal to 1 wt%. PHBV crystallinity ratio is evolving with the process used. Solvent casting the polymer tends to increase its crystallinity. The degrees of crystallinity, calculated with a ΔH m 0 value of 109 J/g for 100% crystalline PHBV, are extremely high. That may be due to an underestimated ΔH m 0 value considering the fact that other higher values are reported in literature (164 J/g [START_REF] Enriquez | Biobased blends of poly(propylene carbonate) and poly(hydroxybutyrateco-hydroxyvalerate): Fabrication and characterization[END_REF], 146 J/g [START_REF] Carli | PHBV nanocomposites based on organomodified montmorillonite and halloysite: The effect of clay type on the morphology and thermal and mechanical properties[END_REF][START_REF] Gogolewski | The effect of melt-processing on the degradation of selected polyhydroxyacids: Polylactides, polyhydroxybutyrate, and polyhydroxybutyrate-covalerates[END_REF]). A slight decrease is observed when oLA is added to PHBV. 5 wt% leads to the highest value of crystallinity ratio for the PHBV/oLA formulations been measured via the relaxation drop method. A typical drop relaxation of PHBV dispersed in a PLA matrix is shown in Figure 2 and the relaxation time from the ellipsoid shape after shear to circular has been used to estimate the interfacial tension for each blend using Equation (1). The results have been summarized in Table 2. From the results, a few remarks can be made. First of all, a decrease of the interfacial tension between PLA and PHBV after addition of the oLA is observed, revealing the compatibilizing effect of the oLA. In fact, when oLA (1 wt% in the PHBV phase) is added to the blend, the interfacial tension decreases from 0.5 to 0.3 mN/m, as presented in Table 2. From 1 to 10 wt% of oLA in PHBV, the interfacial tension remains almost unchanged. This phenomenon can be explained by the fact that the interface may be already saturated with 1 wt% of oLA.

Interfacial tension between PLA and additivated PHBV The interfacial tension between PLA and PHBV additivated with different concentrations of oLA have

Films characterization 3.3.1. Thermal analysis

The DSC thermograms from the first heating run of neat PLA, PHBV, PLA/PHBV and PLA/PHBV/oLA films realized by single-screw extrusion are shown in Figure 3. The curves reveal different thermal events as follows: the glass transition (T g ), the cold-crystallization (characterized by T cc and the enthalpy of crystallization ∆H cc ), the melting point (characterized by T m and the melting enthalpy ∆H m ). All these data are summarized in Table 3.

The T g value for the PLA neat film was about 60 °C, while the exothermic peak is observed at about 118 °C and the multi-step melting endothermic peak is obtained at about 151-157 °C. The DSC thermograms for the PLA/PHBV blend display a slight shift to lower T g values with respect to neat PLA and a similar multi-step melting with the first and second peaks (T m1 and T m2 ) due to the PLA and the third one (T m3 ) corresponding to the PHBV around 170 °C. This suggests no complete miscibility between the two polymers [START_REF] Arrieta | Multifunctional PLA-PHB/cellulose nanocrystal films: Processing, structural and thermal properties[END_REF].

The double melting peak in the range of 147-156 °C can be attributed to the formation of different crystal structures for the PLA [START_REF] Nanda | The effects of process engineering on the performance of PLA and PHBV blends[END_REF][START_REF] Phuong | Compatibilization and property enhancement of poly(lactic acid)/polycarbonate blends through triacetin-mediated interchange reactions in the melt[END_REF]. In the first heating run, characterizing the thermal properties of the films 'as prepared', a slight decrease of the glass transition temperature of PLA from 58 °C for the neat PLA film to 56°C for the PLA/PHBV film and 54 °C for the PLA/PHBV/oLA 89/10/1 is observed. This gradual decrease of the T g of PLA phase can be understood as an increase of compatibility of PLA/PHBV blends due to the presence of the oLA, miscible with both PLA and PHBV. The cold crystallization temperature of PLA is decreased from 119 to 105°C with the addition of PHBV, as already observed by previous authors [START_REF] Ferreira | Films of PLLA/PHBV: Thermal, morphological, and mechanical characterization[END_REF][START_REF] Nanda | The effects of process engineering on the performance of PLA and PHBV blends[END_REF][START_REF] Tri | Crystallization behavior of poly(lactide)/poly(β-hydroxybutyrate)/talc composites[END_REF], and is not influenced by oLA addition. The PLA degree of crystallinity also remains extremely low in every case (<5%). When blended with high molar mass PLA, PHBV loses its crystallinity, with a sharp decrease from 85 to 51%, as confirmed by Ma et al. [START_REF] Ma | Toughening of poly (lactic acid) by poly (β-hydroxybutyrateco-β-hydroxyvalerate) with high β-hydroxyvalerate content[END_REF]. But the degree of crystallinity of PHBV increases when an optimum value of oLA is added (from 51% without any oLA to 74 with 0.5 wt% oLA). This follows the same trend as what has been observed for PHBV/oLA blends, with a maximum of crystallinity for 0.5 wt% oLA and can be seen as a plasticizing effect of the oLA on the dispersed PHBV phase. This effect may then hinder the crystallization at higher oLA amount (1 wt%). Such effect is also observed in the second heating run with an increase from 18 to 31% in the PHBV degree of crystallinity with 0.5 wt% of oLA compared to pure PHBV. Finally, it can be noted that PHBV and oLA seems to have an impact on the crystalline forms of PLA. PHBV almost prevents the appearance of the crystalline form that melts at 150 °C while the addition of oLA has the opposite effect. With 1% of oLA, the 150 °C form is even predominant and in higher concentration than in PLA alone.

Morphology

Figure 4 shows the AFM micrographs of the fracture surface of PLA/PHBV (Figure 4a) and PLA/PHBV/ oLA blends (Figure 4b). AFM images of PLA/PHBV films underlined that the dispersed PHBV phase had relatively small average diameter with a typical nodular morphology dispersed in the PLA polymer matrix. With the addition of 1 wt% of oLA, the size of the PHBV domains is reduced but the shape remains unchanged. The mean area goes from 0.14 µm 2 for the PLA/PHBV film to 0.11 µm 2 for the PLA/ PHBV/oLA 89/10/1 film (Table 4) while the ratio length/width is around 2 in both cases. To support this observation and explain the standard deviation of the domains area measured for PLA/PHBV/oLA 1 wt%, the global distribution of the nodules sizes of PHBV is presented in Figure 5. If the PLA/PHBV film presents a broad single distribution, the additivated film shows a double distribution, with a sharper main distribution shifted to smaller size and the presence of a small portion of remaining large domains of PHBV (>0.4 µm 2 ) that increase the average size of the nodule observed. However a clear tendency to the decrease of PHBV size domain is observed: excluding the few remaining large domains, the average area of the compatibilized nodules is less than 0.09 µm 2 (while it is 0.13 µm 2 for the PLA/PHBV nodules), confirming the compatibilizing effect of oLA.

Mechanical properties

The results of mechanical analysis for PLA/PHBV and PLA/PHBV/oLA with different amount of oLA films are summarized in Table 5 and representative stress strain curves are reported in Figure 6.

The addition of small amounts of oLA reveal a slight improvement in the Young's modulus (E Young ) over PLA/PHBV (around 10%) while no effect is observed for single PLA, which may be linked to the increase in the crystallinity of the PHBV nodules. Yield strength however, remains similar for all the samples tested. The elongation at break (ε b ) of PLA/PHBV films strongly increases with the addition of PLA short chains (Table 3), showing an important enhancement in blends ductility. Films with 1 wt% oLA shows the highest elongation at break. This improved elongation at break can be interpreted by two ways: The first explanation is that the oLA acts as a compatibilizer so it creates an interphase between the PLA and PHBV facilitating stress transmission at the interface. The second explanation is based on the fact that the oLA acts as a plasticizer. A nonsignificant evolution in elongation at break from 3.1% for the neat PLA to 2.8% for PLA/oLA (99/1) blend is observed. This result shows that the plasticizing effect is negligible for the PLA phase and the compatibilizing effect of the oLA led to the improvement of the elongation at break for PLA/PHBV/oLA films. Tan δ evolutions obtained by DMTA are shown in Figure 7 for neat PLA, PLA/PHBV and PLA/PHBV/ oLA films. Neat PLA has a T g equal to 61 °C (in agreement with DSC measurements) and, as expected, all the PLA/PHBV/oLA films display a lower glass transition temperature compared to neat PLA. As already observed in DSC results, the presence of oLA chains in the PLA/PHBV blend does have a slight effect on the glass transition temperature of PLA which can be attributed to an increased compatibility of PLA and PHBV due to the oLA chains. Tan δ also shows a slight increase of T g of PHBV that can be seen from -2 °C for PLA/PHBV films to a value ranged between 2 and 5 °C with oLA addition which confirms the enhancement of compatibility. This shift attributed to a compatibilizing effect of oLA is consistent with what has been observed for the mechanical properties.

Gas barrier properties

Results of helium permeability of the films are presented in Table 6. The helium permeability shows a small decrease (10%) for PLA/PHBV film compared to neat PLA. This can simply be attributed to better barrier properties of PHBV compared to PLA alone. According to the morphology, it is found that PHBV is dispersed in the form of small nodules in the PLA matrix related to an improvement of permeability because PLA crystallinity degrees is always low with or without oLA.

Increasing the amount of oLA in PLA/PHBV film results in a small reduction of the permeability, with a minimum (32% compared to neat PLA) obtained for the film made with 1% of oLA (see Table 6). Indeed, as the oLA amount increases, the PHBV degree of crystallinity increases which counterbalances the decrease in the size of the nodules described above, to impact beneficially the gas barrier properties even if oLA acts also as a plasticizer.

Conclusions

The effect of a small amount of oLA on the morphology, thermomechanical and gas barrier properties of PLA/PHBV extruded films (90-10 wt%) was studied.

The miscibility of those oligomers with both polymers improved the compatibility of PLA with PHBV, which was confirmed by a drop in the interfacial tension from 0.5 to 0.3 mN/m with only 1 wt% of oLA in the PHBV dispersed phase. This also resulted in a decrease in the size of PHBV domains and a shift of both T g s in the blends. An unexpected plastification of the PHBV dispersed phase was also confirmed by an increase in the degree of crystallinity of PHBV. These combined effects resulted in films displaying much higher elongation at break (from 3 to 14% with 1 wt% of oLA) along with a slight increase in the Young's modulus. Moreover, the increase in cristallinity of PHBV also leads to a decrease of 30% of the permeability of the films.

Figure 1 .

 1 Figure 1. Chemical structure of PLA and PHBV

Figure 2 .Table 2 .Figure 3 .

 223 Figure 2. PHBV drop relaxation in PLA matrix, a) 0 s, b) 8 s, c) 16 s, d) 24 s, e) 30 s, f) 60 s

Figure 4 .

 4 Figure 4. AFM micrographs of PLA/PHBV (90/10) (a) and PLA/PHBV/oLA (89/10/1) (b)

Figure 5 .

 5 Figure 5. PHBV nodule size repartition in PLA/PHBV and PLA/PHBV/oLA 1 wt% blends. Shapes of the distributions are eyeguides only

Figure 6 .Figure 7 .

 67 Figure 6. Traction curves for the PLA/PHBV blends as a function of the amount of oLA [wt%]

  18 (m3 •m)/(m 2 •s•Pa)] PLA 136.5±5.9 PLA/PHBV 125.2±0.4 PLA/PHBV/oLA 0.1% 112.3±16.2 PLA/PHBV/oLA 0.5% 101.6±2.3 PLA/PHBV/oLA 1% 92.5±12.8

Table 1 .

 1 Thermal analysis properties of PLA, PHBV, PLA/oLA and PHBV/oLA blends (all obtained by solvent cast) as a function of the amount of oLA [wt%] (2 nd heating run)

		T g	T m PLA	T m PHBV	T cc	∆H cc	∆H m	X	c (PHBV)	X	c (PLA)
		[°C]	[°C]	[°C]	[°C]	[J/g]	[J/g]		[%]		[%]
	PLA (solvent cast)	60	151	-	121	22.6	23.70		-		1.0
	PLA (pellet)	61	153	-	-	5.0	6.10		-		1.0
	PHBV (solvent cast)	-	-	172	-	-	100.80		92.5		-
	PHBV (pellet)	-	-	173		-	92.50		85.0		-
	oLA (pellet)	43	143	-	94	42.1	59.60		-		18.0
	PHBV/oLA 99/1	-	-	172	-	-	78.60		73.0		-
	PHBV/oLA 95/5	-	-	167	-	-	86.40		83.0		-
	PHBV/oLA 90/10	-	-	169	-	-	76.80		70.0		-
	PLA/oLA 99/1	60	150	-	116	23.2	28.20		-		1.5
	PLA/oLA 95/5	59	150	-	115	26.6	26.95		-		0.2

Table 3 .

 3 Thermal analysis of PLA, PLA/PHBV and PLA/PHBV/oLA films as a function of the amount of oLA[wt%]

			T g	T m1	T m2	T m3	T cc	ΔH m3 (PHBV)	ΔH cc	X	c (PLA)	X	c (PHBV)
			[°C]	[°C]	[°C]	[°C]	[°C]	[J/g]	[J/g]		[%]	[%]
	heating run 2 nd	PLA PLA/PHBV PLA/PHBV/oLA 0.1% PLA/PHBV/oLA 0.5% PLA/PHBV/oLA 1%	60 57 57 57 58	151 150 150 151 151	157 156 156 156 -	-169 167 167 168	118 111 114 117 117	-1.96 1.76 3.38 1.47	34 29 29 31 31		0.1 3.2 1.4 0.6 0.0	-18 16 31 13
	heating run 1 st	PLA PLA/PHBV PLA/PHBV/oLA 0.1% PLA/PHBV/oLA 0.5% PLA/PHBV/oLA 1%	58 56 55 55 54	152 147 147 148 148	157 156 156 156 156	-170 170 170 171	119 105 104 105 105	-5.53 5.54 8.08 5.92	34 29 29 28 26		1.1 2.7 1.5 1.2 0.6	-51 51 74 54

Table 4 .

 4 PHVB and PHBV/oLA nodules size in PLA film

	Nodules

Mean area [µm 2 ] Mean width [µm] Mean length [µm] PLA

/PHBV 0.14±0.09 0.27±0.09 0.59±0.06 PLA/PHBV/oLA 1 wt% 0.11±0.10 0.26±0.12 0.50±0.

[START_REF] Yang | Transesterification induced mechanical properties enhancement of PLLA/PHBV bio-alloy[END_REF]

Table 5 .

 5 Results from tensile test for PLA, PLA/oLA and PLA/PHBV films as a function of the amount of oLA[wt%]

		Mechanical properties
	Samples	σ y	ε b	E Young
		[MPa]	[%]	[MPa]
	PLA	45±4	3.1±1.0	2050±481
	PLA/oLA 1%	48±7	2.8±0.5	2107±210
	PHBV *	39	2.0	3500-2800
	PLA/PHBV	43±5	3.0±1.0	2507±252
	PLA/PHBV/oLA 0.1%	47±1	8.0±2.0	2747±452
	PLA/PHBV/oLA 0.5%	45±4	6.0±1.0	2709±382
	PLA/PHBV/oLA 1%	44±2 14.0±3.0	2799±133

* http://www.tiananenmat.com/pdf/TDS_Y1000P_Dec2011.pdf. 2011

Table 6 .

 6 Helium film permeability as a function of the amount of oLA[wt%]

Acknowledgements

The authors wish to thank M. Boujelben and I. Charfeddine who helped for the extrusion and the DSC measurements.