David Pobel

Joël Robin

Jean-François F Humbert

Influence of sampling strategies on the monitoring of cyanobacteria in shallow lakes: lessons from a case study in France

Keywords: sampling strategy, cyanobacteria, spatiotemporal dynamic, Microcystis aeruginosa, Aphanizomenon flos-aquae

Sampling cyanobacteria in freshwater ecosystems is a crucial aspect of monitoring programs in both basic and applied research. Despite this, few papers have dealt with this aspect, and a high proportion of cyanobacteria monitoring programs are still based on monthly or twicemonthly water sampling, usually performed at a single location. In this study, we conducted high frequency spatial and temporal water sampling in a small eutrophic shallow lake that experiences cyanobacterial blooms every year. We demonstrate that the spatial and temporal aspects of the sampling strategy had a considerable impact on the findings of cyanobacteria monitoring in this lake. In particular, two peaks of Aphanizomenon flos-aquae cell abundances were usually not picked up by the various temporal sampling strategies tested. In contrast, sampling once a month was sufficient to provide a good overall estimation of the population dynamics of Microcystis aeruginosa. The spatial frequency of sampling was also important, and the choice in the location of the sampling points around the lake was very important if only two or three sampling points were used. When four or five sampling points were used, this reduced the impact of the choice of the location of the sampling points, and allowed to obtain fairly similar results than when six sampling points were used. These findings demonstrate the importance of the sampling strategy in cyanobacteria monitoring, and the fact that it is impossible to propose a single universal sampling strategy that is appropriate for all freshwater ecosystems and also for all cyanobacteria.

Introduction

Due to eutrophication and, to a lesser extent, to climatic changes [START_REF] Markensten | Simulated lake phytoplankton composition shifts toward cyanobacteria dominance in a future warmer climate[END_REF][START_REF] Paerl | Climate change : a catalyst for global expansion of harmful cyanobacterial blooms[END_REF] cyanobacterial blooms seem to be increasing in freshwater ecosystems worldwide. These blooms severely disrupt the functioning of these ecosystems and potential water use. Furthermore, many cyanobacterial species are able to produce a variety of toxic metabolites, which can be harmful to both human [START_REF] Kuiper-Goodman | Human health aspects[END_REF] and animal [START_REF] Codd | Harmful cyanobacteria[END_REF] health. For these reasons, numerous attempts have been made in the last 20 years to elucidate the factors that control cyanobacterial blooms and toxin production, and thus to make it possible to evaluate better the health risks associated with bloom events. From all these studies, it is clear that the spatial distribution of cyanobacteria in freshwater ecosystems can display marked horizontal and vertical variations [START_REF] Porat | Diel buoyancy changes by the Cyanobacterium Aphanizomenon ovalisporum from a shallow reservoir[END_REF][START_REF] Welker | Toxic Microcystis in shallow lakes Müggelsee (Germany) -dynamics, distribution, diversity[END_REF]. Moreover, by means of a real time PCR analysis of a gene involved in the biosynthesis of microcystins we have shown that considerable fluctuations can also occur in the proportions of potentially microcystin-producing and non-producing cells during the course of Microcystis aeruginosa blooms [START_REF] Briand | Spatiotemporal changes in the genetic diversity of a bloom-forming Microcystis aeruginosa (cyanobacteria) population[END_REF]. Similar results have been found for various different M. aeruginosa populations located in the same geographic area [START_REF] Sabart | Spatiotemporal changes in the genetic diversity in French bloom-forming populations of the toxic cyanobacteria Microcystis aeruginosa[END_REF], which makes it difficult to manage the health risks associated with these events.

All these studies indicate that the sampling strategy used for monitoring cyanobacteria is a critical aspect, both in basic research on cyanobacteria, (e.g. investigation of the factors and processes involved in the development of the blooms), and in applied research, (e.g. implementing monitoring programs of these microorganisms in freshwater ecosystems used to provide drinking water or for recreational activities). In recent years, new tools have been tested with the intention of improving cyanobacterial sampling, for example, remote sensing reconnaissance to determine the horizontal distribution of cyanobacteria in freshwater ecosystems [START_REF] Hunter | Using remote sensing to aid the assessment of Human health risks from blooms of potentially toxic Cyanobacteria[END_REF], or spectrofluorometric probes to reveal the vertical distribution of these cyanobacteria in the water column [START_REF] Leboulanger | Application of a submersible spectrofluorometer for rapid monitoring of freshwater cyanbacterial blooms : a case study[END_REF]. Moreover, these spectrofluorometric probes and other sensors have now been integrated into buoys, to provide real-time monitoring of cyanobacteria in freshwater ecosystems (Le Vu et al., in press).

However, despite the great potential interest of these tools, their cost will remain prohibitive for their routine use in the foreseeable future, and most of the monitoring programs worldwide for the survey of cyanobacteria will continue to be based on more conventional methods for some years to come. Taking discrete samples of various volumes of water taken from the shoreline of ecosystems is probably the method one most often used in studies. Unfortunately, as a result of spatial and temporal differences in the distribution of cyanobacteria, this approach can often provide a very poor estimation of cyanobacterial abundance and, consequently, of the associated health risk. We therefore need to devise simple sampling strategies for the low cost monitoring of cyanobacteria in shallow lakes. In an attempt to do this, we performed intense spatiotemporal monitoring of cyanobacteria in a shallow lake known to experience cyanobacterial blooms every year.

Materials & Methods

Study site

This study was performed in a shallow lake named Place (0.08 km 2 , 2.5 m max depth, 45°43'N, 4°14'E) located in the plain of Forez (Central France), (Fig. 1). This lake is used for extensive fish production and its trophic status is eutrophic to hypereutrophic (OCDE, 1982).

Microcystis aeruginosa blooms occur every summer.

Data acquisition

Sampling strategy and cell counting

In order to assess the variations in the horizontal distribution of cyanobacteria in this pond, we monitored six sampling points located around the lake at one meter from the shore (V1-V6; Fig. 1). The water depth in each of these sampling points was around 1 meter. Samples were taken every two days, between 09:00 and 10:00 a.m., from early June 2008 to early October 2008. The first 40-centimeters of the water column were sampled using a watersampler (Uwitech, Austria). This water sample was shacked and then divided into two 1-L bottles, one liter being stored at room temperature with Lugol's iodine solution, and the other at 4°C.

In order to evaluate the vertical distribution of cyanobacteria in the water column, we performed a 22-hour survey (from 4:30 p.m. August 4, 2009 to 2:30 p.m. August 5, 2009), of the variations of cyanobacterial biomass at five sampling points (A-E; Fig. 1) using a BBE Algaetorch (Moldaenke, Germany). This torch is based on the same principle as the BBE spectrofluorometric probe [START_REF] Beutler | A fluorometric method for the differentiation of algal populations in vivo and in situ[END_REF], but provides only an estimation of the concentrations of cyanobacteria and total chlorophyll in water. Every hour, the torch was immersed to a depth of 20 centimeters at the five sampling points, and triplicate measurements were performed in each point.

The cyanobacterial cell concentrations were estimated using a Nageotte cell and an optical microscope, as described in [START_REF] Brient | A phycocyanin probe as a tool for monitoring cyanobacteria in freshwater bodies[END_REF]. For each rectangular area, we counted at least 400 cells of each cyanobacterial species.

Meteorological data

The speed and direction of wind during our study were obtained from the Metéo France meteorological station at St Etienne-Bouthéon (4°18'E -45°32'N). The wind direction rose for this station is given in Supplemental Figure 1, and shows that the two dominant wind directions were NW and SE. The direction of winds blowing from 240° -60° was classified as NW, and that of winds blowing from 60° -240° as SE.

Data analysis

The spatial distribution of cyanobacteria in the lake was represented using Surfer (v. 7.0, Golden Software Inc.), and statistical analyses (Wilcoxon test, Spearman correlation) were performed using the R package version 2.10 (R development core Team, 2010).

Results

Change over time in the population dynamics of the two dominant cyanobacterial species

Two cyanobacterial species, Microcystis aeruginosa and Aphanizomenon flos-aquae, dominated the phytoplankton community during the summer of 2008. The population dynamics of these two species displayed very contrasting patterns (Fig. 2). The population dynamics of Microcystis aeruginosa was characterized by a steady increase in the cell abundance from June to August, apart from a brief dip in the middle of July. The maximum population was reached on August 21 (264,000 cells/mL), and subsequently the cell concentration remained stable until the end of September, and then decreased in October. In contrast, the population dynamics of Aphanizomenon flos-aquae were much more chaotic, with the cell abundance reaching two very high and short-lived peaks in July (400,000 cells/mL on July 17, and 560,000 cells/mL on July 23).

Influence of sampling frequency on the estimation of the population dynamics

Our assessment of the changing population dynamics of the two cyanobacteria were obtained using a very frequent high temporal sampling regime (every two days), which would not be practicable in the context of normal monitoring programs. In order to evaluate the impact of the sampling frequency, we simulated weekly, twice-monthly and monthly sampling frequencies to our data set. The results of these simulations are shown in Fig. 3 and4. From this figure, we can see that changes in M. aeruginosa cell abundance over time would have been fairly accurately estimated at all these sampling frequencies. Moreover, for all sampling frequencies, the quality of the estimation of the M. aeruginosa population dynamics was not influenced by choice of the first sampling date (Fig. 3). In contrast, the population dynamics of A. flos-aquae would have been badly or even very badly estimated by using weekly, twicemonthly and monthly sampling frequencies (Fig. 4). We would only have detected both A. flos-aquae peaks in one of the three trials testing the weekly sampling strategy, and we would never have detected these peaks with twice-monthly and monthly sampling frequencies.

Evolution of the horizontal distribution of cyanobacteria in the lake during the bloom

As shown in the video (Supplemental Fig. 2), the horizontal distribution of both cyanobacteria displayed marked variations during the course of the study. Moreover, when the spatial distributions of the two species at the same sampling dates were compared, it could be seen that similar or contrasting patterns in the horizontal distribution of M. aeruginosa and A. flosaquae cells would have been found, depending on the dates chosen (some examples are provided in Fig. 5).

In order to obtain a better picture of this spatial variability in the cell concentrations of the two species, we estimated the coefficients of variation in the mean abundance for each sampling date and for each species from the results obtained at the six sampling points (Fig. 6). These coefficients were usually higher for A. flos-aquae than for M. aeruginosa (Wilcoxon test, p=-3.25.10 -05 ), suggesting that the horizontal distribution of A. flos-aquae was more variable. Finally, there was no correlation (Spearman coefficient) between the coefficient of variation and the mean cell abundance for Aphanizomenon, and only a weak correlation was found for Microcystis (Spearman coefficient, p=0.003 r =-0.4; Supplemental Fig. 3).

In order to find out whether wind speed/direction could account for the variations in the horizontal distribution of cyanobacterial cell abundance in the lake, we recorded in a first time, for each species and for each sampling date, the sampling point (out of the six) at which the highest cell abundance was detected. We then constructed a table in which we related these findings to the wind direction and speed in the five hours before the sampling, knowing that only data with wind speed values ≥2.0 m/s were taken into consideration. For M. aeruginosa, the detection of the highest cell abundances in the southernmost sampling points V2 and V3 were associated with winds blowing from the NW (Table 1), whereas those at the V1 and V4 sampling points were more surprisingly associated with winds from the SE. High cell abundances in the northern most sampling points V5 and V6 were equally associated with winds from NW and SE. For A. flos aquae, the results were more complicated, and no obvious link could be seen between the direction of the wind and the distribution of the cyanobacteria (Table 1). The same analyses were performed by taking into account the wind data one and two days before sampling (instead 5-10 hours before sampling), but no obvious relationship was detected (data not shown).

Influence of the number of sampling points on the estimated cyanobacterial cell abundances in the lake

The cyanobacterial cell abundances in the shallow lake were estimated by calculating the average value for the six sampling points (see Fig. 1). In order to determine the number of sampling points required to obtain a good estimation of cyanobacterial cell abundances in the lake, we compared the estimations of cell abundance based on using samples from just one, two, three, four or five sampling points with that based on all six. To do this, we calculated the correlation coefficients (Spearman) between the estimations based on the six sampling points and those based on one to five sampling points for each species (Fig. 7). We considered all possible combinations of points, and the results are classified in the figure on the basis of increasing order of r values within each combination of groups. For both species, we found that the estimations of cell abundances based on only one or two sampling points were generally rather badly correlated with those obtained using all six sampling points. On the other hand, it appeared that good correlations (around or > 0.9) were obtained when at least three sampling points were used, but also that the variations due to the choice of the sampling points was still considerable when only three sampling points were used.

In order to find out which combinations of sampling points provided the best results when only two or three sampling points were used, we classified all the possible combinations of points. To do this, we added the rank of each combination of sampling points obtained for the two species (M. aeruginosa and A. flos-aquae). From Figure 8, we can see that the best estimations obtained using only two or three sampling points were provided by combinations in which the sampling points used were on the shore opposite to the prevailing wind direction over the lake.

Diel variations in the subsurface cyanobacterial biomass in the lake

Finally, we carried out a 24-hour estimation of the variations in the total cyanobacterial biomass in the subsurface water (20 cm depth) of the lake, at five sampling points using the BBE torch (A-E, see Fig. 1). As shown in Fig. 9, there was a steady fall in the cyanobacterial biomass at all sampling points during the afternoon and evening, and conversely an increase late at night and in the morning. Moreover, the differences in biomass between the five sampling points were smaller during the night than during the day, as was the standard error (three measurements per sampling point). A multidimensional scaling analysis performed on the same values confirmed these observations, with all the night sampling times being grouped together, whereas the sampling times during the day were much more scattered (Data not shown).

Discussion

As far as we are aware, this is the first attempt to investigate the influence of sampling strategies on the evaluation of spatial and temporal variations in cyanobacterial abundances in shallow lakes, which constitute unstable and complex ecosystems. These lakes are used by humans for numerous activities, including recreational activities and the supply of drinking water, which makes the monitoring of cyanobacteria in such ecosystems of particular importance, especially as part of the evaluation of the health risks linked to cyanobacterial blooms and their toxins. Sampling strategy is also very important in the context of basic studies, because the quality of sampling has a major impact on the quality of the final results.

In this study, we found that the sampling frequency required to obtain a good estimation of the temporal evolution of the cyanobacterial abundance depends on the blooming species, M. aeruginosa or A. flos-aquae. Twice-monthly or monthly sampling provided good results for M. aeruginosa, whereas this was not often enough to monitor the chaotic population dynamics of A. flos-aquae. These findings are in contradiction with the recommendations of [START_REF] Codd | Toxic Cyanobacteria in water: A guide to their public health consequences, monitoring and management[END_REF], who proposed weekly or a twice-monthly sampling for species that do not form scum (A. flos-aquae for example), and more frequent sampling for scum-forming species (such as M. aeruginosa), because they can display more rapid changes in concentration. On the other hand, in agreement with these authors, our findings also demonstrate that a reactive approach to cyanobacterial sampling is called for, and that appropriate monitoring programs must be devised for each ecosystem based on what is known about how these systems function. It is clear that sampling only once or twice a month can lead to a very considerable under-estimation of cyanobacterial concentrations, and thus of the health risks associated with the bloom. As a result, a weekly sampling frequency seems to be required for cyanobacteria in small freshwater ecosystems.

Our data on the variability of the spatial distribution of cyanobacteria in the lake indicate that at least three sampling points were needed to obtain a good estimation of the abundance, based on a comparison with estimations based on six sampling points. It appeared also that if only three sampling points are used, the choice of the location of these sampling points is very important for the quality of the estimation. The most reliable results were obtained using sampling points located on the opposite side of the lake shore to the main axis of the wind direction, and that adding more sampling points reduces the impact of the choice of the location of the sampling points. Such horizontal variability in the distribution of cyanobacteria has been previously documented for many ecosystems, and also for many cyanobacterial species. For example, in a recent study, [START_REF] Briand | Spatiotemporal changes in the genetic diversity of a bloom-forming Microcystis aeruginosa (cyanobacteria) population[END_REF] showed that the spatial distribution of M. aeruginosa in a large freshwater reservoir on a given date could vary from 7.10 3 cells/mL to 2.10 8 cells/mL, depending on the location of the sampling points in the reservoir. Many factors and processes can influence the horizontal distribution of cyanobacteria in a freshwater ecosystem. Among them, wind and surface currents seem to have the greatest impact. For example, the distribution of Microcystis spp. in lake Taihu (see the review paper of [START_REF] Qin | A drinking water crisis in Lake Taihu, China: Linkage to climatic variability and lake management[END_REF] and in Lake Ontario [START_REF] Hotto | Molecular characterization of potential microcytsin-producing cyanobacteria in lake Ontario embayments and nearshore waters[END_REF] is clearly influenced by both winds and currents. Similarly, [START_REF] Moreno-Ostos | Spatial heterogeneity of Cyanobacteria and Diatoms in a thermally stratified canyon-shaped reservoir[END_REF] have shown that in a Spanish reservoir currents have a marked effect on the distribution of cyanobacteria, and more globally on the phytoplankton community. In this study, we found that the horizontal distribution of M. aeruginosa in the lake was influenced more by wind direction than that of A. flos-aquae. This could be explained by the fact that M. aeruginosa colonies are located at the surface of the lake at the end of the night, and thus are more subjected to the influence of the wind than A. flos-aquae filaments, which are distributed over the entire water column. We found also that two sampling points in the lake (V5 and V6) were less influenced by wind direction that the others. This could be explained by the fact that these two sampling points are protected from the influence of winds blowing from the NW by an embankment located in the North part of the lake. Finally, we also demonstrated that in such a small lake, the impact of wind occurred at the scale of a few hours, in contrast to the previous findings of [START_REF] Welker | Toxic Microcystis in shallow lakes Müggelsee (Germany) -dynamics, distribution, diversity[END_REF] showing that the distribution of cyanobacteria was influenced by winds that had been blowing one or two days earlier.

In addition to this variability in their horizontal distribution; the vertical distribution of cyanobacteria was also variable. Indeed, during the 24 h for which we used the BBE Torch to monitor the concentrations of cyanobacteria, we found that they were lower in the subsurface layer early at night than during the day. The greatest variations in biomass were recorded during the daytime, both at the scale of one sampling point when the three measurements were compared, and at the scale of the five sampling points monitored during this study.

These findings also suggest that several sampling points are necessary to obtain an accurate assessment of the cyanobacterial biomass and that integrated sampling of the first meter of the water column reduces the variability in the estimation of the biomass due to the position of cyanobacteria in the water column. This finding is consistent with data reported by [START_REF] Ahn | Comparison of sampling and analytical methods for monitoring of cyanobacteria-dominated surface waters[END_REF] showing that an integrated method was the most appropriate sampling method for Oscillatoria and Microcystis blooms. The causes of these variations in the position of cyanobacteria in the water column have been studied for different species. Several papers [START_REF] Porat | Diel buoyancy changes by the Cyanobacterium Aphanizomenon ovalisporum from a shallow reservoir[END_REF]Rabouille and Salençon, 2005;Rabouille et al., 2005;[START_REF] Visser | The ecophysiology of the harmful cyanobacterium Microcystis -Features explaining its success and measures for its control[END_REF][START_REF] Walsby | Gas vesicles[END_REF] have shown that migrations of cyanobacteria in the water column are probably due to the dynamics of the carbon-reserve metabolism, and are strongly influenced by light, temperature, and water mixing.

From all these findings, guidelines should be proposed for the monitoring of cyanobacteria in shallow lakes [START_REF] Codd | Toxic Cyanobacteria in water: A guide to their public health consequences, monitoring and management[END_REF] propose that the choice of sampling frequency and the choice of the number and location of the sampling sites should depend on the purpose of monitoring. For example, sampling near public bathing sites was recommended in freshwater ecosystems used for recreational activities. However, this strategy might generate data relevant only to the immediate vicinity of the bathing area, which do not reflect the global distribution of cyanobacteria in the lake. This is especially true when this distribution is very varied, and could make it difficult to prevent or manage blooms. On the basis of our findings, we proposed a different sampling strategy, which does not depend on the purpose of the monitoring. In order to minimize the cost of the cyanobacteria survey, twice-monthly sampling could be the norm for monitoring, but only if it is complemented by regular visual surveys. Changes in the appearance of the water (e.g. its color) between two successive dates would lead to an immediate increase in the sampling frequency. If it is not possible to carry out this visual survey, only a weekly sampling strategy can ensure that a sporadic cyanobacterial bloom is not missed. With regard to the number of sampling points, we found that at least three sampling points were necessary to obtain an accurate assessment of the cyanobacterial biomass (based on comparison with six sampling points). However, even when three sampling points were used, we found that the choice of the location of the sampling points was also very important (Fig. 8), even though the lake was fairly rectangular in shape and its perimeter small (around 1.3 Km). These findings suggest that for large lakes and also for lakes with a more complex shape, a large number of sampling points would be necessary to obtain a good estimation of the cyanobacterial abundance. Clearly such sampling is time consuming and expensive. One way to reduce these costs would be to collect a large number of samples and then pool equal volumes of these samples in the same flask, before carrying out a single analysis. In this study, as in most of the monitoring programs performed in small lakes, all samples were taken from the shoreline of the lake. This kind of sampling is suitable for small lakes, but it has been shown that for large lakes [START_REF] Rogalus | Evaluation of sampling and screening techniques for tiered monitoring of toxic cyanobacteria in lakes[END_REF] shoreline sampling may miss early warning signs of bloom development, and also lead to the overestimation of the concentration of microcystins, when compared to data obtained from offshore samples. For bigger lakes, therefore, the sampling strategy must include offshore samples.

Different programs worldwide are testing alternatives to water sampling for the monitoring of cyanobacteria in freshwater ecosystems. Two main approaches have been investigated. The first one is based on the use of remote sensing, which has long been in use in marine ecosystems (see for example [START_REF] Bracher | Quantitative observation of cyanobacteria and diatoms from space using PhytoDOAS on SCIAMACHY data[END_REF]. In freshwater ecosystems, the paper of [START_REF] Hunter | The spatial dynamics of vertical migration by Microcystis aeruginosa in a eutrophic shallow lake: A case study using high spatial resolution time-series airborne remote sensing[END_REF] has shown the potential of high resolution images for the assessment of the spatial distribution of M. aeruginosa in a shallow eutrophic lake. However, the cost of these images and the impact of meteorological conditions are limiting factors for envisaging the use of this tool in routine cyanobacteria monitoring programs. One alternative, lower-cost solution could be based, in the future, on the use of drones to take aerial photographs of freshwater ecosystems, but these tools are still in development. Moreover, they will be only useful for cyanobacterial species that live in the surface water of lakes.

The second way of monitoring of cyanobacteria without sampling the water being investigated is the use of buoys equipped with a variety of sensors, including, for example, a submersible spectrofluorometer to quantify the biomass of the cyanobacteria. This kind of tool permits the real-time monitoring of phytoplankton, including cyanobacteria, as shown for example in the paper of Le Vu et al. (in press). The two obstacles to their use in routine cyanobacteria monitoring programs are the high price of these systems, and the fact that they only provide estimations for one sampling point. Despite this, the possible use of such buoys, combined with the spatial monitoring of cyanobacteria by water sampling looks very promising for surveying cyanobacteria in freshwater ecosystems.

Conclusion

The sampling of cyanobacteria in freshwater ecosystems is a hot topic, in particular in the context of programs for surveying these toxic microorganisms in ecosystems used for the production of drinking water or for recreational activities. Paradoxically, fewer studies deal with the impact of sampling strategies on the estimation of cyanobacterial cell abundances in freshwater ecosystems. In this study, we demonstrate that the choice of sampling strategy can lead to very different estimations of the cell abundances of two blooming species in a shallow lake and also that, depending on the cyanobacterial species involved, different sampling strategies are required to obtain a good estimation of their population dynamics. All these findings suggested that monthly or twice-monthly sampling strategies at just one sampling point do not allow to provide an accurate estimation of cyanobacterial abundances, and thus of the health risks associated with the presence of toxic species in aquatic ecosystems.

Moreover, although promising new technologies are being developed for monitoring freshwater cyanobacteria, their cost and some other drawbacks mean that at present they cannot replace water sampling, which will remain the basis of most of these monitoring programs for the foreseeable future. (July,9,17 & 23;August,8)

Table 1 :

 1 Fig. 1: Geographical location of the study site in France (left), and of the sampling points in

Fig. 3 :

 3 Fig. 3: Simulation of the change over time of Microcystis aeruginosa cell concentrations

Fig. 4 :

 4 Fig. 4: Simulation of the change over time of the biomass of Aphanizomenon flos-aquae

Fig. 5 :

 5 Fig. 5: Spatial distribution of two cyanobacteria, Microcystis aeruginosa and Aphanizomenon

Acknowledgment

This work was funded by the Région Rhône-Alpes and the Conseil Général de la Loire.

Monika Ghosh is acknowledged for improving the English version of the manuscript. The comments and suggestions of the two anonymous reviewers were greatly appreciated.


P.-E Dumouchel

J.-J Marigo
email: marigo@lmm.jussieu.fr

M Charlotte

Dynamic fracture: an example of convergence towards a discontinuous quasi-static solution

Considering a one-dimensional problem of debonding of a thin film in the context of Griffith's theory, we show that the dynamical solution converges, when the speed of loading goes down to 0, to a quasi-static solution including an unstable phase of propagation. In particular, the jump of the debonding induced by this instability is governed by a principle of conservation of the total quasi-static energy, the kinetic energy being negligible.

Introduction

Griffith's theory of fracture mechanics is universally used to treat the propagation of cracks in brittle elastic bodies. When the load varies slowly with time, the quasi-static approach is preferred to the dynamic approach for its simplicity. But it turns out that, even when the crack path is known in advance, Griffith's law is unable to give the evolution of the crack if the total (quasi-static) energy is not a convex function of the crack length. Indeed, in this case, an unstable phase of propagation takes place and the quasi-static Griffith law, which does not allow discontinuous crack evolutions, cannot be applied any more. It is generally claimed that in a such situation the quasistatic framework must be abandoned because the crack propagation becomes dynamic and that the kinetic energy plays an important role. The goal of the paper is to show that, in a one-dimensional particular context, that claim is not true and that the unstable phase of propagation can also be treated by considering only quasi-static energies. But the price to pay is to replace the Griffith law of propagation based on the concept of critical energy release rate by that of conservation of the total quasi-static energy, at least during the phases of unstable propagation. This result is proved by solving the full dynamical problem and by passing to the limit when the loading speed goes down to zero. To our best knowlegde, a similar argument has been used only by Berry [START_REF] Berry | Some kinetic considerations of the Griffith criterion for fracture -I: Equations of motion at constant force[END_REF][START_REF] Berry | Some kinetic considerations of the Griffith criterion for fracture -II: Equations of motion at constant deformation[END_REF] in his approximate analysis of dynamic fracture (in 2D), which refered to a quasi-static G for evaluating the kinetic energy and the fracture energy at crack arrest in an homogeneous medium. Consequently, that gives a theoretical basis to the variational approach suggested by [START_REF] Mielke | Evolution of rate-independent systems[END_REF], [START_REF] Francfort | Existence results for a class of rate-independent material models with nonconvex elastic energies[END_REF] and [START_REF] Bourdin | The variational approach to fracture[END_REF] where the crack evolution is governed by two criteria: a stability criterion and an energy balance. Indeed, the variational formulation can be applied even if the total quasi-static energy is a non convex function of the crack length, the jumps of the crack are governed by that energy balance whereas the two criteria are equivalent to Griffith's law when the crack length evolves smoothly.

The paper is organized as follows. In Section 2, we first introduce the onedimensional problem which corresponds to a simplified version of the peeling test. Then, we set the dynamic problem and introduce the two cases of toughness repartition that are considered in the sequel, i.e. the so-called Case a and Case b. Follows a discussion where it is shown that we cannot solve the non convex Case b in a quasi-static context by using Griffith's law alone because the debonding evolution is necessarily discontinuous. As an alternative the energy conservation principle is proposed to determine the value of the jump. Section 3 is devoted to the dynamic analysis of each case. We first introduce the general structure of the dynamic response made of a sequence of shock waves which interact with the front of debonding. Then, we obtain the solutions in a closed form. Finally, we pass to the limit when the speed of loading goes down to 0 to prove the convergence to the quasi-static solutions in each case.

Setting of the problem

The dynamical problem

Notation and main assumptions

The problem deals with a simplified version of the peeling test, cf [START_REF] Obreimoff | The splitting strength of mica[END_REF], [START_REF] Lawn | Fracture of Brittle Solids -Second Edition[END_REF], [START_REF] Jaubert | Justification of Paris-type fatigue laws from cohesive forces model via a variational approach[END_REF] and Figure 1. We consider a semi-infinite perfectly flexible and inextensible thin film which is initially perfectly bonded to a rigid substrate with normal e 2 . The end x 1 = 0 of the film is submitted both to a constant tension -N e 1 , N > 0, and to an opening W linearly increasing with time t so that

W = ǫct, c = N ρ . (1)
In [START_REF] Berry | Some kinetic considerations of the Griffith criterion for fracture -I: Equations of motion at constant force[END_REF], ρ is the mass of a unit length of the film, c is the velocity of the transversal waves propagating in the film and ǫ represents the dimensionless speed of loading. In the sequel, ǫ will be considered as a (small) parameter. Therefore, the dependence of the response on ǫ is made explicit by denoting the solution with ǫ as a superscript. The main goal of the paper is to determine the asymptotic behavior of the dynamical response when ǫ goes down to 0. Moreover, in order to simplify the expressions we introduce dimensionless quantities. Specifically, if L denotes a characteristic length of the film (which will correspond to the point where there is a change of toughness in the model developed below), we set

x = x 1 L , T = ǫ ct L (2)
and consider that the displacement field u ǫ of the film is a function of (x, T ) defined on the quadrant Q := (0, +∞) 2 ,

u ǫ (x, T ) = u ǫ (x, T )Le 1 + w ǫ (x, T )Le 2 . (3)
In (2), T can be considered as a dimensionless rescaled time. Assuming that w ǫ is small and using a linearized inextensibility condition, u ǫ can be expressed in terms of w ǫ as follows:

u ǫ (x, T ) = 1 2 ∞ x ∂w ǫ ∂x 2 (s, T )ds. ( 4
)
The debonded part of the film at time t corresponds to the points x such that sup S≤T w ǫ (x, S) > 0. To simplify the presentation, we assume that the debonding grows from x = 0 at T = 0 in such a way that the debonded part of the film at time t corresponds to the interval (0, ℓ ǫ (T )). Thus ℓ ǫ (T )L represents the physical debonded length at time t. Therefore,

u ǫ = w ǫ = 0 in Q ǫ 0 := {(x, T ) ∈ Q : x ≥ ℓ ǫ (T )} (5)
with ℓ ǫ (0) = 0.

The equations of motion

The motion of the debonded part of the film is governed by the (classical) wave equation, cf [START_REF] Achenbach | Wave Propagation in Elastic Solids[END_REF]:

∂ 2 w ǫ ∂x 2 -ǫ 2 ∂ 2 w ǫ ∂T 2 = 0 in Q \ Q ǫ 0 . (6)
Let us note that the wave velocity in the (x, T )-plane is equal to 1/ǫ because of the time rescaling and the use of dimensionless coordinates. Moreover, since the gradient of the displacement is discontinuous across some curves of the (x, T )-plane, it is better to set the problem in terms of the gradient components. Thus, let us denote by ω ǫ and v ǫ the partial derivatives of w ǫ :

ω ǫ = ∂w ǫ ∂x , v ǫ = ∂w ǫ ∂T . (7)
The field ω ǫ represents the infinitesimal rotation of the film and the field v ǫ is the rescaled transverse velocity of the film (the real transverse velocity of the material points of the film is equal to ǫcv ǫ because of the time rescaling).

The rescaled transverse velocity must satisfy the boundary condition

v ǫ (0, T ) = 1, ∀T > 0. ( 8
)
Both fields have to satisfy the Hadamard compatibility condition and the equation of motion

0 = ∂ω ǫ ∂T - ∂v ǫ ∂x , 0 = ∂ω ǫ ∂x -ǫ 2 ∂v ǫ ∂T . (9)
Let S ǫ be the set of points (x, T ) ∈ Q \ Q ǫ 0 where ω ǫ and v ǫ are discontinuous; that corresponds to the shock waves. The equations [START_REF] Griffith | The phenomena of rupture and flow in solids[END_REF] are to be satisfied in Q \ (S ǫ ∪ Q ǫ 0 ). On S ǫ , they have to be replaced by the Hadamard compatibility condition and the Rankine-Hugoniot condition, see [START_REF] Truesdell | The Classical Field Theories[END_REF]:

0 = s[[ω ǫ ]] + [[v ǫ ]], 0 = [[ω ǫ ]] + sǫ 2 [[v ǫ ]]. (10)
In [START_REF] Achenbach | Wave Propagation in Elastic Solids[END_REF], s denotes the local velocity of a shock wave at a point of S ǫ , [[f ]] denotes the jump discontinuity of the function f across S ǫ :

[[f ]] = f + -f -on S ǫ , (11)
f + and f -denoting respectively the limit after and before the passage of the wave. We deduce from (10) that s 2 ǫ 2 = 1, i.e. , the shock waves propagate with a velocity equal to ±1/ǫ (in the (x, T )-plane). Accordingly, we can divide S ǫ into two subsets: the forward shock waves S ǫ + where s = +1/ǫ and the backward shock waves S ǫ -where s = -1/ǫ. Then, the only jump condition which remains to satisfy is the Hadamard compatibility condition which becomes:

0 = [[ω ǫ ]] + ǫ[[v ǫ ]] on S ǫ + , 0 = [[ω ǫ ]] -ǫ[[v ǫ ]] on S ǫ -. (12)

Griffith's law of debonding

Let us denote by Γ ǫ the front of debonding and lǫ its local speed in the (x, T )plane:

Γ ǫ := {(ℓ ǫ (T ), T ) : T > 0}, lǫ = dℓ ǫ dT . ( 13
)
The displacement w ǫ is continuous on Γ ǫ (and equal to 0), but ω ǫ and v ǫ are discontinuous and must satisfy the Hadamard compatibility condition

0 = [[v ǫ ]] + lǫ [[ω ǫ ]] on Γ ǫ . (14)
In ( 14) the double brackets denote, like in [START_REF] Truesdell | The Classical Field Theories[END_REF], the difference between the limit after and the limit before the passage of the front of debonding. The debonding of the film is governed by Griffith's law [START_REF] Griffith | The phenomena of rupture and flow in solids[END_REF] which is is formulated in terms of the dynamic energy release rate G ǫ . In a general two-dimensional context, G ǫ is defined as the limit of a path integral where the path tends to the tip of the crack, see [START_REF] Freund | Dynamic Fracture Mechanics[END_REF]. In our one-dimensional context, G ǫ is given by

G ǫ = N 2 [[ω ǫ2 ]] -ǫ 2 [[v ǫ2 ]] , (15)
where the jumps are defined on Γ ǫ , cf the appendix. Therefore, since the limits before the passage of the debonding front vanish, upon using [START_REF] Evans | Measure theory and fine properties of functions[END_REF], G ǫ reads as

G ǫ = N 2 1 -(ǫ lǫ ) 2 (ω ǫ + ) 2 . ( 16
)
Denoting by x → G c (x) the spatial repartition of the toughness, the Griffith law reads as, cf [START_REF] Griffith | The phenomena of rupture and flow in solids[END_REF][START_REF] Bui | Mécanique de la Rupture Fragile[END_REF]:

lǫ ≥ 0, G ǫ ≤ G c (ℓ ǫ ), G ǫ -G c (ℓ ǫ ) lǫ = 0. ( 17
)
In (17) the first inequality is the irreversibility condition, the second inequality is the Griffith criterion requiring that the dynamic energy release rate must be always less than the local toughness G c (ℓ), and the last equality stipulates that the debonding evolves only when the dynamic energy release rate is equal to the local toughness. According to ( 16) and (17), we obtain

0 < ǫ lǫ < 1, (18)
meaning that the speed of debonding is necessarily less than the wave speed.

Last definitions

The interpretation of the solution in energetic terms is fundamental in the sequel. Let us define the various energies involved. Since the film is perfectly flexible, the bending energy is neglected and the potential energy at time t is simply equal to the opposite of the work produced by the tension N in the displacement u ǫ (0, T ). Therefore, using (4), the potential energy of the film at time t can be written

P ǫ (T ) = N L 2 ∞ 0 ∂w ǫ ∂x 2 (x, T )dx. ( 19
)
The kinetic energy of the film at time t can be read as

K ǫ (T ) = ǫ 2 N L 2 ∞ 0 ∂w ǫ ∂T 2 (x, T )dx, ( 20
)
the factor ǫ 2 being due to the rescaling (2) of the time. Following Griffith's assumption the surface energy of the film at time t reads as

S ǫ (T ) = L ℓ ǫ (T ) 0 G c (x) dx. (21)
To summarize this setting of the dynamic debonding evolution, the fields ω ǫ , v ǫ and the debonding length ℓ ǫ must satisfy ( 8), ( 9), ( 12), ( 14), ( 16)-( 17). We will solve that problem in two particular cases of spatial toughness repartition:

Case a : G c (x) =    γ 1 N if 0 < x < 1 γ 2 N if x ≥ 1 , (22)
Case b : G c (x) =    γ 2 N if 0 < x ≤ 1 γ 1 N if x > 1 , (23)
where in both cases γ 2 > γ 1 > 0.

Remark 1 We construct in each case a (particular) solution. The issue of the uniqueness will not be discussed in this paper.

The quasi-static approach

In the quasi-static approach, the kinetic energy is neglected as well as any inertial effect and the film is assumed to be at equilibrium for every value of the given opening T . Therefore, let T and ℓ be the given opening and the debonding length. By virtue of the theorem of potential energy minimum, the transverse displacement field at equilibrium, w 0 , is that which minimizes the potential energy on all kinematically admissible vertical displacement fields w. Since w is admissible if and only if w ≥ 0, w = 0 in [ℓ, +∞) and w(0) = T , and, since the associated potential energy reads as ℓ 0 N Lw ′ (x) 2 /2 dx, then the potential energy of the film at equilibrium, P(T, ℓ), is given by

P(T, ℓ) = min {w≥0 : w(0)=T } ℓ 0 N L 2 w ′ (x) 2 dx. (24)
In ( 24) and above the prime denotes the derivative with respect to x. Then we easily deduce that the minimizer w 0 is given by

w 0 (x) = T 1 - x ℓ + , (25)
the superscript + denoting the positive part. Therefore the potential energy and the surface energy are given by

P(T, ℓ) = N LT 2 2ℓ , S(ℓ) = L ℓ 0 G c (x) dx (26)
and the quasi-static potential energy release rate reads as

G := - 1 L ∂P ∂ℓ (T, ℓ) = N T 2 2ℓ 2 . ( 27
)
Let us note that P is a strictly convex function of ℓ at given T > 0, while S is a convex function of ℓ if and only if x → G c (x) is monotone increasing. Therefore, the total energy E(T, ℓ) := P(T, ℓ) + S(ℓ) is a strictly convex function of ℓ at given T > 0 in Case a, but not in Case b.

If we still assume that the debonding evolution is governed by Griffith's law, then the function T → ℓ(T ) giving the evolution of the length debonding must be an absolutely continuous function (see [START_REF] Evans | Measure theory and fine properties of functions[END_REF] for a precise definition of absolutely continuous functions) which satisfies, for almost all T > 0:

l(T ) ≥ 0, N T 2 2ℓ(T ) 2 ≤ G c (ℓ(T )), N T 2 2ℓ(T ) 2 -G c (ℓ(T )) l(T ) = 0, (28)
with the initial condition ℓ(0) = 0. We have so the fundamental result Proposition 2.1 In Case a, Griffith's law (28) admits the solution given by

ℓ(T ) =        T / √ 2γ 1 , if 0 ≤ T ≤ √ 2γ 1 1, if √ 2γ 1 ≤ T ≤ √ 2γ 2 T / √ 2γ 2 , if T ≥ √ 2γ 2 . ( 29
)
In Case b there is no (absolutely continuous) solution.

Proof. In Case a let us verify that T → ℓ(T ) given by ( 29) is a solution of (28). Its continuity and its monotonicity are easily checked. We have

G(T ) = γ 1 N = G c (ℓ(T )) for T ∈ (0, √ 2γ 1 ) and G(T ) = γ 2 N = G c (ℓ(T )) for T ≥ √ 2γ 2 . For T ∈ ( √ 2γ 1 , √ 2γ 2 ), we have G(T ) < G c (ℓ(T )) = γ 2 N and l(T ) = 0.
Let us consider Case b and assume that there exists an absolutely continuous solution T → ℓ(T ). Then ℓ(T ) grows from 0 to ∞ when T grows from 0 to +∞. (Indeed, assume the contrary and set ℓ M = lim T →∞ ℓ(T ) < ∞. Then, according to (28) 2 , we should have

T 2 ≤ 2γ 2 ℓ M , ∀T > 0, what is impossible. Hence lim T →∞ ℓ(T ) = ∞.
) Therefore, ℓ(T ) = 1 for some T . Moreover, because of the monotonicity of T → ℓ(T ), there exist T 0 and T 1 with 0

< T 0 ≤ T 1 < +∞ such that ℓ(T ) < 1 for T < T 0 , ℓ(T ) = 1 for T ∈ [T 0 , T 1 ] and ℓ(T ) > 1 for T > T 1 . We should have G(T 0 ) = γ 2 N . (Indeed, if not, G(T 0 ) < γ 2 N
and by continuity G(T ) < γ 2 N in some interval (T 0 -h, T 0 ]. Hence, according to(28) 3 and by the absolute continuity assumption, ℓ(T ) = 1 in that interval, what is in contradiction with the definition of T 0 .) But, we should have also

G(T 1 ) ≤ γ 1 N . (Indeed, G(T ) ≤ G c (ℓ(T )) = γ 1 N for T > T 1
and hence by continuity the inequality holds at T 1 .) Finally, we should have 2γ 2 = T 2 0 ≤ T 2 1 ≤ 2γ 1 , which is in contradiction with the assumption on the toughness. Accordingly, there is no absolutely continuous solution.

✷

In Case a, a more careful analysis could show that (29) is the unique (absolutely continuous) solution of (28) and that ℓ(T ) is, at every T , the global minimizer of l → E(T, l) on {l ≥ 0}. The proof is not reproduced here, but the interested reader is referred to [START_REF] Francfort | Revisiting brittle fracture as an energy minimization problem[END_REF] for similar results relating global minimization to Griffith's law. The graph of T → ℓ(T ) is plotted in Figure 2. The evolution of the potential and surface energies are given by

P(T ) =        √ 2γ 1 T N L/2 if T ≤ √ 2γ 1 T 2 N L/2 if √ 2γ 1 ≤ T ≤ √ 2γ 2 √ 2γ 2 T N L/2 if T ≥ √ 2γ 2 (30)
and

S(T ) =        √ 2γ 1 T N L/2 if T ≤ √ 2γ 1 γ 1 N L if √ 2γ 1 ≤ T ≤ √ 2γ 2 γ 1 -γ 2 + √ 2γ 2 T /2 N L if T ≥ √ 2γ 2 (31)
In Case b, the debonding evolution is necessarily discontinuous. We can expect that the debonding evolves continuously and follows Griffith's law until the debonding reaches the point x = 1 where there is a loss of toughness.

In other words, we can expect that ℓ(T

) = T / √ 2γ 2 for T ∈ [0, √ 2γ
2 ). As it is shown in the proof of Proposition 2.1, at T = √ 2γ 2 , ℓ must jump from 1 to a certain ℓ c > 1. This phase of instability of the debonding can be interpreted from the evolution of the graph of the total energy as a function of the debonding length when the given opening grows, cf Figure 3 and see also [START_REF] Nguyen | Stability and Nonlinear Solid Mechanics[END_REF] for the use of the concepts of instability in Fracture Mechanics. For T < √ 2γ 2 ,

ℓ(T ) = T / √ 2γ 2 is a local minimum of l → E(T, l). But at T = √ 2γ 2 , this local
minimum disappears and the debonding length must jump from l = 1. The question is: "what is the equation giving the value of the jump?". Clearly, the quasi-static Griffith law is unable to give the answer because it is written only for smooth evolutions. We must introduce a new criterion and the main goal of the paper is to deduce the right equation from a dynamic analysis. Considering the graph of l → E( √ 2γ 2 , l) and by analogy with the motion of a ball in a basin, we could propose that the debonding length "falls" in the energy well, i.e. at the point ℓ m such that ∂E/∂l( √ 2γ 2 , ℓ m ) = 0 or in other words to the point such that G = γ 1 N . The dynamic analysis will show that this is not the good value, but that the right principle is the energy conservation: the jump of the debonding is such that the total quasi-static energy just after the jump remains equal to that just before the jump. Let us assume that this rule holds. Then, the length debonding after the jump, ℓ c , is such that ℓ c > 1 and

E( √ 2γ 2 , 1) = E( √ 2γ 2 , ℓ c ). A direct calculation gives ℓ c = γ 2 γ 1 > ℓ m = γ 2 γ 1 . ( 32
)
After this jump, the debonding front will propagate in a medium with a constant toughness and we can assume that the debonding evolution will follow the quasi-static Griffith law again. Since ℓ c > ℓ m , the energy release rate just after the jump is less than γ 1 N and the debonding length will remain at its value ℓ c as long as the opening T is not sufficiently large so that the energy release rate be equal to γ 1 N . That holds as long as T ≤ T c with T c given by

T c = 2γ 2 √ 2γ 1 . (33)
Then, when T > T c , the debonding restarts and the debonding length grows so that the energy release rate remains equal to γ 1 N . Therefore, ℓ(T

) = T / √ 2γ 1 for T > T c .
To summarize this analysis of Case b, if we assume that the debonding evolution follows always Griffith's law except during the unstable phase where the jump is governed by the energy conservation principle, the graph of T → ℓ(T ) should be that plotted in Figure 4,

ℓ(T ) =        T / √ 2γ 2 if 0 ≤ T < √ 2γ 2 γ 2 /γ 1 if √ 2γ 2 < T ≤ 2γ 2 / √ 2γ 1 T / √ 2γ 1 if T ≥ 2γ 2 / √ 2γ 1 . ( 34
)
The main goal of the paper is to justify this conjecture by solving the full dynamic problem and by passing to the limit when the speed of loading goes down to 0. We will also state a result of convergence for the energies. Their values corresponding to the conjectured quasi-static solution are given by

P(T ) =        √ 2γ 2 T N L/2 if T < √ 2γ 2 T 2 N L/2ℓ c if √ 2γ 2 < T ≤ T c √ 2γ 1 T N L/2 if T ≥ T c (35)
and Let us remark that their respective values just before and just after the jump are given by

S(T ) =        √ 2γ 2 T N L/2 if T < √ 2γ 2 (2γ 2 -γ 1 )N L if √ 2γ 2 < T ≤ T c γ 2 -γ 1 + √ 2γ 1 T /2 N L if T ≥ T c . (36)
P -2γ 2 = γ 2 N L, P + 2γ 2 = γ 1 N L, (37)
S -2γ 2 = γ 2 N L, S + 2γ 2 = (2γ 2 -γ 1 )N L. ( 38
)
Of course, since the total energy is conserved, all the released potential energy is so transformed into surface energy.

3 The dynamic analysis

The structure of the dynamic solution

The structure of the dynamic solution is the same for each case. The solution contains a first phase of debonding evolution where the debonding front propagates at a constant speed until it reaches the point x = 1 where the toughness changes. That change of the toughness generates a backward (traveling) shock wave as well as a change of the speed of debonding. The backward shock wave is reflected at x = 0 and is transformed into a forward (traveling) shock wave.

Since the wave speed is necessarily greater than the speed of debonding, the forward shock wave intersects the front of debonding. Then, the speed of debonding changes again and the forward shock wave is transformed into a backward shock wave. This second backward shock wave is reflected itself at x = 0 into a forward shock wave, which in its turn will intersect the front of debonding. Then, the speed of debonding changes, a backward shock wave is generated and so on.

We seek for a solution for which the speed of debonding as well as the rota- tion and the velocity fields are piecewise constant, i.e. constant in each sector delimited by the backward shock wave, the forward shock wave and the front of debonding, see Figure 5. Those sectors are denoted by {Q ǫ i } i∈N . The sector Q ǫ 0 corresponds to the sector ahead of the front of debonding, i.e. where ω ǫ = v ǫ = 0. The sector Q ǫ 1 is delimited by the front of debonding before x = 1, Γ ǫ 0 , the backward shock wave generated at x = 1, S ǫ 1 , and the boundary x = 0. Then the index i of the sector is incremented at each reflection of the shock wave. In the sector Q ǫ i , i ≥ 0, the (constant) values of (ω ǫ , v ǫ ) are denoted (ω ǫ i , v ǫ i ). Thus, ω ǫ 0 = v ǫ 0 = 0. The speed of debonding before x < 1, that is the slope of the line segment Γ ǫ 0 separating Q ǫ 0 and Q ǫ 1 , is denoted by lǫ 0 and the associated dynamic energy release rate by G ǫ 0 . The front of debonding reaches x = 1 at T = T ǫ 1 . Once the debonding propagates in the part x > 1, the speed of debonding and the energy release rate on the line segment Γ ǫ i separating Q ǫ 0 and Q ǫ 2i , i ≥ 1, are denoted by lǫ i and G ǫ i . The values of the opening when the backward shock wave is reflected at x = 0 and is transformed into a forward shock wave are denoted T ǫ 2i , i ≥ 1. The values of the opening when the forward shock wave intersects the front of debonding and is transformed into a backward shock wave are denoted T ǫ 2i-1 , i ≥ 2. The associated position of the debonding front is denoted ℓ ǫ i .

We put T ǫ 0 = 0 and ℓ ǫ 1 = 1. The equation of the first front of propagation reads as

Γ ǫ 0 : l = lǫ 0 T, 0 < T < T ǫ 1 . ( 39
)
The equations of backward shock wave S ǫ 2i-1 between Q ǫ 2i-1 and Q ǫ 2i , and of the forward shock wave S ǫ 2i between Q ǫ 2i and Q ǫ 2i+1 can be written as

S ǫ 2i-1 : T = T ǫ 2i -ǫl, 0 < l < ℓ ǫ i , S ǫ 2i : T = T ǫ 2i + ǫl, 0 < l < ℓ ǫ i+1 , ( 40
)
whereas the equation of Γ ǫ i , for i ≥ 1 reads as

Γ ǫ i : l = ℓ ǫ i + lǫ i (T -T ǫ 2i-1 ), T ǫ 2i-1 ≤ T < T ǫ 2i+1 . (41)
Therefore

S ǫ = i≥1 S ǫ i , Γ ǫ = i∈N Γ ǫ i (42)
and the T ǫ i 's are related to the ℓ ǫ i 's by

T ǫ 2i -T ǫ 2i-1 = ǫℓ ǫ i , T ǫ 2i+1 -T ǫ 2i = ǫℓ ǫ i+1 , ∀i ≥ 1. ( 43
)
Since the fields are piecewise constant, ( 9) is automatically satisfied in Q\ S ǫ ∪ Γ ǫ . It remains to find the sequences ω ǫ i , v ǫ i and lǫ i for i ∈ N. The procedure is the same in both cases and we use in this subsection the following notation for the toughness

G c (x) =    γ -N if x < 1 γ + N if x > 1 . ( 44
)
Let us begin by the first phase of debonding. The boundary condition (8) 1 gives v ǫ 1 = 1 and ( 14) leads to lǫ 0 ω ǫ 1 = -1. Hence lǫ 0 = 0 and, according to (17) 3 , G ǫ 0 = γ -N . Using ( 16), we get lǫ2 0 = (2γ -+ ǫ 2 ) -1 . Since lǫ 0 > 0 by virtue of (17) 1 , we finally have

lǫ 0 = 1 2γ -+ ǫ 2 , ω ǫ 1 = -2γ -+ ǫ 2 , v ǫ 1 = 1. ( 45
)
To find the debonding evolution in the part x > 1, we use the following statement:

Proposition 3.1 Let i ≥ 1. In Q ǫ 2i-1 , v ǫ 2i-1 = 1;
let us assume that ω ǫ 2i-1 and ℓ ǫ i are known. Then ω ǫ 2i , v ǫ 2i and lǫ i are given by

lǫ i = 1 ǫ (ω ǫ 2i-1 -ǫ) 2 -2γ + + (ω ǫ 2i-1 -ǫ) 2 + 2γ + , ω ǫ 2i = ω ǫ 2i-1 -ǫ 1 + ǫ lǫ i , v ǫ 2i = -lǫ i ω ǫ 2i , (46)
the superscript + in the first relation denoting the positive part. Furthermore, ω ǫ 2i+1 , v ǫ 2i+1 and ℓ ǫ i+1 are given by

ω ǫ 2i+1 = 2ω ǫ 2i -ω ǫ 2i-1 , v ǫ 2i+1 = 1, ℓ ǫ i+1 = 1 + ǫ lǫ i 1 -ǫ lǫ i ℓ ǫ i . (47)
Moreover the different energies are piecewise linear function of T . Their values at T ǫ i are given by

P ǫ (T ǫ 2i-1 ) = N L 2 ω ǫ 2i-1 2 ℓ ǫ i , (48)
K ǫ (T ǫ 2i-1 ) = ǫ 2 N L 2 ℓ ǫ i , (49)
S(T ǫ 2i-1 ) = γ -+ γ + (ℓ ǫ i -1) N L, (50)
P ǫ (T ǫ 2i ) = N L 2 ω ǫ 2i 2 (1 + ǫ lǫ i )ℓ ǫ i , (51)
K ǫ (T ǫ 2i ) = ǫ 2 N L 2 v ǫ 2i 2 (1 + ǫ lǫ i )ℓ ǫ i , (52)
S(T ǫ 2i ) = γ -+ γ + (1 + ǫ lǫ i )ℓ ǫ i -1 N L. (53)
Proof. Most of the above relations are direct consequences of the boundary condition (8) 1 , the Hadamard compatibility conditions ( 12) and ( 14) or the geometric properties (40), ( 41) and (43). In particular, (46) 2 is obtained by a simple combination of ( 12) and ( 14). The only part which needs a careful analysis is the formula (46) 1 for lǫ i . Let us first use ( 14) and ( 16) to obtain

v ǫ 2i = -lǫ i ω ǫ 2i , G ǫ i = N 2 1 -ǫ 2 lǫ2 i ω ǫ 2i 2 . ( 54
)
According to (17), if G ǫ i < γ + N then lǫ i = 0. In such a case, because of (54),

|ω ǫ 2i | < √ 2γ + . Otherwise, if G ǫ i = γ + N
, then, because of (54) again, we have

|ω ǫ 2i | ≥ 1 -ǫ 2 lǫ2 i |ω ǫ 2i | = 2γ + .
Therefore, it is sufficient to compare |ω ǫ 2i | to √ 2γ + to determine whether the debonding grows and find lǫ i . But, according to (46

) 2 , |ω ǫ 2i | ≥ √ 2γ + is equiva- lent to ω ǫ 2i-1 + ǫ ≥ √ 2γ + .
Inserting the relationship (46) 2 giving ω ǫ 2i in terms of lǫ i and ω ǫ 2i-1 into the equation G ǫ i = γ + N , we finally obtain lǫ i in terms of ω ǫ 2i-1 .

The energies are piecewise linear functions of T because ω ǫ , v ǫ and lǫ are piecewise constant. By using their definition (19), ( 20) and (21), and by remarking that

ℓ ǫ (T ǫ 2i ) = (1 + ǫ lǫ i )ℓ ǫ i leads to (48)-(53). ✷ 3.2 Case a

The dynamic solution

The speed of loading ǫ is assumed to be small enough to remove some particular cases. We can distinguish three phases in the debonding evolution, cf Figure 6.

(46), if

√ 2γ 1 + ǫ 2 + (2i -1)ǫ ≤ √ 2γ 2
, then lǫ i = 0 and all the other properties follow. Hence, (58) holds as long as i ≤ n ǫ with

n ǫ = 1 2ǫ 2γ 2 -2γ 1 + ǫ 2 + ǫ , (59)
• denoting the integer part. In terms of the opening, this phase of arrest corresponds to the interval (T ǫ 1 , T ǫ 2n ǫ +1 ) with

T ǫ 2n ǫ +1 = T ǫ 1 + 2ǫn ǫ . (60)
Let us note also that

T ǫ i = -ω ǫ i = 2γ 1 + ǫ 2 + (i -1)ǫ, 1 ≤ i ≤ 2n ǫ + 1. (61)
3.2.1.3 Debonding in the zone with the higher toughness. At T = T ǫ 2n ǫ +1 the debonding restarts and propagates inside the zone where the toughness is γ 2 . To find the solution, we use Proposition 3.1 and the following Lemma Lemma 3.2 Let η and X 0 be such that 0 < 1 -2η < X 0 < 1. Then the sequence {X i } i∈N defined by

1 X i+1 = X i + 2η (62)
remains always between 1 -2η and 1, its generic element X i is given by

1 X i + η + √ 1 + η 2 - 1 2 √ 1 + η 2 = (-1) i 1 + η 2 -η 2i 1 X 0 + η + √ 1 + η 2 - 1 2 √ 1 + η 2 (63)
and converges, with oscillations, to √ 1 + η 2 -η.

Proof. Let us verify that the sequence remains between 1 -2η and 1. By assumption, it is true for

X 0 . Let i ≥ 1. If X i ∈ (1 -2η, 1), then 1 > 1 X i + 2η = X i+1 > 1 1 + 2η > 1 -2η, thus X i+1 ∈ (1 -2η, 1
). Let us establish (63). We first make the following change of variable:

Y i = X i + √ 1 + η 2 + η. That allows to change the former sequence into 1 Y i+1 = 1 √ 1 + η 2 + η - √ 1 + η 2 -η 2 Y i ,
from which we easily deduce (63). Since √ 1 + η 2 -η < 1, the term on the right hand side of (63) tends to 0 when i grows to infinity and the result of convergence follows. The oscillations around the limit are due to the change of sign with i in that term. ✷ Let us set, for j ≥ 0,

X ǫ j = √ 2γ 2 ω ǫ 2n ǫ +1+2j + ǫ , η ǫ = ǫ √ 2γ 2 . ( 64
)
By definition of n ǫ , we have

√ 2γ 2 < ω ǫ 2n ǫ +1 + ǫ = √ 2γ 1 + ǫ 2 + (2n ǫ + 1)ǫ ≤ √ 2γ 2 + 2ǫ and hence, for ǫ small enough, 0 < 1 -2η ǫ < (1 + 2η ǫ ) -1 ≤ X ǫ 0 < 1.

Furthermore, since lǫ

n ǫ +1 > 0, we get from ( 46) and (47), 1/X ǫ 1 = X ǫ 0 + 2η ǫ . Since 0 < X ǫ 1 < 1, we have lǫ n ǫ +2 > 0 and we can iterate the procedure to obtain (62). Hence, we can use Lemma 3.2. By (46), ( 47) and (63), we obtain the debonding evolution for T ≥ T ǫ 2n ǫ +1 . In particular, lǫ i can be read as

ǫ lǫ i = 1 -(X ǫ i-n ǫ -1 ) 2 1 + (X ǫ i-n ǫ -1 ) 2 (65)
and we have the results of convergence lim

i→∞ lǫ i = lǫ ∞ := 1 √ 2γ 2 + ǫ 2 , lim i→∞ ω ǫ i = -2γ 2 + ǫ 2 ,
setting that the speed of debonding and the rotation of the film tend to those corresponding to the propagation in a homogeneous medium with toughness γ 2 . The speed of debonding lǫ i oscillates around that limit value and tends exponentially to it (with i). That phenomenon of oscillation can be interpreted as a boundary layer effect generated by the phase of propagation into the first zone and the phase of arrest at the interface which modify the initial conditions corresponding to the problem of propagation into a homogeneous medium.

Convergence to the quasi-static solution when ǫ → 0

When the speed of loading goes down to 0 we obtain the following expected convergence result Proposition 3.3 The dynamic solution T → ℓ ǫ (T ) converges to the quasistatic solution T → ℓ(T ) given by (29), uniformly on any compact, see Figure 7. Moreover, the potential and surface energies converge also uniformly on any compact to their quasi-static homologue given by (30)-(31), the kinetic energy converges to 0.

Proof

Step 1: Convergence of the first two phases. We deduce from (55) and (56) that lim ǫ→0 lǫ 0 = 1/ √ 2γ 1 and lim ǫ→0 T ǫ 1 = √ 2γ 1 .

Furthermore, since (59) gives lim ǫ→0 2ǫn ǫ = √ 2γ 2 -√ 2γ 1 , we also deduce from (60) that lim ǫ→0 T ǫ 2n ǫ +1 = √ 2γ 2 .

Step 2: Convergence of the third phase. It is the most difficult part of the proof, because of the oscillations. We know from Lemma 3.2 that lǫ i oscillates around the limit speed (2γ 2 + ǫ 2 ) -1/2 and converges to it when i → ∞. Let us define the mean value of the slope of the debonding front during one oscillation, i.e.

p ǫ i := ℓ ǫ i+2 -ℓ ǫ i T ǫ 2i+3 -T ǫ 2i-1 , i ≥ 1. (66)
Using ( 47), ( 62) and (65), we obtain after a tedious calculation

lǫ ∞ -p ǫ i = 2 lǫ ∞ X ǫ i-n ǫ -1 + η ǫ -1 + η 2 ǫ 2 (X ǫ i-n ǫ -1 + 2η ǫ ) 2 + (X ǫ i-n ǫ -1 ) 2 + 2 . Since X ǫ i + η ǫ -1 + η 2
ǫ decreases as i increases, at given ǫ, (that quantity tends to 0 when i → ∞,) it is less than X ǫ 0 + η ǫ -1 + η 2 ǫ which is of the order of ǫ. Therefore, denoting generically by C any positive constant independent of ǫ and i, we have Step 3: Convergence of the energies. Since ℓ ǫ (T ) → ℓ(T ), we also have S ǫ (T ) → S(T ), uniformly on any compact.

0 ≤ lǫ ∞ -p ǫ i ≤ Cǫ 2 . Furthermore, since 0 ≤ (2γ 2 ) -1/2 -lǫ ∞ ≤ Cǫ 2 and T ǫ 2n ǫ -1 - √ 2γ 2 ≤ Cǫ, we get |ℓ ǫ (T ) -ℓ(T )| ≤ CǫT for T ≥ √ 2γ 2 .

T QuasiStatic

For the potential and the kinetic energies, we have for 1 ≤ i ≤ 2n ǫ + 1:

P ǫ (T ǫ i ) = N L 2 (T ǫ i ) 2 , K ǫ (T ǫ i ) =    ǫ 2 N L/2, if i odd 0, if i even (67)
with T ǫ i given by (61). Since T ǫ 1 → √ 2γ 1 , since T ǫ 2n ǫ +1 → √ 2γ 2 , since the T ǫ i 's are equi-partitioned between T ǫ 1 and T ǫ 2n ǫ +1 and since the energies are piecewise linear, we obtain the uniform convergence of P ǫ to P given by ( 30)-(31) and of K ǫ to 0, in the interval [0, √ 2γ 2 ].

Since lim ǫ→0 X ǫ j = 1 for all j ≥ 0 (with a uniform convergence with respect to j), we obtain from (46), ( 64) and (65) that lim ǫ→0 ǫ lǫ i = 0, lim ǫ→0 ǫv ǫ 2i-1 = lim ǫ→0 ǫv ǫ 2i = 0 and lim ǫ→0 ω ǫ 2i-1 = lim ǫ→0 ω ǫ 2i = √ 2γ 2 , for all i ≥ n ǫ + 1 (with a uniform convergence with respect to i). Inserting that into (48) and (49), we obtain the uniform convergence of P ǫ to P given by ( 30)-(31) and of K ǫ to 0 for T ≥ √ 2γ 2 . ✷

Case b

The dynamic solution

The speed of loading ǫ is assumed to be small enough to remove some particular cases. We can divide the debonding evolution into four phases, see Figure 8. and of the final values are those of the quasi-static response, because the kinetic energy is then negligible. However, the kinetic energy is not always negligible during this phase because it takes a finite limit at T ǫ 2 . Specifically, let us examine the evolution of the kinetic energy. It is negligible during the first phase. Then, it grows rapidly during the propagation of the first backward shock wave to become maximal when that wave is reflected at x=0. But the kinetic energy decreases rapidly during the propagation of the first forward shock wave to become negligible again when that wave intersects the front of debonding, see Figure 10. During this come back of the wave, all the kinetic energy is transformed into surface energy. Finally, the jump of the debonding satisfies the principle of conservation of the quasi-static energy.

Hence, the phase arrest converges to that of the quasi-static solution. Since v ǫ i oscillates between 0 and 1 during that phase, the kinetic energy converges to 0 as ǫ 2 . We easily deduce that the other energies converge to their quasi-static homologue.

3.3.2.4 Convergence of the last phase. The proof that the dynamic last phase (where the debonding restarts and propagates with a oscillating speed) converges to the quasi-static last phase is quite similar to that given in Case a, see Step 2 of Proposition 3.3. During this phase the kinetic energy remains of the order of ǫ 2 and all the energies converge to their quasi-static homologue. We can summarize the convergence result obtained in Case b by Proposition 3.4 When the speed of loading ǫ goes down to 0, the dynamic solution of Case b converges to the extended quasi-static solution (34). Thus the evolution is given by the quasi-static Griffith's law during the first, third and fourth phases where the debonding evolution is continuous, but the second phase leads to a jump of the debonding which satisfies the principle of conservation of the total quasi-static energy. The role of the kinetic energy is transitory. It takes finite values only during the first back and forth of the shock waves, being negligible before and after the jump.

Conclusion

We can learn two major lessons from this model problem of dynamic fracture:

(1) Contrary to what is generally claimed, we can treat the phase of rapid propagation by considering only quasi-static quantities, because the kinetic energy plays only a transitory role. (2) The right criterion giving the value of the crack jump is not Griffith's criterion formulated in terms of the energy release rate but the conservation of the total energy.

Of course, it remains to generalize those results to two-or three-dimensional heterogeneous bodies. In particular, by considering layered composite materials we could obtain their effective toughness.

Fig. 1 .

 1 Fig. 1. Debonding of a thin inextensible film.

1 2 Γ 1 2 Γ 2 TFig. 2 .

 2122 Fig.2. The quasi-static evolution of the debonding in Case a by using Griffith's law.

Fig. 3 .

 3 Fig. 3. The graph of the total energy as a function of the debonding length ℓ for a given opening T in Case b. left: when T < √ 2γ 2 ; right: when T = √ 2γ 2 and the value ℓ c of the debonding length after the jump given by the energy conservation principle.

Fig. 4 .

 4 Fig.[START_REF] Francfort | Existence results for a class of rate-independent material models with nonconvex elastic energies[END_REF]. The quasi-static evolution of the debonding in Case b by using both Griffith's law for the continuous phases and the principle of energy conservation for the jump.

Fig. 5 .

 5 Fig. 5. Structure of the dynamic solution and numbering of the sectors. Thick line: debonding front; thin lines: backward and forward shock waves.

Fig. 7 .

 7 Fig. 7. Case a : Comparison of the dynamic solution when ǫ = 0.075 with the quasi-static solution for γ 1 = 0.5, γ 2 = 2.

Fig. 8 .Fig. 10 .

 810 Fig. 8. Case b : Dynamic solution with the four phases of the debonding propagation for γ 1 = 0.5, γ 2 = 2. and ǫ = 0.05

3. 3 . 2 . 3

 323 Convergence of the phase of arrest. Passing to the limit in (76) and (77) we get lim ǫ→0 T ǫ 2n ǫ +1 = 2γ 2 √ 2γ 1 = T c .

Acknowledgements

We thank the referee who has carefully read the paper and who has verified each formula and each proof. His remarks have allowed us to improve the

The debonding front reaches x = 1 at T ǫ 1 given by

Then the debonding evolution is governed by Proposition 3.1 with γ + = γ 2 .

3.2.1.2 Arrest at the change of toughness. For ǫ small enough, since |ω ǫ 1 | + ǫ = √ 2γ 1 + ǫ 2 + ǫ < √ 2γ 2 , we deduce from (46) that lǫ 1 = 0. The debonding stops at x = 1 in the interval (T ǫ 1 , T ǫ 1 + 2ǫ) corresponding to the first back and forth of the shock wave. Upon inserting into (46) and (47), we deduce the values of the constants in

Let us show that the debonding does not evolve during a certain number n ǫ of back and forth of the shock wave. Specifically, let us show by induction that, if 1 ≤ i ≤ n ǫ , with n ǫ to be determined, then lǫ i = 0 and

(58) That holds for i = 1. Let us assume that it is true until i -1 ≥ 1 and let us find under which condition that remains true for i. Using the induction assumption, we have

For T > T ǫ 1 , the debonding evolution is governed by Proposition 3.1.

3.3.1.2

Rapid debonding in the zone with the higher toughness. There is major difference between the two cases when the front of debonding reaches the point x = 1 where the toughness changes. In Case a since the toughness increases, the debonding stops. On the contrary, in Case b since the toughness decreases, the debonding accelerates. Specifically, we deduce from (46) that

Hence, ǫ lǫ 1 is less than 1 but of the order of 1, i.e. the debonding propagates at a speed of the same order that the wave speed which remains an unreachable limit speed. The first backward wave reaches x = 0 at T ǫ 2 = T ǫ 1 + ǫ. Then, it is reflected and the generated forward shock wave intersects the front of debonding at (ℓ ǫ 2 , T ǫ 3 ) given by

At T ǫ 1 and T ǫ 3 , the kinetic energy is of the order of ǫ 2 while it is of the order of 1 at T ǫ 2 . Specifically, we have

At T = T ǫ 3 , the debonding stops. Indeed, we have

Therefore, for ǫ small enough, since |ω ǫ 3 | + ǫ < √ 2γ 1 , we deduce from (46) that lǫ 2 = 0. Furthermore, by virtue of ( 46) and (47) we get

With the same procedure as in 3.2.1.2, using (46), we can show that

and

with

In terms of T , the debonding does not evolve in the interval (T ǫ 3 , T ǫ 2n ǫ +1 ) with

(77)

Slow debonding in the zone with the higher toughness. At T = T ǫ 2n ǫ +1 the debonding restarts. Indeed, for j ≥ 0 let us put

By definition of n ǫ , we have

Furthermore, since lǫ n ǫ +1 > 0, we get from ( 46) and (47), 1/X ǫ 1 = X ǫ 0 + 2η ǫ . Since 0 < X ǫ 1 < 1, we have lǫ n ǫ +2 > 0 and by induction we obtain 1/X ǫ j+1 = X ǫ j + 2η ǫ . We can use Lemma 3.2. Owing to (46), (47) and (63), we obtain the debonding evolution for T ≥ T ǫ 2n ǫ +1 . In particular, for i ≥ n ǫ + 1, lǫ i can be read as

and we get the convergence result

The speed of debonding lǫ i oscillates around this limit value and tends exponentially to it (with respect to i).

Convergence to the quasi-static solution when ǫ → 0

We pass to the limit when ǫ goes down to 0 and compare the limits of the dynamic solution to that of the conjectured quasi-static response given by (34).

3.3.2.1

Convergence of the first phase. We deduce from (68) that lim ǫ→0 lǫ 0 = 1/ √ 2γ 2 and lim ǫ→0 T ǫ 1 = √ 2γ 2 . Hence, the convergence to the corresponding quasi-static first phase. Let us note that the kinetic energy converges to 0 like ǫ 2 .

Convergence of the second phase to the predicted jump.

Let us remark first that lim

then that lim

Therefore, the second phase tends to the jump of the debonding length from ℓ = 1 to ℓ = ℓ c at T = √ 2γ 2 as in (34). This convergence result can be interpreted in terms of the energies. Thus, at T ǫ 1 , T ǫ 2 and T ǫ 3 , the kinetic energy takes the values

Therefore, at the limit

For the surface energy, we have

and at the limit

The potential energy takes the following values

and at the limit

While comparing to (37) and (38), we see that the limits of the initial values

quality of the paper.

A Definition of the dynamic potential energy release rate

Let us consider a two-dimensional homogeneous elastic medium in which a crack propagates in the direction 1. In dynamics, the potential energy release rate G can be defined by the limit of an integral over a path which tends to the tip of the crack. Let Γ r be a such path -for instance a circle of radius r and of center the tip of the crack -G can be read as (cf [START_REF] Freund | Dynamic Fracture Mechanics[END_REF])

where ρ is the mass density of the material, u the velocity vector field, u the displacement vector field, ε(u) the strain tensor field, W the elastic potential, σ the stress tensor field and n the output unit normal to the path. The dynamic energy release rate differs from its static homologue by the first term corresponding to the kinetic energy.

This relationship of G can be used for the peeling test as follows. Let ℓ be the position of the debonding tip. The path Γ r corresponds to the points {ℓ-r, ℓ+r}, the normal to the path at ℓ±r is equal to ±1, all the indices take the value 1, the displacement u corresponds to the opening w and the strain reads as ε = w ′ . Assuming that the elastic potential is quadratic, W (w ′ ) = N w ′2 /2, the stress becomes σ = N w ′ and the dynamic energy release rate reads as

where the dot denotes the time derivative and the prime the spatial derivative. Furthermore, if we assume that the debonding grows in the direction 1 and if the front of debonding is oriented in the space-time so that the + side corresponds to the points after the passage of the tip debonding, we get

remaining that the double brackets denote the jump.