Mohammad Saied Aghighi
email: ms.aghighi@gmail.com

Amine Ammar
email: amine.ammar@ensam.eu

Christel Metivier
email: christel.metivier@univ-lorraine.fr

M Gharagozlu

Rayleigh-Bénard convection of Casson fluids

Keywords:

come

Introduction

Viscoplastic fluids, characterized by a yield stress τ y , are known to present a complex transition between a solid-like and a fluid-like behavior. If the material is not stressed enough, i.e. less than the yield stress, the material does not flow and behaves as a solid. Above the yield stress, it flows with a shear-thinning behavior. The most used models describing viscoplastic fluids are the inelastic Bingham [START_REF] Bingham | The behavior of plastic materials[END_REF], Herschel-Bulkley [START_REF] Herschel | Measurement of consistency as applied to rubber-benzene solutions[END_REF] and Casson [START_REF] Casson | Rheology of disperse systems[END_REF] models. The Bingham model is the simplest model which deals with a yield stress. Due to its relative simplicity, the Bingham model is the most frequently used model in theoretical and numerical studies. However, the Bingham model is not realistic. Most of real materials behave as shear-thinning yield stress fluids such as the Herschel-Bulkley or Casson fluids. The Casson model is currently used in the field of food stuff [START_REF] Steffe | Rheological methods in food process engineering[END_REF], in particular, it has been chosen by the International Office of Cocoa and Chocolate to model the chocolate rheological behavior. Furthermore, the Casson model is often adopted in medicine to fit the blood rheology [START_REF] Scott Blair | The success of Casson's equation[END_REF][START_REF] Iida | Influence of plasma layer on steady blood flow in microvessels[END_REF][START_REF] Chaturani | Pulsatile flow of Casson's fluid through stenosed arteries with applications to blood flow[END_REF][START_REF] Dash | Estimation of increased flow resistance in a narrow catheterized artery-a theoretical model[END_REF][START_REF] Chakravarty | Unsteady flow of a two-layer blood stream past a tapered flexible artery under stenotic conditions[END_REF][START_REF] Sankar | Two-fluid Casson model for pulsatile blood flow through stenosed arteries: a theoretical model[END_REF]. The flow of some particulate suspensions can also be modeled by Casson fluids [START_REF] Macosko | Rheology : principles, measurements, and applications[END_REF], [START_REF] Chhabra | Non-Newtonian flow and applied rheology[END_REF].

The Rayleigh-Bénard Convection (RBC) is a buoyancy-driven instability, i.e. it occurs in a fluid layer heated from below and submitted to a vertical temperature gradient. This configuration is widely encountered in nature and also in industrial processes, justifying the impressive volume of work devoted to its understanding since more than a century especially in the case of Newtonian fluids, reviews are provided by Refs. [START_REF] Chandrasekhar | Hydrodynamic and hydromagnetic stability[END_REF][START_REF] Koschmieder | Bénard cells and taylor vortices[END_REF][START_REF] Bodenschatz | Recent developments in Rayleigh-bénard convection[END_REF] and some recent studies are done by Refs. [START_REF] Park | Rayleigh-Bénard convection of nanofluids based on the pseudo-singlephase continuum mode[END_REF][START_REF] Aghighi | Parametric solution of the Rayleigh-Benard convection model by using the PGD Application to nanofluids[END_REF][START_REF] Sanmiguel | On the onset of horizontal convection[END_REF][START_REF] Yada | Different routes to chaos in low Prandtl-number Rayleigh-Bénard convection[END_REF][START_REF] Kaur | Heat transfer in thermally modulated two-dimensional Rayleigh Bénard convection[END_REF]. Compared to the Newtonian case, the RBC in viscoplastic fluids has been little studied. However, the RBC in viscoplastic fluids is experiencing a rising interest these last decades leading to theoretical and numerical studies [START_REF] Zhang | Yield stress effects on Rayleigh-Bénard convection[END_REF][START_REF] Vikhansky | Thermal convection of a viscoplastic liquid with high Rayleigh and Bingham numbers[END_REF][START_REF] Turan | Laminar Rayleigh-Bénard convection of yield stress fluids in a square enclosure[END_REF][START_REF] Yigit | Effects of aspect ratio on natural convection of Bingham fluids in rectangular enclosures with differentially heated horizontal walls heated from below[END_REF][START_REF] Hassan | Rayleigh-Benard convection in Herschel-Bulkley fluid[END_REF][START_REF] Aghighi | Aspect ratio effects in Rayleigh-Bénard convection of Herschel-Bulkley fluids[END_REF][START_REF] Li | Natural convection in shear-thinning yield stress fluids in a square enclosure[END_REF][START_REF] Yigit | Numerical investigation of laminar Rayleigh-Bénard convection of Bingham fluids in square cross-sectioned cylindrical enclosures[END_REF], as well as experimental investigations [START_REF] Darbouli | Rayleigh-Bénard convection for viscoplastic fluids[END_REF], [START_REF] Kebiche | Experimental investigation of the Rayleigh-Bénard convection in a yield stress fluid[END_REF]. By considering either Bingham or Herschel-Bulkley fluids, all these studies show that the yield stress has a stabilizing effect, meaning that the conductive solution (motionless solution) remains stable for larger Rayleigh values when τ y is increased. It is also shown (i) that the motionless solution is linearly stable in the case of an inelastic viscoplastic fluid [START_REF] Zhang | Yield stress effects on Rayleigh-Bénard convection[END_REF] [START_REF] Vikhansky | Thermal convection of a viscoplastic liquid with high Rayleigh and Bingham numbers[END_REF], (ii) that convection can occur if the amplitude perturbation is large enough and (iii) that the transition from conductive to convective regime is subcritical [START_REF] Li | Natural convection in shear-thinning yield stress fluids in a square enclosure[END_REF], [START_REF] Vikhansky | On the stopping of thermal convection in viscoplastic liquid[END_REF]. Another result of great interest is the structure of the convective flow, via the velocity and temperature profiles, and the topology of the flow, i.e. the existence, size and location of the yielded and unyielded regions because it can have consequences on heat transfer for instance. Indeed, in the case of Bingham and Herschel-Bulkley fluids, the increase in τ y involves a decrease in the Nusselt number as well as the development and the increase in unyielded regions [START_REF] Turan | Laminar Rayleigh-Bénard convection of yield stress fluids in a square enclosure[END_REF] [26], [START_REF] Li | Natural convection in shear-thinning yield stress fluids in a square enclosure[END_REF]. Concerning the experimental investigations, the first to study the RBC in yield stress fluids is Darbouli et al. [START_REF] Darbouli | Rayleigh-Bénard convection for viscoplastic fluids[END_REF], who focus on the onset of convection in Carbopol gels in a cylindrical cavity. Onset of the Rayleigh-Bénard convection in Carbopol gels has also been studied experimentally in a rectangular cavity by Kebiche et al. [START_REF] Kebiche | Experimental investigation of the Rayleigh-Bénard convection in a yield stress fluid[END_REF]. In these studies, it is shown that the governing parameter is the yield number Y which represents the ratio between the yield stress and the buoyant effects. However, the studies do not obtain similar critical conditions in terms of Y c at the onset of convection, this is certainly due to different boundary conditions. A more recent study given by Métivier et al. [START_REF] Métivier | Origin of the onset of Rayleigh-Bénard convection in a concentrated suspension of microgels with a yield stress behavior[END_REF] highlight similar critical conditions than that given by Ref. [START_REF] Darbouli | Rayleigh-Bénard convection for viscoplastic fluids[END_REF]. Furthermore, by considering the development of thermal plumes due to a localized heat source in Carbopol gels, Davaille et al. [START_REF] Davaille | Thermal instabilities in a yield stress fluid: existence and morphology[END_REF] obtain the same order of Y c than [START_REF] Darbouli | Rayleigh-Bénard convection for viscoplastic fluids[END_REF] and [START_REF] Métivier | Origin of the onset of Rayleigh-Bénard convection in a concentrated suspension of microgels with a yield stress behavior[END_REF].

The aim of our paper is to consider numerically the RBC in Casson fluids since, to our best knowledge, it does not exist any result while the Casson model is widely used to describe viscoplastic fluids. For this respect, we extend the numerical study given by Turan et al. [START_REF] Turan | Laminar Rayleigh-Bénard convection of yield stress fluids in a square enclosure[END_REF] and Aghighi & Amar [START_REF] Aghighi | Aspect ratio effects in Rayleigh-Bénard convection of Herschel-Bulkley fluids[END_REF] to the Casson fluid case and we consider the bench case of a two-dimensional square cavity heated from below.

First, this paper set the mathematical formulation of our configuration. The numerical method and its validation are then described. In section 4, results are presented and discussed. Correlations for the Nusselt number are also proposed in this section. Our paper finally ends with a concluded section.

Mathematical formulation

Casson model

The stress-deformation behavior of viscoplastic materials based on Casson model is given by:

⎜ ⎟ = ⎛ ⎝ ⎜ + ⎛ ⎝ ⎞ ⎠ ⎞ ⎠ ⎟ > = < τ Bn γ γ for τ τ γ f o rτ τ 1 ˙ȧnd ˙0 ij ij y y 1 2 2 (1)
For = i j , 1,2 with = x x x y ( , ) ( , ) 1 2
are defined in Fig. 2 (2)

The Prandtl number, Pr, and Rayleigh number, Ra, are defined as follows:

= = Pr Ra μC k gβ TH αν Δ p 3 (3)
The various dimensional quantities above are defined as follows: τ y is the yield stress, μ is the constant plastic viscosity, C p is specific heat capacity, k is the thermal conductivity, g is the acceleration due to gravity, β is the coefficient of thermal expansion, α is the thermal diffusivity, T Δ is the temperature difference between hot and cool walls and ν is the kinematic viscosity.

γ˙and τ are the second invariant of the respective rate-of-strain and shear stress tensors

= γ γγ ˙˙i j ij 1 2 and = τ ττ ij ij 1 2
,Where γ˙i j is defined by:

= ∂ ∂ + ∂ ∂ γ u x u x ˙ij i j j i (4)
Papanastasiou modifications method [START_REF] Papanastasiou | Flows of materials with yield[END_REF] can be applied to Casson model to avoid discontinuity between yielded and unyielded regions. This method consists in an exponential regularization of the viscosity, leading to the expression of the stress [START_REF] Mitsoulis | Flows of viscoplastic materials: models and computations[END_REF]:

⎜ ⎟ = ⎛ ⎝ ⎜ ⎜ + ⎛ ⎝ ⎞ ⎠ ⎛ ⎝ ⎜ - ⎛ ⎝ ⎜ - ⎞ ⎠ ⎟ ⎞ ⎠ ⎟ ⎟ τ Bn γ mγ 1 ˙1 exp ˙γ ij ij 1 2 2 ( 5
)
where m is the regularization parameter which is a large value. It is worth noting that, the effective viscosity, = μ e τ γ ˙ij ij , tends to mBn when γṫ ends to zero (Fig. 1). This is a large but finite value which depends on the Papanastasiou parameter and the Bingham number. On the other hand, viscosity μ e tends to 1 when γ˙tends to infinity. Furthermore, one notices that at fixed value of m, the decrease in Bn involves a decrease in the effective viscosity (whatever the value of γ˙) but also a decrease in the shear-thinning degree which could be defined by dμ dγ / ė .

Rayleigh-Bénard convection

The configuration under analysis is depicted in Fig. 2. The model is defined in the unit square domain fulfilled by a yield stress fluid modeled by the Casson law. Because one knows that the transient phenomena do not modify the steady state result (see for instance [START_REF] Turan | Laminar Rayleigh-Bénard convection of yield stress fluids in a square enclosure[END_REF] [24] [START_REF] Ibrahim | Magnetohydrodynamic stagnation point flow and heat transfer of casson nanofluid past a stretching sheet with slip and convective boundary condition[END_REF]), in this study we consider the steady state set of equations corresponding to the two-dimensional natural convection case. Assuming the Boussinesq approximation and using the characteristic Where u v θ , , and p are dimensionless horizontal velocity, vertical velocity, temperature and pressure respectively.

+ = + =-+ + + =-+ + + + = + ∂ ∂ ∂ ∂ ∂ ∂ ∂ ∂ ∂ ∂ - ∂ ∂ ∂ ∂ ∂ ∂ ∂ ∂ ∂ ∂ - ∂ ∂ ∂ ∂ ∂ ∂ ∂ ∂ - ∂ ∂ ∂ ∂ ( ) ( ) ( ) u v
No-slip conditions are considered at walls. = = u v 0 , Adiabatic conditions are written at vertical walls and imposed temperatures at horizontal walls:

∂ ∂ = = = θ x at x and x 0 0 1 (7) = = θ a ty 0.5 0 = - = θ a ty 0.5 1
In these equations, the dimensionless temperature θ is defined by:

= - - θ T T T T r H C ( 8
)
Where

Numerical method

In order to solve the system of equation ( 6) with their corresponding boundary conditions [START_REF] Chaturani | Pulsatile flow of Casson's fluid through stenosed arteries with applications to blood flow[END_REF], by using Matlab we have developed a numerical code that is presented in this section. Similarly to Aghighi et al. [START_REF] Aghighi | Non-incremental transient solution of the Rayleigh-Bénard convection model by using the PGD[END_REF] in what follows we consider a penalty formulation of the incompressibility constraint. In this way, equation ( 6) can be written as follows: Where λ is a sufficiently large constant which allows to ensure the in- compressibility condition.

+ - = + + = + + + -= + + - + = ∂ ∂ ∂ ∂ ∂ ∂ ∂ ∂ ∂ ∂ - ∂ ∂ ∂ ∂ ∂ ∂ ∂ ∂ ∂ ∂ - ∂ ∂ ∂ ∂ ∂ ∂ ∂ ∂ - ∂ ∂ ∂ ∂ ( ) ( ) ( )
The matrix form of the system of equation [START_REF] Sankar | Two-fluid Casson model for pulsatile blood flow through stenosed arteries: a theoretical model[END_REF] writes:

+ = L T f [ ]{ } { } 0 ( 11
)
where

(13) = ⎡ ⎣ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ - - ⎤ ⎦ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ∂ ∂ ∂ ∂ ∂ ∂ ∂ ∂ L l l l [ ] 0 0 0 1 0 0 0 0 x y x y λ 1 1 1 2 (14) = ∂ ∂ + ∂ ∂ l u x v y 1 ( 15
) ⎜ ⎟ = ∂ ∂ + ∂ ∂ - ⎛ ⎝ ∂ ∂ + ∂ ∂ ⎞ ⎠ - l u x v y Ra Pr x y ( . ) 2 1 2 2 2 2 2
The numerical solutions of system [START_REF] Macosko | Rheology : principles, measurements, and applications[END_REF] with the corresponding boundary conditions are obtained by developing a numerical code based on the Galerkin weighted residual finite element method with quadrilateral, eight nodes elements. The Galerkin weak statement is obtained by minimizing the residual error over the discretized domain [START_REF] Dhatt | Finite element method[END_REF].

Space discretization

Using the set of shape function

= N { } k k n
1 the velocity components u v ( , ), pressure p and temperature θ can be expanded in the domain A x y ( , ) for ≤ ≤ x y (0 , 1) as bellow:

≈ ∑ ≈ ∑ ≈ ∑ ≈ ∑ = = = = u uN x y v vN x y p pN x y θ θN x y ( , ), ( , ), ( , ), ( , )
k n k k k n k k k n k k k n k k 1 1 1 1 ( 16
)
where n is the total number of nodes.

Using the quadrilateral, eight nodes elements to discrete the domain and isoparametric method for mapping from x y ( , ) to a square reference domain ξ η ( , ) one can write:

∑ ∑ = = = = x xN ξ η y yN ξ η ( , ), ( ,
) i i i i i i 1 8 1 8 (17)

Weak formulation

The integral form of equation ( 11) can be written as:

∫ = 〈 〉 + = W ψ L T f dA {[ ]{ } { }} 0 A ( 18
)
where ψ is a weighting(test) function.

〈 〉 = 〈 〉 ψ ψ ψ ψ ψ u v p θ ( 19
)
Based on the Galerkin weighted residual finite element method the weighting functions ψ comprise all the variations δT of the functions T:

= =〈 〉 δ N δ T ψ T { } (20)
Substituting Eq. ( 19) into Eq. ( 16) leads to:

∫ = + = W δ fd A T T ( ( ) ) 0 A ( 21
) So L ∫ = 〈 〉 〈 〉 + = W δT N N T f dA { }( ( { }) ) 0 A ( 22
)
The expansion of Eq. ( 22) given W W W and W , ,

i u i v i p i θ is detailed in Appendix A.
The integrations are evaluate numerically by using a direct integration formula on a square [START_REF] Mitsoulis | Flows of viscoplastic materials: models and computations[END_REF] in ξ η ( , ) domain using the fol- lowing relationships:

⎧ ⎨ ⎩ ⎫ ⎬ ⎭ = ⎧ ⎨ ⎩ ⎫ ⎬ ⎭ = ∂ ∂ ∂ ∂ - ∂ ∂ ∂ ∂ J dxdy J dξdη [ ] N x N y N ξ N η 1 i i i i ( 23
)
where J is the Jacobian matrix defined by: After applying the numerical integration, a system of algebraic equations can be written based on elementary matrices as follows:

= ⎡ ⎣ ⎢ ⎢ ⎤ ⎦ ⎥ ⎥ ∂ ∂ ∂ ∂ ∂ ∂ ∂ ∂ J [ ] x ξ y ξ x η y η (24)
= K T F [ ]{ } { } (25)
Where K [ ] is the global left hand side matrix, F { } is the global load vector and T { } is the global unknown vector:

= ⎡ ⎣ ⎢ ⎢ ⎢ ⎢ ⎢ - - - - - ⎤ ⎦ ⎥ ⎥ ⎥ ⎥ ⎥ - K k wk k k wk k k k k k Ra Pr k [ ] 0 0 0 0 0 0 0 ( . ) k λ 1 2 3 1 2 4 5 3 4 1 2 5 1 2 (26)
where

= + = -∇ ∇ = = = = ⎡ ⎣ ⎢ ⎤ ⎦ ⎥ ∂ ∂ ∂ ∂ ∂ ∂ ∂ ∂ k u N v N k N H N k N k N k N N H { } { } ({ }[ ] ) { } { } { } [ ] 1 0 0 1 p N x p N y N x N y 1 2 3 4 5 (27)
and

= ⎧ ⎨ ⎩ ⎫ ⎬ ⎭ = - ⎧ ⎨ ⎪ ⎩ ⎪ ⎫ ⎬ ⎪ ⎭ ⎪ T u v p θ F f f { } ; { } 0 0 nx ny (28) = - ⎡ ⎣ ⎤ ⎦ + ⎡ ⎣ ⎤ ⎦ + = - ⎡ ⎣ ⎤ ⎦ + ⎡ ⎣ ⎤ ⎦ + - ∂ ∂ ∂ ∂ - ∂ ∂ ∂ ∂ ( ) ( ) f P r R a N τ N τ w k u f P r R a N τ N τ w k v { } { } { } { } ([ ] { }) { } { } { } { } ([ ] { }) nx N x xx N y yx p ny N y yy N x xy p 1 2 1 2 2 1 2 1 2 2 (29)
Here u p and v p are the previous values of velocity and w is a con- stant. At the end of each iteration, the previous values of velocity are replaced by the new computed values:

= = u u and v v p p
. This procedure is used in order to linearize the non-linear terms which come from the advective and convective terms in the momentum and energy equations. The linear system of equation ( 25) is solved by using the Matlab function 'mldivide'. Using the new computed values, the stress terms are evaluated and the process is repeated until the convergence is reached. The convergence is obtained when the norm 2 of the relative difference between the new and the previous values of the vector T { } is less than - 10 5 . The velocity distribution of Newtonian fluids is considered as the initial values of u p and v p . It is worth mentioning that the numerical analysis of Eq. ( 10) leads to an unstable system of equations. In order to circumvent this difficulty, the Laplacian velocity components were added in Eqs. ( 26) and [START_REF] Darbouli | Rayleigh-Bénard convection for viscoplastic fluids[END_REF]. Because these terms are added on both sides of the momentum equations, they vanish by considering the whole system. It means that these terms do not affect the results but allow to make the equations stable to introduce diffusive terms which make the equation stable.

Numerical method validation

The grid independent solution has been tested by considering two cases = = = Ra Pr Bn ( 10 , 10, 0) We found that the mesh consisting in 1825 nodes guarantees a grid independent solution within the relative tolerance level of - 10 5 . In order to ensure convergence of results in the case of viscoplastic fluid, the convergence of the solutions was checked by varying the penalty and regularization parameters. Results show that Nu converges within 0.1% by varying λ from 10 3 to 10 4 and also by varying m from 10 3 to 10 4 . In the following, all results are obtained for = λ 10 4 and = m 10 4 . The numerical method has been validated by considering the Newtonian case = Bn ( 0). Our results have been compared with that ones given in Refs. [START_REF] Aghighi | Non-incremental transient solution of the Rayleigh-Bénard convection model by using the PGD[END_REF], [START_REF] Ouertatani | Numerical simulation of twodimensional Rayleigh-Bénard convection in an enclosure[END_REF]. A very good agreement has been obtained within a mean relative difference of 1%.

The numerical code was also validated by testing the Bingham case. Results of mean Nusselt number Nu for = * Ra 5 10 4 and = Pr 10 are shown in Fig. 3, one observes a good agreement between our results (continuous line) and the results obtained by Turan et al. [START_REF] Turan | Laminar Rayleigh-Bénard convection of yield stress fluids in a square enclosure[END_REF] (circles). Bn at fixed Pr. It highlights a stabilizing effect of the yield stress via Bn. Similarly, the increase in Pr stabilizes the motionless solution. In order to compare our results to experimental studies given by Refs. [29] [30], and [START_REF] Métivier | Origin of the onset of Rayleigh-Bénard convection in a concentrated suspension of microgels with a yield stress behavior[END_REF], Fig. 5 worth noting that criticality is sensitive to boundary conditions and the value of aspect ratio. In Refs. [START_REF] Darbouli | Rayleigh-Bénard convection for viscoplastic fluids[END_REF] and [START_REF] Métivier | Origin of the onset of Rayleigh-Bénard convection in a concentrated suspension of microgels with a yield stress behavior[END_REF], several values of aspect ratio are investigated between 6 and 18, while in the present study, the aspect ratio has been chosen equal to 1 corresponding to the benchmark case of the Rayleigh-Bénard configuration. When considering slip boundary conditions, the value of Y 1/ c decreases as highlighted by Ref. [START_REF] Darbouli | Rayleigh-Bénard convection for viscoplastic fluids[END_REF]. In Ref. [START_REF] Kebiche | Experimental investigation of the Rayleigh-Bénard convection in a yield stress fluid[END_REF], results provide that ≈ Y 1/ 5 c for slip conditions and insulated horizontal walls conditions which are quite far from the present study.

Velocity and temperature

When the convection occurs, the velocity vectors are no longer equal to zero as well as the linear profile of the temperature is modified by the motion. Figs. 6 and7 depict the variations of non-dimensional horizontal component of velocity u and non-dimensional vertical component of velocity v along the vertical and horizontal mid-planes of the cavity when = Pr 100 and for different values of Ra and Bn. Results show that the magnitude of velocity increase with increasing Ra due to stronger buoyancy force. In other words, the increase in Ra strengthens the intensity of convection. On the other hand, one can see that the magnitudes of velocity and temperature, displayed in Figs. 678, decrease with increasing Bingham number indicating the stabilizing effect of the yield stress, which corresponds to an additional flow resistance. This result agrees with those obtained in the Bingham and Herschel-Bulkley cases [START_REF] Turan | Laminar Rayleigh-Bénard convection of yield stress fluids in a square enclosure[END_REF] [26], [START_REF] Li | Natural convection in shear-thinning yield stress fluids in a square enclosure[END_REF]. In these articles, it is shown that the increase in the Bingham number (or Yield number) leads to stop the convective motion. It means that for fixed Ra and Pr values, there exists a value of Bn above which the convection does not occur, say Bn max . Figs. [START_REF] Iida | Influence of plasma layer on steady blood flow in microvessels[END_REF][START_REF] Chaturani | Pulsatile flow of Casson's fluid through stenosed arteries with applications to blood flow[END_REF][START_REF] Dash | Estimation of increased flow resistance in a narrow catheterized artery-a theoretical model[END_REF] show that the value of Bn max increases with increasing Ra due to the stronger convection force which can overcome the flow resistance up to greater values of the Bingham number.

Similarly, variations of non-dimensional temperature θ with Ra and Bn along the vertical mid-plane of the cavity have been shown in Fig. 8. As expected, the non-linearity of temperature distribution increases with increasing Ra due to stronger convection force. On the other hand, the variation of the temperature θ tend to a linear profile when the value of Bingham number increases and reaches Bn max , meaning that the heat transfer is only due to thermal conduction and there is no flow within the enclosure.

Contours of temperature and velocity vectors

Fig. 9 shows contours of non-dimensional temperature θ and velo- city vectors for different values of Bingham number at = Ra 10 5 and = Pr 100. The Shaded areas in these figures show the unyielded regions where < τ τ y . One can observe different types of regions: the Appar- ently Unyielded Regions (AUR) and the Truly Unyielded Regions (TUR) [START_REF] Mitsoulis | Flows of viscoplastic materials: models and computations[END_REF]. At the corners of the cavity, the unyielded regions correspond to AUR since they are dead zones where the velocities are very small and the fluid is stagnant. In the center of the cavity which coincides with the center of the convective cell, TUR appear when Bn is increased at nominal Ra and Pr values as one can see in Fig. 9. The unyielded re- gions occur where the shear rate is minimal (equal to zero). By contrast, maximal values of shear rate are located at walls where no slip conditions are considered. The increase in Bn is accompanied by an increase in the TUR area as well as a decrease in the convective intensity. Finally, unyielded regions invade the whole cavity at = Bn Bn max . Ob- viously, the increase in TUR area highlights the decrease in the convection intensity as already discussed in the paragraph 4.1. These results are in agreement with results proposed by Turan et al. [START_REF] Turan | Laminar Rayleigh-Bénard convection of yield stress fluids in a square enclosure[END_REF] for Bingham fluids and by Refs. [START_REF] Aghighi | Aspect ratio effects in Rayleigh-Bénard convection of Herschel-Bulkley fluids[END_REF] [START_REF] Li | Natural convection in shear-thinning yield stress fluids in a square enclosure[END_REF],for Herschel-Bulkley fluids. = Nu 1, meaning that convection stops leading to a mo- tionless state. The transition between the convective to conductive regimes is abrupt in terms of Nusselt values, highlighting a subcritical bifurcation which is in agreement with [START_REF] Vikhansky | On the stopping of thermal convection in viscoplastic liquid[END_REF] and [START_REF] Li | Natural convection in shear-thinning yield stress fluids in a square enclosure[END_REF]. As expected, at fixed value of Bn, the increase in the Rayleigh number involves larger heat transfer (Nu values) via larger buoyancy effects. It also means that Bn max increases with increasing Ra.

Concerning the influence of the Pr number, one notices that results are insensitive to Pr variations in the Newtonian case. However, for viscoplastic fluids ( ≠ Bn 0), the increase in Pr values decreases the heat transfer as well as the value of Bn max as displayed in Figs. 10 and11.

Our results obtained for Casson fluids can be compared with the results obtained for Bingham [START_REF] Turan | Laminar Rayleigh-Bénard convection of yield stress fluids in a square enclosure[END_REF] and the shear-thinning Herschel-Bulkley fluids [START_REF] Aghighi | Parametric solution of the Rayleigh-Benard convection model by using the PGD Application to nanofluids[END_REF]. At fixed Ra and Bn values, = Pr 10, one can observe that Casson fluids lead to smaller Nu values than the two other models which involves a larger stabilizing effect of the yield stress in the case of Casson fluids. In other words, the heat transfer is the smallest one for the Casson fluids while the largest heat transfer is obtained for Herschel-Bulkley model. This result is correlated with the non-linear (shear-thinning) variations of the viscosity [START_REF] Aghighi | Non-incremental transient solution of the Rayleigh-Bénard convection model by using the PGD[END_REF], [START_REF] Turan | Laminar natural convection of power-law fluids in a square enclosure with differentially heated side walls subjected to constant temperatures[END_REF].

Correlation for the mean Nusselt number

A correlation for the mean Nusselt number Nu can be proposed in a general form as below:

⎜ ⎟ = = + -⎡ ⎣ ⎢ -⎛ ⎝ ⎞ ⎠ ⎤ ⎦ ⎥ Nu f Ra Pr Bn Nu Bn Bn ( , , ) 1 ( 1
The correlation we propose for Nu 0 in Eq. ( 31) implies a maximal relative difference with our numerical results which is less than 2% and a mean relative difference smaller than 1% for < < Ra 5.10 10 3 5 and = Pr 10, 100, 1000 as one can observe in Figs. 10 and 11. These cor- relations are also very useful in the case of viscoplastic fluids, since they provide the heat transfer via Eq. ( 30) and also Bn max (Eq. ( 32)), the value above which the flow cannot occur. This last value is determined as a function of Ra and Pr , with a maximal relative difference with our numerical results of 3.2% and a mean relative difference close to 1%.

Conclusions

In this study, a finite element numerical code has been developed in order to investigate the Rayleigh-Bénard convection in viscoplastic fluids modeled by a Casson law. A two-dimensional square enclosure is considered. Results obtained for a Casson fluid are qualitatively similar to that obtained for Bingham and Herschel-Bulkley fluids. It means that the stabilizing effect of the yield stress via the Bingham number is highlighting since the increase in Bn reduces the convection intensity. From a threshold value, Bn max , the convective motion stops leading to a conductive state characterized by a linear temperature profile and = Nu 1. One recovers that the transition between the conductive to the convective state ( < Bn Bn max ) is abrupt corresponding to a subcritical bifurcation. Furthermore, we also obtain that the truly unyielded regions appear in the center of the convective cell where the shear rate is minimal. By contrast, maximal values of shear rate are reached at walls. The areas of unyielded regions increase with the increase in Bn, un- yielded regions invade the whole cavity from Bn max . Finally, the heat transfer is evaluated and correlations for the mean Nusselt number as well as Bn max are provided as a function of Ra Pr , and Bn. By com- paring with the Bingham and Herschell-Bulkley models, the Casson model leads to the smallest heat transfer. (A.4)

  the steady dimensionless continuity, momentum and energy equations are:

Fig. 1 .

 1 Fig. 1. Variations of viscosity μ with γ˙at the constant Bingham number (left) and constant regularization parameter (right).

  T C and T H ( > T T H C ) are the temperatures enforced at the top and bottom cavity boundaries, respectively, and T r is a reference tem- perature: flux averaged over the hot wall is defined via the Nusselt number:

Fig. 2 .

 2 Fig. 2. Schematic diagram of the physical model and co-ordinate system.

Fig. 3 .

 3 Fig. 3. Mean Nusselt number Nu for Bingham fluids at * 5 10 4 and = Pr 10 (line: present study, points: Turan study [23]).

Fig. 4 .

 4 Fig. 4. The variations of critical Rayleigh number Ra c with Bingham number Bn for = Pr 10, 100, 1000.

Fig. 5 .Fig. 6 .

 56 Fig. 5. The variations of critical Yield number Y c with Rayleigh number Ra for = Pr 100.

Fig. 7 .

 7 Fig. 7. Variations of non-dimensional velocity v of Casson fluids for different Bn values along the horizontal mid-plane for different values of Ra at = Pr 100.

Fig. 8 .

 8 Fig. 8. Variations of non-dimensional Temperature θ of Casson fluids for different Bn values along the vertical mid-plane for different values of Ra at = Pr 100.

5 .

 5 The convergence is tested on the values of the average Nusselt number.

Fig. 9 .

 9 Fig. 9. Contours of non-dimensional temperature θ and velocity vectors with plug regions of Casson fluids for six different values of Bingham number at = Ra 10 5 and = Pr 100.

Fig. 10 .

 10 Fig. 10. Evolution of the mean Nusselt number of Casson fluids with the Bingham number for three different values of Pr at = Ra 5.10 3 (left) and = Ra 10 4 (right). [Numerical results (-) and prediction of the correlation (o)given by Eq. (30)].

4 5 10 c 2 )

 4102 represents the onset of convection in terms of the critical Yield number Y c (where the Yield number is defined as = Y τ ρβg TH Δ y ) as a function of the Rayleigh number for = Pr 100. In the range of the Rayleigh values ), we obtain realistic values of Y c since in Refs. [29] and [32], the studies highlight ≈ Y 1/ 80 c (i.e. ≈ - Y for no-slip conditions and (perfect) conductive horizontal walls. Furthermore, it is

Fig. 11 .

 11 Fig. 11. Evolution of the mean Nusselt number of Casson fluids with the Bingham number for three different values of Pr at = Ra 5.10 4 (left) and = Ra 10 5 (right). [numerical results (-) and prediction of the correlation (o)given by Eq. (30)].

4. 3 . 5 .

 35 Heat transfer (Nusselt number) The variations of mean Nusselt number Nu with Bingham number Bn are shown in Figs. 10 and 11 for different values of Prandtl number ( = Pr 10, 100, 1000) It can be seen that the maximal value of the mean Nusselt number occurs for Newtonian fluids where = Bn 0. For visco- plastic fluids ( > Bn 0), the mean Nusselt number decreases with in- creasing Bingham number, underlying once again that the yield stress decreases the convective motion intensity. As already discussed, for a given Ra and Pr set of values, results show that the mean Nusselt number decreases with increasing Bingham number and for large value of Bingham number ≥ Bn Bn max the value of mean Nusselt number fall abruptly to


Céline Bellard

Göran Englund

Bernard Hugueny

Biotic and abiotic drivers of

Keywords: aquatic ecosystems, isolation, piscivores, age, fragmentation, Words: 6, 128

HAL is

Introduction

Today, we are facing the sixth wave of extinctions of our planet [START_REF] Ceballos | Biological annihilation via the ongoing sixth mass extinction signaled by vertebrate population losses and declines[END_REF][START_REF] De Vos | Estimating the normal background rate of species extinction[END_REF]. In most of deep time, extinctions occurred at a fairly constant rate, called background extinction rate. Today, human activities have raised extinction rates several orders of magnitude above this level. For example, the current rate of extinctions for vertebrates is more than 100 times the background rate, with literally ca. 4 vertebrates going extinct every year since 1900 [START_REF] Ceballos | Accelerated modern human-induced species losses: Entering the sixth mass extinction[END_REF]. Land use changes are responsible for most of the modern extinctions [START_REF] Collen | Global patterns of freshwater species diversity, threat and endemism[END_REF][START_REF] Leclerc | Insular threat associations within taxa worldwide[END_REF][START_REF] Pereira | Scenarios for global biodiversity in the 21st century[END_REF][START_REF] Sala | Global Biodiversity Scenarios for the Year 2100[END_REF]. Habitat destruction and altered land use have reduced previous continuous habitats to isolated fragments. Local and regional species extinctions may follow but not immediately, leading to decreasing species richness over time at several spatial scales. This means that the ecological cost of current habitat fragmentations could be delayed, a phenomenon termed extinction debt. Studies of extinction debts in a wide range of taxa and ecosystems [START_REF] Harding | Stream biodiversity: The ghost of land use past[END_REF][START_REF] Piessens | Does the heathland flora in north-western Belgium show an extinction debt?[END_REF][START_REF] Piha | Amphibian occurrence is influenced by current and historic landscape characteristics[END_REF] has shown that the time until extinctions take place may vary from weeks and months to hundreds or thousands years [START_REF] Shaw | Decadal changes and delayed avian species losses due to deforestation in the northern Neotropics[END_REF][START_REF] Uezu | Time-Lag in Responses of Birds to Atlantic Forest Fragmentation: Restoration Opportunity and Urgency[END_REF]. This gradual loss of species (also known as faunal relaxation to refer, in physics, to the return of perturbed system into equilibriumhas first studied by [START_REF] Diamond | Biogeographic kinetics: estimation of relaxation times for avifaunas of Southwest Pacific islands[END_REF] and have since been studied in many different systems [START_REF] Cousins | Detection of extinction debt depends on scale and specialisation[END_REF][START_REF] Hugueny | Habitat fragmentation and extinction rates within freshwater fish communities: a faunal relaxation approach[END_REF][START_REF] Hylander | The mechanisms causing extinction debts[END_REF][START_REF] Kuussaari | Extinction debt: a challenge for biodiversity conservation[END_REF]. This observation is particularly important for conservation, as the time from impact to extinction provides a time window within which remedial actions can be taken. Moreover, understanding mechanisms of delayed extinction is essential for predicting the impact of global changes on biodiversity. A range of intrinsic and extrinsic mechanisms may influence the timing and sequence of species losses, including the increased risk of stochastic extinction faced by small populations size and species have narrow habitat requirement (e.g., [START_REF] Almeida | Species abundance plays a dominant role in local and regional extinctions in freshwater communities and allows the identification of selective extinctions[END_REF]. At the regional scale, habitat modifications may result in disequilibrium meta-population dynamics (i.e., extinction and colonization rates are no more balanced), which also can lead to delayed extinctions [START_REF] Helm | Slow response of plant species richness to habitat loss and fragmentation[END_REF]. Species loss rates may be accelerated by inter-specific competition, predation, and anthropogenic impacts such as overexploitation and species invasions. Life history traits could also influence species' responses to extinctions [START_REF] Krauss | Habitat fragmentation causes immediate and time-delayed biodiversity loss at different trophic levels[END_REF], e.g., long lived species are expected to persist longer than short lived species [START_REF] Cronk | Plant extinctions take time[END_REF][START_REF] Lindborg | Historical landscape connectivity affects present plant species diversity[END_REF].

Recent studies emphasize that species loss rates in isolated habitat fragments can be accurately modelled with a few predictors, like age and area of the isolate, for a large array of terrestrial taxa and for different types of fragments (e.g. true islands, forest remnants, mountaintops) [START_REF] Halley | Dynamics of extinction debt across five taxonomic groups[END_REF][START_REF] Hugueny | Age-area scaling of extinction debt within isolated terrestrial vertebrate assemblages[END_REF]. The negative relationship between species loss rate and area is well documented in the literature [START_REF] Ferraz | Rates of species loss from Amazonian forest fragments[END_REF][START_REF] Halley | Neutral theory as a predictor of avifaunal extinctions after habitat loss[END_REF][START_REF] He | Species-area relationships always overestimate extinction rates from habitat loss[END_REF][START_REF] Hugueny | Age-area scaling of extinction debt within isolated terrestrial vertebrate assemblages[END_REF]. Some of these studies also demonstrate that species decay curves often differ from the exponential decay expected for equivalent and non-interacting species. Instead of the constant species loss rate per unit time expected under an exponential decay, empirical data show that species loss rates decrease with increasing age of the isolate. This age effect could arise because the loss of species leads to reduced risk of extinctions due to species interactions, i.e., diversity-dependent dynamics (e.g., [START_REF] Halley | Dynamics of extinction debt across five taxonomic groups[END_REF]). Alternatively, it can be the result of interspecific variability in extinction proneness. As species with high probability of extinctions disappear soon after isolation, the remaining community will consist of species with comparably low extinction risk [START_REF] Ferraz | Rates of species loss from Amazonian forest fragments[END_REF][START_REF] Hugueny | Age-area scaling of extinction debt within isolated terrestrial vertebrate assemblages[END_REF].

Most studies of delayed extinctions have focused on terrestrial assemblages, and little is known about natural extinction rates of aquatic organisms in the absence of catastrophic events. Fortunately, freshwater fish populations frequently occur in isolated habitats such as lakes [START_REF] Englund | Predation leads to assembly rules in fragmented fish communities[END_REF] and coastal rivers [START_REF] Hugueny | Habitat fragmentation and extinction rates within freshwater fish communities: a faunal relaxation approach[END_REF] that receive no or few colonists. Recently, [START_REF] Hugueny | Age-area scaling of extinction debt within isolated terrestrial vertebrate assemblages[END_REF] built a new theoretical model to test the effect of age and area on species loss rates in isolated terrestrial vertebrate communities. We believe the modeling framework built by [START_REF] Hugueny | Age-area scaling of extinction debt within isolated terrestrial vertebrate assemblages[END_REF] should prove useful in explaining species extinctions in aquatic systems and allow us to test for the contribution of aquatic variables to explain species loss rates. For instance, small and/or young (recently isolated) lakes and rivers should experience a higher annual rate of extinctions per species than larger and/or older ones.

A former study of fish communities in coastal rivers isolated since the beginning of the Holocene confirmed the important contribution of the area of the isolate to extinction rates, but as the rivers became isolated at approximatively the same time the effect of age could not be studied [START_REF] Hugueny | Habitat fragmentation and extinction rates within freshwater fish communities: a faunal relaxation approach[END_REF]. In fact, few fragmented systems in both terrestrial and aquatic environments allow studying species losses over long and continuous temporal scales.

For instance, the fragments analyzed in recent meta-analyses are either young (a few hundred years for the older) and the result of human driven fragmentation, or very old (about 10,000 years) resulting from environmental changes that occurred at the end of the Pleistocene, with no values in between [START_REF] Dias | Anthropogenic stressors and riverine fish extinctions[END_REF][START_REF] Halley | Dynamics of extinction debt across five taxonomic groups[END_REF][START_REF] Hugueny | Habitat fragmentation and extinction rates within freshwater fish communities: a faunal relaxation approach[END_REF].

Here, we will study extinction rates of freshwater fishes in a natural experiment -bays transformed into lakes by land uplift caused by post-glacial isostatic rebound in Northern Sweden. Because the uplift rate was approximately constant over the last 4000 years, we can estimate the age of isolated lakes. This process has created lakes of various ages covering 0 to 5,000 years that provides an opportunity to study the time from isolation to extinction.

Besides factors acting across taxon and types of fragment, we expect that specific mechanisms of lacustrine ecosystems are also important in shaping species decay through time. For instance, acidification is known to cause extinctions of fish species [START_REF] Holmgren | Challenges in assessing biological recovery from acidification in Swedish lakes[END_REF], and depth is known to mediate the effects of piscivorous fishes on their prey in this system [START_REF] Englund | Predation leads to assembly rules in fragmented fish communities[END_REF]. Indeed, shallow lakes are known to serve as a refuge for prey species, while deep lakes are particularly favorable to the presence of two dominant piscivores, Eurasian perch and northern pike. Yet, the relative importance of area, age, depth, pH, and piscivory for extinction patterns have not been unexplored in aquatic ecosystems (justification of used variables are provided in the material and method section). In this study, we adapt the modelling approach introduced by Hugueny (2017) to examine the effects of these factors on extinction rates.

Material and Methods

Study sites:

Lakes are classified as either (i) connected if they are linked to the Baltic Sea by a permanent stream with a gentle gradient or (ii) isolated if they lack permanent stream connections or if their outlet streams have steep sections that act as a dispersal barrier (see [START_REF] Englund | Predation leads to assembly rules in fragmented fish communities[END_REF]. The isolated lakes are not connected to upstream lakes with fish. Thus, all the isolated lakes are fully disconnected from other sources and no event of colonization has been recorded. The study sites also included closed bays that were in the process of forming a lake, but were still well connected to the Baltic Sea with the water surface being in level with the surface of the Baltic Sea. The land uplift caused by the post-ice-age isostatic recovery is thus transforming bays in the Baltic sea to lake environments. Because of the low salinity of the Baltic Sea, it is inhabited by numerous freshwater fish species that could successfully invade bays and maintain self-sustaining populations in newly isolated lakes.

The dataset includes 71 isolated lakes, 41 connected lakes, and 12 closed bays located along the coast of the northern Baltic Sea, between latitudes and longitudes of 61.78°-64.48° and 17.30°-21.42°.

Environmental variables:

Detailed physical and chemical properties for a subset of the analyzed lakes are provided in Öhman et al. (2006). In summary, these lakes are small, shallow, and oligotrophic. Primary productivity in the area is limited by light rather than nutrients [START_REF] Seekell | The influence of dissolved organic carbon on primary production in northern lakes[END_REF]. The abiotic factors identified by previous studies as being important for fish species composition include lake age, lake area, isolation from the Baltic Sea, pH and oxygen [START_REF] Öhman | Associations between water chemistry and fish community composition: a comparison between isolated and connected lakes in northern Sweden[END_REF][START_REF] Englund | Predation leads to assembly rules in fragmented fish communities[END_REF]. We do not have data on oxygen levels for most of the lakes, but previous studies have shown that shallow depth is a good predictor of low oxygen levels [START_REF] Englund | Predation leads to assembly rules in fragmented fish communities[END_REF][START_REF] Öhman | Associations between water chemistry and fish community composition: a comparison between isolated and connected lakes in northern Sweden[END_REF]. Anthropogenic activities that may also affect species composition include treatment with the piscicide rotenone, liming, and species introductions.

Thus, when selecting lakes, we interviewed fishery rights owners and excluded lakes where the interviewees claimed that such activities had occurred. Human activities in the catchments are minor and mostly concerns forestry and ditching for forestry production. The effects on fish extinction rates of these activities are largely unknown.

Consequently, we focus our analysis on the following environmental variables for bays, isolated lakes, and connected lakes were extracted from [START_REF] Englund | Predation leads to assembly rules in fragmented fish communities[END_REF]. Digital maps were used to estimate each lake's area, catchment area, and distance to the sea while other variables were measured directly in the field. Elevations were estimated using either a highresolution GPS instrument or an automatic level. The maximum depth was measured from a boat with an echo-sounder and pH was measured with a portable pH-meter.

For isolated lakes at low elevations (0-30 m) the time since isolation could be approximated by a linear equation, age = 100*elevation (Renberg & Segerstrom 1981) as the uplift rate was approximately constant the last 4,000 years (Berglund et al., 2004). The age of isolated lakes at higher elevations (>30m) was predicted with a log-log model based on data in Påsse et al., (1996) (log (Age) = 0.57* log (elevation) + 6.048; see Appendix S1 for further details, with log as natural logarithms).

Species richness

Data on fish species richness (Sr) were extracted from [START_REF] Englund | Predation leads to assembly rules in fragmented fish communities[END_REF]. Each lake was sampled summer time (June-Aug) on two different occasions during the period 1999-2007, with the aim to detect all fish species present. The methods used in each lake included minnow traps (5-10), spinning rods (90 min), and 2-4 multi-mesh gill nets. The sampling effort for gill nets and minnow traps was adapted to lake size following [START_REF] Appelberg | Development and intercalibration of methods in nordic freshwater fish monitoring[END_REF]. The nets each had twelve 2.5 m sections, with knot-to-knot mesh sizes of 5, 6.25, 8, 10, 12.5, 15.5, 19.5, 24, 29, 35, 43, and 55 mm. To ensure detection of stickleback species, which were sometimes missed by the other methods, we also surveyed the littoral zone of each lake using either electricity (45 min) or small detonations (n = 10 per lake; (Öhman et al., 2006). The salinity of the northern Baltic sea is low (0.1-0.4%) so most (>90%) of the fish species living in lakes are also present in the sea. Among isolated lakes, 18 lakes were fishless and for the other lakes species richness ranges from 1 to 4. Species occurring in these lakes are nine-spined stickleback (Pungitius pungitius L.), three-spined stickleback (Gasterosteus aculeatus L.), roach (Rutilus rutilus L.), crucian carp (Carassius carassius L.), and two piscivores: perch (Perca fluviatilis L.) and pike (Esox Lucius L.). In bays and connected lakes we infrequently encountered eight additional freshwater species: bleak (Alburnus alburnus L.), bream (Abramis brama L.), ide (Leuciscus idus L.), ruffe (Gymnocephalus cernuus L.), rudd (Scardinius erythrophthalmus L.), dace (Leuciscus leuciscus L.), white bream (Blicca bjoerkna L.) and European minnow (Phoxinus phoxinus L.), and two marine species (herring: Clupea harengus and sprat: Sprattus sprattus). The dataset of the isolated lakes used for this study is available from the Dryad Digital Repository [START_REF] Bellard | Data from: Biotic and abiotic drivers of species loss rate in isolated lakes[END_REF].

Theoretical models to explore drivers of extinction rates

Because no colonization events have been recorded in our isolated lakes and speciation are unlikely to occur in such short time scales, we assumed that the isolation has caused species richness to decrease over time. Under these conditions, the observed species loss rate equals species extinction rate. To model species loss rate, we adapted the modelling approach described in Hugueny (2017, Appendix S2, eq. S9) to deal with more explanatory variables.

Using this approach, we modeled the effects of time since isolation (t), lake area (A), pH (P), depth (D), and presence/absence of piscivores (pisc) on the proportion of extant species at t years after isolation:

Sr t /Sr 0 =1/[1+ (c +vb*pisc)*m*(A b )*(P x )*(D y )*t ] 1/m (1)
It is expected that, all else being equal, Sr t /Sr 0 decreases with age (as extinction events accumulate). The parameter m, in conjunction with the variable t, sets the shape of this species decay curve. Eqn ( 1) is undefined when m=0, but as m→0, the dynamics approaches an exponential decay (constant species loss rate per species per year through time) with ln(Sr t /Sr 0 ) being a linear function of t. Setting m>>0 results in a convex decreasing relationship between ln(Sr t /Sr 0 ) and t, meaning that species loss rate per species per year is no more constant but slows down as time since isolation elapses (hereafter "Age effect").

The other parameters set the contribution of the selected variables to species loss, and can be interpreted as coefficients in a multiple regression. For instance, if x is negative it means that, all else being equal, there is a negative relationship between pH and species loss rate. Because pH varies with age, we did not model the direct effect of pH. Instead, we included as predictor the residuals of the model (log(pH) ~ log(Age)) to account for the pH variation that is not related to age.

Note that the parameter vb, which describes the effect of piscivory, enters the equation differently to account for the binary nature of this variable (0: piscivores absent; 1: piscivores present). Note also that the parameter b has been constrained to be null or negative to be consistent with theoretical expectations. From eqn (1) different models could be fitted to the data by setting some parameters to zero (Table 1).

Estimating species richness at the time of isolation

Because of the constant rise of the land, the coastal bays are becoming isolated lakes with time. Moreover, as salinity in the bays are low (<0.3 %), all the species occurring in coastal lakes are also present in bays. Thus, we can expect that before isolation, species richness within isolated lakes were similar to either bays or connected lakes near the coastline. A common method to estimate species richness before isolation is to use the information about species richness from reference regions within a continuous piece of habitat but otherwise similar to the isolate of interest. Hence, we made the assumption that initial species richness for isolated lakes was similar to the species richness observed in connected lakes and bays and estimated Sr 0 for isolated lakes based on stepwise regression models obtained. Based on connected lakes and bays dataset, we found that only area and elevation significantly influenced species richness when piscivores are included, while elevation and pH were the only significant predictors when piscivores were excluded (see Appendix 2 for details). Based on these minimum adequate models (Appendix S2 for details), we predicted Sr 0 of isolated lakes, by setting elevation and distance to the sea equal to 0 and using modelled pH-values and area. Thus, Sr 0 was also calculated (i) when including all species that occurred in the connected lakes and bays, and (ii) when excluding the two piscivores, pike and perch. This resulted in two different estimates of Sr 0 for each lake. Piscivores were excluded from Sr 0 to allow assessment of their impact on the extinction rate of potential prey species. In all cases, if Sr t > Sr 0 we assumed that Sr t =Sr 0 and no extinction occurred. Table 2 presents the best models and their AICc values, and Appendix S4 and S5 show the parameters values.

Estimations of parameter values and comparison between models

For the sake of graphical comparison of observed and predicted values, we computed an average annual extinction rate per species per year over the period 0-t as follows:

k=-ln(Sr t /Sr 0 )/t (3)

Observed and predicted k are computed using the observed and predicted Sr t , respectively, with the constraint 0< Sr t /Sr 0 <1 in the case of observed values.

Parameter values of eqn [START_REF] Whittaker | A general dynamic theory of oceanic island biogeography[END_REF] were estimated by maximizing the likelihood of observing Sr t using the function 'nlminb' in the R environment (version 3.4.3, R Core Team, 2017). A binomial distribution for Sr t was assumed, using the integer part of Sr 0 to set the number of trials and the right side of eqn [START_REF] Whittaker | A general dynamic theory of oceanic island biogeography[END_REF] to compute the probability of success. In this way, and in contrast with [START_REF] Hugueny | Age-area scaling of extinction debt within isolated terrestrial vertebrate assemblages[END_REF], communities with 0 and 100% of extinct species can be considered. Then, we compared all the models using the Akaike information criterion (AICc)

with AIC = 2k -2log(L) with k the number of estimated parameters and L the log likelihood and AICc criteria as follows AICc = AIC + (2k 2 + 2k) / (n-k-1) with n the sample number, defined for small samples (n/k <40) [START_REF] Burnham | Model selection and multimodel inference: a practical information-theoretic approach[END_REF]. The sensitivity of parameter estimates to initial species richness estimations (i.e., Sr 0 ) were examined. In this aim, we modified Sr 0 values by 10% (in other words, we first increase and then decrease Sr 0 by 10%) and assess how parameter values from equation ( 1) vary, compared to results from Table 2 and Appendix 4 and 5.

Results

Extinction rates at community level

The average of the initial species richness (Sr 0 ) was 4.05 when including piscivore species, and 2.57 when excluding the piscivores. The observed extinction rate (k in eq. 3) ranged from 0.0001 to 0.008 (mean 0.0015) when including piscivores, and from 0 to 0.006 (mean 0.0013) when the piscivores were excluded (Appendix S6). From Figure 1, it is clear that observed extinction rates decrease with the time elapsed since isolation. If extinction rates are constant among species and through time, a flat relationship with age would have been expected. The relationship with area is negative, as expected, but relatively weak. We also observed a clear negative relationship between species loss rates and pH residuals (from a regression of pH on age).

Eq. 1 was used to examine the relationships between extinction rate in isolated lakes and the following variables: area, age, piscivores, pH, and depth. In all models the parameter m, related with age, is positive and high, indicating that species loss rate departs strongly from an exponential decay, as shown by figure 1. For both Sr 0 estimates (with and without piscivores), the best model is the ApH model, which includes negative effects of area and pH on extinction rates (see Table 2 and Appendix S4 and S5 for parameter values). This model slightly outperforms the full model (∆AICc 1-3.3) and is much better than all other models.

The full model shows that extinction rate decreases with increasing area, pH and depth when piscivores are included. We also found a positive effect of piscivores on extinction rates (when they are excluded of the calculation of Sr 0 ), meaning that extinction rate increases with the presence of piscivores (i.e., perch and pike), as well as depth. Overall, the ApH model provides a good fit to observed data (Figure 2). Note that models have been fitted with all the lakes, while only lakes having 0<Sr t <Sr 0 are plotted in figure 2. The apparent overestimation of observed extinction rates occurs because fishless lakes (with high extinction rates) are not shown. Sensitivity analyses show parameter estimates are robust to variation in initial species richness, Sr 0 (Appendix S6).

Discussion

Extinctions are rarely observable, particularly when they occur over centuries and millennia.

Habitats that became isolated many years ago, such as land-bridge islands created by sea level rise at the end of the Pleistocene, are valuable natural experiments to study extinction at work in natural communities. Such natural experiments, also known as relaxation, have greatly contributed to the estimation of extinction rates in terrestrial communities and to the identification of underlying ecological mechanisms. In contrast, freshwater ecosystems have been largely ignored in this regard and to our knowledge our study is the first to estimate and study species loss rate in isolated lakes following a relaxation approach.

Area effect

If area is a surrogate for population size then, on theoretical grounds [START_REF] Hanski | A practical model of metapopulation dynamics[END_REF][START_REF] Lande | Risks of population extinction from demographic and environmental stochasticity and random catastrophes[END_REF], an inverse relationship is predicted between extinction rate and area. Our finding confirms this expectation and is consistent with earlier empirical studies [START_REF] Diamond | Biogeographic kinetics: estimation of relaxation times for avifaunas of Southwest Pacific islands[END_REF][START_REF] Halley | Neutral theory as a predictor of avifaunal extinctions after habitat loss[END_REF][START_REF] Hugueny | Age-area scaling of extinction debt within isolated terrestrial vertebrate assemblages[END_REF]. The estimated strength of the extinction-area relationship varies depending on the method used to estimate initial species richness. In the absence of co-variables and for models accounting for an age effect (A models with and without piscivores), the mean estimate is -0.67, which is similar to the values estimated with the same modelling framework for terrestrial vertebrates, -0.53 [START_REF] Hugueny | Age-area scaling of extinction debt within isolated terrestrial vertebrate assemblages[END_REF], and river fishes, -0.56 [START_REF] Hugueny | Habitat fragmentation and extinction rates within freshwater fish communities: a faunal relaxation approach[END_REF]. However, the range of our estimates is rather wide (-1.25 to -0.08), and when co-variables are added to the model, the estimated parameter value is even more variable, and sometimes null (the highest value that we considered realistic are returned). We suspect that the narrow range of lake sizes in this data set results in a rather weak area signal (see also Öhman et al., 2006) that can be swamped by the stochasticity inherent to species poor communities or erased when the effects of co-varying variables are accounted for. On this basis, it is not surprising that the area effect was weaker in the analyses conducted by excluding the piscivores as the initial species richness is particularly low in this case. Indeed, for the analyses that exclude the piscivores the average estimated value was low, -0.05, suggesting that area had a weak effect at best. In light of our results, it is therefore difficult to assess whether or not the relationship between extinction rate and area is lacustrine-specific or if a general pattern exists for freshwater fish or all vertebrates.

pH effect

Natural and anthropogenic acidity is a chemical factor known to structure lacustrine fish communities in many regions [START_REF] Bergquist | Extinction and natural colonization of fish in acidified and limed lakes[END_REF][START_REF] Haines | Acidic Precipitation and Its Consequences for Aquatic Ecosystems: A Review[END_REF][START_REF] Rahel | Factors Structuring Fish Assemblages Along a Bog Lake Successional Gradient[END_REF]. In our study system acidity is due primarily to natural processes (Öhman et al., 2006). The negative effects of acidity on fish populations are both direct, through impaired gill function and ionoregulation for instance [START_REF] Fromm | A review of some physiological and toxicological responses of freshwater fish to acid stress[END_REF], and indirect such as increasing aluminium toxicity [START_REF] Wigington | Episodic Acidification of Small Streams in the Northeastern United States: Ionic Controls of Episodes[END_REF]. The inclusion of pH in our models largely increased AICc values compared to models without it. Overall, we observed that acid-stressed lakes experienced higher rates of extinctions than more neutral ones, as observed for anthropogenic acidification in other parts of Sweden [START_REF] Bergquist | Extinction and natural colonization of fish in acidified and limed lakes[END_REF]. This result is coherent with a previous study showing that, in the same region, the occurrence of most fish species was negatively correlated to lake acidity (Öhman et al., 2006), pike and perch being the more tolerant (see also [START_REF] Bergquist | Extinction and natural colonization of fish in acidified and limed lakes[END_REF]. For instance, survey of many lakes showed small probability of finding roach at pH lower than six [START_REF] Holmgren | Critical values of different acidity indicesas shown by fish communities in Swedish lakes[END_REF]. Studies of anthropogenic acidification suggest that extinction of populations is due to recruitment failure [START_REF] Bergquist | Extinction and natural colonization of fish in acidified and limed lakes[END_REF] because the young stages (eggs , fry and alevins) are more vulnerable to acidity.

Age effect

We observed that the average extinction rate per species per year was negatively related to the time elapsed since isolation of the lake, as already noted for terrestrial communities isolated in habitat fragments or islands [START_REF] Hugueny | Age-area scaling of extinction debt within isolated terrestrial vertebrate assemblages[END_REF]. When extinction rates are estimated considering all the species and when area is the only explanatory variable, estimated values for parameter m setting the strength of the age effect is estimated at 5.55. It is likely that those estimates are inflated because of the action of environmental variables that are correlated to altitude, and hence to age of the lakes. When depth and/or pH are added as co-variables, the estimates for the age effect are lower, about 4.4 on average, and comparable to the value of 3.9 estimated for terrestrial vertebrates [START_REF] Hugueny | Age-area scaling of extinction debt within isolated terrestrial vertebrate assemblages[END_REF].

A plausible explanation for an age effect is interspecific variability in extinction proneness [START_REF] Ferraz | Rates of species loss from Amazonian forest fragments[END_REF]. Indeed, species that are more sensitive to go extinct (e.g., due to small population size or long generation time) are more likely to disappear rapidly after the isolation events, while species more resistant to extinctions are more likely to persist. This mechanism is analogous to the term "extinction filter" which was used by [START_REF] Balmford | Extinction filters and current resilience: the significance of past selection pressures for conservation biology[END_REF] to explain why communities may appear more resilient to particular threats if they have faced similar challenges in the past. Restricting the analyses to non-piscivorous species led to a considerably lower age effect (ca. 1) then when all the species are included, giving strong support to the extinction filter hypothesis. Predator species generally have lower population sizes than other species because of the higher energetic demand imposed by their large body size and top position in the food web and, as a result, are extinction-prone species. The fact that the two focal piscivorous species (pike and perch) are particularly vulnerable to winter kill arising from low oxygen levels under the ice, is also worth considering in this regard. In this line of reasoning, excluding pike and perch from analyses is likely to result in lake assemblages with lower interspecific variability in extinction proneness and to weaken the strength of the age effect on species loss rate.

Another hypothesis for explaining that extinction rate per species is lower in older lakes is ecological release resulting from the absence of competitors or predators. If some pairs of coexisting species are involved in biotic interactions sufficiently strong to induce population declines or to reduce population densities, then extinction rate per species should be higher in species rich, young, assemblages than in old and species-poor ones. Some diversity dependent models (e. g. [START_REF] Halley | Dynamics of extinction debt across five taxonomic groups[END_REF] allow integrating the number of coexisting species (more specifically the number of species present at the time of isolation) as a factor controlling extinction rate per species. Unfortunately, they are of limited use in our case because there is not enough variability in the estimated species richness at the time of isolation to assess its contribution to species loss rate. In addition, because our initial species richness is estimated from lake features and not directly extracted from empirical data, it is not possible to separate the effect of the variable entered into the equation used to predict the initial species richness from the direct effect of species richness itself.

Piscivores effect

In contrast with competition, the effects of predators on the loss rates of prey species can be assessed in our study system. The best model including the presence of piscivores as a covariable is outperformed by the best model without a predation effect, however the difference in AICc values is small (∆AICc =3.3). The presence of northern pike and perch positively influenced extinction rates showing that predation by northern pike and perch can cause extinctions at community levels in isolated lakes. We also found that the best models with a predation effect were those with pH and depth as additional co-variables. Interestingly in those models we observed that extinction rates and depth are positively correlated while in models fitted with all species depth always contributed negatively to extinction rates. The contribution of pike and perch to species loss rate is concordant with results from a previous study of the same system of isolated lakes in which negative co-occurrence between predators and prey were observed and interpreted as the result of prey extinction in the presence of predators [START_REF] Englund | Predation leads to assembly rules in fragmented fish communities[END_REF]. There are other strong evidences in the literature that northern pike influence the abundance and size distributions of crucian carp [START_REF] Brönmark | Predation as a determinant of size structure in populations of crucian carp (Carassius carassius) and tench (Tinea tinea)[END_REF]. The presence of large piscivorous perch is also known to hinder the re-establishment of roach, unless coexistence of perch and roach is mediated by the presence of pike [START_REF] Persson | State-dependent invasion windows for prey in size-structured predator-prey systems: whole lake experiments[END_REF]. Monitoring small lakes after the introduction of perch or pike revealed rapid extinction (within less than 50 years) of nine-spined and three-spined sticklebacks [START_REF] Englund | Predation leads to assembly rules in fragmented fish communities[END_REF] and small cyprinidae [START_REF] Nicholson | Apparent extirpation of prey fish communities following the introduction of northern pike (Esox Lucius)[END_REF]. The fact that piscivores enter the best models only in combination with depth clearly suggest that their impact on prey species is sufficiently strong to induce extinction in deep lakes but this outcome is less likely in shallow lakes confirming the complex relationships between piscivores, preys and depth already emphasized in a previous study of the very same system [START_REF] Englund | Predation leads to assembly rules in fragmented fish communities[END_REF]. Perch and northern pike are more likely to occur in deep lakes [START_REF] Englund | Predation leads to assembly rules in fragmented fish communities[END_REF][START_REF] Macdougall | Context-dependent interactions and the regulation of species richness in freshwater fish[END_REF][START_REF] Magnuson | Isolation vs. extinction in the assembly of fishes in small northern lakes[END_REF] which are buffered against winter kill (i.e., arising from low oxygen levels under the ice) which induce fish to suffocate [START_REF] Petrosky | Behavioral Responses of Northern Pike, Yellow Perch and Bluegill to Oxygen Concentrations under Simulated Winterkill Conditions[END_REF].

Therefore, when these two species are taken into account to estimate extinction rates, depth and species loss rate are negatively related as they are more likely to persist in deep lakes. In contrast, when those two piscivores are excluded of the estimation of extinction rates, there is a positive relationship between depth of lakes and extinction rates as in deep lakes piscivores are more likely to occur and are responsible for extinctions of prey species. When predators occur in shallow lakes they are probably present at low density and coexistence with the more vulnerable preys could be facilitated. Shallow lakes are also known to offer refugee from piscivores, such as abundant macrophytes (Sand Jensen & Borum 1991, Snickars 2009).

Finally, pike and perch are generalist predators that can potentially drive one or several prey species to extinction without undermining their own persistence. For instance, in our database there are lakes where pike is the only recorded fish species, which is plausible as pike may feed on insects, ducklings and small pike [START_REF] Dessborn | Pike predation affects breeding success and habitat selection of ducks[END_REF][START_REF] Venturelli | Diet and Growth of Northern Pike in the Absence of Prey Fishes: Initial Consequences for Persisting in Disturbance-Prone Lakes[END_REF], even rather large pike can have insects as a major part of the diet. However, it is likely that extinction is not the outcome of all prey-predator interactions taking place in the studied lakes. First, as discussed before, predator impact is dependent upon lake environmental features such as depth or pH that may influence trophic dynamics but also population sizes.

Second adults of the largest prey species such as carp and roach may be difficult to handle for a perch and even a pike. Third both preys and predators feed on zooplankton in their youngest stages and a prey may therefore have a negative impact on predator population dynamics through competition [START_REF] Byström | Competing predators and prey: juvenile bottlenecks in whole-lake experiments[END_REF]. For these reasons it is likely that the contribution of predation to species loss rate in isolated lakes is limited in our system mostly to the two stickleback species as suggested by the fact that these two species are rapidly extirpated after piscivore introductions and possibly to carp which rarely coexist with predatory fishes.

Limitations

Some limitations inherent to our study suggest that some of our findings should be treated with caution. For instance, species richness at time of isolation was modelled, rather than observed, which likely introduces some errors in our calculations (though such data are rarely available, especially at such time scale). The sensitivity analyses of initial species richness had little effect on the estimation of species loss rate per lake on average (Appendix 6). This robustness suggests that, as a whole, our estimates are trustworthy and could be compared to those from other systems. The comparison of species loss rate among lakes with the aim to identify underlying biotic or abiotic factors is much more dependent on the way initial species richness was estimated. However, even in this case, robust results emerged such as a strong and consistent effect of age and pH and the synergetic action of piscivory and depth. Based on geomorphological arguments we assumed that no natural colonization events have taken place in those lakes but we cannot discard with certainty that some (but few) populations have been introduced by humans. Indeed, high quality historical data is a key limiting factor for studying extinction at the scale of centuries and millenia. Based on interviews of fishery rights owners, only lakes where the interviewees claimed that species introductions had not occurred were included. This does not preclude very old introductions or introductions made without informing the fishery rights owner. Such introductions would lead to underestimates of the extinction rate. Lack of data for population abundances also prevented us from exploring extinction rates at species level and more generally to rank species according to their extinction proneness. Moreover, because of the low prevalence of some species (Nine-spined stickleback, 90% of absences; three-spined stickleback: 97% of absences), it was not possible to conduct statistical analyses to investigate factors that may influence their extinction rate, and in particular the role of predators. We also investigate extinction rates considering each lake as an independent unit, but it is most likely that nearby lakes share environmental features not considered in our analyses. Finally, other drivers (not taking into account in our analyses due to data availability) might influence species loss rates in this system. For instance, the effects of anthropogenic activities (e.g., forestry and ditching for forestry production) on fish extinction rates in this area are largely unknown. Moreover, productivity in those lakes is mostly limited by light rather than nutrients. Light limitation is related to water color, which in turn is well correlated with the size of the catchment. Thus, had there been strong effects of productivity on extinction rates, we would expect catchment area to represent this effect.

Implications

Our study revealed that to some extent isolated lakes are not so different from terrestrial isolates (islands, habitat patches, …) with regard to the dynamics of species loss. In both cases species loss rate decreases as time since isolation elapses and extinction rates are higher in small isolates. While based on different processes, age is also a strong component of the unified theory of island biogeography [START_REF] Borregaard | Oceanic island biogeography through the lens of the general dynamic model: assessment and prospect[END_REF][START_REF] Whittaker | A general dynamic theory of oceanic island biogeography[END_REF] stressing that this factor should not be overlooked in future works conducted on both true islands and island-like habitats. We also observed peculiarities specific to lacustrine habitat regarding extinction patterns. In particular, we observed that prey extinction rates seem to be regulated by the interplay between depth and predation, which is in line with a conceptual framework for lacustrine community dynamics elaborated by Tonn W. M. [START_REF] Tonn | Patterns in the Species Composition and Richness of Fish Assemblages in Northern Wisconsin Lakes[END_REF]. Those observations suggest that our study has a good level of representativeness in light of the available knowledge about vertebrate extinction patterns and lacustrine community dynamics.

Conclusion

Our results add to recent studies [START_REF] Halley | Dynamics of extinction debt across five taxonomic groups[END_REF][START_REF] Hugueny | Age-area scaling of extinction debt within isolated terrestrial vertebrate assemblages[END_REF] that emphasize the implications of unpaid extinction debts in isolated habitats for conservation and environmental impact assessments. Many studies modeling future effects of climate change, land use changes, or biological invasions on biodiversity have neglected extinction debts and the implied community and population processes. Knowing at which rate those anticipated extinctions will occur is paramount when prioritizing conservation actions. This is particularly true for freshwaters which cover less than 1% of the earth's surface while supplying 12% of the fish consumed by humans and are increasingly fragmented by human activities. Recently it has been shown that delayed extinctions in terrestrial vertebrates could be analyzed within a common framework. Unfortunately, such a synthesis is not at hand for freshwater fishes because too few relevant studies have been conducted. A more promising avenue is to test whether generalities observed for terrestrial vertebrate hold true for freshwater ones.

Acknowledgements

This work has been supported by the BiodivERsA project ODYSSEUS funded from the European Union's Horizon 2020 research and innovation program. BH is supported by the French Laboratory of Excellence project "TULIP" (ANR-10-LABX-41 ; ANR-11-IDEX-0002-02).

Data Accessibility

Data of isolated lakes used in this manuscript are available in a dryad deposit (https://doi:10.5061/dryad.5sv75bd/1).

Author's contributions

GE provided the data, CB and BH designed the study and conducted the analyses. All the authors interpret the results and CB wrote the first draft of the manuscript. Table 1:

Tables