Sami Kaidi

Margny Cerema-Dtecemf

Emmanuel Lefrançois

NUMERICAL MODELLING OF THE MUDDY LAYER EFFECT ON SHIP'S RE-SISTANCE AND SQUAT

), the World Association for Waterborne Transport Infrastructure (PIANC) established the concept of the nautical bottom. Using this concept, ships can navigate with both small and negative under keel clearance (UKC) relative to the water-mud interface. Hence, the aim of this work, is to conduct a numerical investigation in order to study the influence of the muddy seabed on the ship's manoeuvrability especially on the ship's resistance and squat. Accordingly, a 3D Computational Fluid Dynamics (CFD) model based on the Volume of Fluid (VoF) method was used to simulate the multiphase flow for various setups. Four parameters were tested: the mud properties, the ship's speed, the mud thickness and the UKC value relative to the water-mud interface. The numerical results of this investigation were in reasonable agreement with experimental data. Through this investigation it was also shown the performances of the CFD method to simulate setups difficult to achieve in towing tank.

Any ship navigating through confined and shallow waters is strongly affected by hydrodynamic effects, as opposed to in open seas. Major effects of the limited navigating width and water depth (h) are the squat effect, and the increase in the ship's resistance. Water in front of the bow is pushed away, and flows down to the sides and under the hull of the ship with an increased velocity (See Figure 2) due to the reduced section. According to Bernoulli's principle, increasing velocity under the hull indicates a vertical pressure drop, and consequently the ship's sinkage increases. In addition, the ship generally trims forward or aft, as the bow or stern may experience more or less pressure drop, depending on the ship type. The effects of sinkage and trim are known as the ship's squat. This has a significant influence on the ship's resistance and can lead to serious safety issues, such as grounding, loss of steering, or collision.

In estuaries, the presence of the mud layer can significantly modify the ship's behaviour, especially when the ship is navigating in negative UKC relative to the mud/water interface. Note that the concept of the nautical bottom was established by the PIANC MarCom Working Group 30 in their 2014 report. This concept enables ships with larger drafts whose physical properties do not exceed the critical limit (whereby contact with the ship's keel causes damage or unacceptable effects on controllability and manoeuvrability) to navigate in the mud layer. The same report also noted that it is difficult to give the critical limit value, hence, different density limits were set for different ports. That said, a critical limit was still provided, based only on the mud density, where the nautical bottom is the level from where the mud density is more than 1200 kg/m 3 . The viscosity of the mud could not be used as a parameter to define the nautical bottom because it changes under shear rates change. In some ports, such as the Port of Emden, Germany, the critical limit is given as a yield point that has been fixed to 100 Pa (Wurpts, 2005). Using this criterion, it was observed that the corresponding bottom density (approximately 1300 kg/m 3 ) considerably exceeded the limit given by PIANC.

In the Gironde estuary, the squat is an essential parameter for the traffic management of ships, where the water level in the estuary depends on the tide. Accordingly, to accommodate larger ships it is necessary to wait until the tide is high.

Ships have to sail at the same speed as the propagation of the tide wave, which is in the order of 10 kn. However, this is not always the case, because in some situations ships can no longer keep up with the speed of the tidal wave for various reasons (mainly related to the significant increase in the ship's resistance caused by the ship's squat in the mud, which slows its speed considerably). In other situations, ships are equipped with a power limiter that stops the operation of the propulsive system if the ship meets a strong resistance. In such a situation, ships will be moored to wait for the next tide.

To manage the estuarine network better, and to ensure safe navigation, it is thus essential to study the phenomenon of ship's resistance, its origins, and any consequences for navigation.

To predict a ship's squat, several empirical formulas have previously been proposed. Barrass and Derrett (1999) concluded some important factors of the squat effect, as follows:

• The main factor is the ship speed relative to the water, and the squat is approximately proportional to the square of this velocity;

• The decrease of water depth will increase the ship's squat;

• The block coefficient of the ship (that is the ratio of the ship's underwater volume to the volume of box surrounding it) is proportional to the squat;

• Similarly, the blockage factor (a ratio of the ship's immersed cross section to that of the canal) has a direct impact on the squat. Over time, numerical efforts have been focused on estimating ship's resistance and squat. The slender-body method assumes that a ship's beam, free surface wave amplitude, and water depth are small compared to its length. This allows the simplification of the flow simulation in two dimensions using the slender-body theory [START_REF] Gourlay | Slender-body methods for predicting ship squat[END_REF][START_REF] Tuck | A systematic asymptotic expansion procedure for slender ships[END_REF].

To take into account the dynamic coupling of a ship's motion with flow, the potential flow theory can be applied, which only assumes the flow to be irrotational. This has been widely used for squat prediction and very good results have been obtained [START_REF] Debaillon | Numerical investigation to predict ship squat[END_REF][START_REF] Ma | Practical estimation of sinkage and trim for common generic monohull ships[END_REF][START_REF] Sergent | Virtual bottom for ships sailing in restricted waterways (unsteady squat)[END_REF], whereas it is difficult to apply the potential flow model to resistance prediction, because it neglects the viscous stresses, crucial for evaluating the ship's resistance.

Modern CFD techniques based on solving the fully viscous Navier-Stokes equations have been extensively applied to ship hydrodynamics with fruitful results, as they consider the important features of the actual flow, such as viscous effects and turbulence. Hence, they are more reliable for predicting ship's resistance and motions. [START_REF] Stern | Computational ship hydrodynamics: nowadays and way forward[END_REF] summarised the achievements made regarding ship hydrodynamics using CFD in the last decade. Further, recent progress in modern computational ship hydrodynamics with respect to shallow and confined water has been made. [START_REF] Eloot | Validation of ship manoeuvring in shallow water through free-running tests[END_REF] performed a turning circle and a zigzag test on a KVLCC2 hull model to determine the manoeuvring performance in a shallow water zone. To study the scale effect, [START_REF] Tezdogan | Full-scale unsteady RANS simulations of vertical ship motions in shallow water[END_REF] performed unsteady Reynolds-averaged Navier-Stokes (RANS) simulations at full-scale for the squat in shallow water with STAR-CCM+ commercial software.

They compared the results to the 3D potential flow theory and the experimental data of [START_REF] Mucha | Technical Note: PreSquat -Workshop on Numerical Prediction of Ship Squat in Restricted Waters[END_REF]. They reported an underestimation of the ship squat and pointed indicated that the ship's resistance is sensitive to sinkage. [START_REF] Linde | Three-dimensional numerical simulation of ship resistance in restricted waterways: Effect of ship sinkage and channel restriction[END_REF] In ports, flow stratification might occur as the non-saline, light river water flows into the colder and heavier saline seawater, leading to large horizontal or vertical fluid density variations. Highly density-stratified waters are known to pose particular challenges to navigating ships. When a ship's keel (bottom) is travelling just above the interface of the water layers, the ship experiences large wave resistance. This resistance occurs particularly if the ship is travelling close to the speed of the fastest internal waves, due to the generation of large internal waves. This phenomenon, known as `dead water', affects the ability of ships to move through stratified water. Accurately assessing the effects of stratified flow on ship navigation requires a detailed knowledge of the flow field, including turbulent mixing and in particular, the generation of internal waves on the interface between the two layers of water. [START_REF] Crapper | Ship waves in a stratified ocean[END_REF] and [START_REF] Hudimac | Ship waves in a stratified ocean[END_REF] presented analytical approaches to study the internal wave modes caused by a moving body in a two-layered ocean. It follows from their work that-just as for surface waves-at ship speeds sufficiently larger than the internal wave speed, only divergent waves travel downstream from the ship, while both divergent and transverse waves are present for slower ships. [START_REF] Tulin | The generation and propagation of ship internal waves in a generally stratified ocean at high densimetric Froude numbers, including nonlinear effects[END_REF] suggested a nonlinear theory to capture internal wave behaviour at high Froude number (Fn) in weakly stratified flow, which compared satisfactorily with available experimental results for a semi-submerged spheroid. [START_REF] Delefortrie | Linear manoeuvring derivatives in muddy navigation areas[END_REF][START_REF] Delefortrie | Modelling navigation in muddy areas, Hydronav'05 -Manoeuvring '05[END_REF] conducted a large number of experiments on towing tanks. They studied the mud layer effect on the ship manoeuvring by considering several parameters. They also developed a mathematical model to take into account the mud effect. (Chang et al., 2006) presented one of the few available examples of the use of CFD for a ship in a stratified medium. [START_REF] Esmaeilpour | Near-field flow of submarines and ships advancing in a stable stratified fluid[END_REF] studied the evolution of the stratified flow in the near field of a surface ship in detail. They demonstrated that the generation of internal waves requires energy, which results in an increase in resistance.

In this paper, we present an overview of a numerical study of the mud layer effect on ship's resistance and sinkage by using a multi-phased CFD method. Note that, the mud is supposed stratified, hence, the water and the mud were modelled as separate layers with average values of density and viscosity. The main objective of this work is to test the ability of the CFD method to simulate and assess the influence of the mud on the hydrodynamic forces acting on the ship's hull.

Four parameters were tested: the mud properties, the ship speed (Vs), the mud thickness and the UKC. A preliminary study was conducted to first show the influence of the non-Newtonian behaviour of mud on the ship's resistance, and the internal waves at the mud-water interface. Based on this preliminary study, the Newtonian model was selected for this investigation. The UKC level was referenced to the water-mud interface; hence, it can take both positive and negative values. The limits of the CFD method are discussed in Section 5.

MATHEMATICAL FORMULATION AND NUMERICAL METHODS

The fluid flow is governed by the incompressible viscous Navier-Stokes equations completed with the continuity equation, as follows:

∇. = 0 (1)

+ ∇. ( ⨂ ) = -∇p + ∇² (2)
where u and p represent the velocity vector and pressure, respectively, and ρ and ν are the fluid properties of density and kinematic viscosity.

TURBULENCE MODELLING

To model the turbulence effect, the Reynolds averaging was computed on the flow variable in time, which gave rise to ∇. # = 0

(3)

+ ∇. (#⨂#) = -∇p + η ∇ % # -∇. ( ⨂ ) (4)
where = # + and p = P + p . The last term in the RANS momentum equation is the Reynolds stress, which is often approximated by turbulence models. In this research, we employed the SST k-ω turbulence model, which is actually a combination of the k-ω and k-ε models while a shifting function is used to switch one from another.

MULTIPHASE APPROACH

The volume of fluid (VoF) method was used to simulate three-phase interactions (the interface of air/water and water/mud). Using this approach, both interfaces can be captured in a fixed grid by solving the continuity equation of the volume fraction (Eq. 5), as follows:

∂α ) ∂t + ∇α ) = 0 (p = 1, 2) (5)
α ) denotes the volume fraction of the p th fluid, and:

. α )

/ )0 = 1 (n = 2 or 3) (6)
Equations presented thus far (Eq. 1-6) are solved using the commercial code Ansys-Fluent 13.0 based on the finite volume method. The pressure-velocity coupling was ensured by using a steady pressure-based coupled algorithm, and the interpolation method selected to compute the cell-face pressure was the PREssure STaggering Option (PRESTO). The second order was set for the VoF's special discretisation.

NON-NEWTONIAN BEHAVIOUR LAW FOR MUD

It should be noted that mud behaviour is often considered non-Newtonian, which means the viscosity depends on the shear rate. Hence, the Herschel-Bulkley model was selected to reproduce this behaviour. The Herschel-Bulkley model is represented by the following equations:

= + 5 67 8 (7)
= 0 67 9 [START_REF] Biography | Hi is responsible for Ship Hydrodynamic and Manoeuvring researches. His previous experience includes: CFD method; Fluid-structure interaction, numerical methods[END_REF] where, and are shear and yield stress, respectively. K is the consistency factor, n is the power law exponent, and is the shear rate.

The non-Newtonian viscosity is computed using one of the following formulas:

= : ; < + 5= 7>? 8 (9) = : ; @%= A A B C < + 5= D(2 -E) + (E -1) < < B F 7>? G (10)
is the critical shear rate.

Figure 3 shows the variation of the shear stress with shear rate according to the Herschel-Bulkley model. works [START_REF] Kaidi | Numerical estimation of bank-propeller-hull interaction effect on ship manoeuvring using CFD method[END_REF](Kaidi et al., , 2018;;[START_REF] Razgallah | The impact of free surface modelling on hydrodynamic forces for ship navigating in inland waterways: water depth, drift angle, and ship speed effect[END_REF]Ali et al., 2018), and were carried out in accordance with the ITTC recommendations. Verification consisted of tests and analysis of the results of several mesh qualities, whilst validation was performed by taking the numerical results of the ship's resistance, the profile of generated waves, and the ship's squat and comparing them to the measurements carried out in a towing tank of the University of Liège and the Central School of Nantes. Note that no mud layer was considered in these works: only the water-air interface was modelled and validated.

Based on the obtained results, it was concluded that the CFD model provided a good estimate of the hydrodynamic forces around the ship's hull, and correctly captured the air-water interface. It was also concluded that the CFD model could simulate the sediment suspension and transport accurately (Kaidi et al., 2018).

STUDIED SHIP, CHANNEL CONFIGURATION, BOUNDARY CONDITIONS, AND MUD PROPER-

TIES

To conduct this investigation, we used a container cargo-hull form (see Figure 4). This kind of ship was selected because it is one of the most common ships to sail in the Gironde estuary. Table 1 provides the main characteristics of the hull, where, LPP is the length between perpendicular, LOA is the length over all, B is the ship beam, T is the ship draft, and CB is the block coefficient. Note that the same reference frame is used in this investigation. The origin of this reference is fixed on the ship, where the x = 0 corresponds to the ship's bow plane, while, the z = 0 corresponds to the ship's keel plane. The X-axis is oriented right while the Z-axis is oriented up as is shown in Figure 4. In the present study, we considered only the confinement effect using the UKC. To prevent large body motions, the reference frame was fixed on the ship; hence, the fluid and other parts moved relative to the hull. The computational domain was chosen with a rectangular section, large enough such that there was little influence of the position of the inlet and outlet: 1-2 Lpp for the inlet and 3-5 Lpp for the outlet are recommended by ITTC [START_REF] Ittc | ITTC -Recommended Procedures and Guidelines -Verification and validation of linear and weakly nonlinear seakeeping computer codes[END_REF]. Half of the computational domain was used to reduce computational time.

For the boundary conditions, at the inlet the flow velocity was imposed, and at the outlet the outflow condition was used.

A symmetrical condition was applied at the top, at the mid-plane, and at the side boundaries. At atmosphere, total pres-sure was applied; at the bottom moving wall condition was employed to take into account the relative motion, and at the hull's surface no-slip wall condition was used.

Four combinations of mud properties were selected to conduct this investigation. These properties represent the average values measured at different zones in the Gironde estuary and some ports. Table 2 presents the combination of the density and viscosity of the mud.

RESULTS AND DISCUSSION

COMPARISON BETWEEN THE NEWTONIAN AND NON-NEWTONIAN MODELS

The effect of the use of non-Newtonian viscosity to simulate the mud behaviour is discussed in this section. Using the CFD method, the mud layer could be modelled in different ways. The first was to suppose that the mud properties were constant and slightly influenced by the shear stress induced by the ship's passage. Hence, we used average values for the density and the viscosity. The second was to consider that the mud was significantly affected by the shear stress, from where a non-Newtonian model was used to estimate the dynamic viscosity.

It should be highlighted that in the estuary environment, the flow was modelled using hydraulic models, which estimate the average turbulent viscosity at the channel bottom that can be used by manoeuvring simulators. This provides indispensable results for assessing the difference between the results obtained by both approaches.

To carry out this study, the ship's draft was set to 10 m (0.125 m in the model) and the ship's speed was set to 10 kn (0.575 m/s). The mud thickness used was 3 m in full scale (0.0375 m in the scaled model), which corresponds to the average thickness in the estuaries. The UKC value was set to +10%*T with respect to the mud/water interface, which corresponded in our case to 1 m in full scale (0.0125 m in the scaled model). To assess the ship's resistance, two UKCs were tested (+10% and -10%). For the two types of mud, the difference between the Newtonian and non-Newtonian models was small, as presented in Table 3 (approximately 2% for both mud A and mud C at positive UKC, and less than 5% at negative UKC. Note that, the verification and validation procedure shows that the uncertainty of the total resistance was about 7% in very confined water. This uncertainty corresponds to a monotonic convergence condition with an order of accuracy of 1.86. Based on these values, the only conclusion that can be drawn from this comparative analysis is that the Newtonian model gives an acceptable estimation of the ship's resistance despite the overestimation of the mud/water undulation. Hence, the Newtonian model can be used to carry out this investigation. Figure 8 illustrates the profile of the mud layer deformation caused by the ship's passage. As can be seen, the deformation is composed of a principal undulation and secondary undulations. The principal undulation is similar to the free surface deformation with a small shift, where a stern divergent wave is observed (Figure 9). The divergence angle, the wave height, and the wavelength of this wave depend on the mud properties, as shown in Figure 10. In this study, we only focused on the principal undulation, which had an impact on the ship manoeuvring, principally on the ship's re-sistance and squat. This undulation was characterised by a maximum trough and crest, where generally the trough is located at the mid hull, whilst the crest is located at the hull's stern. For all tested mud properties, the mud layer trough started from the same position (the ship's bow). However, compared to the initial mud setup, the trough level and length increased by decreasing the mud viscosity. We note that the origin of this trough was principally the pressure variation along the ship hull caused by the return flow, which was influenced by the mud properties (see Figure 11). From the same figure, it can be seen that the relative increase of the mud trough shows a linear variation for viscosities varying between 0.025 and 0.12 Pa.s. It can also be observed that the physical properties of the mud played an important role on the mud crest, the location of this crest, and in some situations the hull/mud contact area. When the density and viscosity of the mud were smaller, the mud was considered more fluid; hence, the later behaved as a denser fluid and followed the water flow. When the viscosity of the mud was greater, the mud layer was more solid and its behaviour was more rigid. From this, we noted a maximum uprising value for the mud D ~20% less than for mud B, whilst an insignificant variation was noted between mud samples B and A. For all tested properties, the mud uprising position varied as the mud properties varied. From simulated cases, it was noted that the lower the viscosity, the more the mud uprising moved backwards. The same observations were also noted by [START_REF] Delefortrie | Modelling navigation in muddy areas, Hydronav'05 -Manoeuvring '05[END_REF]. Contact between the hull and the mud was also observed for mud A and B and the contact area was slightly larger in the case of mud A. The internal waves corresponding to the viscosity and the density variation were plotted in Figure 12. As it can be seen at the left of this figure, the viscosity has a very important influence on the internal waves pattern, and for the used ship's speed the internal bow waves length is longer than the ship's length which is a characteristic of shallow water navigation. For the smallest viscosity, a stern waves pattern was observed, where the first is convergent, while the second is divergent which resembles to the free surface behaviour in a shallow water. Reflected waves were also observed far behind the ship's stern. By increasing the viscosity value to 0.01 Pa.s, the same waves pattern were observed, however, transverses waves appeared behind the ship's stern. For viscosity value of 0.05 Pa.s, the waves pattern becomes more apparent, while the angle of the diverging waves increases and tends to be perpendicular to the ship's heading direction. The same observations were noted for the viscosity value of 0.1 Pa.s, The converging waves remains unchanged and apparent, whereas the diverging waves becomes completely transversal and slightly less apparent. For the highest mud viscosity (0.2 Pa.s), only the transverses waves pattern were affected and becomes even less apparent and tend to disappear. This behaviour, is in accordance with experimental findings, however, viscosity limits defining waves patterns were slightly different.

The density variation effect on the internal waves pattern is shown at the right of Figure 12. The only observation noted from this figure is that the density has any effect on the transverses waves, while the converging waves increase slowly by increasing the mud density.

Basing on these results it can be concluded that the internal Froude number Figure 13 shows the free surface attenuation caused by the mud layer. For this, we compared the free surface deformation of a channel (with a muddy layer) to the free surface deformation of a channel without a muddy layer (rigid bottom). The same total depth was maintained for both tests. Note that for the present study, the total depth of the channel is the sum of the ship's draft, the UKC, and the mud thickness. It may be observed that the free surface elevations closely resemble the mud layer undulation. It was also observed that the trough and crest of the free surface was approximately the same between the rigid seabed and mud samples A, B, and C. However, for mud D, we noted a lowering of the free surface. This lowering was essentially caused by the effect of shallow water. As mentioned above, the higher the viscosity, the more solid the mud; hence, the seabed can be considered solid. The shear stress due to the high viscosity of the mud slowed down the flow velocity of the mud/water interface under the ship's hull, inducing an acceleration of the water flow and consequently a pressure drop.

SHIP'S SPEED AND MUD THICKNESS EFFECT ON THE SEABED UNDULATION (INTERNAL WAVES

PATTERN)

The study of the influence of the ship's speed and the mud's thickness on the seabed undulation is presented in this section. The influence of the ship's speed was performed by setting the mud layer thickness to 3 m in full scale, and the UKC to +10%*T. Mud A was selected for this simulation. The ship speeds that were tested were 6, 8, and 10 kn in full scale (0.345, 0.46, and 0.575 m/s in the scaled model). It can be seen from Figure 14 that for the selected mud, the ship's speed has an influence on the position of the maximum rise of the undulation. The higher the speed, the more the undulation crest moved backwards. It was also seen that the crest width increased as the ship's speed increased. This behaviour was due first to the mud type. Here, the mud density was small and it was thus was more fluid than solid. Second, the return current amplification caused by the increase in the ship's speed affected the behaviour. A slight influence of the ship's speed on the trough depth and the crest height of the mud/water interface undulation was noted. thickness (1m in full scale), the internal waves appear behind the ship's stern and its pattern is very divergent and nearly transversal to the ship's heading direction. For the medium thickness (2 m in full scale), the undulation is principally transversal, however, a small Kelvin pattern appears at the ship's stern, which leads us to consider this thickness as a transition thickness. For the larger mud thickness (3m in full scale), the internal waves have a Kelvin pattern form.

These patterns were also noted in the experimental work done by Delefortrie(2016) (see Figure 17). The authors of this work related these patterns only to the mud viscosity and the ship's speed. According to their findings: the Kelvin pattern appeared for larger ship's speed and lower viscosities ( < 0.12 Pa.s ), and the transversal pattern appeared in the case of lower ship's speed. It was also noted that the transversal pattern were observed for larger ship's speed and larger mud viscosities ( > 0.12 Pa.s). In the present work the both undulation patterns (transversal and Kelvin pattern) were observed by varying only the mud thickness. Where, the ship's speed can be considered larger, however, the mud viscosity is lower (0.025 Pa.s).

To separate the effect of the mud thickness and the effect of the confinement, the total depth was considered constant and the mud thickness was varied. The obtained results were depicted in Figure 16-b. The simulation setup for the smaller mud thickness (1m in full scale) is the same as in the case of variable total depth hence the results of the previous simulation were kept. For the medium mud thickness (2m in full scale) the UKC is 0%*T, the similar waves pattern as for the mud thickness of 1m was observed, however, the undulation amplitude was amplified. For larger mud thickness (UKC = -10%*T) the waves change pattern where a convergence and divergence waves appear. Transversal waves were also observed far behind the ship. From these observations, it can be concluded that the mud thickness has a significant influence on the internal waves pattern which change as a function of the total depth and the UKC, however, it is difficult to specify exactly the influence of the mud thickness. In fact, the mud thickness effect is always coupled with another parameter, either the total depth or the UKC. From this study it can also be concluded that the internal waves pattern are independent of the Frh contrary to the generated waves pattern at the water -air interface.

Basing on the results obtained in this section the pattern of internal waves are not only depends on the mud viscosity and ship's speed as is given in the literature but by the combined effect of the viscosity, the ship's speed, the mud thickness, the water depth and the UKC. It is known that ship's resistance is greatly affected by channel configuration (such as confinement and restrictions). In shallow water, the ship's resistance increases significantly due to the accelerated water around the hull, as explained previously. The presence of the mud layer in turn affects the flow under the ship's hull, inducing a variation in the ship's resistance. The effect of the latter can be considerably amplified if the UKC is negative.

In this section, findings are presented from studying the impact of the mud layer, first by testing the mud properties effect for a given mud thickness (3 m in full scale). Hence, the four mud properties were tested for an UKC of +10%*T with respect to the mud/water interface. The ship's speed was set to 10 kn (0.575 m/s). No squat was considered in these simulations.

Figure 18 shows the ship's resistance variation caused by variations in the mud properties. From this figure, it is evident that the ship's resistance increased with mud viscosity, although there was no contact between the hull and the mud with mud samples C and D. This leads us to conclude that this increase essentially concerns the frictional component of the resistance. In fact, when the mud is consolidated it seems as though the total depth of the water is reduced, which makes the navigation environment more confined, and consequently, the return current velocity increases. From these figures, it can be seen that the ship's resistance increases with increasing ship's speed. However, this increase is different according to the mud layer thickness as well as the type of the mud. In fact, by decreasing the mud layer thickness, the ship's resistance increases considerably. This increase concerns both types of forces: pressure and friction, however, the wave-making resistance remains the most dominant and the most impacted by this decrease in thickness.

The wave-making resistance is amplified about 9 times when the ship is sailing at a speed of 10 kn on a mud thickness of 1 m than 3 m, and about 4 times when the ship's speed is 6 kn. While the amplification of the friction resistance is of the order of 2 times for the two ship's speeds. It was also observed that the wave-making resistance was the most dominant compared to the friction resistance except in the case of the largest thickness (3 m), where both types of forces were approximately in the same range. This leads us to conclude that Note that for navigation in channels with solid seabed the confinement is often defined by the ratio of water depth to ship's draft (hw/T). This ratio is an essential element for the calculation of the ship's resistance. In the case of navigation in turbid water with a muddy bottom, this ratio may be valid however, the definition of the depth hw must be modified by including the thickness of the mud layer. A proposal has already been made by [START_REF] Delefortrie | Ship Manoeuvring Behaviour in Muddy Navigation Areas[END_REF] proposing the use of the hydraulic depth (h*) given by the following formula:

ℎ * = ℎV + ∅ ℎX (11)
Where, ∅ is the fluidization parameter which represents the mud type. So far, all works studying the influence of the mud quality on the ship's resistance were performed by varying both values of viscosity and density. In order to distinguish the influence of each physical property we varied separately the viscosity and the density of the mud. The same process as in the sub-section 4.3 was used: first, the density is fixed at 1100 kg / m and the viscosity is varied, then the viscosity is fixed at 0.1 Pa.s and the density is varied. For a better presentation of the influence of these properties on the ship's resistance, we considered the highest speed of navigation (10 kn) and assume that the ship's UKC is -10% * T.

Three mud layer thicknesses were tested (2 m, 3 m and 4 m) corresponding to a distance between the ship's keel and the solid seabed of 1 m, 2 m and 3 m as is presented in Figure 20. The computed ship's resistances as a function of the mud viscosity and density were presented in Figure 21by its two components: wave-making and friction. Basing on these results it was noted that for smaller and medium mud thicknesses the wave-making resistance is dominant compared to the friction resistance, however, for the larger mud thickness, the friction resistance is dominant. This behavior is physical, because it depends on the depth Froude number.

Where for high values of depth Froude number the wave-making resistance is dominant while the friction resistance is dominant for very lower values. By analyzing the numerical results of the ship's resistance under the variation of the viscosity, it can be seen that the effect of the later begins to be visible when the viscosity is more than 0.01 Pa.s in smaller and medium mud thicknesses. For the larger mud thickness the effect is visible only when the viscosity of the mud is greater than 0.05 Pa.s. It can also be seen that despite the dominance of the wave-making resistance , its variation at a given mud thickness is insignificant compared to the variation of the friction resistance. Where, the maximum variation of the wave-making resistance is about 30% computed between the largest and the smallest value of the viscosity at the larger mud thickness. While, the friction resistance was amplified by 2.8, 1.9 and 1.3 times at smaller, medium and larger mud thickness respectively. These behavior is completely realistic for both types of resistances: the friction resistance depends principally on the shear stress on the ship's hull which increases with confinement. However, the wave-making resistance depends on the length and the amplitude of the generated waves which are important when the mud thickness is larger (see figure 16).

Through Figure 21-b, it can be noted that the effect of the density on the ship's resistance is unimportant for the tested thicknesses. The wave-making resistance variation is insignificant at a given mud thickness, whereas, the friction resistance is slightly affected by mud density, where an increase of 23% was computed between the largest and smallest density at the larger mud thickness. This insensitivity to density is related on one side to the non-variation of the generated waves as is illustrated in Figure 12 and on the other hand to the phenomenon of fluidization of the interface of the mud layer caused by the ship passage especially at high speed. This point will be discussed in the next paragraph.

density variation

The second part of this section presents the results for the ship's resistance, studied as a function of the UKC. Mud A was used for 6 values of UKC, as follows: +10%, +5%, 0%, -5 %, -10%, and -15% of the ship's draft. The ship's draft here was 10 m, the mud thickness was 3 m, and the speed of the ship was 10 kn.

From Figure 22, we observe an increase of the total resistance with the decrease of the UKC. The pressure resistance dominates for UKC, varying between +10 and -5%*T. Less than that, the frictional force dominates considerably. We also observe that ship's resistance increase is very slight for UKC values between +5% and -5% and less than -10%. The ship's resistance increase is significant only between +10% and +5%, and between -5% and -10%. By analysing Figure 23-a plotting the area of contact between the hull and the mud, we note that for UKC range of +10%-+5%, the contact area is almost the same and the total resistance increase is principally due to the shallow water effect. For the UKC range -5%---10%, the resistance increase is principally due to the hull/mud contact. As is shown in Figure 22, the frictional resistance dominates while the pressure resistance increases slightly.

It can also be observed from Figure 23-a that when the ship is sailing inside the mud layer (negative UKC); the keel of the hull is not fully covered by the mud. This is one of the relevant phenomena observed in this work. In fact, when the ship is sailing in the mud, we observe two different behaviours. When the ship's speed is low, the ship's keel is fully covered by the mud. However, when the ship's speed increases, the top boundary of the mud layer tends to be more liquefied, especially at the ship's bow, and a film of very liquefied mud (or turbid water) is created between the mud and , 2016). The standard dynamic Newton algorithm was not used because of several numerical complications encountered. The origin of these complications is essentially the bad estimation of the added mass due to the high blockage coefficient, which considerably affects the stability and convergence of the numerical solution.

Because of the large computation time, only one mud layer thickness was considered (2 m) for an UKC of +10%*T. The effect of the mud properties on the squat was simulated for the four types of mud and for three ship speeds (6, 8, and 10 kn), which correspond to a Froude depth number (Frh) of 0.297, 0.396, and 0.495, respectively.

The ship's sinkage as a function of the mud type was plotted in Figure 24 and compared to the experimental and numerical sinkage for a rigid bottom.

Concerning the ship's sinkage, similar observations given by Delefortrie (2016) were noted. First, we observed that the sinkage obtained numerically for a rigid seabed was in accordance with measurements. We also observed that the sinkage increased by increasing the ship's speed in all configurations with or without the mud layer. However, the sinkage values decreased slightly with the change in mud properties. For larger viscosities (Mud C and Mud D), we observed an insignificant decrease, whilst a moderate decrease was observed for Mud A and Mud B. This decrease augmented in turn with the increase in the ship's speed. In fact, this increase was due to the added buoyancy generated by the contact between the hull and the mud. This contact, as mentioned previously, was located at the ship's stern when the undulation crest was larger, as for Mud A and Mud B.

-3.5

CONCLUSIONS

In this paper, an overview of a numerical investigations on the impact of muddy seabed on a ship's resistance and squat was presented. A multi-phase CFD model was used to estimate the ship's resistance and squat as a function of several parameters: different configurations.

Based on observations noted in the present work, it was concluded the following:

• The obtained numerical results are in agreement with physical models results;

• The internal waves crests depend strongly on the mud properties;

• The internal waves patterns depend strongly on several parameters: the viscosity, the total depth, the UKC and the thickness of the mud layer and the Frh is not adapted to characterize the waves pattern.

• The internal waves influence the ship's resistance and squat especially when the UKC is negative;

• The effect of the mud layer on the ship's sinkage is significant only when the UKC is negative;

• The effect of the mud on the ship's resistance can be felt even when the UKC is positive, and this depends on the mud properties;

• The ship's speed tends to move the mud/water interface undulation in the backwards.

The ability of the CFD method to simulate multiphasic flow and interaction between the fluid flow and the structure has been demonstrated. Some difficulties were encountered in the modelling of the depth-dependent density and viscosity of the mud. Some difficulties were also encountered in the simulation of the dynamic ship squat, especially when the mud contacted the ship's hull. An improvement can be made in future works by integrating a new numerical algorithm that stabilises calculations and takes into account the real vertical profile of the density.

6

Figure 2 .

 2 Figure 2. Nautical bottom representation with flow around a ship.

  validated this observation with FLUENT, by simulating the ship's resistance both with and without consideration of ship sinkage. The predicted value of resistance with sinkage was closer to the experimental data.(Kaidi et al., studied the ship manoeuvring with FLUENT under the effect of bank-propeller hull interaction in shallow water.

Figure 3 .

 3 Figure 3. Shear stress as a function of the shear rate basing on the Herschel-Bulkley model (This figure was taken from the Fluent software manual)

Figure 4 .

 4 Figure 4. Container carrier-hull form.

Figures 5 Figure 5 .Figure 6 .

 556 Figures 5 and 6 present a comparison between the undulations of the mud/water interface obtained using both methods for mud types A and C at two different plans (see Figure7). The first plane is located at Y/B of 0.621 while the second is located at Y/B = 1.243, where Y is the lateral distance from the ship's mid-plan and B is the ship's beam. The x-axis of these figures represented the dimensionless distance X/LPP, where X is the longitudinal distance from the ship's bow and LPP is the length of the ship, and the z-axis represented the undulation elevation. Generally, we remark that the non-Newtonian model tends to underestimate the mud layer crests. As is evident, the difference between both models shows that mud A has a smaller viscosity and density compared to mud C. The height difference at the crest and trough is almost similar; approximately 11% at plane 1 and 18% at plane 2 for mud A, and approximately 16% at plane 1 and 30% at plane 2 for mud C. The difference is larger at plane 2, because, the area affected by the ship's passage is reduced when the mud is considered non-Newtonian, especially at high density and viscosity, as shown in Figure7.

Figure 7 .

 7 Figure 7. Iso-surface of internal waves for mud A and C by considering a Newtonian and non-Newtonian model (UKC = +10%).

  TIONHere, the mud layer thickness was set to 3 m (0.0375 m in the scaled model). The ship's draft and speed were set to 10 m (0.125 m in the scaled model) and 10 kn (0.575 m/s in the scaled model), respectively. The value chosen for the UKC with respect to the mud/water interface was +10%*T.

Figure 8 .

 8 Figure 8. Profile of the mud layer undulation at the ship symmetry plane.

Figure 9 .

 9 Figure 9. Iso-surface of the internal waves as a function of mud properties.

Figure 10 .Figure 11 .

 1011 Figure 10. Mud layer undulation along the channel as a function of mud properties (cut at ship's mid-plan).

Figure 12 .

 12 Figure 12. Internal waves: a) viscosity variation and b) density variation

Figure 13 .

 13 Figure 13. Comparison between the free surface elevation of seabed with and without the muddy layer.

Figure 14 .Figure 15 .

 1415 Figure 14. Mud layer undulation along the channel as a function of the ship speed (cut at the ship's mid-plan). X is the longitudinal position.In the second part of this section the influence of the mud thickness on the internal waves pattern was studied. Three thicknesses of the mud A were tested: 1 m, 2 m and 3 m in full scale. Note that the mud thickness effect on the internal waves patterns can be tested using a constant UKC (10%*T ) and variable total depth as is depicted in Figure15-a or using a variable UKC (+10%*T, 0%*T and -10%*T ) and constant total depth as is depicted in Figure15-b. The ship's speed was set to 10 kn (0.575 m/s in the scaled model). The Frh values were the same for both type of test 1.64, 1.16 and 0.94 (using the mud thickness as a characteristic length) for mud thickness of 1m, 2m and 3m respectively. These values correspond to critical and supercritical regimes. The numerical results of the waves patterns were shown in Figure16.

Figure 16 .

 16 Figure 16. Internal waves patterns as a function of mud thickness : (a) constant UKC and variable total depth and (b) variable UKC and constant total depth.

Figure 16 -

 16 Figure 16-a shows the influence of the mud thickness on the internal waves patterns for a constant UKC and variable total depth. From this figure it can be shown three different waves patterns of the seabed. Where for the smaller mud

Figure 17 .

 17 Figure 17. Measured mud layer undulation carried out by Delefortrie (2016)

Figure 18 .

 18 Figure 18. Ship's resistance (half of ship) as a function of mud properties (for a ship speed of 10 kn and UKC of +10%*T).

Figure 20 .

 20 Figure 20. Mud thickness variation for an UKC value of -10%*T.

  's resistance variation (half of ship) as a function of mud thickness: a) viscosity variation and b)

  Figure 22. Ship's resistance (half of ship) as a function of UCK (mud A).

Figure 24 .

 24 Figure 24. Ship sinkage as a function of the mud properties for an UKC of +10%. ZG is the sinkage at the mid ship

Figure 25 .

 25 Figure 25. Ship trim as a function of the mud properties for an UKC of +10%.

Figure 1. Ship sailing in the Gironde estuary. NOMENCLATURE

			Volume fraction
	B	Ship's beam (m)
	CB	Ship's block coefficient
	Fn	Froude number
	Fni	Internal Froude number
	h	Total depth (m) (water + mud)
	hw	Water depth (m)
	hm	Mud thickness (m)
			Consistency factor (Kg.s^n.2/m)
	LOA	Ship's Length over all (m)
	LPP	Ship's length between perpendiculars (m)
	n	Power law exponent
	T	Ship's draft (m)
	t		Time (s)
			Fluid pressure (Pa)
			Fluctuation of the fluid pressure (Pa)
	UKC	Under Keel Clearance (m)
			Fluid velocity vector (m/s)
			Fluctuation of the fluid velocity vector (m/s)
	Vs	Ship speed (m/s)
			Dynamic viscosity (Pa.s)
	w	Water density (kg/m 3 )
	m	Mud density (kg/m 3 )
			Shear rate (1/s)
			Critical shear rate (1/s)
			Shear stress (Pa)
			Yield stress (Pa)
	1	INTRODUCTION

Table 1 . Ship dimensions and scales

 1

	Ship	Container cargo
	Scale	1/1	1/80
	LPP (m)	230.0 2.875
	LOA (m)	232.5 2.906
	B (m)	32.2	0.402
	T (m)	10.0	0.125
	CB	0.681 0.681

Table 2 . Physical properties of tested mud

 2

	Mud type	density (kg/m L ) viscosity (Pa.s)
	Mud A	1085	0.025
	Mud B	1160	0.068
	Mud C	1210	0.128
	Mud D	1230	0.260

Table 3 . Computed ship's resistance (half of ship) using the Newtonian and non-Newtonian models.

 3

			Ship's resistance	∆
			Newtonian	Non-Newtonian
			UKC = +10%
		Mud A	1.125±0.078 N	1.149±0.080 N	2.08%
		Mud C	1.328±0.093 N	1.356±0.095 N	2.06%
			UKC = -10%
		Mud A	2.302±0.160 N	2.415±0.170 N	4.76%
		Mud C	2.652±0.185 N	2.623±0.182 N	1.10%
	4.2	INFLUENCE OF MUD PROPERTIES ON SEABED UNDULATION AND FREE SURFACE ATTENUA-

ACKNOWLEDGEMENTS

The work presented in this paper was performed within the project "Gironde-XL" funded by "Le Grand Port Maritime de Bordeaux" and the European commission.


Alexandre Génin
email: alexandre.genin@umontpellier.fr.

Steven R Lee

Eric L Berlow

Steven M Ostoja

Sonia Kéfi

Mapping hotspots of potential ecosystem fragility using commonly available spatial data

Keywords: global changes, management, subalpine ecosystems, alternative states, perturbations, spatial data

Introduction

While many ecosystems respond gradually to environmental changes, others can exhibit ecological shifts, where disproportional changes in an ecosystem state occur following gradual changes in environmental conditions [START_REF] Scheffer | Catastrophic shifts in ecosystems[END_REF][START_REF] Suding | Alternative states and positive feedbacks in restoration ecology[END_REF]). This represents a global conservation concern, as a wide array of terrestrial and marine ecosystems have shown such large degradation events, including coral reefs [START_REF] Hughes | Catastrophes, Phase Shifts, and Large-Scale Degradation of a Caribbean Coral Reef[END_REF], fish stocks [START_REF] Sguotti | Catastrophic dynamics limit Atlantic cod recovery[END_REF], arid ecosystems [START_REF] Kéfi | Spatial vegetation patterns and imminent desertification in Mediterranean arid ecosystems[END_REF]), shallow lakes [START_REF] Scheffer | Alternative equilibria in shallow lakes[END_REF]) and kelp forests [START_REF] Steneck | Ecosystem Flips, Locks, and Feedbacks: the Lasting Effects of Fisheries on Maine's Kelp Forest Ecosystem[END_REF]. In general, because these systems have attracted a large research effort, the ecological mechanisms underpinning these shifts have been identified.

Yet, this knowledge does not always translate into actionable conservation directives because local factors can alter ecological dynamics in seemingly idiosyncratic ways [START_REF] Suding | Threshold models in restoration and conservation: a developing framework[END_REF]. For example, in rangelands, erosion and fire regimes can tip ecosystems between grass and shrub-dominated states [START_REF] Suding | Alternative states and positive feedbacks in restoration ecology[END_REF]). Yet, because such dynamics are altered by local factors (such as topography or grazing pressure), it remains difficult to pinpoint the ecological conditions at which a state shift will occur for a given piece of land [START_REF] Bestelmeyer | Threshold concepts and their use in rangeland management and restoration: the good, the bad, and the insidious[END_REF]. Such uncertainty can severely hinder effective conservation efforts, because the potential of some areas to undergo large changes and become degraded can be left overlooked. Conversely, degraded areas that would respond favorably to restoration may be left unattended. For informed land management in the current context of global changes, it is important to identify, a priori, how different zones in a given area of interest (e.g. a conservation area) may respond to perturbations. Rigorously and precisely documenting large changes after they happen is too late.

Approximately forecasting its potential to occur can help prioritize resources for proactive decisions.

To improve the practical forecast of large changes in ecosystems, ecologists have sought to build more operational frameworks and tools to identify from empirical data if a given ecosystem can exhibit shifts, and if so, for which environmental conditions [START_REF] Scheffer | Catastrophic regime shifts in ecosystems: linking theory to observation[END_REF][START_REF] Andersen | Ecological thresholds and regime shifts: approaches to identification[END_REF][START_REF] Suding | Threshold models in restoration and conservation: a developing framework[END_REF][START_REF] Ratajczak | Abrupt Change in Ecological Systems: Inference and Diagnosis[END_REF]). However, not all methods can always be applied to a focal ecosystem. Ideally, efforts to identify upcoming shifts would involve documenting temporal responses to changes in environmental conditions, or perturbations [START_REF] Bestelmeyer | Analysis of abrupt transitions in ecological systems[END_REF][START_REF] Petraitis | Multiple stable states in natural ecosystems[END_REF], and then applying predictions based on those observations to new locations where environmental changes have yet to occur. Such an approach relies on time-series from either observations or experiments [START_REF] Scheffer | Catastrophic regime shifts in ecosystems: linking theory to observation[END_REF]. However, for many locations, obtaining such datasets, or setting up experiments at the appropriate scales can be costly and require perturbation experiments that are often incompatible with management policies. Conservation areas more commonly have inventory data, where descriptive attributes (e.g. abundance of species, soil characteristics, etc.) are surveyed along known environmental gradients (effectively replacing temporal sampling by spatial sampling --a 'space-for-time' approach). These surveys are often not repeated, resulting in data collected from multiple locations along environmental gradients, but at a single or several close points in time (what we name hereafter 'snapshot' data). As an example, such data could correspond to plant surveys along altitudinal or aridity gradients during a single year. We show here how such data can help land managers prioritize their conservation efforts by predicting where perturbations may cause potentially large changes in an ecosystem.

Our work focuses on ecologically important subalpine meadows in the Sierra Nevada (California, USA).

Many of these meadows are located in two major national parks, Yosemite and Sequoia, and exemplify conservation areas where, due to their remoteness, it is often impractical to regularly monitor for temporal changes, and the use of snapshot data to gain insight on the potential effect of perturbations can be particularly useful. The small-scale species composition of meadow plant communities is known to be tightly linked to water availability [START_REF] Benedict | Plant Associations of Subalpine Meadows, Sequoia National Park, California[END_REF][START_REF] Allen-Diaz | Water Table and Plant Species Relationships in Sierra Nevada Meadows[END_REF][START_REF] Lowry | Groundwater controls on vegetation composition and patterning in mountain meadows[END_REF]). However, this general response can be altered, for example because of local changes in soil organic matter content [START_REF] Ankenbauer | The effects of soil organic matter on soil water retention and plant water use in a meadow of the Sierra Nevada, CA: Soil organic matter affects plant water use[END_REF], or because of a locally-increased importance of erosion processes [START_REF] Klikoff | Microenvironmental Influence on Vegetational Pattern near Timberline in the Central Sierra Nevada[END_REF]. Due to such local effects, some meadow areas may be more fragile than others to similar perturbations. We show in this study how a novel approach to the analysis of snapshot datasets may help better prioritize conservation and research efforts to protect these habitats.

Theoretically, three broad types of changes in a variable of interest (e.g. productivity) can be observed along an environmental gradient: linear, continuous non-linear or discontinuous non-linear with discrete alternative states (Box 1 panels a1, b1, c1) [START_REF] Ratajczak | Abrupt Change in Ecological Systems: Inference and Diagnosis[END_REF]. Such variable can be a population-level (e.g. abundance of a focal species), community-level (e.g. species composition) or a more ecosystemic attribute (e.g. productivity), depending on what is most relevant for management. When a linear response is observed, the observed state changes proportionally to changes in environmental conditions along the gradient (panel a1 in Box 1). In contrast, with a continuous non-linear response, small changes in conditions result in an apparent large shift as an ecological threshold is crossed (panel b1 in Box 1). A specific case arises when the response exhibits alternative states (panel c1 in Box 1), a case in which several discrete states can be observed for a range of environmental conditions along the gradient. These states can arise from differences in an environmental variable not captured by the main gradient (e.g. different local soil types), or because ecological processes make the ecosystem persist in one of several contrasted states despite similar environmental conditions. In this latter case, the discrete states are effectively alternative stable states, which can arise due to a variety of ecological mechanisms, such as priority effects, or reinforcing feedbacks [START_REF] Beisner | Alternative stable states in ecology[END_REF][START_REF] Petraitis | Multiple stable states in natural ecosystems[END_REF]. Regardless of the ecological origin of non-linear responses along a gradient, for practical conservation purposes, the areas where such responses are observed can be considered more fragile, either because they are susceptible to larger shifts than elsewhere, or because the community types in those areas strongly depend on a specific, narrower set of environmental conditions not reflected in the main gradient. Finding ways to identify where they occur is thus essential for informed management.

Previous studies have suggested approaches to identify the ranges of external conditions where discrete states are observed, based on detecting distinct ecosystem states in time series (Livina et al. 2010, Vasilakopoulos and[START_REF] Vasilakopoulos | Resilience and tipping points of an exploited fish population over six decades[END_REF] or ecosystem data along spatial gradients [START_REF] Hirota | Global Resilience of Tropical Forest and Savanna to Critical Transitions[END_REF], Scheffer et al. 2012). These approaches have mostly been used for cases where a single, well-known index (e.g., tree cover), characterizes the observed ecosystem state. When the latter cannot be adequately characterized by a single index, it is more challenging to visualize discontinuous changes in data [START_REF] Magurran | Long-term datasets in biodiversity research and monitoring: assessing change in ecological communities through time[END_REF], and fewer statistical tools are available to identify discrete states. Yet many cases of interest to both conservation and ecological theory require identifying discontinuous changes in complex multi-species communities. For example, when focusing on changes in species compositions, the observed state is accurately characterized by a set of species abundances, which cannot be reduced into a single index without losing information [START_REF] Barros | N-dimensional hypervolumes to study stability of complex ecosystems[END_REF]. We show here how non-linear responses can be detected along environmental gradients by applying such previous approaches to inventory data. We apply our work to detect areas in subalpine meadows that could be more fragile to upcoming perturbations, and discuss the implications of such fragility for the local conservation of habitats and species.

Box1: Identifying non-linear responses along a gradient from the modality of observations

Proof of concept model: to illustrate the link between multimodality in observations and the type of response of an ecosystem along a gradient, we can consider a simple model of ecosystem dynamics [START_REF] Noy-Meir | Stability of Grazing Systems: An Application of Predator-Prey Graphs[END_REF][START_REF] May | Thresholds and breakpoints in ecosystems with a multiplicity of stable states[END_REF], which can produce different types of deterministic ecosystem responses (red lines) along an environmental gradient. When a linear response is observed (a1), the average ecosystem state (which could here represent the vegetation biomass; red line, y-axis) changes proportionally to the changes in environmental conditions (x-axis). This is not the case for non-linear responses (b1, c1), where large changes in the observed state can arise following small changes in the environmental conditions captured by the gradient (blue arrows), or between alternative states after a perturbation (green arrow).

We modified the model to have similar dynamics, but subjected to typical sources of ecological noise. We retained the final state of a set of independent simulations with varying parameters to simulate a sampling scheme that would lead to snapshot data (i.e. sampling independent instances of the same system in different environmental conditions, but at a single point in time). This dataset is represented by the black points in panels a1, b1, c1, each point corresponding to one observation. Using a moving-window approach, we can investigate how the distribution of these observations changes along the grazing gradient, and thereby characterize how their 'modality' relates to the type of ecosystem response (panels a2-c2). Wherever the observations tend to exhibit two modes along the gradient, their distribution is better described by a density with two modes, and thus the frequency of counts for two modes increases (see Methods for details about this approach). As this arises in the vicinity of non-linear responses (panels b2 and c2), modality can be used as a criterion to identify ranges of environmental conditions where apparent non-linear responses occur.

It is important to note that a non-linear response along a gradient does not mean that intrinsic ecological processes drive the response: external variables not captured by the main gradient could also explain such a pattern (e.g. different soil types producing alternative states). More information about the model is provided in Supplementary Material S2.

Methods

Data source

We used a dataset resulting from a study documenting meadow plant communities in 47 meadows located in Yosemite and Sequoia National Parks in 2011 and 2012 [START_REF] Lee | A multiscale evaluation of pack stock effects on subalpine meadow plant communities in the Sierra Nevada[END_REF]. Meadows were situated between 2517m and 3355m a.s.l. and had an area between 0.42ha and 22.1ha. Data were collected for each meadow from 2x2m plots, regularly spaced (20 m intervals) along belt transects that were oriented perpendicularly to the main meadow drainage, spaced 40 m apart and ran from meadow edge to meadow edge. Volumetric water content of the top 12 cm of soil was obtained for all plots using a handheld TDR soil moisture probe. In every third plot, ocular cover estimates of all plant species were recorded in 8 25x25cm sub-plots, and averaged to estimate the relative abundance of each species in the plot. In addition, ocular estimates of the cover of 28 variables describing the local water regime, upper soil characteristics and disturbance by herbivores were measured (e.g. cover of bare ground, rock bed, etc., see Supplementary Material S1). Overall, sampling resulted in 2860 plots with soil moisture and local environmental conditions; of these, 1287 also had species composition recorded.

We focused on the variations of a community-level attribute, the species composition of vegetation communities, along the natural moisture gradient in the meadows. Decades of work show that moisture is a key factor in determining meadow species composition, as within a meadow, plant species sort along the gradient (Allen-Diaz 1991, [START_REF] Lowry | Groundwater controls on vegetation composition and patterning in mountain meadows[END_REF]. We checked that this aspect of meadows was also captured by our dataset (Supplementary Material S2, section 2), and used water content as the main gradient. We attenuated the effect of year-to-year variations in precipitation on the measurements of the moisture gradient by standardizing the measured water content within year. As a result, in what follows, water content values are given in standard deviations (s.d.) to the annual mean.

Analyzing changes in modality along a gradient

When a set of observations is unimodal, observations tend to cluster around a single mean. For example, a unimodal set of observations of a given species' abundance can be effectively characterized by a mean abundance, from which most observations would deviate only slightly (Figure 1a). When a set of observations is multi-modal, observations cannot be well-summarized by a single mean and instead, several, different clusters are present in the data. This could for instance be the case when a species is either locally very rare or very abundant (these abundances being modes, Figure 1b), but rarely observed with intermediate abundances.

Characterizing the multimodality for a set of observations can be done using the density of observations. The density is a function that captures the relative frequency of different observed values: those that are more frequent in a given set of observations will have higher density values. Conversely, values that are scarcely observed correspond to lower density values. This results in local peaks ('modes') in the density around values that are frequently observed in a given set of observations. A unimodal set of observations will tend to have a single peak in its density, i.e. most observations cluster around a single, average value (Figure 1 a,c). A multi-modal set of observations will have a density with multiple well-defined peaks (Figure 1 b,d). Computing the density of a given set of observations, and estimating the number of wellseparated peaks in it is thus a natural way to characterize the number of modes in a set of observations (its modality, [START_REF] Silverman | Using Kernel Density Estimates to Investigate Multimodality[END_REF], [START_REF] Müller | Excess Mass Estimates and Tests for Multimodality[END_REF]).

For each subset of observations, we computed its density, and considered n cut thresholds regularly spaced between zero and the maximum value of the density. For each threshold, we considered a "horizontal slice" of the density and counted the number of observed disjunct parts of the density above the threshold (λ-clusters sensu [START_REF] Müller | Excess Mass Estimates and Tests for Multimodality[END_REF], see Figure 1). This number of clusters provides an estimate of the number of well-separated modes in a given density above a certain threshold [START_REF] Müller | Excess Mass Estimates and Tests for Multimodality[END_REF][START_REF] Fisher | Mode testing via the excess mass estimate[END_REF][START_REF] Scott | Multivariate density estimation: theory, practice, and visualization[END_REF]. For a given density, i.e. a set of observations, this process yields a set of n cut integer values, each corresponding to an estimated number of modes. How these counts vary along an environmental gradient informs about changes in the modality of the distribution of states. A typical example could be that of a density switching from unimodal to bimodal around a certain value of the environmental gradient, which would be reflected in a higher number of counts for two well-separated modes.

This non-parametric method can be applied to univariate densities (Figure 1, top; Box 1), which corresponds to the case where the state of the system is defined by a single value. It can also be applied to multivariate densities, obtained through multi-dimensional kernel density estimation (Figure 1b), which makes it applicable to cases where each observation is multi-dimensional (e.g. when each observation consists of abundances of individual species). It only depends on a single parameter (the density bandwidth) and it is more robust to outliers than previous approaches for the detection of multimodality [START_REF] Fisher | Mode testing via the excess mass estimate[END_REF]. n cut only affects the precision of computations, and was set to a large value to identify differences in height of the density (we used 512 here, which was enough to detect changes in modality along the gradient).

Given a fixed set of observations, this procedure allows estimating the number of modes in the density, but it needs to be adapted to investigate change in this number along a gradient. To do so, we adopted a moving-window approach. We chose a window size w, here corresponding to a range of water content values, and estimated the number of modes for each subset of observations falling within this range, each time 'moving' the window by a small amount along the gradient. This allows determining how the number of modes (i.e. the modality) of observations change along an environmental gradient. that are described by a single variable such as the abundance of a single species), the density can be represented as a one-dimensional function of that variable (a1, b1). This function may have one peak, and thus be unimodal (a1) or two and be bimodal (b1). When an observation is multivariate (c1, d1), i.e. when two or more variables are used for each observation, such as abundances of several species, the density is represented as a two-dimensional plot (c1, d1) along these variables, where red values represent high density values. Again, multivariate densities can be unimodal (c1) or multi-modal (d1). To estimate the number of modes, we can consider the number of disjunct parts of a given density with values above a certain threshold (dashed lines, a1, b1 and continuous lines in c1 and d1). Doing so for a range of thresholds between zero and the maximum value of the density yields a majority of counts for one mode for unimodal densities (red area in panels a2, c2). For bimodal densities, the same analysis yields an increased number of counts for two modes (yellow areas, b2, d2). Note that very low thresholds identify spurious clusters because of outliers in the density, but this effect only concerns a reduced number of threshold values, so the influence of outliers on this estimate of modality is minimal. Data was generated by drawing samples out of a univariate or multivariate Gaussian distribution (a,c), or a mixture of two Gaussian distributions (b,d; see Supplementary Material S2 section 5 for the distribution parameters).

Empirical data analyses

Classical univariate synthetic indices (e.g. total species richness or total cover) fail to capture the response of meadow vegetation communities to environmental parameters [START_REF] Lee | A multiscale evaluation of pack stock effects on subalpine meadow plant communities in the Sierra Nevada[END_REF]. We thus defined the state of meadow vegetation communities as the set of relative abundances of all species in a given plot [START_REF] Barros | N-dimensional hypervolumes to study stability of complex ecosystems[END_REF]. Each single observation had thus a dimension equal to the total number of species in the dataset (257). Because carrying kernel density estimation cannot be done on data with this many dimensions (typically, computational costs become too high after 5-10 dimensions [START_REF] Blonder | The n-dimensional hypervolume[END_REF], we summarized the variations of vegetation species compositions to a few axes using an ordination. In principle, using Principal Coordinate Analysis (also named metric Multi-Dimensional Scaling) is the method of choice as it preserves dissimilarities between sites in the resulting axes. However, empirically it often produces strong arch artefacts on vegetation data, in particular when species turnover is high along a gradient [START_REF] Legendre | Numerical ecology[END_REF]. Non-metric multidimensional scaling could be considered as an alternative but, because it preserves only the ranks of the original dissimilarities between sites, the resulting axes summarizing species turnover cannot be used quantitatively. We therefore chose a Correspondence-analysis-based method -Detrended-Correspondence Analysis -which (i) produces axes that can be used as quantitative variables and (ii) is less subject to arch-like artifacts. We carried out sensitivity analyses to confirm that our results were robust to the choice of dimensionality reduction method (Supplementary Material S2, section 4). We used the two first axes of the DCA and measured changes in the modality of observations along the moisture gradient using the previously-described approach. We set the window size to 1 and used the bivariate diagonal density estimation kernel [0.70, 0; 0, 0.62], as obtained from automatic bandwidth selection (function bw.nrd in R). Because the sample size affects the shape of the density (and thus the counts of modes), we always computed densities using a fixed number of observations (150). Because there were often more samples falling within a window of water content, we repeated the analysis 100 times for each window, using a random subsample of 150 observations. Based on results from the previous analysis, patterns identified two ranges of water content with increased support for two well-separated modes, dry areas and wet areas. We extracted these subsets from the dataset and investigated soil attributes that could underpin their more bimodal distributions. We refer to dry areas for sites where water content was below -1 s.d. and wet areas for those where water content was above 1 s.d. Because the multimodality analysis suggested a dominance of two discrete plant community types, we clustered each subset into two groups based on species composition, using hierarchical clustering with a Ward distance on a Bray-Curtis distance matrix. We contrasted the environmental covariates of each group to discuss possible underlying mechanisms for this apparent non-linear response along the water-content gradient. We used two-sample Mann-Whitney tests (wilcox.test in R) to contrast the distribution of each plot's characteristics across groups, adjusting p-values for multiple tests using the Benjamini-Hochberg correction.

Because the sampling was regularly spaced in meadows, we computed for each meadow an estimate of the proportion of area with increased multimodality as the proportion of plots laid in that meadow with a measured standardized water content above 1 or below -1. To map in space where this increase of modality occurred, we interpolated the standardized soil water content within meadows over a regular grid of points using ordinary kriging (using a spherical variogram for each meadow separately). We then modelled the empirical relationship between water content and the proportion of counts for more than one cluster as a one-dimensional generalized additive model. We used this empirical relationship to infer from the interpolated water content the proportion of counts for more than one cluster, and used it as an index of meadow areas with apparent non-linear responses (Supplementary Material S2, section 6). Analyses were conducted in R (version 3.6.1, R Core Team 2019), along with the R package 'vegan' v2.5 [START_REF] Oksanen | vegan: Community Ecology Package[END_REF], 'gstat' v2.0.3 [START_REF] Pebesma | Multivariable geostatistics in S: the gstat package[END_REF]) and 'mgcv' v1.8.28 [START_REF] Wood | Stable and efficient multiple smoothing parameter estimation for generalized additive models[END_REF]).

Results

We used a proof-of-concept model (Box 1) to investigate how the modality of observations, i.e. how many peaks in the distribution are identified, changes along gradients for different types of ecosystem responses: linear, non-linear continuous and discontinuous. Simulations show that, in the case of a single ecosystem state along a gradual/linear ecosystem response, the distribution is predominantly characterized by a single mode, while multimodality arises when the state changes non-linearly with environmental conditions (Box 1, panels b2, c2).

The vegetation composition of the sites of our data set was summarized to the two first axis of a Detrended-Correspondence Analysis (Figure 2). The analysis of changes in modality in DCA scores along the moisture gradient highlighted ranges of the moisture gradient with an increased modality in the species composition of the meadows (Figure 2): dry areas and wet areas. In these two ranges of the water gradient, we found an increased count for two well-defined modes (Figure 3b), highlighting an apparent non-linear response of species composition to changes in water content. In wet areas (plots with water content above 1 s.d.), the two groups defined based on species composition were contrasted. One was dominated by Carex vesicaria (CARVES), while the other was a more mesic community that included Oreostemma alpigenum (Figure 4a, b-top). The C. vesicaria-dominated group represented 25% of all the plots in wet areas and had different water regime-related attributes: a more saturated and less moist soil, as well as a higher observed cover of silt and bare ground (Figure 4b). The covers of dry and flooded soils were significantly different between the two groups but the estimated effect size was close to zero.

In dry areas (plots with water content below -1 s.d.), one of the two groups was dominated by Carex filifolia (CARFIL) (10% of plots in dry areas), while the other was dominated by more mesic plants (e.g.

Vaccinium caespitosum (VACCAE); Calmagrostis brewerii (CALBRE)

; Figure 5a). The C. filifolia-dominated group had a higher ocular cover of coarse particles (sand, pebble) and bare ground (Figure 5b). It also showed a significant difference in water-regime attributes (cover of dry and moist soil), and disturbance attributes (stock prints, rodent mounds), although they had a very small estimated cover difference (below 10 -4 ).

Overall, the proportion of area with increased multimodality in meadows represented on average 22% of a given meadow area of (across surveyed meadows, this number reached from 0 to 77%) and could be identified on maps based on the spatial interpolation of soil moisture (Figure 6).

Discussion

Our analysis revealed that in wet and dry areas of meadows, there was no smooth response of vegetation composition to the moisture gradient but instead, discrete vegetation types were present. These types were associated with contrasted local soil characteristics. As a result, current and upcoming perturbations may have an unexpected effect in these areas because general changes in moisture regime may alter the specific environmental factors on which vegetation states depend. Reviewing the factors that co-varied with species composition (Figure 4 and 5) helps identify specific environmental variables and potential ecological processes that are critical to the conservation of these habitats, and provide informed possible scenarios for their response to current perturbations in the Sierra Nevada.

Flagging potentially fragile meadow regions.

In wet areas of the meadows, a Carex-vesicaria-dominated community contrasted with a more diverse mesic-wet meadow community. While both communities co-occurred in regions of similar average surface soil moisture (as measured by 12 cm TDR probes); the C. vesicaria-dominated community had a larger amount of flooded or saturated soil and more silt, which are typical characteristics of streamside flooded pools where Carex spp. dominate [START_REF] Baldwin | The Jepson manual: vascular plants of California[END_REF]. The mesic-wet community had a more diverse species composition, and occurred more frequently outside of areas with standing water.

Here, the discrete vegetation states seemed to depend on the absence or presence of standing water, despite the fact that average, empirically measured surface soil moisture was similar. Local changes in water regimes could thus have a strong effect on species composition -or vice versa. These local changes can arise following changes in regional climate (see below), but may also be caused by local factors.

Grazing for example has been shown to increase run-off through compaction of the soil [START_REF] Kauffman | Livestock Impacts on Riparian Ecosystems and Streamside Management Implications[END_REF]Krueger 1984, Ostoja et al. 2014), which could diminish the probability of creation of areas of standing water in meadows. However, given the overall low levels of grazing, and the relatively small importance of stock disturbance in our dataset (Figure 4, 5), their potential impacts are most likely minor within parks compared to regional-scale perturbations [START_REF] Holmquist | Patch-Scale Effects of Equine Disturbance on Arthropod Assemblages and Vegetation Structure in Subalpine Wetlands[END_REF][START_REF] Lee | A multiscale evaluation of pack stock effects on subalpine meadow plant communities in the Sierra Nevada[END_REF].

In dry meadow areas, our analysis identified a Carex filifolia-dominated community contrasting against a more diverse mesic community which both occur in areas of similar empirically measured soil moisture.

The C. filifolia community was observed to have coarser upper soil particles (sand, pebble), and a higher amount of bare ground. The two vegetation states differed significantly in visible disturbance by herbivores and water regime. This difference is consistent with the fact that dry communities are known to respond to external parameters such as rodent disturbance, erosion processes and soil nutrients [START_REF] Klikoff | Microenvironmental Influence on Vegetational Pattern near Timberline in the Central Sierra Nevada[END_REF]. Differences in local soil surface conditions that are not captured by the main gradient could thus explain the presence of different discrete states. However, these local soil conditions are likely tied to the vegetation state as plant-soil feedbacks have been shown to exist in these areas. For example, sods formed by mesic meadow species can improve local microhabitat conditions by trapping smaller soil particles which hold nutrients and water during the dry season (Wood 1975, Ankenbauer and[START_REF] Ankenbauer | The effects of soil organic matter on soil water retention and plant water use in a meadow of the Sierra Nevada, CA: Soil organic matter affects plant water use[END_REF]. These plant-soil feedbacks could be strong enough to impair a rapid regrowth of the vegetation after soil disturbance, so that an acute disturbance of dry areas could have much longer-term effects than expected (e.g. through trampling). It is unknown whether such feedbacks could be strong enough to impair the vegetation regrowth completely following disturbance (effectively underpinning alternative stable states).

However, the evidence as a whole highlights the need to adopt a conservation perspective that takes into account the possibility of long-term, poorly-reversible degradation: based on the available evidence, mesic communities in dry areas are likely to be very fragile to upcoming changes in conditions and recovering from acute degradation following an extreme event would likely take decades, if possible at all without restoration efforts [START_REF] Ratliff | Restoring plant cover on high-elevation gravel areas, Sequoia National Park, California[END_REF].

More broadly, in both wet and dry areas, soil water-related parameters stood out as possible drivers of vegetation shifts. With the increase in frequency and intensity of drought events associated with climate change, snow to rain ratios and total spring snowpack are decreasing across the Sierra Nevada [START_REF] Knowles | Trends in Snowfall versus Rainfall in the Western United States[END_REF][START_REF] Barnett | Human-Induced Changes in the Hydrology of the Western United States[END_REF]). These long-term changes in climate dynamics have the potential to extensively alter meadow hydrology, to which the wet and dry plant communities appear to be the most sensitive. As these communities represent on average 20% of meadow areas (up to 77%), this could constitute a significant concern for subalpine meadows across the central and southern Sierra Nevada. The Sierra Nevada has recently experienced an intense drought event [START_REF] Belmecheri | Multi-century evaluation of Sierra Nevada snowpack[END_REF], which included the four driest years (2012)(2013)(2014)(2015) of the last 2000 years [START_REF] Adams | Annually resolved late Holocene paleohydrology of the southern Sierra Nevada and Tulare Lake, California[END_REF]. While the repercussions of the drought have clearly manifested across Sierra Nevada forest ecosystems [START_REF] Potter | Landsat Image Analysis of Tree Mortality in the Southern Sierra Nevada Region of California during the 2013-2015 Drought[END_REF], the impacts on meadow ecosystems have yet to be thoroughly evaluated. Monitoring changes and setting up experimental research in potentially fragile areas occur could inform us greatly about such events. Our approach allows planning this next step by creating maps identifying these areas.

A workflow to identify and map fragile areas

With this work, we showed that investigating snapshot data for changes in modality of ecosystem states can provide information on the potential fragility of specific meadow areas. Importantly, our results show that we may identify these areas more accurately by taking into account the full state of the community (here, species composition) as we would not have detected changes in modality based on simple indices such as richness or total plant cover (Supplementary Material S2). It is a relatively simple analysis that can be carried out on commonly-available inventory data to investigate the possible dynamics of the ecological system under focus. An example workflow could be the following:

1. Gather a dataset of observations describing how the state of an ecosystem varies along known influential environmental gradients. This can be done over a relatively short period of time, but the sample size of observations must be sufficiently large (at least 10-50 points per level of the environmental gradient). If no prior knowledge is available to identify main environmental gradients, this can be done based on redundancy analysis or canonical correspondence analysis (Supplementary Material S2, section 2).

2. When the state is described using a high number of variables (e.g. species abundances), reduce the dimensionality of the dataset (e.g. using Detrended Correspondence Analysis) to 1-5 dimensions (variables). It is noteworthy that reducing the dataset to 2 dimensions or less drastically speeds up the numerical computations.

3. Analyze the dataset for changes in modality along the main gradient, to identify areas along the gradient with increased modality, i.e. where apparent non-linear responses are present 4. For these areas, use local environmental data to identify the source of increased multimodality.

Such increase may arise as a combination of (1) effects of secondary environmental drivers of plant community composition (e.g. different soil types) and (2) intrinsic ecological processes driving non-linear responses (e.g. reinforcing feedbacks, [START_REF] Wilson | Positive-feedback Switches in Plant Communities[END_REF]. For conservation purposes, in both cases, these areas may be considered fragile because the local state (e.g. species composition) depends on a narrower set of environmental conditions and/or the response to perturbations may be non-linear (possibly exhibiting ecological shifts).

5. Use spatial data to map where these areas are present, i.e. where the effect of perturbations may be larger. It is important to note that identifying the ranges of environmental conditions where non-linear responses occur only requires a reduced subset of data. For example, in our work, only 45% of the plots had a recorded species composition. However, the measurement of water content was carried out for a much higher number of plots, which allowed for higher resolution in our maps of community fragility.

The measurement of modality is an open area of statistical research [START_REF] Scott | Multivariate density estimation: theory, practice, and visualization[END_REF], and other approaches, not based on the use of a density, could be useful to evaluate multimodality in a set of observations [START_REF] Clarke | Testing of null hypotheses in exploratory community analyses: similarity profiles and biota-environment linkage[END_REF]) and determine the statistical significance of patterns. Computing multivariate densities can be computation intensive, although approximate numerical methods may help [START_REF] Blonder | The n-dimensional hypervolume[END_REF].

Comparing different methods is beyond the scope of this article, but it is important to carry out sensitivity analyses to test for the robustness of the result (e.g. to the bandwidth used, or the method of ordination to reduce dimensionality, Supplementary Information S2, section 4).

Detecting large ecological changes for conservation

The observation of discrete community states does not mean that a given ecosystem necessarily exhibits irreversible shifts between alternative stable states. For example, discrete states can be also due to changes in an environmental parameter not captured by the ones measured. This is probably the case for the different community states observed on the wetter end of the gradient in the data set studied here.

Nonetheless, complementing this approach with other sources of evidence can suggest where non-linear and/or poorly-reversible degradation may be possible in a given ecosystem of interest. Such additional knowledge can come from previously-known ecological mechanisms which drive irreversible shifts (e.g. known reinforcing feedbacks between vegetation and soil quality like here), historical information about shifts that occurred in the past and appear irreversible, or evidence from independent approaches (e.g. experiments, or based on indicators of ecosystem shifts; [START_REF] Dakos | Methods for Detecting Early Warnings of Critical Transitions in Time Series Illustrated Using Simulated Ecological Data[END_REF][START_REF] Kéfi | Early Warning Signals of Ecological Transitions: Methods for Spatial Patterns[END_REF][START_REF] Nijp | Spatial early warning signals for impending regime shifts: A practical framework for application in real-world landscapes[END_REF]. While no approach by itself is likely to prove the existence of irreversible shifts, from a conservation perspective it is particularly important to consider this possibility when the alternative community is considered to be degraded (e.g. because it is poorer in species or because it represents a loss of conservation-critical habitat). It is better to imprecisely predict a poorly-reversible shift than to precisely document it after it happens.

As the conditions for the emergence of alternative stable states are quite restrictive, ecosystems probably respond more often to gradual changes in a gradual way than in an abrupt way (e.g. because of spatial heterogeneity that favors gradual responses (van Nes and Scheffer 2005)). Cases of irreversible shifts have however been shown to be possible in a broad range of ecosystems and in a response to a number of drivers (lake eutrophication due to increase nutrient loading, coral reef degradation due to warming and/or water pollution, dryland desertification due to increased grazing and/or droughts). Because abrupt ecosystem responses can happen in about any ecosystems and their occurrence depends on a joint set of conditions which are not straightforward to monitor, it is extremely useful to find approaches that allow rapid predictions of where they are most likely to happen. The approach proposed here is promising because it provides a way to do so for species-rich communities (where most theoretical developments so far focus on species-poor models), and because it can be based on data typically gathered in conservation/protected areas.

Figure 1 :

 1 Figure 1: Principle of the measure of modality based on synthetic data. With univariate observations (i.e.

Figure 2 :

 2 Figure 2: Detrended Correspondence Analysis results on the subalpine dataset (a). The distribution of each point in the two axes is based on similarity of species composition (two close points tend to have a more similar composition). Species labels indicate points (plots) where a given species is the most abundant, with font size being proportional to the overall abundance of the species in the dataset (see supplementary materialS2 for a list of species codes). The bottom panels (b1-4) show the same graph, but only with points comprised within different ranges of water content. Increased modality (i.e. multimodality) can be observed for water contents below -1 s.d. and above +1 s.d.

Figure 3 :

 3 Figure 3: (a) Moving averages (window width: 0.4 s.d.) of the relative covers of the ten most abundant species in the subalpine meadow dataset along the gradient of soil water content. (b) results of modality analysis along the same gradient. (c) Changes in the first axis values of the Detrended Correspondence Analysis (DCA) along the water moisture gradient (each point corresponds to a summarized species composition), and density contour highlighting local modes.

Figure 4 :

 4 Figure 4: (a) Detrended Correspondence Analysis results on the wet meadow subset. Species labels indicate where individual species are dominant, with label size proportional to the total abundance of the species in the subset (see supplementary material S2). A density estimate is added on top of the points. (b) Differences in species composition, number of species S, and in plot attributes between the two groups (only attributes with significant differences are shown). Stars represent the significance of a Mann-Whitney test (one star for P < 0.05 and two for P < 0.01). The grey bars show the estimated difference in medians between the two groups.

Figure 5 :

 5 Figure 5: (a) Detrended Correspondence Analysis results on the dry meadow subset. Species labels indicate where individual species are dominant, with label size proportional to the total abundance of the species in the subset. (b) Differences in species composition, number of species S, and in plot attributes between the two groups (only attributes with significant differences are shown), stars represent the significance of a Mann-Whitney test (one star for P < 0.05 and two for P < 0.01). The grey bars show the estimated difference in medians between the two groups.

Figure 6 :

 6 Figure 6: Example of spatial use of a multimodality index (the empirical support, between zero and one, for more than one community type for a given soil moisture level, see Methods). Zones in red represent areas where non-linear responses to changes in soil water content are present, and could thus be more fragile to perturbations.