Benoit Lebreton
email: benoit.lebreton@univ-lr.fr

Jennifer A Pollack
email: jennifer.pollack@tamucc.edu

Brittany A Blomberg
email: brittany.blomberg@tamucc.edu

Terence A Palmer
email: terry.palmer@tamucc.edu

Leslie A Adams
email: leslie.adams@tamucc.edu

Gaël Guillou
email: gguill01@univ-lr.fr

Paul A Montagna
email: paul.montagna@tamucc.edu

Jennifer Beseres Pollack

Origin, composition and quality of suspended particulate organic matter in relation to freshwater inflow in a South Texas estuary

Keywords: suspended particulate organic matter, Chl a: chlorophyll a, GoM: Gulf of Mexico, POC: Particulate organic matter Suspended particulate organic matter, Freshwater inflow, Stable isotope ratio, Chlorophyll a, Gulf of Mexico, Texas estuary

des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Introduction

Estuaries are among the most productive ecosystems on Earth, with production equal to that of tropical rain forests [START_REF] Mclusky | The Estuarine Ecosystem[END_REF]. The quantity and timing of freshwater delivery to the mixing zone is essential to their functioning [START_REF] Montagna | Freshwater inflow: science, policy, management[END_REF]. Freshwater inflows play key roles in carrying continental organic matter from the watershed to the estuary [START_REF] Hedges | What happens to terrestrial organic matter in the coean?[END_REF], and in balancing effects of tidal inputs of saltwater and of evaporation, consequently affecting the structure of habitats. As a result, freshwater inflows affect inputs of continental and marine organic matter and autochthonous pelagic and benthic primary production [START_REF] Longley | Estuaries[END_REF]. Increasing human populations will result in an increasing demand for freshwater, which will affect freshwater inflow into estuaries and have consequences on the functioning of estuarine ecosystems [START_REF] Montagna | Freshwater inflow: science, policy, management[END_REF].

One consequence is the potential change of composition of suspended particulate organic matter (SPOM), an integral component of estuarine food webs [START_REF] Longley | Freshwater inflows to Texas bays and estuaries: ecological relationships and methods for determination of needs[END_REF]. Organic matter locally produced in estuaries or provided by rivers or the sea varies in terms of quantity and quality depending on its origin and its state of decay [START_REF] Dunton | Response of estuarine marsh vegetation to interannual variations in precipitation[END_REF][START_REF] Raymond | Riverine export of aged terrestrial organic matter to the North Atlantic Ocean[END_REF]. This diversity of primary productions can be used by many different consumers, making highly complex estuarine food webs. Consequently, estuaries provide many ecological functions, including breeding grounds, nurseries and feeding sites for a large diversity of pelagic and benthic species [START_REF] Beck | The identification, conservation, and management of estuarine and marine nurseries for fish and invertebrates[END_REF]. For example, SPOM is generally the main food source of suspension feeders like oysters. As a result, a change of abundance or composition of SPOM may severely affect the oyster diet [START_REF] Riera | Temporal variation of δ 13 C in particulate organic matter and oyster Crassostrea gigas in Marennes-Oléron bay (France): effect of freshwater inflow[END_REF], and therefore the general functioning of oyster reefs, which provide many ecological functions [START_REF] Beck | Oyster reefs at risk and recommendations for conservation, restoration, and management[END_REF][START_REF] Grabowski | Economic valuation of ecosystem services provided by oyster reefs[END_REF].

The origins, quantity, and quality of SPOM can be determined by measuring several variables, including stable isotope ratios, chlorophyll a (chl a) content and C/chl a ratios. Stable isotope ratios of carbon and of sulfur provide information about the origin of organic matter. Stable isotopes of carbon allow for discrimination between the different estuarine primary producers [START_REF] Fry | δ 13 C measurements as indicators of carbon flow in marine and freshwater ecosystems[END_REF][START_REF] Peterson | Stable isotopes in ecosystem studies[END_REF] and are often used to determine composition of SPOM in estuaries [START_REF] Cifuentes | Stable carbon and nitrogen isotope biogeochemistry in the Delaware estuary[END_REF][START_REF] Riera | Temporal variation of δ 13 C in particulate organic matter and oyster Crassostrea gigas in Marennes-Oléron bay (France): effect of freshwater inflow[END_REF][START_REF] Savoye | Dynamics of particulate organic matter δ 15 N and δ 13 C during spring phytoplankton blooms in a macrotidal ecosystem (Bay of Seine, France)[END_REF][START_REF] Survey | National Water Information System data available on the World Wide Web (USGS Water Data for the Nation)[END_REF][START_REF] Modéran | Organic matter exploitation in a highly turbid environment: Planktonic food web in the Charente estuary, France[END_REF].

Due to potential overlaps of carbon stable isotope composition between different primary producers and the effects of decay processes on isotopic composition [START_REF] Benner | Depletion of 13 C in lignin and its implication for stable carbon isotope studies[END_REF], stable isotopes of sulfur can be used to discriminate between the different food sources. Primary producers from continental areas are generally much more 34 S-depleted than primary producers from oceanic environments [START_REF] Peterson | Multiple stable isotopes used to trace the flow of organic matter in estuarine food webs[END_REF][START_REF] Peterson | Sulfur, carbon, and nitrogen isotopes used to trace organic matter flow in the salt-Marsh estuaries of Sapelo Island, Georgia[END_REF][START_REF] Thode | Sulphur isotopes in nature and the environment: an overview[END_REF]. Stable isotope composition of nitrogen is generally used to determine the trophic levels of consumers [START_REF] Minagawa | Stepwise enrichment of 15 N along food chains: further evidence and the relation between δ 15 N and animal age[END_REF][START_REF] Vander Zanden | Primary consumer 13 C and N and the trophic position of aquatic consumers[END_REF], but can also be an indicator of anthropogenic influence (i.e., release of wastewater from treatment plants) due to increase of δ 15 N values with bacterial activity [START_REF] Mcclelland | Linking nitrogen in estuarine producers to land-derived sources[END_REF][START_REF] Mooney | Watershed export events and ecosystem responses in the Mission-Aransas National Estuarine Research Reserve, South Texas[END_REF]. The measurement of these three stable isotopes enables discrimination between the different pools of organic matter in order to determine its origin. The measurement of total organic carbon, chl a and C/chl a ratios along with these stable isotope ratios help determine the quantity and the quality of SPOM along a salinity gradient. The aim of this study was to determine how the different primary producers and inputs of allocthonous organic matter (i.e., from continental and marine origins) affect the composition, quantity and quality of SPOM, and to determine the evolution of these parameters in relation to changes in freshwater inflows. This study was conducted in the Mission-Aransas estuary, Texas, a semi-arid estuary that is vulnerable to changes in natural and anthropogenic changes in freshwater inflow.

Material and methods

Study area

The Mission-Aransas estuary is one of seven major estuarine systems located along the Texas coast (Fig. 1). It is a shallow, bar-built estuary composed of Aransas bay, the primary bay located closest to the Gulf of Mexico inlet, and Copano Bay, the secondary bay located closest to the two main rivers: the Mission and Aransas rivers. The Mission-Aransas estuary is 2 m deep on average, its surface is 463 km 2 and its volume is 0.93 km 3 [START_REF] Armstrong | The Ecology of Open-bay Bottoms of Texas: a Community Profile[END_REF]. The Copano bay watershed is 4851 km 2 [START_REF] Borel | Estimating E. coli and Enterococcus loads in a coastal Texas watershed[END_REF]; 49% of the freshwater comes from the Mission river and 15% from the Aransas river [START_REF] Orlando | Salinity Characteristics of Gulf of Mexico Estuaries[END_REF].

The Mission-Aransas estuary is located in a sub-tropical humid climate. Rainfalls generally occur as brief and intense rain showers and as a result, freshwater inflows occur as episodic pulses.

Rainfall, and consequently freshwater inflows, varies from year to year with alternating dry and wet periods [START_REF] Evans | Hydrography and oceanography[END_REF]. From 2001-2010, freshwater inflows from both rivers ranged from 17 to 430 million m 3 y -1 (U.S. Geological [START_REF] Survey | National Water Information System data available on the World Wide Web (USGS Water Data for the Nation)[END_REF]. On an annual basis, Mission-Aransas estuary is losing water due to higher evaporation (1513 mm y -1 ) than precipitation (886 mm y -1 ) [START_REF] Armstrong | The Ecology of Open-bay Bottoms of Texas: a Community Profile[END_REF]. The estuary has low mixing efficiency (< 0.05) and long residence times (~360 days, [START_REF] Solis | Hydrography, mixing characteristics, and residence times of Gulf of Mexico estuaries[END_REF]. The salinity structure is driven by episodic freshwater pulses that depress salinities and then maintain low salinities for a prolonged period [START_REF] Orlando | Salinity Characteristics of Gulf of Mexico Estuaries[END_REF].

The Mission-Aransas estuary hosts diverse primary producers and is fueled by different allocthonous sources i.e., continental organic matter from Mission and Aransas rivers in Copano bay, pelagic organic matter from the Gulf of Mexico in Aransas bay, organic matter from San Antonio bay transported southward across the Mesquite bay via the Intracoastal Waterway [START_REF] Solis | Hydrography, mixing characteristics, and residence times of Gulf of Mexico estuaries[END_REF].

Sampling stations

Five stations were sampled along a salinity gradient in the Mission-Aransas estuary, with station 1 (St. 1) the closest to the Aransas river mouth and station 5 (St. 5) the closest to the Gulf of Mexico (Fig. 1). Three other stations were sampled to collect freshwater from the two main rivers (Mission river (MIS) and Aransas river (ARA), and sea water from the Gulf of Mexico (GoM). Freshwater inputs into Copano bay were determined by summing flows of Mission and Aransas rivers using data from USGS gauge stations (Fig. 1) located at Refugio (Mission river, Site ID: 08189500) and Skidmore (Aransas river, Site ID: 08189700) (U.S. Geological Survey, 2012).

Collection of water samples and sample processing

At each station, three samples of surface water were collected on a monthly basis, from November 2010 (relatively high freshwater inflow) to August 2011 (relatively low freshwater inflow).

All water samples were collected far enough from the shore to be representative of the river or bay water mass. Due to restricted access, it was impossible to sample water at the mouth of the Mission river.

Salinity was measured with an YSI 6920 multiprobe sonde, using practical salinity scale. Water samples were pre-filtered on a 250-µm sieve to eliminate large zooplankton and detritus. For each water sample, SPOM was collected on four different filters (precombusted Whatman GF/F glass fiber filters, 0.7 µm porosity) in order to measure stable isotope composition (δ 13 C, δ 15 N, δ 34 S) and chl a concentration. A volume ranging from 5 to 325 ml was filtered under moderate vacuum; filters were then frozen at -20°C and freeze-dried. Carbonates were removed from filters for δ 13 C and %C analyses by contact with HCl fumes in a vacuum-enclosed system.

Collection of primary producers for stable isotope analyses

To determine the composition of SPOM, the main primary producers (C 3 and C 4 salt marsh plants, seagrass, benthic microalgae) were sampled at different locations along the salinity gradient and during different seasons to examine possible spatio-temporal changes in isotopic composition.

Vascular plants (leaves) and seagrass (leaves and roots) were collected, placed in plastic bags and stored in ice until return at the lab. They were then rinsed with tap water to remove detritus and sediment, frozen at -20°C, freeze-dried and ground to a fine and homogeneous powder using a ball mill.

Benthic microalgae were collected by scraping surface sediments in the field and then by extracting the microalgae in the laboratory following the method of [START_REF] Riera | Isotopic determination of food sources of Crassostrea gigas along a trophic gradient in the estuarine bay of Marennes-Oléron[END_REF], slightly modified by [START_REF] Herlory | Methodology of light response curves: application of chlorophyll fluorescence to microphytobenthic biofilms[END_REF]. Extracted samples were checked under a microscope for purity, then benthic microalgae were collected on three different filters (precombusted Whatman GF/F glass fiber filters) in order to measure stable isotope composition (δ 13 C, δ 15 N, δ 34 S). Samples were frozen at -20°C then freeze-dried.

Chl a and stable isotope ratio measurements

Chl a was extracted overnight from filters and read on a Turner Trilogy fluorometer (Turner Designs, Sunnyvale, USA) using a non-acidification technique [START_REF] Welschmeyer | Fluorometric analysis of chlorophyll a in the presence of chlorophyll b and pheopigments[END_REF]. Detection limit of the fluorometer was 0.025 µg.L -1 with a linear range of 0-300 µg.L -1 . Calibration was performed once a year using commercially available fluorometric chlorophyll calibration standards for chlorophyll a (Turner Designs).

Samples for stable isotope analyses and particulate organic carbon percentages were analyzed using an elemental analyzer (Flash EA 1112, Thermo Scientific, Milan, Italy) coupled with an isotope ratio mass spectrometer (Delta V Advantage with a Conflo IV interface, Thermo Scientific, Bremen, Germany) at the LIENSs stable isotope facility of the University of La Rochelle, France. Results are expressed in the δ notation as deviations from standards (Vienna Pee Dee Belemnite for δ 13 C, N 2 in air for δ 15 N and Canyon Diablo Troilite for δ 34 S) following the formula: δ 13 C, δ 15 N or δ 34 S = [(R sample /R standard ) -1] x 10 3 , where R is 13 C/ 12 C, 15 N/ 14 N or 34 S/ 32 S, respectively. Calibration was done using reference materials (USGS-24, IAEA-CH6, -600 for carbon; IAEA-N2, -NO-3, -600 for nitrogen; IAEA-S-1, -S-2 for sulfur). Analytical precision based on the analyses of acetanilide (Thermo Scientific) used as laboratory internal standard was <0.06 and <0.1‰ for carbon and nitrogen, respectively. Analytical precision based on the analyses of sulfanilamide (Thermo Scientific) used as laboratory internal standard was <0.15‰. Vanadium pentoxide was added in each SPOM sample as a catalyzer during sample combustion for sulfur elemental analysis.

Data analyses and statistics

Mean daily freshwater discharges were determined from September 1 st , 2010 to September, 30 th 2011, using data from USGS gauge stations. Loads of POC into the Copano bay from the Aransas and the Mission rivers were computed by multiplying the mean daily freshwater discharge with the concentrations of POC (µg.ml -1 ). C/chl a ratios were computed using concentrations of particulate organic carbon (µg l -1 ) from elemental analyses and of chlorophyll a from fluorimetry measurements (µg l -1 ). Comparisons between stable isotope values of primary producers and SPOM samples were conducted using nonparametric procedures (i.e., Kruskal-Wallis rank sum test and Wilcoxon signed rank test), to determine any differences between locations, sampling dates and tissues (roots and leaves). Kruskal-Wallis tests were followed by multiple comparisons of means using the pgirmess package [START_REF] Giraudoux | pgirmess: data analysis in ecology[END_REF] in R statistical software (R Development Core Team, 2010).

Results and discussion

River discharges and hydrological features of the Mission-Aransas estuary in 2010-2011

The year 2011 was a very dry year-the second driest (after 2008) of the past 10 years [START_REF] Mooney | Watershed export events and ecosystem responses in the Mission-Aransas National Estuarine Research Reserve, South Texas[END_REF]. Inputs of freshwater from the Aransas and Mission rivers to Copano bay were very low, with 21 million m 3 , compared to an average of 170 million m 3 y -1 over the past 10 years (U.S. Geological [START_REF] Survey | National Water Information System data available on the World Wide Web (USGS Water Data for the Nation)[END_REF]. Freshwater inflows in the Mission-Aransas estuary decreased by more than 2 orders of magnitude from September 2010 (2.30 10 6 m 3 d -1 ) to September 2011 (8.05 10 3 m 3 d -1 ). A large pulse of freshwater occurred in September 2010-before the start of the experiment-with the highest input of 1.92 10 7 m 3 d -1 on September 20 th , leading to a strong salinity gradient in the estuary (Fig. 2 ; [START_REF] Mooney | Watershed export events and ecosystem responses in the Mission-Aransas National Estuarine Research Reserve, South Texas[END_REF][START_REF] Bruesewitz | Estuarine ecosystem function response to flood and drought in a shallow, semiarid estuary: Nitrogen cycling and ecosystem metabolism[END_REF]. A minor pulse occurred in January 2011, with the highest input of 5.15 10 6 m 3 d -1 on January 17 th , 2011.

During the January pulse (January 17 th and 18 th , 2011), 96 and 4 % of the freshwater inflow were issued from Mission and Aransas rivers, respectively. Salinity increased along a gradient from Aransas river mouth (Station 1) to the Gulf of Mexico (GoM) (Fig. 2). This gradient was more evenly distributed along the length of the sampling transect in November 2010 (salinity from 5.2 to 24.3) than in August 2011 (salinity from 30.5 to 37.3) due to the overall increase of salinities in Copano and Aransas bays. The salinities of the river station ARA were close to zero throughout the study period although the salinities at the MIS river station rose to 5 by August 2011. Salinities measured at Station 5 and GoM were similar throughout the sampling period.

On average, evaporation largely exceeds precipitation in the Mission-Aransas estuary and south east winds generally predominates [START_REF] Armstrong | The Ecology of Open-bay Bottoms of Texas: a Community Profile[END_REF], which leads to a flow of water from Aransas bay into Copano bay, this flow being twice more than the volume of the freshwater issued from rivers [START_REF] Johnson | A general method for modeling coastal water pollutant loadings[END_REF]. Minor rain events in the Mission-Aransas watershed, like the one observed in January 2011, do not have any effect on the general hydrology of the bay, as evidenced by the general increase of salinity in the bay throughout the study period [START_REF] Bruesewitz | Estuarine ecosystem function response to flood and drought in a shallow, semiarid estuary: Nitrogen cycling and ecosystem metabolism[END_REF].

Nevertheless, a slight decrease in salinity was observed at Station 4 in January 2011, indicating a potential flow of freshwater originating from San Antonio bay, which is connected to Aransas bay via the Intracoastal Waterway in Mesquite bay [START_REF] Solis | Hydrography, mixing characteristics, and residence times of Gulf of Mexico estuaries[END_REF]. As a result, from November 2010 to August 2011, a combination of low freshwater inflow, high evaporation and flow of saline water from Aransas bay into Copano bay led to a general increase in salinity throughout the estuary and compression of the salinity gradient, which was observed in November 2010. Thus, from November 2010 to August 2011, the Mission-Aransas estuary went from a system dominated by precipitation (i.e., positive estuary) to a system dominated by evaporation (i.e., negative estuary) [START_REF] Pritchard | Estuarine hydrography[END_REF][START_REF] Montagna | Conceptual model of estuary ecosystems[END_REF].

During the study period, particulate organic carbon (POC) concentrations ranged from 0.48 (December 10) to 2.8 µg ml -1 (August 11) in Aransas river and from 0.53 (March 11) to 2.12 µg ml -1 (August 11) in Mission river. POC concentration data in Mission river in November and December are missing due to technical issues. Except in January, computed load of POC into the Mission-Aransas estuary ranged from 334 to 939 kg month -1 and from 15 to 1241 kg month -1 for Aransas and Mission rivers, respectively. In January, POC load into the Mission-Aransas estuary was much higher, particularly for Mission river (13097 kg month -1 ) while Aransas river POC load was 3292 kg month -1 (Fig. 2).

Discrimination of the different benthic primary producers using stable isotopes

Salt marsh plants had typical isotopic compositions for carbon as C 3 and C 4 plants were relatively 13 C-depleted (-29.2 to -26.8‰) and 13 C-enriched (-15.5 to -12.5 ‰), respectively (Table 1, Fig. 3). This disparity is mostly related to different enzymatic processes occurring during photosynthesis and differential utilization of HCO 3 and CO 2 [START_REF] Smith | Two categories of 13 C/ 12 C ratios for higher plants[END_REF][START_REF] Fry | δ 13 C measurements as indicators of carbon flow in marine and freshwater ecosystems[END_REF][START_REF] O'leary | Physical and chemical basis of carbon isotope fractionation in plants[END_REF]. Seagrass material was the most 13 C-enriched, with δ 13 C values ranging from -14.7 to -10.7‰, typical for seagrasses [START_REF] Hemminga | Stable carbon isotopes in seagrasses: variability in ratios and use in ecological studies[END_REF] and slightly higher than C 4 vascular plants (Table 1, Fig. 3). Detrital matter from seagrass had a very low δ 13 C value at station 3 (-18.1‰) whereas δ 13 C value at station 4 (-11.9‰) was in the range of fresh material. This large difference is probably related with the degradation process which leads to higher levels of lignin in detrital matter, which is 13 C-depleted [START_REF] Benner | Depletion of 13 C in lignin and its implication for stable carbon isotope studies[END_REF]. As a result, degradation of seagrass may lead to a larger overlap between δ 13 C values of seagrass detritus and those of C 4 salt marsh plants, either as fresh or as detritus. Benthic microalgae had a wide range of δ 13 C values, ranging from -21.7 to -12.0‰, largely overlapping the ranges of values observed for seagrass and C 4 vascular plants.

No pattern was observed for δ 13 C values of C 3 and C 4 plants between species, sampling dates and locations in the bay. Leaves and roots of seagrass were more 13 C-enriched in downstream sections than in upstream sections of the Mission-Aransas estuary (Table 1, Kruskal-Wallis rank sum test, pvalues: Ruppia maritima: leaves: 0.013, roots: 0.016, Halodule wrigthii: leaves: 0.001, roots: 0.004).

A similar pattern was observed for microphytobenthos (Kruskal-Wallis rank sum test, p-value: 0.035) and was clearly related to salinity: microphytobenthos had lower δ 13 C values (from -21.7 to -16.4‰) at low salinities (Stations 1 and 3 in November, station 1 in March) than at high salinities (Fig. 4). This is probably due to the higher influence of dissolved inorganic carbonates from freshwater, which are 13 C-depleted [START_REF] Fry | δ 13 C measurements as indicators of carbon flow in marine and freshwater ecosystems[END_REF]. It is thus possible to discriminate these benthic algae from the other primary producers. The δ 13 C values of benthic algae collected in marine influenced environments, ranging from -14.2 to -12.0‰, were in the range of those from the literature [START_REF] Currin | The role of standing dead Spartina alterniflora and benthic microalgae in salt marsh food webs: considerations based on multiple stable isotope analysis[END_REF][START_REF] Riera | Isotopic determination of food sources of Crassostrea gigas along a trophic gradient in the estuarine bay of Marennes-Oléron[END_REF][START_REF] Riera | Determination of food sources for benthic invertebrates in a salt marsh (Aiguillon bay, France) by carbon and nitrogen stable isotopes: importance of locally produced sources[END_REF][START_REF] Kharlamenko | Determination of food sources of marine invertebrates from a subtidal sand community using analyses of fatty acids and stable isotopes[END_REF][START_REF] Lebreton | Trophic importance of diatoms in an intertidal Zostera noltii seagrass bed: Evidence from stable isotope and fatty acid analyses[END_REF].

Seagrass, C 3 and C 4 vascular plants had a large range of δ 15 N values, ranging from 2.0 to 9.2‰, from 4.5 to 10.0‰ and from -1.1 to 7.7‰, respectively. Very different δ 15 N values were observed between the two samples of seagrass detrital matter, with values of 7.5‰ at station 3 and 1.8‰ at station 4 in November 2010. Halodule wrightii and Spartina alterniflora clearly had higher δ 15 N values in upstream parts of the estuary than in downstream parts (Table 1, Wilcoxon signed rank test, p-values: Halodule wrigthii: leaves: 0.003, roots: 0.004, Spartina alterniflora: <0.001). δ 15 N values of microphytobenthos ranged from -0.6 to 5.3‰, with slightly higher δ 15 N values at lower salinities.

Sulfur stable isotope composition of primary producers was highly variable, even at the species level (e.g., Spartina alterniflora, Borrichia frutescens, Ruppia maritima) (Table 1, Fig. 3). Among primary producers, C 3 vascular plants were the most enriched in 34 S, with δ 34 S values ranging from 7.1 to 18.4‰. C 4 vascular plants were less enriched in 34 S, with values ranging from -5.8 to 7.1‰.

However, these values were higher than those previously observed during other studies [START_REF] Peterson | Sulfur and carbon isotopes as tracers of salt-marsh organic matter flow[END_REF][START_REF] Peterson | Sulfur, carbon, and nitrogen isotopes used to trace organic matter flow in the salt-Marsh estuaries of Sapelo Island, Georgia[END_REF][START_REF] Sullivan | Edaphic algae are an important component of salt marsh foodwebs: evidence from multiple stable isotope analyses[END_REF][START_REF] Currin | The role of standing dead Spartina alterniflora and benthic microalgae in salt marsh food webs: considerations based on multiple stable isotope analysis[END_REF]. This large range of δ 34 S values is linked with differential assimilation of sediment sulfides, which are 34 Sdepleted, and of sulfates, which are 34 S-enriched, depending on seasonal and spatial conditions [START_REF] Carlson | Uptake of dissolved sulphide by Spartina alterniflora. Evidence from natural sulphur isotope abundance ratios[END_REF][START_REF] Fry | Sulphur uptake by salt grasses, mangroves, and seagrasses in anaerobic sediments[END_REF][START_REF] Currin | The role of standing dead Spartina alterniflora and benthic microalgae in salt marsh food webs: considerations based on multiple stable isotope analysis[END_REF][START_REF] Stribling | Variability of stable sulfur isotopic ratios in Spartina alterniflora[END_REF]. The much higher δ 34 S values measured on C 3 plants in comparison with those of C 4 plants indicate that C 3 plants assimilate sulfates originating from precipitation (0 to 15‰) and seawater (21‰) [START_REF] Peterson | Sulfur and carbon isotopes as tracers of salt-marsh organic matter flow[END_REF][START_REF] Peterson | Sulfur, carbon, and nitrogen isotopes used to trace organic matter flow in the salt-Marsh estuaries of Sapelo Island, Georgia[END_REF][START_REF] Krouse | Pedosphere and biosphere[END_REF][START_REF] Connolly | Sulfur stable isotopes separate producers in marine food-web analysis[END_REF]. No clear pattern was observed between δ 34 S values and salinity, location or sampling time for C 3 and C 4 vascular plants.

Microphytobenthos had intermediate δ 34 S values, going from 4.1 to 11.8‰, in the range of those reviewed by [START_REF] Currin | The role of standing dead Spartina alterniflora and benthic microalgae in salt marsh food webs: considerations based on multiple stable isotope analysis[END_REF]. No relation was observed between δ 34 S values of microphytobenthos and salinities, in contrast to the pattern observed in carbon isotopic composition, as these samples were collected at salinities (>5.2) for which isotopic composition of sulfates is already relatively high. In the Mississippi river estuary, [START_REF] Fry | Sulfur stable isotope indicators of residency in estuarine fish[END_REF] highlighted that the stable isotope composition of sulfates mostly varies at very low salinities (<0.6) and then remains relatively stable (21‰), close to the isotopic composition of sulfates in seawater [START_REF] Nriagu | Hydrosphere[END_REF].

Seagrass also had a wide range of δ 34 S values, with roots generally much more 34 S-depleted than leaves: -11.1 to 7.8 ‰ vs. -1.5 to 6.1‰, respectively. The differences observed between leaves and roots are probably related with a differential uptake of sulfur by these two tissues. Roots mostly use sediment sulfur which contains high quantities of sulfides resulting from anaerobic bacterial activity, and are generally very 34 S-depleted (-10 to -35 ‰), whereas leaves use dissolved sulfates from the water column (21‰) [START_REF] Fry | Sulphur uptake by salt grasses, mangroves, and seagrasses in anaerobic sediments[END_REF][START_REF] Nriagu | Hydrosphere[END_REF][START_REF] Moncreiff | Trophic importance of epiphytic algae in subtropical seagrass beds: evidence from multiple stable isotope analyses[END_REF]. As a result, there is good discrimination between seagrass roots and benthic algae or C 4 plants thanks to δ 34 S values. However, the large range of δ 34 S values observed on seagrass leaves and C 4 plants do not allow for discrimination between these two primary producers. The large variations described for δ 34 S values, depending on species, tissues, locations and seasons [START_REF] Fry | Sulphur uptake by salt grasses, mangroves, and seagrasses in anaerobic sediments[END_REF][START_REF] Moncreiff | Trophic importance of epiphytic algae in subtropical seagrass beds: evidence from multiple stable isotope analyses[END_REF][START_REF] Oakes | Causes of sulfur isotope variability in the seagrass, Zostera capricorni[END_REF][START_REF] Belicka | Stable isotope and fatty acid biomarkers of seagrass, epiphytic, and algal organic matter to consumers in a pristine seagrass ecosystem[END_REF], signify that information issued from sulfur stable isotope composition should be used with caution in large scale ecosystem studies [START_REF] Connolly | Sulfur stable isotopes separate producers in marine food-web analysis[END_REF].

As a result, C 3 vascular plants are very well differentiated from all other primary producers and some microphytobenthos samples are clearly separated from seagrass and C 4 plants. δ 34 S measurements allow for good discrimination between microphytobenthos and seagrass detritus (Table 1, Fig. 3). The combination of three different stable isotope markers (δ 13 C, δ 15 N and δ 34 S) thus allowed the discrimination of three different groups of benthic primary producers from the Mission-Aransas estuary: 1) C 3 salt marsh plants, 2) benthic microalgae, and 3) C 4 salt marsh plants and seagrass [START_REF] Peterson | Sulfur, carbon, and nitrogen isotopes used to trace organic matter flow in the salt-Marsh estuaries of Sapelo Island, Georgia[END_REF][START_REF] Sullivan | Edaphic algae are an important component of salt marsh foodwebs: evidence from multiple stable isotope analyses[END_REF][START_REF] Currin | The role of standing dead Spartina alterniflora and benthic microalgae in salt marsh food webs: considerations based on multiple stable isotope analysis[END_REF]Connolly et al., 2004).

Origin and quality of allocthonous inputs from the Aransas and Mission rivers and from the

Gulf of Mexico

The Mission-Aransas estuary receives inputs of organic matter from its different watersheds, mostly the Aransas and Mission rivers, which drain the highest quantity of freshwater into this estuary. Some allocthonous inputs may also originate from the Gulf of Mexico, via the pass at Port Aransas, due to movements of water masses due to tides, winds and evaporation of water into Aransas bay [START_REF] Solis | Hydrography, mixing characteristics, and residence times of Gulf of Mexico estuaries[END_REF]. During this survey, the origin and quality of allocthonous inputs of organic matter may have been influenced by two patterns: 1) From November 2010 to August 2011, a general decrease of freshwater inflow, and 2) A minor rain event in January 2011. Chl a concentrations were relatively low all year long, with values ranging from 2.1±0.2 (May 2011) to 6.7±0.3 µg.l -1 (February 2011) (Fig. 6). Combined to the low C/chl a values (ranging from 154 to 484, Fig. 6), indicating that this organic matter was fresh [START_REF] Cifuentes | Stable carbon and nitrogen isotope biogeochemistry in the Delaware estuary[END_REF][START_REF] Savoye | Origin and composition of particulate organic matter in a macrotidal turbid estuary: The Gironde Estuary, France[END_REF], stable isotope values indicated that the composition of organic matter from the Gulf of Mexico varied little all year long, being mostly composed of fresh marine phytoplankton. Seawater SPOM δ 13 C values range indeed typically from -18 to -24‰ (Fry andSherr, 1984, Riera et al., 2000) and δ 34 S values are close to 21‰ [START_REF] Fry | Sulphur uptake by salt grasses, mangroves, and seagrasses in anaerobic sediments[END_REF][START_REF] Nriagu | Hydrosphere[END_REF].

Composition of SPOM from Aransas and Mission rivers was generally stable and much more 13 Cdepleted (i.e., ranging from -35.5 to -27.7‰) than SPOM at GoM station, except in January 2011, when δ 13 C values of SPOM from Aransas (-24.0‰) and Mission rivers (-25.1‰) were similar to that from GoM (Fig. 5, Tables 2 and3, Appendix A). δ 15 N values of SPOM from rivers were much more variable (Table 2, Fig. 5). SPOM from Mission river typically had δ 15 N values similar to those observed in the Mission-Aransas estuary, but low δ 15 N values in winter (3.6‰ in January 2011) (Fig. 5, Table 3). δ 15 N values of SPOM from Aransas river were the most variable, with a large enrichment in 15 N between January and February 2011 (17.5‰) and then a gradual decrease in 15 N, until reaching δ 15 N values observed in the Mission-Aransas estuary in July 2011 (4.8‰) (Fig. 5). SPOM sulfur stable isotope composition was generally relatively depleted in the Aransas and Mission rivers (from 3.3 to 9.6‰), except in November and December 2010 in the Mission river, when δ 34 S values of SPOM were slightly higher (12.1 and 10.7‰, respectively) (Fig. 5).

Due to sampling complications, chl a concentrations and C/chl a ratios are not available at station MIS in November and in December 2010. From November 2010 to August 2011, chl a concentrations were highly variable in Aransas and Mission rivers, ranging from 0.7 (MIS, January) to 43.2 µg.l -1 (ARA, February) (Fig. 6, Table 2). Chl a concentrations were particularly high in Aransas river in February and March. From May to August, chl a concentrations were frequently higher in Aransas and Mission rivers than at station 5 or GoM (Table 3). C/chl a ratios were high in January in Mission (1638) and Aransas (1385) rivers (Fig. 6).

Except in January and February 2011, the composition and quality of organic matter issued from the Aransas and Mission rivers was pretty constant. The combination of C/chl a ratios and of δ 13 C and δ 34 S values demonstrates indeed that the organic matter from the Aransas and Mission rivers was mostly composed of freshwater phytoplankton. The δ 13 C values of SPOM, lower than -27 ‰, are typical of freshwater phytoplankton [START_REF] Fry | δ 13 C measurements as indicators of carbon flow in marine and freshwater ecosystems[END_REF][START_REF] Cifuentes | Stable carbon and nitrogen isotope biogeochemistry in the Delaware estuary[END_REF][START_REF] Boutton | Stable carbon isotope ratios of natural materials: II. Atmospheric, terrestrial, marine, and freshwater environments[END_REF] and δ 34 S values (from 3.3 to 12.2‰) indicated that these primary producers were relying on sulfates from riverine waters [START_REF] Nriagu | Isotopic composition of sulphur in atmospheric precipitation around Sudbury, Ontario[END_REF][START_REF] Chukhrov | The isotopic composition of plant sulfur[END_REF][START_REF] Nriagu | Hydrosphere[END_REF]. Carbon to chl a ratios, with most of the values lower than 200, confirmed that SPOM was mostly composed of fresh material [START_REF] Cifuentes | Stable carbon and nitrogen isotope biogeochemistry in the Delaware estuary[END_REF][START_REF] Savoye | Origin and composition of particulate organic matter in a macrotidal turbid estuary: The Gironde Estuary, France[END_REF]. Blooms of freshwater phytoplankton were observed in June and July in the Mission river and in July and August in the Aransas river, when Chl a concentrations were high (>15 µg l -1 ) and C/chl a ratios were lower than 100, demonstrating that SPOM was only composed of fresh phytoplankton.

In January 2011, a minor rain event strongly affected the composition of the riverine SPOM.

Much higher C/chl a values (> 1300) in January demonstrated that SPOM was mostly composed by refractory material (i.e., detrital matter) [START_REF] Cifuentes | Stable carbon and nitrogen isotope biogeochemistry in the Delaware estuary[END_REF][START_REF] Savoye | Origin and composition of particulate organic matter in a macrotidal turbid estuary: The Gironde Estuary, France[END_REF]. Higher δ 13 C values (-24.0 and -25.1‰ in Aransas and Mission rivers, respectively) indicated a higher influence of C 4 plant material in the riverine SPOM composition than during previous months. The Mission-Aransas estuary watershed has a high number of C 4 primary producers, since it is composed of 30% of grassland [START_REF] Mooney | Watershed export events and ecosystem responses in the Mission-Aransas National Estuarine Research Reserve, South Texas[END_REF], with 68% of the grasses being C 4 plants [START_REF] Teeri | Climatic patterns and the distribution of grasses in North America[END_REF]. Moreover, main crops in the Mission-Aransas estuary, particularly in the Aransas watershed are C 4 plants (i.e., sorghum and corn) [START_REF] Mooney | Watershed export events and ecosystem responses in the Mission-Aransas National Estuarine Research Reserve, South Texas[END_REF]. The 13 C enrichment indicates rainfall drained high quantities of detrital matter from continental origin within the watershed into the rivers. In February 2011, the peak of δ 15 N values in the Aransas river is very likely related with the release of highly 15 N-enriched nitrogen originating from wastewater treatment plants [START_REF] Mcclelland | Linking nitrogen in estuarine producers to land-derived sources[END_REF], probably due to the overwhelming of their treatment capacities during the January rain event [START_REF] Mooney | Watershed export events and ecosystem responses in the Mission-Aransas National Estuarine Research Reserve, South Texas[END_REF]. This higher release of enriched nitrogen only occurred in the Aransas river, which is subject to a much higher discharge of wastewater than the Mission river [START_REF] Mooney | Watershed export events and ecosystem responses in the Mission-Aransas National Estuarine Research Reserve, South Texas[END_REF]. The release of enriched nitrogen, in higher quantities than during base flow [START_REF] Bruesewitz | Estuarine ecosystem function response to flood and drought in a shallow, semiarid estuary: Nitrogen cycling and ecosystem metabolism[END_REF], stimulated the production of freshwater phytoplankton in February and March as illustrated by the higher chl a concentrations. The gradual decrease of δ 15 N values after this event indicates that there was a slow dilution of this enriched nitrogen during the following months, due to the very low flows of fresh water. This highlights that the different management and land use land cover of these two riversone with numerous discharge points of wastewater and dominated by a watershed dominated by agricultural fields, and the other one with low wastewater release and a watershed dominated by shrub-have a different effect on the riverine SPOM composition. Rain events also have a clear effect on the composition of riverine organic matter: In periods of normal flow, riverine waters contain low concentrations of high quality SPOM whereas they contain high concentrations of low quality SPOM during rain events.

As evidenced by the decrease in salinity observed at St. 4 in January, there are also exchanges of organic matter between Copano and Aransas bays, San Antonio bay and Corpus Christi bay through the Intracoastal Waterway. Unfortunately, our sampling design did not take into account these potential flows of organic matter, and it is therefore not possible to determine the composition and the quantity of organic matter carried by the Intracoastal Waterway.

Autochthonous phytoplankton dominates SPOM composition during droughts

δ 13 C values of SPOM in the Mission Aransas estuary (i.e., stations 1 to 5) ranged from -26.9 (January 2011, Station 2) to -20.8‰ (April 2011, station 4) and were relatively stable throughout the year. No clear temporal pattern was observed even though some δ 13 C values of SPOM were occasionally significantly lower during some fall and winter months than during some spring and summer months (Fig. 5, Table 2). No clear spatial difference was observed between δ 13 C values in the estuary, except in January 2011, when SPOM at station 2 was significantly more depleted than at station 5 (Fig. 5, Table 3). δ 15 N values of SPOM at stations 1 to 5 were also relatively stable and similar between stations, ranging from 4.7 (August 2011, station 1) to 8.8‰ (November 2010, station 2) (Fig. 5, Tables 2 and3). δ 34 S values of SPOM were rather relatively stable and similar throughout the year from stations 1 to 5, ranging from 11.9 (January 2011, station 1) to 20.0‰ (April 2011, station 4) (Fig. 5, Tables 2 and3). Chl a concentrations were less variable in the Mission-Aransas estuary (St. 1 to St. 5) than in the rivers, ranging from 0.7 (July 2011, station 5) to 20.1 µg.l -1 (April 2011, station 1). Two peaks of chl a concentrations in November and April were observed at station 1 (Fig. 6, Table 2). C/chl a ratios observed in the Mission-Aransas estuary ranged from 114 (January 2011, station 1) to 803 (February 2011, station 3) (Fig. 6).

As a result, the composition of SPOM was relatively constant in the Mission-Aransas estuary during the survey, even if some slight spatio-temporal differences were observed. On a yearly basis, δ 13 C and δ 34 S values of estuarine SPOM were very different from the carbon and sulfur isotopic compositions of benthic primary producers (Fig. 3), demonstrating that these primary producers had a negligible influence on the composition of estuarine SPOM, either as fresh (i.e., resuspended benthic microalgae) or as detrital matter (i.e., seagrass, C 3 or C 4 plant) [START_REF] Qian | Suspended particulate organic matter (SPOM) in Gulf of Mexico estuaries: compound-specific isotope analysis and plant pigment compositions[END_REF]. The carbon and sulfur isotopic composition of estuarine SPOM demonstrated that this organic matter was, all year long, mostly composed by locally produced phytoplankton. This indicates minimal benthicpelagic coupling in the Mission-Aransas estuary, which is relatively surprising because windy conditions can increase turbidity [START_REF] Beseres Pollack | A Restoration suitability index model for the eastern oyster (Crassostrea virginica) in the Mission-Aransas estuary, TX, USA[END_REF], allowing for resuspension of sediment organic matter. Organic matter of benthic origin is thus either exported as large particles (e.g., seagrass leaves) or buried in the sediment.

The carbon and sulfur isotopic compositions of estuarine SPOM were very different from those of riverine SPOM (Fig. 3). From March to July, there was also no influence of the highly 15 N-enriched SPOM discharge from the Aransas river in Copano bay (Fig. 5). This indicates that organic matter from rivers had a very low influence on the general functioning of the Mission-Aransas estuary during droughts, probably because inputs of riverine SPOM were in very low quantity, due to the very low freshwater inflows, whereas autochthonous phytoplanktonic production is high in the Mission-Aransas estuary, like in other Texas estuaries [START_REF] Pennock | Nutrient behavior and phytoplankton production in Gulf of Mexico estuaries[END_REF].

Effect of freshwater inflow into composition of SPOM

A slight influence of riverine SPOM was occasionally and very locally observed during the minor rain event -which occurred in January -at stations close to the Aransas and Mission river mouths and at station 4, probably due to an influence of riverine organic matter from San Antonio bay drained through the Intracoastal Waterway. We used data from our survey (i.e. during drought) and from the survey carried out by [START_REF] Mooney | Watershed export events and ecosystem responses in the Mission-Aransas National Estuarine Research Reserve, South Texas[END_REF] (i.e. during drought and flood) for Copano bay, the secondary bay of the estuary, to determine how freshwater inflows affect the composition and the quality of the SPOM. We determined that an exponential relationship exists between SPOM δ 13 C values and monthly freshwater inflows (moving average of the present and the two previous monthly freshwater inflows) from the Aransas and Mission rivers (Fig. 7).

We then considered two scenarios: droughts and rain events. In our study, we observed that the riverine SPOM is dominated by freshwater phytoplankton (δ 13 C = -31.4‰) during droughts and by detrital matter of C 3 -plants (δ 13 C =-24.7‰) during rain events. We considered that the SPOM δ 13 C value in Copano bay without any freshwater inflow (i.e., autochthonous production only) would be equal to the δ 13 C value of SPOM from the Gulf of Mexico (-22.1 ‰). For both scenarios, we computed the average δ 13 C value which would result from an equal mixture of SPOM from rivers and SPOM from autochthonous production. We then determined the threshold of monthly freshwater inflow between a system more influenced by riverine inputs vs. a system more influenced by autochthonous production. In other words, a monthly freshwater inflow higher than the threshold indicates that SPOM in the secondary bay is dominated by riverine inputs whereas a monthly freshwater inflow lower than the threshold indicates that SPOM in the secondary bay is dominated by autochthonous production.

The scenario based on drought conditions gave an inflow threshold of 6.8 million m 3 month -1 . This inflow of freshwater is much higher than the base flow conditions, and corresponds to a very high discharge of freshwater, which only occurs during large rain events. As a result, this computation demonstrates that the required inflow of freshwater for a significant effect of riverine inputs into the SPOM composition of Copano bay corresponds to a large rise in river level (i.e., rain event, which does not correspond to this scenario anymore). This means that riverine inputs during droughts (composed by riverine fresh phytoplankton) cannot have any major influence in the SPOM composition of this secondary bay. As a result, during droughts, SPOM will always be dominated by locally produced phytoplankton.

The scenario based on rain event conditions gave an inflow threshold of 910,000 m 3 month -1 . Consequently, when freshwater inflow is higher than this threshold the SPOM in Copano bay shifts gradually from fresh autochthonous phytoplankton to highly degraded terrestrial material. Indeed, the very high C/chl a ratios of riverine SPOM drained during rain events indicate that this material is highly degraded, and is probably of very low nutritive quality for potential consumers (e.g., suspension feeders) [START_REF] Ittekkot | Global trends in the nature of organic matter in river suspensions[END_REF]. The quality of the SPOM thus decreases into Copano bay, which may affect its fate and very likely the functioning of the estuarine food web. In 2007, the large floods, which occurred in the Mission Aransas estuary watershed affected oyster filtration rates, due to low salinities [START_REF] Beseres Pollack | Role of flood disturbance in natural oyster (Crassostrea virginica) population maintenance in an estuary in South Texas, USA[END_REF]. We suggest that the poor quality of the organic matter provided to the oysters during rain events may also affect the filtration rates, and thus the growth of oysters.

Minor rain events may also affect local production of phytoplankton, as demonstrated by the peak of chlorophyll a in January [START_REF] Pennock | Nutrient behavior and phytoplankton production in Gulf of Mexico estuaries[END_REF][START_REF] Mooney | Watershed export events and ecosystem responses in the Mission-Aransas National Estuarine Research Reserve, South Texas[END_REF].

Phytoplankton production increase due to a higher discharge of nutrients and remineralization of terrestrial organic matter drained into the estuary [START_REF] Mooney | Watershed export events and ecosystem responses in the Mission-Aransas National Estuarine Research Reserve, South Texas[END_REF][START_REF] Bruesewitz | Estuarine ecosystem function response to flood and drought in a shallow, semiarid estuary: Nitrogen cycling and ecosystem metabolism[END_REF]. As a result, even if pulses of freshwater do not drain SPOM of sufficient nutritional quality for its assimilation by consumers, the release of nutrients allows for local increase of autochthonous phytoplankton production, which may locally support the functioning of food webs.

Conclusion: Functioning of a South Texas estuary

Overall, autochthonous phytoplankton was the dominant component of SPOM in the Mission-Aransas estuary, during droughts but also during the minor rain event. Benthic organic matter did not influence SPOM composition, either as fresh material or as detritus. The influence of riverine inputs was very low and restricted to areas close to river mouths, due to the high residence time in the Mission-Aransas estuary [START_REF] Solis | Hydrography, mixing characteristics, and residence times of Gulf of Mexico estuaries[END_REF].

A comparison of the functioning of the Mission-Aransas estuary with the Sabine-Neches estuary (i.e., neutral vs. highly positive estuary) reveals that freshwater inflow clearly affects composition of SPOM in these two estuaries. On a global basis, the Sabine-Neches receives high loads of low quality organic matter which are relatively quickly flushed outside of the system due to the large freshwater inflow whereas the biomass of autochthonous phytoplankton is low due to the high turbidity [START_REF] Cifuentes | Particulate organic matter in Gulf of Mexico estuaries -Implications for net heterotrophy[END_REF]. The general functioning of the Mission-Aransas estuary is different, with high phytoplanktonic production and very low inputs of refractory continental organic matter. The Mission-Aransas estuary is thus dominated by autochthonous production. Because of predicted global change, which will result in a drier climate along Gulf of Mexico coastline (Nielsen-Gammon, 2011 ; IPCC, 2013) and the increasing demand for fresh water for human activities [START_REF] Montagna | Freshwater inflow: science, policy, management[END_REF], it is therefore probable that general functioning of Texas estuaries will shift towards systems less based on allocthonous inputs of SPOM and much more influenced by SPOM from autochthonous pelagic production. These changes in organic matter quality will probably have a large effect on the functioning of food webs in Texas estuaries.

A potential shift from estuaries where SPOM is dominated by poor quality allocthonous inputs to estuaries dominated by high quality autochthonous production would also lead to issues related with nutrients, due to the importance of nutrients for phytoplankton production. A better understanding of the relationships between nutrient inputs and freshwater inflows, such as that carried out by [START_REF] Mooney | Watershed export events and ecosystem responses in the Mission-Aransas National Estuarine Research Reserve, South Texas[END_REF] in Copano bay, may be necessary at a much larger regional scale. A monitoring of these changes along the Texas coastline will require a combination of studies based on SPOM composition and nutrient loads along the coastal shoreline.

Acknowledgments

This study was partially funded thanks to a postdoctoral research grant from the Conseil Général of Charente Maritime in France. B. Lebreton is grateful to R. Kalke, M. Besonen, S. Arismendez, D. Yoskowitz, P. Richard and G. Blanchard for their very helpful advices and their support. Authors C/chl a of SPOM between sampling months at each station. Stations: ARA: Aransas river, MIS: Mission river, GoM: Gulf of Mexico. survey and from the survey carried out by [START_REF] Mooney | Watershed export events and ecosystem responses in the Mission-Aransas National Estuarine Research Reserve, South Texas[END_REF].

upplementary material for online publication only lick here to download Supplementary material for online publication only: Supplementary_data_Lebreton_et_al.docx

δ

  13 C, δ 15 N and δ 34 S values of SPOM at GoM station were stable throughout the year, ranging from to -23.2 (January 2011) to -20.8‰ (April 2011), from 6.5 (December 2010) to 8.6‰ (May 2010) and from 16.5 (May 2011) to 18.5‰ (July and January 2011), respectively (Fig.5, Table2, Appendix A).

Figure captions Figure 1 :

 captions1 Figure captions

Figure 2 :

 2 Figure 2: Salinity (full and dotted lines), monthly inputs of particulate organic carbon from Aransas

Figure 3 :

 3 Figure 3: δ 13 C and δ 34 S values of primary producers and SPOM in the Mission-Aransas estuary. SPOM

Figure 4 :

 4 Figure 4: δ 13 C values of benthic microalgae and salinities of sampling stations in the Mission-Aransas

Figure 5 :

 5 Figure 5: δ13C, δ15N, δ34S values (full and dotted lines) and total freshwater inflow (i. e., Mission

Figure 6 :

 6 Figure 6: Chlorophyll a concentration (µg.l-1), C/chlorophyll a ratios (full and dotted lines) and total

Figure 7 :

 7 Figure 7: Relationship between freshwater inflow and stable isotope composition of the SPOM in the

Table 3 :

 3 Summary of Kruskal-Wallis tests and comparisons of means for δ 13 C, δ 15 N, δ 34 S, Chl a and C/chl a of SPOM between sampling stations for each month. Stations: ARA: Aransas river, MIS: Mission river, St. 1 to 5: Stations 1 to 5, GoM: Gulf of Mexico.

			δ 13 C		δ 15 N		δ 34 S		Chl a		C/chl a
	Station	p-values	Comparisons of means	p-values	Comparisons of means	p-values	Comparisons of means	p-values	Comparisons of means	p-values	Comparisons of means
	November 10	0.007	MIS < St. 5 MIS < GoM	0.005	ARA > St. 1	0.039	-	0.004	St. 1 > St. 5	0.004	ARA < St.5 St.2 < St. 5
	December 10	0.003	MIS < St. 5 MIS < GoM	0.096	-	0.017	-	0.030	-	0.030	St. 2 > St. 4
	January 11	0.003	St. 2 < St. 5 St. 2 < GoM	0.003	MIS < St. 4	0.013	-	0.003	MIS < St. 1	0.012	-
	February 11	0.002	ARA < GoM	0.007	ARA < MIS	0.004	ARA < St. 5	0.002	ARA > MIS ARA > St. 5	0.002	ARA < St. 3
	March 11	0.003	ARA < St. 5	0.006	ARA < MIS	0.009	MIS < GoM	0.003	ARA > MIS ARA > St. 5	0.004	ARA < St. 1 ARA < St. 2
	April 11	0.004	ARA < St. 4	0.008	ARA > St. 1	0.013	-	0.003	MIS < St. 1	0.012	ARA < St. 1
	May 11	0.002	ARA < St. 5 ARA < GoM	0.003	ARA > MIS	0.025	-	0.007	ARA > St. 5 ARA > GoM	0.074	ARA < St. 3
	June 11	0.005	-	0.010	ARA > MIS	0.041	-	0.006	MIS > St. 5	0.015	-
	July 11	0.030	-	0.006	ARA < MIS ARA < GoM	0.008	-	0.003	ARA > St. 5 MIS > St. 5	0.009	MIS < St. 5
	August 11	0.007	-	0.002	ARA < St. 3 St.1 < St. 3	0.174	-	0.002	ARA > St. 5 ARA > GoM	0.070	-

thank E. Hill and K. Withers for their help in benthic algae identification and G. Wise who allowed the collection of water samples on his property. B. Lebreton also thanks E. Morgan, A. Lord, C. Santos, L. Hutchison, K. Welder, B. Paudel, L. Hyde, D. Mangold, E. Chalkley, C. Taylor, M. Proffitt, E. Buschfort, C. Hampton, C. Carollo, M. Barajas for their invaluable help during collection and preparation of samples, and G. Sutton, L. Eckert and A. Knight for their support in administrative tasks. Authors also thank the two reviewers who helped to improve the manuscript thanks to their comments.

Tables

 Table 1:Stable isotope compositions (δ 13 C, δ 15 N, δ 34 S, mean ± standard deviation, ‰) of the different benthic primary producers in the Mission-Aransas estuary. Sample size is displayed in parentheses.

Species


J.-Y Faou

Sergey Grachev

Etienne Barthel

Guillaume Parry

From telephone cords to branched buckles : a phase diagram

Keywords: Buckling, Delamination, Finite element method, Cohesive zone model, Telephone cord, mode-mixity

Thin films with low adhesion and large residual stresses may buckle and delaminate from their substrates. This delamination often results in oscillating patterns known as 'telephone cords'. Other configurations can be observed as well, such as entangled networks of blisters. The present study aims at elucidating how these networks are generated through branching. A model coupling a geometrically non-linear plate model and a cohesive zone with mode dependent interfacial toughness is used. Previous work on the early stage delamination of initially circular blisters during interfacial crack front propagation are revisited and extended to the remote post-critical regime to capture front branching. We also build a partial phase diagram of buckling driven delamination morphologies. An intrinsic critical stress parameter is proposed which can be used to predict delamination morphology and to define 'safe' conditions where buckle propagation can be avoided. This critical stress also defines an intrinsic lengthscale which plays the role of a minimum nucleation size.

List of notations

E, ν

Young's modulus and Poisson's ratio of the film h Film thickness δ ij Kronecker symbol σ ij = -σ 0 δ ij Equibiaxial stress tensor in the film prior to buckling (σ 0 > 0) E * = E/(1 -ν) For an equibiaxial loading σ 0 = E * 0 σ c Critical buckling stress ψ Mode mixity parameter G 0

Strain energy per unit area stored in the film G c

Toughness of the film/substrate interface G Ic

Pure mode I interface toughness η

Parameter defining the interface sensitivity to mode mixity T Interface traction vector (T n , T t )

Normal and tangential components of T (T 0 n , T 0 t ) Maximum normal and tangential components allowed for T (T i n , T i t ) Normal and tangential components of T actually reached at peak traction

Introduction -Phases of thin film blisters

Thin film blisters form because as the coating buckles it breaks the film/substrate interface. The shape of the buckle is eventually controlled by the coupling between these two complex problems. Buckling has been intensively studied in the framework of the Föppl-Von Kàrmàn thin plate equilibrium equations [START_REF] Landau | Theory of Elasticity[END_REF][START_REF] Hutchinson | Mixed Mode Cracking in Layered Materials[END_REF]. A coating with large compressive stresses can buckle above a critical stress σ c . In thin films the internal stress incurred during deposition can easily exceed this critical buckling stress. The post-critical regime is characterized by large out-of-plane displacements, so the geometrically nonlinear character of the plate equations plays a significant role, and analytical solutions exist only in the simplest cases. The difficulty is even greater if, instead of fixed boundaries, the process of delamination itself is taken into account. The only analytical results available in the post-critical regime are the straight-sided blister (i.e. the Euler column) and the circular blister [START_REF] Hutchinson | Mixed Mode Cracking in Layered Materials[END_REF], characterized by their simple geometries. With the help of numerical calculations, some of the mysteries surrounding the ubiquitous 'telephone cord buckle' (TCB) [START_REF] Gille | Buckling instability and adhesion of carbon layers[END_REF][START_REF] Wang | Effects of strain cycling on buckling, cracking and spalling of a thermally grown alumina on a nickel-based bond coat[END_REF][START_REF] Yu | Spatial and kinetic evolutions of telephone cord buckles[END_REF] have been unraveled. Various studies [START_REF] Parry | Stability diagram of unilateral buckling patterns of strip-delaminated films[END_REF][START_REF] Audoly | Stability of straight delamination blisters[END_REF][START_REF] Audoly | Secondary buckling patterns of a thin plate under in-plane compression[END_REF][START_REF] Moon | The characterization of telephone cord buckling of compressed thin films on substrates[END_REF][START_REF] Moon | Buckle delamination on patterned substrates[END_REF] identified the TCB as resulting from the secondary buckling of a straight-sided blister. Recently, significant progress have been made in understanding the TCB propagation mechanism by coupling geometrical nonlinearity for the plate with a cohesive zone for the interface [START_REF] Faou | How does adhesion induce the formation of telephone cord buckles?[END_REF]. In this respect, the mode mixity dependence of interfacial toughness proved to be a key point. In particular, the model allowed for a very good understanding of the relationship between adhesion and blister morphology with a good correlation with experiments [START_REF] Faou | Telephone cord buckles: A relation between wavelength and adhesion[END_REF].

The morphology of a telephone cord blister on a rigid substrate can now be anticipated given full details about elastic properties, coating thickness, residual stresses and adhesion, including mode mixity dependence. However, the question of blister bifurcation (or "branching") remains open. Indeed, a large variety of blister configurations can be observed experimentally. For example, in Fig. 1a, we show isolated telephone cords obtained with Mo thin films (see below). Decreasing film thickness and increasing residual stresses we observe entangled networks of telephone cords (Fig. 1-b) with separation of the delamination front into several branches at particular steps of the delamination process. The aim of this paper is to determine conditions under which blisters branch. The case of the initial destabilization of a circular blister, originally studied in [START_REF] Hutchinson | Growth and configurational stability of circular, buckling-driven film delaminations[END_REF], is first investigated. It contains the key features for the analysis of the delamination front destabilization and branching and serves as a reference case for our numerical study. Then we build a partial phase diagram for blister morphologies which we also illustrate with experimental observations.

Methods

The numerical scheme and experimental methods are as described in our previous papers [START_REF] Faou | How does adhesion induce the formation of telephone cord buckles?[END_REF][START_REF] Faou | Telephone cord buckles: A relation between wavelength and adhesion[END_REF].

In brief, the FEM model (Fig. 2) consists of plate elements with geometrical nonlinearity bonded to a rigid substrate.

To generate compressive residual stresses in the simulation, we start from a relaxed film and apply a uniform thermal expansion to the film with components xx = yy = 0 and xy = 0. This expansion (the so called eigenstrain, which is here positive, with 0 > 0) is in fact restrained by the rigid substrate, so that the actual total film strain is zero. At equilibrium, the eigenstrain is therefore compensated by an exactly opposite elastic contraction, which generates an equibiaxial compressive stress state

σ xx = σ yy = -E/(1 -ν) 0 = -σ 0 , σ xy = 0.
Interface debonding is described by mixed-mode cohesive elements inserted between the plate elements and the substrate. The traction vs. separation law is linear/softening. It is linear and reversible (T = Kδ where K is the stiffness) before peak traction is reached. One of the difficulties in the linear/softening model is to determine a criterion for peak traction which duly takes mode mixity into account. For that purpose we use a simple quadratic criterion:

T n T n 0 2 + T t T t 0 2 = 1 (1)
where T n and T t are the normal and tangential components of the traction T , T 0 n is the maximum traction in mode I and T 0 t is the maximum traction in mode II. Thus, for any combination of traction components (T n , T t ), peak values (T n i , T t i ) are reached when the quadratic criterion is satisfied. Once peak traction is attained, softening sets in with linearly decreasing traction. The separation energy per unit area (i.e. the work of the traction) G c (ψ), depends on the mode mixity ψ. This relation is the main component of the model. It is actually independent from the peak traction criterion. Following [START_REF] Hutchinson | Growth and configurational stability of circular, buckling-driven film delaminations[END_REF], we use a phenomenological expression for the interface toughness G c :

G c (ψ) = G Ic (1 + tan 2 (ηψ)) (2)
with G Ic the mode I separation energy (or toughness), and a parameter η used to tune the mode II dependence. Values of η close to zero imply a weak dependence on the mode mixity, while values close to unity imply a strong dependence. Note that the measure of the mode mixity ψ in our cohesive law is defined in terms of the combination of tractions attained at the peak, i.e. tan(ψ) = T i t /T i n , for each cohesive element meeting the initiation criterion. This definition differs from that employed in linear fracture mechanics in that it is defined pointwise on the interface in the vicinity of crack front, whereas the conventional definition of mode mixity, ψ, is defined in terms of the modes I and II stress intensity factors at the crack front as tan(ψ) = K II /K I (e.g. [START_REF] Hutchinson | Growth and configurational stability of circular, buckling-driven film delaminations[END_REF]). In both descriptions, a mode I situation corresponds to high bending and low stretching at the edge of the plate, whereas a mode II situation is one with high stretching and low bending. We have previously shown that there is good agreement between the present numerical model (especially the definition of mode mixity angle) and the standard results for plates ( [START_REF] Faou | How does adhesion induce the formation of telephone cord buckles?[END_REF], Supplementary Material).

The numerical implementation of the model is carried out in the framework of the finite element method using the FEM code ABAQUS [START_REF]ABAQUS Manuals Collection[END_REF].

For the experiments, we monitored buckles generated by molybdenum coatings of thickness ranging between 50 nm and 200 nm, deposited by magnetron sputtering on 300 µm thick silicon wafers with a native oxide layer. A thin Ag layer (about 10 nm thick) is deposited onto the substrate prior to Mo deposition in order to ensure a relatively low and well controlled adhesion. Delamination occurs at the Ag/SiOx interface.

Initial destabilization of a circular blister

Growth of a circular defect

Consider a film submitted to an equibiaxial compressive stress σ 0 . A circular blister of initial radius R 0 will expand as long as the normalized driving force for delamination G/G c (ψ) is larger than 1. The mode mixity dependence of the interface toughness, given here by Eq. 2, plays a key role in the existence of an equilibrium radius. The energy release rate G increases monotonically with increasing radius, but so does G c (ψ) because |ψ| increases with radius as well. As a result of this competition, the energy release rate G needed to advance the interface crack increases as the blister spreads. Crack arrest is obtained if interface toughness increases sharply enough with mode II.

Actually, circular growth may eventually become unstable as the radius increases as shown in [START_REF] Hutchinson | Growth and configurational stability of circular, buckling-driven film delaminations[END_REF] where the sensitivity of the circular crack front to harmonic perturbation and subsequent destabilization are studied using linear stability analysis. It is shown that the loading ratio, defined as the ratio of the residual stress to the critical buckling stress σ 0 /σ c , selects the number of lobes n formed along the circular front. The formation of lobes is associated with a non-uniform distribution of crack driving force along the crack front and a modulation of mode mixity. Such modifications of the crack front enable further growth of the lobes. We denote

Incipient destabilization of the circular blister

G 0 = 1 -ν E hσ 2 0 ( 3
)
the elastic strain energy per unit surface stored in the film. E is Young's modulus, ν Poisson's ratio and h the film thickness. Given a set of adhesion parameters (G Ic , η) (Eq. 2), there is a maximum ψ that can be reached at the delamination front during debonding, ψ lim . We then have G c = G 0 or:

ψ lim = 1 η arctan G 0 G Ic -1 (4)
In our simulations, a value of η = 0.95 is usually taken, and ψ lim ranges between 36 • and 82 • depending on the ratio G 0 /G Ic .

An initial debonded zone is set as a circular defect of radius R 0 . Upon compressive loading, the initial defect buckles into a circular blister when σ 0 = σ c . Here the critical buckling stress σ c is given by:

σ c = 1.2235 E 1 -ν 2 h R 0 2 (5)
An example of evolution of the blister is shown in Fig. 3, where the out-ofplane displacement is plotted at various growth stages. In the first stage of growth, the blister remains circular (Fig. 3-b). For higher values of the loading ratio, destabilization of the circular front is observed, and the blister evolves into a non-axisymmetric morphology (Fig. 3-c). In the given example we find four lobes, with a convex local geometry, separated by concave areas where the film has sagged.

In order to better understand this phenomenon, we plot the local mode mixity along the crack front (Fig. 4) as a function of the polar angle θ (see Fig. 3). Two states are reported: a circular shape before stabilisation and a 4 lobes shape after destabilization.

3). In the non-axisymmetric stage, we find alternating areas with lower and higher values of mode mixity. The mode mixity at the crack front is close to ψ lim in the concave (collapsed) regions, where the tensile contribution (mode II) is maximum; it is close to |ψ| ≈ 40 • in the convex (lobe) regions where mode I is enhanced. As interfacial toughness is mode dependent, a decrease in mode mixity promotes crack growth whereas high mode mixity keeps the front pinned.

Note that, at this early stage of destabilization, there are four identical lobes, in agreement with the n lobe perturbation used in [START_REF] Hutchinson | Growth and configurational stability of circular, buckling-driven film delaminations[END_REF].

As previously stated, the loading ratio σ 0 /σ c is the critical parameter for destabilization of the circular buckle and directly affects the number of lobes upon destabilization. In Fig. 5 we plot the number of lobes obtained as a function of σ 0 /σ c for various G 0 /G Ic and σ 0 . The analytical results of [START_REF] Hutchinson | Growth and configurational stability of circular, buckling-driven film delaminations[END_REF] are also reported. We observe that the number of lobes formed upon destabilization follows the trend found analytically ( [START_REF] Hutchinson | Growth and configurational stability of circular, buckling-driven film delaminations[END_REF]) and increases gradually with σ 0 /σ c . However, we observe that with a cohesive zone, the same number of lobes n is formed at a larger loading ratio, especially for small n values. We believe this deviation is due to a finite size effect of the crack process zone. The process zone stabilizes the edge of the buckle as the stress is less released over its width. In the numerical simulations, the process zone width Λ is about 6h for the range of interface properties (i.e. peak traction components and mode I toughness) chosen. For large loading ratios, or large radii, the ratio Λ/R diminishes and the effect of the process zone vanishes as seen in Fig. 5.

The loading ratio is a convenient parameter to characterize the destabilization of circular buckles. However it is computed from the radius of the circular Results from numerical simulation with a cohesive zone are reported for η = 0.95 and ν f = 0.31 (labelled CZM) and from [START_REF] Hutchinson | Growth and configurational stability of circular, buckling-driven film delaminations[END_REF] (labelled LEFM) for η = 1 with α = 0 and ν f = 1/3 buckle (i.e. initial defect) so that it is not an intrinsic parameter of the problem. Therefore, it cannot provide a criterion for growth of TC and branching buckles. We will propose an alternative critical parameter in section 5.

Branching in the propagation regime

In the previous section, our investigations have been restricted to the study of the very first stage of the destabilization process, where the n lobes are equivalent, following [START_REF] Hutchinson | Growth and configurational stability of circular, buckling-driven film delaminations[END_REF][START_REF] Nilsson | A finite element analysis of configurational stability and finite growth of buckling driven delamination[END_REF]. In order to understand the subsequent evolution of the blister morphologies out of this initial unstable defect, we have to explore the conditions under which lobes further propagate or come to a halt, eventually leading to the experimentally observed morphologies (Fig. 1). This involves computation of G along a rather complex crack front, which is handled by numerical calculations.

Case of the single TC blister

The TC morphology is characterized by a width 2b and a wavelength λ (Fig. 6). It was shown in a previous study [START_REF] Faou | How does adhesion induce the formation of telephone cord buckles?[END_REF] that during propagation of the buckle the crack front is pinned in the sagged areas due to the mixed-mode dependence of interfacial toughness. To investigate the evolution of the mode mixity along the buckle edge during propagation and branching, we plot ψ along the crack front (Fig. 6) between two collapsed points A and B taken on each side of the propagation, the distance AB is normalized to 1. The first curve corresponds to the state of the crack front after formation of point A. We observe the high mode mixity, |ψ| ≈ ψ lim in the sagged areas around A and B. Between A and B however, |ψ| is lower. Note that the mode mixity is not constant in these areas, unlike the case of the circular blister destabilization for which a plateau is observed. This observation suggests that the propagating front of the TC buckle is not quite equivalent to a circular blister. After propagating further, we observe the formation of another area with |ψ| ≈ ψ lim around point C, which becomes a new pinning point. Interestingly, point C does not appear at middistance between A and B. The position of C varies during propagation but is definitely closer to B. In addition, the average mode mixity over AC is lower (in absolute value) than between BC. We believe this difference is the source of arrest of the front BC whereas front AC will keep propagating and form a new undulation. We have found a similar kinematics for all the TC buckles we have simulated.

A critical value for the formation of telephone cord buckles seems to be correlated to the value of the loading ratio σ 0 /σ c upon destabilization of the circular buckle. In fact, no TC buckle was observed for values below σ 0 /σ c ≈ 16 nor above σ 0 /σ c ≈ 22, a rather limited range.

Case of branching blisters

Above σ 0 /σ c ≈ 22, the growth of lobes still occurs but in this case the morphology obtained, while still periodic, develops a branching structure. Separation of the crack front into two branches is visible on both experimental propagation and numerical simulations on Fig. 8. It is interesting to note that the transition between isolated TC buckles and branching structures is not sharp, as it is possible to simultaneously observe the two structures under the same conditions, both numerically and experimentally (Fig. 1). To better understand branching, we plot ψ along the propagating crack front of a buckle (Fig. 9) for the three different states depicted in Fig. 8. Locations of points A, B, C and D are indicated in Fig. 8. The distance between AB is normalized to 1. In the first frame, captured after the formation of pinning point A, the mode mixity is rather constant, except in the neighbourhood of pinning points A and B, and similar to the plateau state observed for the lobes upon destabilization of a circular buckle (Fig. 4). After propagation, we observe the formation of pinning point C between A and B similar to the TC case. However the position of point C differs as it is right at mid-distance between A and B. The consequence of the position is that the average |ψ| on both sides of C is lower than in the TC case. As a result none of the crack fronts arrest: they both propagate. Later on, the symmetry between the two propagating branches is broken as another pinning point D appears between C and B. The consequence is a lower |ψ| along CD than DB. This difference causes arrest of the DB front while both AC and CD fronts keep propagating.

A phase diagram for thin film blisters

The loading ratio impacts both the wavelength of the TC buckle and the ability of the front to experience branching. This combination of phenomena gives rise to the wide variety of morphologies observed experimentally. Starting from a rectangular defect, and by varying the parameters of the model, we have simulated a variety of buckling morphologies shown in Fig. 10 under the form of a phase diagram. Each point indicates a calculated morphology for different h, σ 0 and G Ic .

We observe that two simple reduced parameters provide an adequate description of this phase diagram. The fist parameter is G Ic /G 0 which is the ratio of the interface toughness to the total elastic energy stored by unit surface in the film. The second parameter is directly the normalized residual stress σ 0 /E * , with E * = E/(1 -ν). With these two parameters, we find a clear separation between three regions (Fig. 10). For large values of mode I interfacial toughness G Ic , the initial defect is stable: no propagation is observed. By decreasing G Ic while keeping the stress σ 0 constant, we first find the TC blister followed by branching blister and complete delamination when interfacial toughness decreases further.

In this diagram, the boundary between the 'no delamination' and the 'telephone cord' areas can be approximated by a straight line. The condition for buckle delamination to occur can be expressed as:

σ 0 E * ≥ c G Ic G 0 ( 6
)
where c is a positive constant. The limit between TC buckles and branched structures is slightly more difficult to define. Note that the boundary may move a little depending on the mesh size adopted in the simulation. Here the element sizes were scaled with h and Ē/σ 0 so that the number of elements along R for an equivalent circular blister would remain constant. Also, two boundaries may be considered: 1) first appearance of a branch in the morphology or 2) complete disappearance of the TC morphology. One possible boundary is plotted in Fig. 10. We find that this boundary is also linear in the present parameter space. The value η = 0.95 that has been chosen here characterizes an interface toughness that is highly sensitive to mode mixity. In this regime, the choice of η is not expected to qualitatively influence the phase diagram. To evaluate quantitative modifications, the two cases η = 0.95 and η = 0.8 are plotted in Fig. 11. The linear boundaries are scaled up by roughly the same factor (about 1.25) when η decreases from 0.95 to 0.8.

Discussion

In view of the previous results, inequality 6 can be rewritten as σ 0 ≥ σ * c with

σ * c = c 1 (E * ) 2/3 G Ic h 1/3 (7)
which defines a new critical stress σ * c . For η = 0.95 and ν = 0.31, c 1 ≈ 0.23. It is interesting that, unlike the critical buckling stress σ c (Eq. 5), this new critical stress is intrinsic since it does not require the introduction of an extrinsic reference length (e.g. radius of a defect). In fact, Eqs. 5 and 7 define a characteristic size L * with

L * h = E * σ * c E * h G Ic 1/6 (8)
Rough estimates suggest that in our systems, L * is of the order of 10 times the film thickness. The prefactor varies weakly with the model parameters. This intrinsic lengthscale is quite different from the telephone cord wavelength we have defined previously [START_REF] Faou | Telephone cord buckles: A relation between wavelength and adhesion[END_REF]. Here a natural interpretation is that for defect sizes larger than L * , buckles will propagate while smaller defects will remain stable. It does not mean that the initial defect can be taken arbitrarily though.

If the initial defect is larger than the resulting blister size, the initial mode mixity will be too large for propagation (G c > G 0 ). But provided that the initial defect is taken small compared to the resulting buckle size after propagation, its size does not influence the final result. We observe that for σ 0 /σ * c = 1.07 the first branching events appear in the phase diagram. TC blister are still observed up to σ 0 /σ * c = 1.15, above which almost only branching blisters are observed. We point out that for identical thickness and interfacial toughness, a small variation of stress may lead to a strong variation of the observed morphology. Such variations are common in the Mo thin films presented above.

Comparing the positions of the phase boundaries for the two values of η we have investigated, we observe that for a given value of G Ic /G 0 , σ * c seems to decrease with η. In other words, smaller values of η promote destabilization. This observation is consistent with earlier results [START_REF] Hutchinson | Growth and configurational stability of circular, buckling-driven film delaminations[END_REF] which show that harmonic destabilization occurs at lower compression for smaller η.

Interestingly, we did not observe any propagation of straight sided blister in the numerical simulations for equi-biaxial loaded films with the mode II dependence chosen (η = 0.95). This result is in accordance with Jensen [START_REF] Jensen | Effects of residual stresses in the blister test[END_REF] who predicted almost immediate instability of a straight front for such loading. We expect the straight sided morphology to appear for smaller vales of η for low values of G I c/G 0 and σ 0 /E * . It seems from the diagram that the domain of stability of the TC buckle is narrow. However, experimentally, TC buckles are observed over a wide range of thicknesses and stress levels. This is not contradictory with the results of or analysis, but is actually related to the way stresses build up in thin films during deposition. For example if we assume that the stress is kept constant and variation of G 0 arises solely from the thickness, TC buckles may be observed over a range of thicknesses. If h 0 is the thickness of the film at the onset of propagation, TC buckle will be observed up to 1.5h 0 .

To illustrate this effect, in Fig 12 we plot the evolution of G Ic /G 0 as a function of σ 0 /E * for two experimental molybdenum thin films taken from a previous study [START_REF] Faou | Stress tuning in sputter-deposited MoOx films[END_REF]. The superlayer is deposited on top of an interface with an interfacial toughness in mode I of G Ic = 0.4 J.m -2 and it is assumed that no relaxation occurs in the thin film after deposition.

-For MoO x thin film, the range over which stable TC buckles are predicted spans from 195 nm to 283 nm, which is consistent with a rapid destabilization of the TC buckle.

-For pure Mo, TC buckles are expected from 230 nm up to 670 nm. This larger range of stability is due to a strong decrease of the average residual stress through the thickness and a small increase in G 0 , which is common in Mo thin films past the first nanometers.

Figure 12 also shows that in the case of the Mo film, given the evolution of the stress and the thickness during the deposition process, it is actually possible to enter and exit the telephone cord buckles stability area to end up inside the 'no delamination' area. It is hence possible to get thick unbuckled films. Indeed, it is usually observed that the buckles develop once the sample have been taken out of the deposition chamber ; so the final state reached by the sample is likely to be the most critical. These observations can be related to the nucleation process of buckle delamination, which is still poorly understood. 10). The thicknesses at the crossing of the boundaries between the stability domains are displayed.

Conclusions

A model coupling non-linear plate buckling and cohesive zone has been used to study buckling driven delamination under equi-biaxial compressive stress and mode-dependent interfacial toughness. The results obtained numerically are well correlated with experimental observations. The critical role of plate instability resulting in interface pinning is identified not only for initial circular defect, but also latter on in the process of buckle propagation for both undulating and branching mechanisms of the crack front.

Both undulating and branching mechanisms are very similar and differ only by the different position of the pining point C when the instability develops. In the case of the standard TC, the locus of the point C breaks the symmetry and leads to selection between an arrested branch and a propagating one. In the case of branching, the pinning point C appears right in the middle of the two candidate branches, which both propagate further. Those two mechanisms are responsible for the generation of the patterns observed experimentally, in particular the telephone-cord like blister.

It has been possible to obtain a phase diagram of the buckling driven delamination morphology based on those two mechanisms. Two boundaries have been identified from the phase diagram: propagation of a crack front from an initial defect and branching of the crack front. With the help of the diagram, we have been able to define a critical stress that is intrinsic to the system, i.e. only depending on the film elasticity modulus, thickness and stresses as well as adhesion. This critical stress also defines an intrinsic size L * which may play the role of a minimum size for propagation. The boundaries on the diagram, which can be determined from the intrinsic data of the system, are of critical importance, whether to avoid propagation of buckles on experimental samples or to control the buckling induced patterning in some way (branching or not). A practical use of the analysis results has been illustrated in the specific case of Mo and MoO x films.

Figure 1 :

 1 Figure 1: Buckles formed with molybdenum coatings deposited on a silicon wafer with a thin Ag interlayer. Two buckling delamination patterns observed experimentally: (a) isolated telephone cords blisters; (b) entangled networks of telephone cords (branching blisters) for reduced film thickness.

Figure 2 :

 2 Figure 2: Model for thin film buckling and delamination: a nonlinear plate on a rigid substrate with a mixed mode cohesive zone model for the adhesive interface. w indicated the out-ofplane deflection of the film due to buckling.

  σ 0 /σ c = 8.7 σ 0 /σ c = 16.7 σ 0 /σ c = 33.5

Figure 3 :

 3 Figure 3: Various stages of growth of a circular defect submitted to an increasing σ/σc loading ratio:(a) Circular growth. (b) Maximum radius before destabilization. (c) Post-destabilization morphology of the circular defect with four lobes. The polar angle θ used for the description of the mode mixity along the crack front is also reported.

Figure 4 :

 4 Figure 4: Local mode mixity ψ along the crack front of an initially circular blister, as a function of the polar angle θ (see Fig.3). Two states are reported: a circular shape before stabilisation and a 4 lobes shape after destabilization.

Figure 5 :

 5 Figure5: Number of lobes n resulting from the unstable growth of an initially circular blister. Results from numerical simulation with a cohesive zone are reported for η = 0.95 and ν f = 0.31 (labelled CZM) and from[START_REF] Hutchinson | Growth and configurational stability of circular, buckling-driven film delaminations[END_REF] (labelled LEFM) for η = 1 with α = 0 and ν f = 1/3

Figure 6 :

 6 Figure 6: Out-of-plane displacement during propagation of a telephone cord buckle. Three points are introduced to support discussion about the mode mixity evolution along the crack front. The front originally developing between the two pinning points A and B progressively extends. A new pining point C eventually forms (G Ic /G 0 = 0.34 and σ 0 /E * = 6.1 × 10 -3 ).

  TC and forming half a wavelength of a TC buckle. Point B is located behind A along the propagation trail of the TC. As the crack front moves forward during

Figure 7 :

 7 Figure 7: Evolution of the mode mixity parameter ψ along the propagating front settling between the two pinning points A and B (see Fig. 6). The crack front length between A and B is normalized. Two states are depicted: before and after the development of a new pinning point C.

Figure 8 :

 8 Figure 8: Detailled snapshot of the out of plane displacement of a blister during branching at the front both numerically (top) and experimentally (bottom). In the experiment, the fracture occurs at the SiOx/Ag interface, with 15 nm of Ag loaded by 132 nm of MoOx (σ 0 = 1.84GP a). For the simulation, the parameters are chosen so that G Ic /G 0 = 0.34 and σ 0 /E * = 6.1 × 10 -3 . Pinning points are introduced in the top picture to support discussion about the evolution of mode mixity along the crack front during propagation.

Figure 9 :

 9 Figure9: Evolution of the mode mixity parameter ψ along the propagating front settling between the two pinning points A and B (see Fig.8). The crack front length between A and B is normalized. Three states are depicted: a uniform plateau of ψ between A and B, the development of a third pinning point C in the middle and finally the development of a fourth pinning point D breaking the previous symmetry.

Figure 10 :

 10 Figure 10: Phase diagram of the delamination morphology with η = 0.95 and ν = 0.31. The morphology of delamination is indicated (symbols) and the calculated limit for propagation from an initial (solid line) defect and branching of a delamination front (dashed line) are figured. ND : No delamination, TC : Telephone cord, B : Branching, CD : Complete delamination.

Figure 11 :

 11 Figure 11: Phase diagram of the delamination morphology with η = 0.8. The lines defining the boundaries between the stability regions have been inserted. The boundary lines obtained in the case η = 0.95 have been recalled.

Figure 12 :

 12 Figure 12: Possible energy range available with two different superlayer, Mo and MoOx at -75 V assuming an interface with G Ic = 0.4 J.m -2 . The experimental data are superimposed to the phase diagram (Fig.10). The thicknesses at the crossing of the boundaries between the stability domains are displayed.

Corresponding author. Email address: guillaume.parry@simap.grenoble-inp.fr Preprint submitted to Acta Materialia December 12,