On biotic and abiotic drivers of the microphytobenthos seasonal cycle in a temperate intertidal mudflat: a modelling study

Microphytobenthos (MPB) from intertidal mudflats are key primary producers at the land-ocean interface. MPB can be more productive than phytoplankton and sustain both benthic and pelagic higher trophic levels. The objective of this study is to assess the contribution of light, mud temperature, and gastropod Peringia ulvae grazing pressure in shaping the seasonal MPB dynamics on the Brouage mudflat (NW France). We use a physical-biological coupled model applied to the sediment 5 first centimeter for the year 2008. The simulated data compare to observations including time-coincident remotely sensed and in situ data. The model suggests a MPB annual cycle characterized by a main spring bloom, a biomass depression in summer, and a moderate fall bloom. In early spring, high simulated photosynthetic rates due to mud surface temperature (MST) values close to the MPB temperature optimum for photosynthesis and to increasing solar irradiance trigger the onset of the MPB spring bloom. After the bloom, high MST values lead to synoptic events when MPB thermo-inhibition (39.5 % of summer) 10 and limitation by P. ulvae grazing (8.7 % of summer) superimpose. During these synoptic events, 14 % of the simulated annual MPB primary production is channeled towards the P. ulvae secondary production through ingestion. The model suggests that such a combined effect is highly linked to the MPB biomass depression in summer. The model ability to infer on biotic and abiotic mechanisms driving the seasonal MPB dynamics could open the door to a new assessment of the export flux of biogenic matter at the land-ocean interface and, more generally, of the contribution of productive intertidal biofilms to the coastal carbon 15


Introduction
Coastal and nearshore waters receive large amounts of organic matter and inorganic nutrients from land that support high biological productivity (Mann, 1982;Admiraal, 1984;Hopkinson and Smith, 2005). However, the high turbidity of estuarine-influenced coastal waters limits the penetration of downward solar irradiance in the water column and, as such, phytoplankton production (Cloern, 1987;Struski and Bacher, 2006). In subtidal and intertidal zones, primary production (PP) sustained by benthic microalgae, or microphytobenthos (MPB), can exceed that of phytoplankton (Underwood and Kromkamp, 1999;Struski and Bacher, 2006). MPB are mostly composed of free motile epipelic diatoms and of epipsammic diatoms that live in close association (attached or free-living) with sediment grains (Round, 1971). Epipelic MPB are associated with fine cohesive intertidal sediments and develop within the top few millimetres . During daytime exposure, they migrate toward the sediment surface, constituting a dense biofilm of a few hundred micrometres Published by Copernicus Publications on behalf of the European Geosciences Union. (Herlory et al., 2004). They are fully exposed to solar irradiance at low tide, promoting PP that can reach values as high as 1.9 g C m −2 d −1 (Underwood and Kromkamp, 1999). During the flood, epipelic MPB move downward within the sediment but can be resuspended into the water column (Demers et al., 1987;van Beusekom, 1992, 1995;Lucas et al., 2001;Orvain et al., 2004;Ubertini et al., 2012). Both epipelic and epipsammic MPB are a key resource for higher trophic levels from benthic fauna to birds on bare mudflats (Herman et al., 2000;Kang et al., 2006;Jardine et al., 2015), but also for pelagic organisms such as zooplankton and planktivorous fishes (Perissinotto et al., 2003;Krumme et al., 2008).
On intertidal mudflats, MPB PP rates are mainly constrained by solar irradiance and temperature . The MPB biofilm faces strong daily and seasonal variations of mud surface temperature (MST) caused by heating through solar irradiance during low-tide emersion periods (Harrison and Phizacklea, 1985;Harrison, 1985;Guarini et al., 1997) and develops phenological adaptations. Blanchard and Cariou-Le Gall (1994), Barranguet et al. (1998), and Pniewski et al. (2015) showed a light-related seasonal adjustment of photosynthetic parameters (the photosynthetic capacity P b max and the light saturation parameter E k ) from production-irradiance (P -E) curves fitted to the model of Platt and Jassby (1976). Photo-inhibition was rarely observed in the field since epipelic diatoms can achieve "micro-migrations", i.e. a negative phototaxic shortterm change of position in the sediment Perkins et al., 2001;Cartaxana et al., 2011). With respect to mud temperature, Blanchard et al. (1996) related P b max to temperature mathematically. Using this relationship, Blanchard et al. (1997) showed that P b max varies according to seasons, suggesting a thermo-inhibition process in response to high mud temperature (> 25 • C). de Jonge (1980) also showed seasonal variations of the carbon (C) to chlorophyll a (Chl a) ratio, which is a proxy for the physiological state of autotrophic cells as a function of air temperature (de Jonge et al., 2012). Regarding nutrients, their limiting role on MPB growth and photosynthetic rate is not evidenced in fine cohesive sediments naturally enriched both from within the sediment and the water column Cadée and Hegeman, 1974;Admiraal, 1984). Vieira et al. (2016) suggested a likely in vitro limitation by dissolved inorganic carbon within biofilms. Benthic diatoms were shown to store ammonium and phosphate within the intracellular matrix (García-Robledo et al., 2010;Yamaguchi et al., 2015), potentially usable for assimilation and growth (García-Robledo et al., 2016). The nutrient limitation of MPB is still in debate.
Along the French Atlantic coast, the spring bloom and summer depression observed in the Brouage mudflat in the Marennes-Oléron Bay are explained by optimal temperature conditions and thermo-inhibition, respectively . Reported differences in the observed MPB seasonal cycles are also attributed to the benthic diatom assemblage (Underwood, 1994). In terms of biomass, epipelic diatoms associated with muddy sediments show a higher seasonality caused by marked exposure to stressful environmental conditions (e.g. cycle of deposition-erosion, desiccation, grazing) than less motile epipsammic species in coarser sandy sediments (Underwood, 1994). In summer, thermo-inhibition and high grazing pressure by deposit feeders are suggested to dampen the MPB biomass (Cadée and Hegeman, 1974;Cariou-Le Gall and Blanchard, 1995;Sahan et al., 2007). On intertidal mudflats, the prosobranch gastropod Peringia ulvae can reach densities up to 30 000 snails m −2 (Sauriau et al., 1989) with a reported maximal ingestion rate of 385 ng Chl a snail −1 h −1 (Coelho et al., 2011). Such grazing activity may translate into a theoretical uptake of 12 g C m −2 d −1 for a C : Chl a ratio of 45 g C g Chl a −1 (Guarini, 1998), which is 6-fold more than the daily maximum MPB PP rate reported for MPB (Underwood and Kromkamp, 1999).
The role of each individual abiotic or biotic factor involved in the MPB short-term dynamics is well documented (e.g. Admiraal, 1977;Admiraal et al., 1983;Blanchard and Cariou-Le Gall, 1994;Montagna et al., 1995;Blanchard et al., 1997;Feuillet-Girard et al., 1997;Barranguet et al., 1998;Light and Beardall, 2001;Pinckney et al., 2003;Coelho et al., 2009;Weerman et al., 2011;Dupuy et al., 2014;Pniewski et al., 2015;Barnett et al., 2015;Cartaxana et al., 2015;Vieira et al., 2016). However, and in light of current knowledge, the quantitative contribution of combined factors in the seasonal MPB dynamics remains uncertain. This impedes any future assessment of how global change might impact the MPB dynamics and carbon cycle in the land-ocean continuum. The goal of this study is to quantify the relative contribution of light, temperature, and grazing to the MPB seasonal cycle and production on an intertidal mudflat (Marennes-Oléron Bay) of the French Atlantic coast. For this purpose, we use a two-layer physicalbiological model representing the MPB and P. ulvae compartments to assess the contribution of the three drivers over an annual cycle. In the paper, we describe first the physicalbiological coupled model and the in situ and remotely sensed data used to investigate the MPB seasonal cycle. Second, we assess the relative contribution of light, MST, and P. ulvae grazing on MPB dynamics and PP, and we analyse the model sensitivity to key biological constants. Finally, we discuss the role of light, temperature, and grazing in the MPB seasonal cycle and the future challenges of modelling the MPB contribution to the carbon cycle at the land-ocean continuum.

Material and methods
The study area is the Pertuis Charentais sea on the French Atlantic coast. It is a shallow semi-enclosed sea characterised by semi-diurnal tides and a macrotidal regime. The tidal range is ∼ 6 m during spring tides. The intertidal zone has two main mudflats composed of fine cohesive sediments, i.e. the Brouage mudflat (42 km 2 ) and the Aiguillon mudflat (28.7 km 2 ) (Fig. 1). The study site (45 • 54 50 N, 01 • 05 25 W) is located on the Brouage mudflat ( Fig. 1). It is composed of fine cohesive sediments (median grain size 17 µm and 85 % of grain with a diameter lower than 63 µm; Bocher et al., 2007) and sheltered from Atlantic swells by the Ile d'Oléron (Pascal et al., 2009).

Observations
A large multiparametric dataset of physical and biological measurements collected in the Pertuis Charentais was used to constrain the model and to compare with the model outputs. We provide here a summary of the data used along with their respective references, within which a detailed methodology of each set of measurements can be found.

In situ data
Atmospheric and hydrological forcings were required to set the temperature and light environment that constrained the physical-biological model. Atmospheric forcings (Fig. 2ae) consisted of meteorological observations (shortwave radiation, air temperature in the shade, atmospheric pressure above the sea, wind speed, and relative humidity) acquired at the Meteo France weather station located near the airport of La Rochelle (46 • 10 36 N, 1 • 11 3 W; data available online at https://publitheque.meteo.fr, last access: 4 December 2018; Fig. 1). Hydrology was represented by the absence or presence of seawater at the study site of the Brouage mudflat. Emersion-immersion periods were determined by the observed water height at the tide gauge of La Rochelle-La Pallice (46 • 9 30 N, 1 • 13 14 W; data Service Hydrographique et Océanographique de la Marine -SHOM-Grand Port Maritime La Rochelle-La Pallice; data available online at http://data.shom.fr/, last access: 4 December 2018) corrected by the bathymetry at the study site. The bathymetry (3.204 m above chart datum) was extracted from a digital elevation model (Litto3D ® 2010; Charente Maritime by the Institut National de l'Information Géographique et Forestière (IGN) and the SHOM) at pixels corresponding to the study site (Fig. 1). The weather and tide gauge stations were located ∼ 30 km away from the study site. Atmospheric and hydrological forcings were 1 h frequency from 1 January 2008 (00:00 UTC) to 31 December 2008 (23:00 UTC). They were linearly interpolated at the time step of the model (6 min).
In order to validate the model, we used daily measurements of MST (first centimetre of sediment), Chl a concentration (first centimetre of sediment), and Peringia ulvae biomass and density from a multiparametric dataset collected on 16-24 February and 13-26 July 2008 at the study site where the model was run (45 • 54 50 N, 01 • 05 25 W; Fig. 1). The sampling protocol is fully detailed in Orvain et al. (2014). In addition to the 2008 dataset, we used data of in situ MPB Chl a concentration collected within the first centimetre of sediment at the same station on 19-22 April 2012, 5 July 2012, 14 November 2012, 11 February 2013, and 10 April 2013. The sampling protocol is fully detailed in Lavergne et al. (2017). Monthly data of P. ulvae abundance and biomass sampled monthly from April 2014 to July 2015 over the Aiguillon mudflat were used to estimate a monthly averaged individual weight. The monthly averaged individual weight was used to convert the simulated biomass per unit of surface into density per unit of surface. The sampling protocol is given in Bocher et al. (2007). We spatially averaged the P. ulvae abundance and biomass data to obtain a monthly mean value for the entire mudflat. Ash-free dry mass (AFDM) was converted to carbon using the relationship derived from Jansson and Wulff (1977) and Remmert (2013) and used by Asmus (1994) for benthic deposit feeders (1 g AFDM = 0.58 g C). When the individual weight was not available, the individual height was used to estimate the AFDM (mg) using the formulation of Santos et al. (2005): where H is the total individual height (mm). L2G Global 250 m SIN Grid product (MOD09GQ) contains 250 m surface reflectance in a red band (620-670 nm, band centre at 645 nm) and a near-infrared band (841-876 nm, band centre at 859 nm). Terra data were used because the morning pass (10-11 h Universal Time) is better adapted than Aqua MODIS data to observe spring low tides at our study site. The data were corrected for atmospheric effects (aerosol, water vapour) and each image was checked for clouds-cirrus and cloud shadows. Cloud-free low-tide scenes were selected to apply a vegetation index. Images were reprojected to the UTM/WGS84 coordinate system. The normalised difference vegetation index (NDVI; Tucker, 1979) was calculated with the reflectance (ρ) in the red (R) and near-infrared (NIR) bands: The NDVI thresholds proposed by Méléder et al. (2003) to identify MPB with SPOT images was adapted for MODIS data and a range of 0 to 0.35 was used in this study. Negative NDVI values were associated with water and null values to bare sediment, while values higher than 0.35 corresponded to macrophytes (macroalgae and seagrass). For the present study, an NDVI time series was extracted for 2008 (47 scene images) at pixels corresponding to the study site ( Fig. 1). Scene images were processed with the ENVI ® software.

The coupled physical-biological model
The coupled model consisted of a mud temperature model coupled to a three-compartment biological model. The mud temperature model was a thermodynamic model developed by Guarini et al. (1997) resolving heat fluxes at the surface in a 1 cm thick sediment layer. Equations are given in Appendix A and Table A1. It was calibrated and validated on the Brouage mudflat by Guarini et al. (1997). During exposure periods, the simulated MST resulted from heat exchanges between the sun, the atmosphere, the sediment surface, from the conduction between mud and air, and from evaporation (Fig. 3). The MST was set to the temperature of the overlying seawater during immersion periods. The seawater temperature was simulated according to heat fluxes resulting from thermal conduction between air and seawater, from upward seawater radiation, and from downward solar and atmospheric radiation. The simulated mud temperature was considered homogeneous at the horizontal scale. The heat fluxes were determined according to equations given in Table A1 (Appendix A). The MST differential equation (Eq. A1 in Appendix A) was solved with an Euler-Cauchy algorithm at a 30 s time step. The mud temperature model constrained a threecompartment biological model, which was modified from Guarini (1998) and Guarini et al. (2000). It is fully detailed in Appendix B. MPB was represented by two compartments, including the Chl a concentration in the first centimetre of R. Savelli et al.: On biotic and abiotic drivers of the microphytobenthos seasonal cycle 7247 Figure 3. Conceptual scheme of heat exchange at the mud surface in the intertidal zone. Fluxes contributing to heat energy balance are represented by arrows during emersion and immersion periods. Modified from Guarini et al. (1997). sediment (F , mg Chl a m −2 ) and the Chl a concentration within the surface biofilm (S, mg Chl a m −2 ). The variable S * represented the S compartment that incorporated the S instantaneous production of biomass (mg Chl a m −2 ), which is directly transferred to F . The model assumed no sediment erosion or deposition and no horizontal movement of MPB within the sediment. It included a scheme of MPB vertical migration between the S and F compartments (Guarini, 1998;Guarini et al., 2000). The migration scheme is summarised in Table 1. The MPB growth rate was constrained by the photosynthetically active radiation (PAR) intensity, the simulated MST, and the grazing pressure. The grazing pressure was represented through a new scalar, Z, representing the P. ulvae biomass (mg C m −2 ). P. ulvae is a very abundant MPB grazer on the Pertuis Charentais intertidal mudflats (Sauriau et al., 1989). The P. ulvae growth rate was constrained by the simulated MPB biomass and the MST. The fourth-order Runge-Kutta method was used to solve the biological differential equations with a 6 min time step.
The coupled physical-biological model was run at the study site ( Fig. 1) from 1 January to 30 December 2008. Initial conditions were 100 mg Chl a m −2 for F and 1000 mg C m −2 for Z. No biomass was set for S at the beginning of the simulation as it started at midnight (i.e. no light). The MST was initialised at the seawater temperature (see Eqs. A5-A9 in Appendix A) at the first period of immersion. A 2008 10-year spin-up was performed before the analysis of the model outputs. The spin-ups and initial biomass conditions allowed for convergence towards similar values of biomass at the end of each run.
We performed a sensitivity analysis to quantify how simultaneous variations of key biological constants might im- pact the simulated MPB production. A Monte Carlo fixed sampling method (Hammersley and Handscomb, 1964) was used to randomly select values of the temperature optimum for photosynthesis (T opt ), the temperature maximum for photosynthesis (T max ), the optimal temperature for grazing (T opt Z ), the shape parameter of temperature-related grazing (α Z ), the light saturation parameter (E k ), and the halfsaturation constant for light use (K E ) within observed ranges (Table 2). A total of 10 000 model runs was performed with the same previous initial conditions. Statistical metrics on simulated annual PP according to parameter values and variations (Spearman's correlation coefficient and parameter average, normalised standard deviation, minimum and maximum) were computed. In addition to the simultaneous variations of parameters, the effect of a gradual variation of each single parameter on the MPB production was investigated. Each single parameter varied, while the others were fixed at the value set by default in the model (T opt = 18 • C, T max = 38 • C, E k = 100 W m −2 , K E = 20 Ein m −2 d −1 , T opt Z = 20 • C, α Z = 15).

Mud surface temperature
The simulated MST followed the seasonal cycle of air temperature (Pearson's r = 0.85, p value < 0.05; Figs. 2d and 4). From November to April, the simulated mud temperature was 9.7 ± 2.6 • C on average. The simulated average temperature from May to October reached 18.3 ± 3 • C. The amplitude (i.e. the difference between the seasonal maximum and the minimum value) of the simulated mud temperature was higher from May to October (32.1 • C) than from November to April (18.1 • C). At the synoptic scale, the model reasonably simulated the high-frequency (1 min) variations of MST measured at the study site in February and July 2008 (RMSE = 2.7 and 1.7 • C, respectively; Fig. 4). Table 1. Conceptual schemes and differential equations of the biological model, including the MPB biomass within the sediment first centimetre (F ), the MPB biomass within the biofilm (S), and the biomass of P. ulvae (Z). The upper case corresponds to daytime emersion periods when MPB cells migrate at the sediment surface () to produce and transfer biomass to the sediment first centimetre (). The middle case corresponds to the day or night-time emersion period when MPB cells migrate down to the sediment first centimetre (). The lower case corresponds to immersion periods when MPB cells are chronically resuspended from the first centimetre to the water column () and the remaining MPB cells within the biofilm finish their downward migration (). P. ulvae grazing is only active during emersion periods (upside down in immersion scheme) (modified from Guarini et al., 2008).

Scheme
Cases Equations  Fig. 5). Then the Chl a concentration decreased to reach a seasonal minimum in July (48-191 mg Chl a m −2 ; Fig. 5). The total MPB biomass (S + F ) simulated by the model within the first centimetre of sediment was the lowest in January and September (∼ 30 and 40 mg Chl a m −2 , respectively) and reached a seasonal maximum in March (∼ 266 mg Chl a m −2 ; Fig. 6a). The simulated seasonal maximum and minimum of MPB biomass during spring and summer were consistent with the observations of 2008 and 2012-2013 (Fig. 5). The model reproduced the fortnightly tidal cycle with maximum values of MPB biomass simulated in spring tides (Fig. 6a). The simulated values of MPB biomass were compared to 2008 time-coincident observations ( Fig. 6a). In February 2008, the simulated biomass was about 140.7 ± 27.7 mg Chl a m −2 , which was close but significantly higher compared to the measured total MPB biomass (106.5 ± 11.3 mg Chl a m −2 ; Mann-Whitney test: p value < 0.05). In July 2008, the model also overestimated (68.1 ± 4.5 mg Chl a m −2 ) the observed (58.6 ± 10.3 mg Chl a m −2 ) MPB biomass (Mann-Whitney test: p value < 0.05). Nevertheless, the simulated values reasonably compared, on average, with the match-up measurements gathered. The simulated daily mass-specific photosynthetic rate followed a seasonal pattern similar to that of the simulated Chl a, with values higher in late winterspring (0.56 ± 0.1 mg C (mg Chl a) −1 h −1 ) than in summer (0.41 ± 0.06 mg C (mg Chl a) −1 h −1 ) and fall-early winter (0.29 ± 0.14 mg C (mg Chl a) −1 h −1 ) (Fig. 6b).
The observed seasonal cycle of MPB retrieved from the NDVI time series was compared to the biomass simulated in the biofilm (S * ). The daily maximum values of S * simulated by the model for 2008 were subsampled to match the 2008 NDVI time series data (Fig. 7). Three distinct seasonal phases were identified in both time series using the amplitude of the sign change of the S * and NDVI second-order time derivatives (Fig. 7). Phase 1 corresponded to the spring bloom during which the biomass in the biofilm and the NDVI data reached their seasonal maximum value (day 1 to 144 and day 1 to 158 in the NDVI and model data, respectively). Phase 2 coincided with a summer depression in the simulated MPB biomass and NDVI data (day 145 to 270 and day 159 to 263 in the NDVI and model data, respectively). Finally, phase 3 showed an increase in both the simulated biomass and NDVI values, suggesting a fall bloom (day 271 to 365 and day 264 to 365 in the NDVI and model data, respectively). With respect to the NDVI data, the model showed a 14-day and 7-day longer spring and fall bloom, respectively, and a 21-day shorter summer depression ( Fig. 7). Overall, the seasonal cycle of the simulated MPB biofilm compared reasonably to that depicted by the remotely sensed NDVI data.
Biological parameters simulated by the model were compared to observed ranges reported in the literature ( Table 3). The yearly averaged value of S * simulated by the model (27.2 ± 3.6 mg Chl a m −2 ) was in agreement with the value given by Herlory et al. (2004, 24 ± 5 mg Chl a m −2 ;). The yearly averaged MPB gross growth rate (µ) simulated within the biofilm was 0.25 ± 0.07 d −1 with values ranging between 0.05 and 0.41 d −1 , which compared to the observed growth rate (0.035-0.86 d −1 ; Table 3). In the model, the MPB growth rate was related to the C : Chl a ratio (see Eq. B8 in Appendix B2). The simulated C : Chl a ratio (16 and 75.5 g C g Chl a −1 ) varied among the observed ranges (18.7-80 g C g Chl a −1 ; Table 3). The simulated annual and daily MPB PP rates (127 g C m −2 yr −1 and 369 ± 281 mg C m −2 d −1 , respectively) were also consistent with the reported in situ estimates (142 ± 82 g C m −2 yr −1 and 690 ± 682 mg C m −2 d −1 , respectively).
In the model, a linear loss term representing the resuspension process was applied to the MPB biomass simulated within the first centimetre of sediment (F compartment; see Appendix B1). In average over a high tide, 1.7 ± 0.3 % of the simulated MPB biomass was resuspended. With respect to primary production, 25 % of the MPB primary production simulated during low tides was resuspended, which corresponded in the model to a total annual resuspension of 31.6 g C m −2 .

P. ulvae dynamics
The MPB biomass simulated by the model was also constrained by the grazing pressure from the gastropod P. ulvae. The simulated density and biomass of P. ulvae increased in late winter with a first seasonal peak of ingestion on 22 February (Fig. 8c). A seasonal maximum of simulated density (25 135 ind m −2 ) and biomass (4 g C m −2 ) was reached on 2 May ( Fig. 8a, b). The simulated density and biomass of P. ulvae were compared to 2008 time-coincident observations ( Fig. 8a, b). In February 2008 the simulated density (2616 ± 371 ind m −2 ) was significantly lower than the measured density (5766 ± 2985 ind m −2 ; Mann-Whitney test: p value < 0.05). In July 2008 an average density of 9020 ± 227 ind m −2 was simulated by the model, while a significantly higher average density of 17 191 ± 7084 ind m −2 was measured (Mann-Whitney test: p value < 0.05). In February 2008 the simulated biomass of P. ulvae was 303.8 ± 40 mg C m −2 , which was significantly lower (Mann-Whitney test: p value < 0.05) than the observed biomass (749.5±388 mg C m −2 ). In July 2008 the model underestimated biomass (2157.2 ± 85 mg C m −2 ), whereas the measured biomass was 4469.8 ± 1841.9 mg C m −2 (Mann-Whitney test: p value < 0.05). The P. ulvae gross secondary production simulated by the model was 27 g C m −2 yr −1 .   Overall, the model reasonably captured the seasonal features depicted by the match-up observations.

Contribution of light, temperature, and grazing to the MPB seasonal cycle
In the model, bottom-up (MST and solar irradiance) and topdown (grazing by P. ulvae) processes constrained the simulated MPB growth rate. Light and temperature limitation terms (see Eqs. B6 and B7 in Appendix B2) varied between 0 and 1. At each time step, the lowest value was set as the most limiting term constraining the computation of the MPB photosynthetic rate. Over each daytime exposure period, the most limiting bottom-up factor was defined as the factor whose limitation was the longest. In phase 1, MST and light limited MPB growth 30 % and 70 % of the time, respectively, because PAR and simulated MST values were lower than the light saturation parameter (E k , 100 W m −2 ) and the temperature optimum for photosynthesis (T opt , 18 • C), respectively (Table 4). In phase 2, light was the most limiting factor (60 %, Table 4). The increasing daytime duration allowed MPB to grow in two daytime emersion periods at the beginning and at the end of the daytime period during neap tides (Fig. 9). However, the simulated MPB was exposed to relatively low light levels during dawn and dusk compared to spring tides conditions,  when the emersion periods occurred in the middle of the day and at relatively high light levels ( Fig. 9). With respect to temperature, the MPB growth was more limited by MST in phase 2 (40 %) than in phase 1 (30 %, Table 4). The high summer air temperature and solar irradiance heated the mud surface (Figs. 2c, d and 4), especially when daytime exposure periods occurred in the middle of the day (10:00-16:00) in spring tides (Fig. 9), with simulated MST consequently higher on average than the MPB T opt value (Fig. 10a). In phase 3, the MPB growth rate was almost limited only by downward irradiance (99 %, Table 4). In fall, the average solar irradiance in daytime exposure periods decreased faster Figure 10a shows the daily occurrence of MPB limitation by the simulated MST over 2008. In phase 1, the simulated MST increased towards T opt and, combined with increasing irradiance, led to a seasonal maximum of the mass-specific photosynthetic rate (Fig. 6b). It resulted in a seasonal maximum of MPB biomass in late March (Fig. 10a). In May (phase 1), the mass-specific photosynthetic rate started to decrease due to thermo-inhibition as soon as the MST exceeded T opt (Figs. 6b and 10a). In phase 2, the simulated  MST was always higher than T opt when temperature limitation occurred (Fig. 10a).
With respect to grazing, the simulated biomass grazed by P. ulvae was compared to the simulated MPB biomass produced over the daytime emersion period (Fig. 10b). During phase 1, the ingested MPB biomass exceeded the MPB PP during 11 days (Fig. 10b). The simulated peaks of ingestion rate during these days varied between ∼ 20 and 90 ng Chl a ind −1 h −1 (Fig. 8c), which was consistent with the reported values from laboratory measurements (0.75-385 ng Chl a ind −1 h −1 ; Table 3). The daily averaged P. ulvae ingestion : MPB production ratio was lower but more variable in phase 1 (0.31 ± 0.45) than in phase 2 (0.47 ± 0.18) (Fig. 10b). Phase 1 was characterised by a marked and synoptic impact of grazing at high MPB biomass levels. By contrast, grazing was moderate but more sustained in phase 2. Grazing contributed with thermo-inhibition to maintain relatively low levels of MPB biomass (Fig. 10). As the ingestion rate of P. ulvae was related to the MPB biomass and to the MST, the peaks of grazing simulated in spring resulted from both the high MPB biomass accumulated during the bloom and the MST close to the temperature optimum for grazing by P. ulvae (T opt Z ).
In the model, the occurrence of temperature or light limitation resulted from the coupling of the fortnightly tidal cycle with the seasonal solar irradiance and air temperature cycles. Over 2008, light was the most limiting factor because of low light levels in fall-winter and the occurrence of early and late daytime exposure periods during neap tides in springsummer. During summer spring tides, the exposure periods occurred in the middle of the day and led to high simulated MST values (> 20 • C), hence limiting the MPB growth rate (T opt = 18 • C). Consequently, the high grazing by P. ulvae in spring driven by the high MPB biomass simulated during the bloom was followed by a low MPB PP due to thermoinhibition along with a moderate but sustained grazing by P. ulvae in summer. It resulted in a marked depression of the simulated MPB biomass in summer.

Annual MPB production sensitivity
A total of 10 000 model runs (N) were performed, in which a set of biological constants (T opt , T max , T opt Z , α Z , E k , and K E ) was randomly selected within the reported observed ranges (Table 2). These biological constants were chosen because they were direct inputs in the mathematical functions used in the calculation of the simulated MPB production rate and P. ulvae ingestion rate. The sensitivity analysis resulted in two kinds of model runs according to the sustainability of the MPB PP over the year. Model runs in which PP was sustained (SPP runs, PP > 40 g C m −2 yr −1 , N = 1632) were distinguished from runs characterised by vanishing PP (VPP runs, PP ≤ 40 g C m −2 yr −1 , N = 8368) according to a graphical representation of the annual PP as a function of the number of runs (Fig. 11). Figure 12 shows the 10 000 parameter combinations and the resulting MPB annual PP. The VPP runs are represented by the dark blueish lines (PP < 40 g C m −2 yr −1 ), while the light blueish to reddish colour gradient represents the SPP runs (PP > 40 g C m −2 yr −1 ). In addition to SPP and VPP runs in which all six biological constants varied simultaneously, simulations were run for which only one of the six constants varied at a time (Fig. 13). Figures 12 and 13a and b show that either a T opt value greater than 24 • C or an MPB temperature maximum (T max ) lower than 26 • C induced the reduction of the annual MPB PP. The annual MPB PP was significantly negatively and positively correlated with T opt and T max , respectively (Fig. 13a,  b). In SPP runs, the annual PP was negatively but not significantly correlated with T opt (Spearman's r = −0.04, p value > 0.05; Table 5) because T opt slightly varied within a range  (days) of the simulated temperature limitation term when daily averaged mud surface temperature during emersion periods was lower (grey vertical bars) or higher (black vertical bars) than the optimal temperature for MPB growth (T opt ). (b) Simulated daily primary production rate (mg C m −2 d −1 ) and P. ulvae ingestion rate (mg C m −2 d −1 ). The dashed vertical lines delimit the three phases shown in Fig. 7. Figure 11. Frequency histogram of the annual primary production (g C m −2 yr −1 ) simulated in the Monte Carlo sensitivity analysis.
(18 ± 2.34 • C) corresponding to the T opt threshold shown in Fig. 13a. Moreover, the annual PP simulated in the SPP runs reflected the combined effect of the variation of T opt with the other biological constants (Fig. 12). The annual PP simulated in SPP runs was positively and significantly correlated with T max (Spearman's r = 0.15, p value < 0.05; Table 5). In SPP runs, the correlation between the annual PP and the MPB temperature amplitude (T amp , the difference between T opt and T max ) was even higher than the correlation between PP and T opt and T max (Spearman's r = 0.21, p value < 0.05; Table 5). Figure 12 indeed shows that an increase in T amp (i.e. a decrease in T opt concomitant with an increase in T max ) led to an increase in PP. The positive effect of an increase in T amp on the annual PP is also shown in Fig. 13a and b as either T opt or T max varied while the other was fixed (T opt = 18 • C, T max = 38 • C). The mean values of T amp , T opt , and T max were 15, 18, and 34 • C, respectively, with relatively low variations of T opt and T max (σ norm ≈ 0.13) in SPP runs (Table 5). With respect to temperature, the use of such a set of values promoted PP in the model. In contrast, runs with combinations which included a T opt above 27 • C or a T max below 20 • C resulted in the vanishing of PP over 2008 (dark blueish lines in Fig. 12). In VPP runs, the mean value of T amp was 10.1 • C lower than in SPP runs because the mean T opt value (29 • C) was higher than in SPP runs (18 • C). The maximum value of T opt was 13 • C higher in VPP runs than in SPP runs. The resulting wider range of T opt values led to higher variations in T amp in VPP runs (σ norm = 0.73). However, SPP runs were also characterised by a T amp minimum of 4.5 • C, which was ∼ 3-fold lower than the T amp mean value (15 • C).
PP was negatively correlated with E k in SPP runs (Spearman's r = −0.71, p value < 0.05) and induced large variations of annual MPB PP (Fig. 13c). Runs in which the simulated annual PP was high were characterised by E k values in the lower part (from 2.5 to 100 W m −2 ) of the full tested range (Fig. 12). In SPP runs, the mean value of E k (77 W m −2 ) was lower than in VPP runs (94 W m −2 ). However, E k variations were comparable (0.55 < σ norm < 0.64) and the minimum (2.5 W m −2 ) and maximum values (180 W m −2 ) were the same in both the SPP and VPP runs. Consequently, annual PP is less sensitive to variations of E k than to variations of T opt and T max , and in SPP runs, a low value of E k could sustain PP if T amp was lower than 15 • C. Annual PP was sensitive to the half-saturation constant for light use (K E ) but to a lesser extent, as a high annual PP was simulated using a K E value within the full tested range (Fig. 12). The annual PP in SPP runs was positively correlated with the half-saturation constant for light use (K E ; Spearman's r = 0.2, p value < 0.05; Table 5 and Fig. 13d). Parameter combinations including high K E values (15-20 Ein m −2 d −1 ) resulted in the highest annual PP.
When either T opt Z or the shape parameter of the temperature grazing function (α Z ) varied individually in the model, it induced only small variations of the simulated annual PP (Fig. 13e, f). In SPP runs, PP showed a low but signifi- Table 5. Metrics obtained from the Monte Carlo sensitivity analysis on the simulated annual primary production of MPB.

Sustainable primary production runs
Vanishing primary production runs r is the Spearman's correlation coefficient between annual production values from the different runs with the associated parameter values (the asterisk indicates p value < 0.05).
T amp corresponds to the difference between T max and T opt . σ norm is the normalised standard deviation, i.e. the standard deviation divided by the mean. Figure 12. Parallel coordinates of the MPB annual primary production (g C m −2 yr −1 ) according to the temperature optimum for MPB growth (T opt ), the temperature maximum for MPB growth (T max ), the light saturation parameter (E k ), the half-saturation constant for light use (K E ), the temperature optimum for grazing by P. ulvae (T opt Z ), and the shape parameter of the temperature grazing function (α Z ) for 10 000 combinations tested in the Monte Carlo sensitivity analysis.
cant correlation with T opt Z (Spearman's r = 0.17, p value < 0.05), suggesting that high T opt Z values resulted in high levels of annual PP. PP was not correlated with α Z in SPP runs (Spearman's r = −0.03, p value > 0.05). However, T opt Z and α Z variations were high and of the same extent in both the SPP and VPP runs (σ norm = 0.21 and σ norm ≈ 0.57, respectively). The mean, maximum, and minimum values of T opt Z and α Z were also very similar in both SPP and VPP runs (Table 5). Compared to other parameters, annual PP was less sensitive to the P. ulvae grazing parameters as SPP runs took place in all the regions of the tested ranges (Fig. 12).
Overall, the simulated annual PP was most sensitive to the MPB light-and temperature-related constants. The specific set of biological constants used in the study promoted realistic levels of MPB primary production. A specific set of these temperature-and light-related parameters allowed for a sustainable level of MPB production and biomass, which resulted in a significant effect of grazing on MPB annual production.

The MPB seasonal cycle
Our study suggests an MPB seasonal cycle on the Brouage mudflat characterised by three phases in 2008, i.e. a bloom in winter-spring, low biomass levels in summer, and a peak of moderate intensity in fall. Cariou-Le Gall and  sampled the MPB Chl a concentration within the top 0.5 cm of sediment on the Brouage mudflat monthly from March 1992 to February 1993. Their measurements suggest a bloom in winter-spring and low Chl a concentrations in summer, which is consistent with the 2008 NDVI data, the observed MPB biomass (2008,(2012)(2013), and MPB biomass simulated by the model. Cariou-Le Gall and Blanchard (1995) did not report any peak of MPB biomass in fall, which may be modulated by the interannual variability driven by meteorological conditions. In northern (De Jong and Sahan et al., 2007) and southern (Brito et al., 2013) European mudflats, MPB spring blooms are also observed. However, the contribution of underlying abiotic (e.g. air temperature, irradiance, rain, wind) and biotic (e.g. autotrophic species community, predators) factors is likely to be different in shaping the seasonal MPB cycle at such contrasting latitudes.
In the Brouage mudflat, the simulated seasonal cycle of MPB at the sediment surface compares to that depicted by the remotely sensed NDVI data and measurements made in 2008 and 2012-2013. The simulated MPB biomass in the biofilm and its instantaneous PP are close to maximum values of biomass previously measured in biofilms developing at the surface of very fine sediments of the Brouage mudflat (Herlory et al., 2004). Once at the surface, the simulated MPB growth is regulated by the mass-specific photosynthetic rate in µg C (µg Chl a) −1 h −1 converted into a growth rate (h −1 ) using a variable C : Chl a ratio. The resulting MPB growth rates simulated by the model were consistent with observations made on epipelic diatoms (Gould and Gallagher, 1990;Underwood and Smith, 1998;Scholz and Liebezeit, 2012). With respect to the simulated C : Chl a ratio, it varies within the range of observed values in mudflats (Guarini, 1998;Gould and Gallagher, 1990;de Jonge et al., 2012).
Contrary to Chl a measurements, there were no PP measurements made in 2008 on the Brouage mudflat. For comparison, we use averages of mass-specific photosynthetic rates computed from previous measurements at different locations on the Brouage mudflat for different years (using CO 2 flux data measured in benthic chambers). Despite the year-to-year variability, the mean mass-specific photosynthetic rates simulated by the model during spring tides (0.66 ± 0.04 mg C (mg Chl a) −1 h −1 in April, 0.52 ± 0.03 mg C (mg Chl a) −1 h −1 in May, and 0.44 ± 0.04 mg C (mg Chl a) −1 h −1 in July) were in the range of measurements for the same months (1.6 ± 1.1 mg C (mg Chl a) −1 h −1 in April 2012, 0.28 ± 0.11 mg C (mg Chl a) −1 h −1 in May 2015, and 0.32 ± 0.13 mg C (mg Chl a) −1 h −1 in July 2015; Johann Lavaud, personal communication, 2018). Moreover, simulated daily and yearly PP rates compared to measurements made across other European intertidal mudflats (Underwood and Kromkamp, 1999). The model-data comparison suggests that the model can resolve with confidence the main patterns of the MPB seasonal cycle.
The relative contribution of light, MST, and grazing to the simulated MPB seasonal cycle resulted from the coupling of the fortnightly tidal cycle and seasonal solar irradiance and air temperature cycles. Such a coupling is reported in intertidal sediments in the Tagus estuary, Portugal (Serodio and Catarino, 1999). In the model, an emersion period takes place in the middle of the day during spring tides, exposing the mud surface to a daily solar irradiance and temperature maximum. In summer, when the seasonal maximum of daily solar irradiance and temperature is reached, the high simulated MST values translate into an enhanced thermoinhibition of MPB growth and P. ulvae grazing pressure. The highest MPB thermo-inhibition in summer spring tides was also highlighted by Guarini et al. (1997) in the Brouage mudflat. During neap tides, light limits MPB growth when exposure periods occur early in the morning and late in the afternoon at low daily light levels. The reduced PP of MPB at low light levels and MST values during neap tides compared to spring tides was also observed by Kwon et al. (2014) on the Hwaseong mudflat, South Korea.
In the model, we do not consider any MPB limitation by inorganic nutrients. In the Brouage mudflat, Feuillet-Girard et al. (1997) highlighted the greater affinity of MPB for ammonium compared to nitrate. They suggest a higher availability of ammonium released from the sediment in summer, making it unlikely that nutrient limitation is responsible for the summer depression of MPB biomass. The high nutrient availability in the sediment in summer can be attributed to faunal activities (bioturbation, bio-irrigation, excretion; Feuillet-Girard et al., 1997;Heilskov et al., 2006;Laverock et al., 2011).
The short-term daily dynamics of MPB are also regulated by resuspension events . The intensity of resuspension of MPB into the water column can be either chronic or catastrophic according to the flow velocity and the sediment stabilisation (Mariotti and Fagherazzi, 2012). Catastrophic events can locally resuspend all the MPB biomass as the resuspended sediment layer is thicker than the vertical distribution of MPB biomass (Mariotti and Fagherazzi, 2012). The repeated occurrence of such events over sev-eral days could contribute to shaping the seasonal cycle of MPB by lowering the biomass of photosynthetically competent MPB. In their model, Guarini et al. (2008) introduced a chronic resuspension of all the MPB biomass remaining in the biofilm when tidal floods occurred. In their parameterisation, the MPB biomass remains at the sediment surface according to the mean time spent at the surface (equivalent to τ s in our study). In our study, the chronic resuspension of MPB biomass is formulated by a linear loss term of the MPB biomass within the first centimetre (0.002 h −1 ). In the absence of MPB biomass deposition, the total simulated MPB biomass that is resuspended into the water column represents 25 % of the simulated benthic MPB annual production. Such a value suggests that the benthic MPB production contributes significantly to the pelagic food web (Perissinotto et al., 2003;Krumme et al., 2008). In light of the work of Mariotti and Fagherazzi (2012), resuspension and deposition are key mechanisms that need to be related to fauna bioturbation, sediment characteristics (e.g. nature and stabilisation), and hydrodynamics (Mariotti and Fagherazzi, 2012). Such an approach requires the availability of wave and current data to estimate the bed shear stress and modulate the intensity of resuspension (from chronic to catastrophic events), which are not available at our study site for 2008.

Role of mud surface temperature on the MPB and P. ulvae activity
On the Brouage mudflat, the simulated MST plays a major role in the MPB seasonal cycle. In spring, the simulated MST increases towards the MPB temperature optimum for photosynthesis. Along with increasing light levels, it contributes to increasing the mass-specific photosynthetic rate and triggers the onset of the MPB spring bloom. As soon as the simulated MST exceeds the MPB temperature optimum for photosynthesis, the MPB PP starts to decrease due to thermo-inhibition, particularly during spring tides. In fall, the average solar irradiance during daytime exposure periods decreased faster than the simulated MST. The simulated MST departs slower from the temperature optimum for photosynthesis than the downward irradiance from the light saturation parameter. Despite decreasing solar irradiance in fall, the simulated MPB PP increases until November, when the simulated MPB growth rate is limited by low light levels and MST values with respect to the MPB light saturation parameter (100 W m −2 ) and temperature optimum for photosynthesis (18 • C), respectively. Using the production-temperature (P -T ) model from Blanchard et al. (1996), Blanchard et al. (1997) and Guarini et al. (2006) also suggested that the MPB PP was temperature-limited in summer on the Brouage mudflat. On a southern intertidal mudflat (Tagus Estuary, Portugal), Brito et al. (2013) suggested that thermo-inhibition was responsible for the summer MPB depression observed in NDVI time series in conditions of high sediment temperature (30 • C). In addition, the detrimental effect of MST ranging between 18 and 24 • C was shown in microcosms using fluorescence (Cartaxana et al., 2015).
In the model, production is related to temperature according the P -T relationship of Blanchard et al. (1996). As a result, the occurrence and intensity of MPB thermo-inhibition depends on the MPB temperature optimum and maximum for photosynthesis used in the relationship. The set of parameters determines the thermal threshold and interval at which thermo-inhibition occurs. The sensitivity analysis shows that annual PP is very sensitive to the temperature amplitude between the two parameters. The annual PP increases as the amplitude increases. On the Brouage mudflat, the MPB temperature optimum and maximum for photosynthesis were estimated to be 25 and 38 • C, respectively, and assumed to be constant over the year . In our study, a lower MPB temperature optimum for a photosynthesis value of 18 • C is required to simulate a spring bloom that compares to the NDVI time series. Such a temperature optimum also implies a more rapid onset and a higher MPB thermoinhibition as the simulated MST increases in summer. Values of both the MPB temperature optimum and maximum for photosynthesis are reported to vary by up to 10 • C ( Table 6). In that respect, the MPB temperature optimum for photosynthesis is a key parameter in the model because it constrains the onset of the MPB spring bloom and the thermo-inhibition span and intensity.
In addition, the strong heating and wind exposure of the mud surface is accompanied by pore water evaporation that results in desiccation and increased salinity (Coelho et al., 2009). A decrease in pore water content can induce even more detrimental effects within the cells through the production of reactive oxygen species (Rijstenbil, 2003;Roncarati et al., 2008), potentially leading to the oxidation of the photosynthetic unit (Nishiyama et al., 2006). The motility of epipelic diatoms is supposed to be a strategy to avoid harmful conditions at the surface of cohesive sediments (Admiraal, 1984). However, Juneau et al. (2015) showed no significant negative effect of salt stress on the photosynthesis of immobile epipelic diatoms. Coelho et al. (2009) highlighted the role of the rate of pore water content decrease in the field. While slow desiccation (reduction by 40 % of the pore water content in 4.5 h) had no significant negative effect on the photosynthesis of microphytobenthic cells within the biofilm, fast desiccation (reduction by 40 % of the pore water con-tent in 2 h) resulted in desiccation and decreased the photosynthetic activity of MPB (Coelho et al., 2009). In addition to micro-migrations, epipelic diatoms produce extracellular polymeric substances (EPSs) to temper the effect of desiccation and high salinity (Steele et al., 2014). High sediment temperature (> 35 • C) is also known to reduce the motility of MPB diatoms and thus their capacity to avoid harmful conditions at the sediment surface (Cohn et al., 2003;Laviale et al., 2015). The detrimental effects of high salinity levels are not explicitly accounted for in the model. The underlying processes could be accounted for in the model in an implicit way by adjusting the MPB temperature-related growth parameters to accentuate the PP reduction in simulated conditions when high evaporation is associated with high MST. The detrimental effects of desiccation on MPB cell motility could also be implicitly represented in the model through more photo-inhibition.
The simulated MST also governs the ingestion rate of MPB by the grazer P. ulvae in the model. The simulated PP rates increase as the value of the optimal temperature for grazing increases because the grazing optimum is not often reached in the model. In the model, the ingestion rate increases when the MST tends towards the optimal temperature for grazing (fixed at 20 • C; Pascual and Drake, 2008). A high metabolism of benthic grazers promoted by high temperature conditions (up to 22 • C) and the resulting increase in the grazing pressure on benthic diatoms was observed by Sahan et al. (2007) on a mudflat in the Netherlands.

Effect of light on MPB photosynthesis
In the model, light is the most limiting factor throughout the year. The low irradiance during fall and winter limits the MPB photosynthesis as irradiance is on average lower than the light saturation parameter. In spring, the increasing irradiance and MST translate into higher mass-specific photosynthetic rates than in fall-winter, leading to the onset of the simulated MPB spring bloom. In summer, photo-inhibition is not accounted for in the model as the simulated mean time spent by an MPB cell at the surface is lower than the time required to induce photo-inhibition at saturating light levels . As a consequence, light limits the simulated MPB growth only during neap tides, when the sediment exposure occurs at low light levels early and late in the day.
Photosynthesis is represented in the model by the production-irradiance (P -E) model of Platt and Jassby (1976). It relies on photosynthetic capacity (P b max ), the light saturation parameter (E k ), and the maximum light utilisation Talling, 1957). Irradiance has no influence on the photosynthetic capacity and maximum light utilisation coefficient (MacIntyre et al., 2002) in our study. Based on the work of Blanchard et al. (1996), the photosynthetic capacity and maximum light utilisation coefficient vary in the model with the simulated MST. Therefore, the seasonal adjustment of photosynthesis to irradiance depends mainly on the photoacclimation status of MPB cells, which can be related to the light saturation parameter (Sakshaug et al., 1997). The light saturation parameter corresponds to the irradiance at which photosynthesis switches from light reactions (light absorption and photochemical energy conversion) to dark reactions (reductant utilisation) (Sakshaug et al., 1997). It has been reported to vary seasonally in benthic microalgae (Blanchard and Cariou-Le Gall, 1994;Barranguet et al., 1998;Light and Beardall, 2001;Pniewski et al., 2015;Barnett et al., 2015). Cells increase their light saturation parameter at high irradiance (summer) and reduce it with decreasing light levels (Sakshaug et al., 1997). In our study, as the light saturation parameter is set as constant throughout the year (100 W m −2 ), photoacclimation is simulated by a variable C : Chl a ratio.
During winter, low-light-acclimated cells have a lower C : Chl a ratio due to an increase in Chl a content (MacIntyre et al., 2002;Brunet et al., 2011). In summer, with the increasing irradiance and day length, high-light-acclimated cells reduce their Chl a content, leading to a higher C : Chl a ratio (MacIntyre et al., 2002;Brunet et al., 2011). In the model, solar irradiance shapes the simulated C : Chl a ratio (Eq. B8 in Appendix B2). The C : Chl a ratio reaches a seasonal maximum value (75.5 g C g Chl a −1 ) in summer when solar irradiance is the highest. Such a result is consistent with estimates (80 g C g Chl a −1 ) reported in summer by de Jonge et al. (2012). In the model, given that the mass-specific photosynthetic rate (µg C (µg Chl a) −1 h −1 ) and the C : Chl a ratio are related to the growth rate (h −1 ), the growth rate increases as the C : Chl a ratio decreases (low-light-acclimated cells). The seasonal variation of the simulated growth rate results from a combination of the variation of the photosynthetic capacity and maximum light utilisation coefficient driven by the simulated MST and the variation of the C : Chl a ratio with irradiance.
Finally, photo-inhibition at high irradiance is not accounted for in the P -E model of Platt and Jassby (1976) used in the model. Epipelic diatoms achieve "micro-migrations" within the sediment to avoid harmful light conditions prevailing at the sediment surface Perkins et al., 2001;Cartaxana et al., 2011). However, combined with high temperature conditions (> 35 • C) at the sediment surface potentially leading to reduced cell motility (Cohn et al., 2003), epipelic diatoms can be photo-inhibited (Laviale et al., 2015). In temperate intertidal mudflats, high light and temperature conditions occur during summer and their combined effect on MPB photosynthetic rate may explain the depression of the MPB biomass observed in summer.

Top-down regulation of MPB dynamics
Grazing by meiobenthos and macrobenthos is often suggested as the main driver of the MPB biomass depression 7258 R. Savelli et al.: On biotic and abiotic drivers of the microphytobenthos seasonal cycle observed in summer on intertidal mudflats (Cadée and Hegeman, 1974;Cariou-Le Gall and Blanchard, 1995;Sahan et al., 2007;Orvain et al., 2014). Weerman et al. (2011) experimentally showed a strong decrease in MPB biomass in the presence of macrofauna driven by direct grazing and by the absence of surface mud stabilisation due to bioturbation by deposit feeders.
In the model, P. ulvae grazing exceeds MPB PP mainly in spring (11 days of MPB biomass removal). P. ulvae depletes a substantial part of the MPB biomass accumulated during the spring bloom. After the bloom, a moderate but sustained grazing by P. ulvae adds to the effect of thermo-inhibition on the MPB dynamics. The simulated gain terms promoting the growth rate of MPB limited by thermo-inhibition do not compensate for the loss terms dominated by the grazing pressure, which leads to a decrease in the MPB biomass. In a conceptual model, Thompson et al. (2000) showed such a seasonal uncoupling between the grazing intensity by intertidal grazing molluscs and the microalgae abundance from observations made on a rocky shore of the Isle of Man (United Kingdom). The authors conceptualised the role played by the light and temperature stress on the microalgae productivity and by the temperature-promoted grazing in the depression of the microalgal standing stocks in summer.
The simulated annual P. ulvae gross secondary production is 27 g C m −2 yr −1 , which represents 21 % of the simulated annual MPB PP (127 g C m −2 yr −1 ). This fraction of PP transferred to P. ulvae secondary production is consistent with the average fraction reported by Asmus and Asmus (15 ± 12 %;1985) on intertidal sand bottom communities of the island of Sylt in the North Sea. In July, the simulated density of P. ulvae lies in the lower range of time-coincident measurements. As the simulated MST fairly agrees with time-coincident measurements, other factors may explain the likely underestimation by the model of the density and ingestion of P. ulvae. First, there may be a bias resulting from the monthly averaged weight estimates used to simulate the P. ulvae density (see Appendix B3). The monthly averaged weights are based on samples gathered in 2014-2015 on the Aiguillon mudflat, in the vicinity of the Brouage mudflat (Fig. 1). Nevertheless, the seasonality of the P. ulvae density is similar on the two mudflats with a peak density in late summer (Haubois et al., 2002), which suggests that such a bias is likely limited. In addition, the simulated ingestion rates (20-90 ng Chl a ind −1 h −1 ) are consistent with ingestion rates measured in experiments with P. ulvae and benthic diatoms collected in our study area and performed at a temperature close to the optimal temperature for grazing in the model (15-20 • C; 0.75-52 ng Chl a ind −1 h −1 ; Blanchard et al., 2000;Haubois et al., 2005;Pascal et al., 2008). Second, the P. ulvae density on the mudflat can change horizontally as a result of the foraging activity of the individuals and transport mediated by the wave-and tidal-induced shear stress on the bottom sediment. Such a process is not accounted for in the model and may lead to an underestima-tion of the P. ulvae biomass and density. Finally, potential MPB grazing by fauna other than P. ulvae is represented in a simple way by a linear and generic loss term in the model, whereas it might be a non-linear process that can vary seasonally (Pinckney et al., 2003). This closure term may be underestimated in the model.
With respect to meiofauna, Pinckney et al. (2003) suggested a more intense grazing by meiofauna in summer than in winter in the Terrebonne Bay estuary (USA). Admiraal et al. (1983) estimated the meiofauna grazing at 300 mg C m −2 d −1 on a mudflat of the Ems Dollard estuary (Netherlands). Comparable rates of meiofauna ingestion (58-189 mg C m −2 d −1 ) are reported for the Brouage mudflat (Montagna et al., 1995). Admiraal et al. (1983) observed a non-significant effect of meiofauna grazing relative to MPB production rates. Nevertheless, their estimated grazing rate exceeds our simulated daily MPB production rates for 36 days in summer, i.e. 34 % of the time in the second phase in the model, suggesting that meiofauna grazing could impact MPB. In addition, Pascal et al. (2008) compared ingestion rates by P. ulvae and a nematode community from the Brouage mudflat in experimental conditions. According to the abundance of organisms selected for the experiment of Pascal et al. (2008) and a constant C : Chl a ratio of 45 g C g Chl a −1 (Guarini, 1998), the amount of Chl a ingested by nematodes per hour was only 1.5 % of the Chl a ingested by P. ulvae per hour in their experiment. However, in regard to the observed abundances on the field and without a density-dependant effect on grazing rates, this theoretical amount of Chl a ingested by nematodes increases to almost 50 % of the Chl a ingested by P. ulvae in the study of Pascal et al. (2008). According to the measured biomass uptake by meiofauna (Montagna et al., 1995) and nematodes (Pascal et al., 2008) for the Brouage mudflat, an explicit representation of meiofauna ingestion in the model might magnify the simulated depletion of MPB biomass in summer months. The representation of grazing in the model can be improved. Nevertheless, the fair agreement between the simulated P. ulvae densities and biomass levels with time-limited but timecoincident observations suggests that overall the model simulates with some confidence the grazing pressure on MPB.

Physical setting of the coupled model
The predictive ability of the physical-biological coupled model depends on the accuracy of the oceanic and meteorological forcings. The frequency of the water height and meteorological time series used to constrain the model is hourly, while the model time step 6 min. The lower frequency of the model forcings over a day partly explains the model-data discrepancies. In addition, the weather station where meteorological data were acquired is located 30 km away from the Brouage study site. Local weather conditions may differ between the two sites, especially the global irradiance and wind speed used to simulate the MST and MPB growth rate.
Global irradiance can be impacted by local cloud cover and the wind regime can be different due to local thermal winds. In the model, the timing of the emersion-immersion cycle is constrained by the observed water heights and bathymetric level. The bathymetric level used to compute the water height above the Brouage study site originates from a digital elevation model with a 1 m horizontal resolution and a 15 cm vertical precision. Even if the Brouage mudflat is relatively flat (1 : 1000), ridges and runnels are present near the study site  and the topography is highly variable at a metre scale. Inaccuracies in the bathymetric level relative to the study site may translate into model-data discrepancies in terms of the timing of the emersion-immersion cycle in the model. Given that the mud temperature model is constrained by the water height and meteorological data, it is sensitive to possible inaccuracies in the forcings that may impact the simulated hourly dynamics of MPB and P. ulvae. Nevertheless, at the seasonal scale, the impact on the biological compartments of such inaccuracies in the forcings may be limited.

Conclusion and perspectives
This study is a first attempt to simulate the MPB seasonal cycle observed on a temperate intertidal mudflat and to quantify the relative contribution of both biotic and abiotic factors to seasonal MPB dynamics. The physical-biological coupled model fairly compares to time-coincident remotely sensed and in situ data and provides key findings about the seasonality of MPB on the Brouage mudflat (French Atlantic coast).
-The 2008 MPB seasonal cycle consists of three phases: a spring bloom, a summer depression of the biomass levels, and a moderate peak of biomass in fall.
-In winter and early spring, the seasonal mass-specific maximum photosynthetic rate mainly driven by the simulated MST and the seasonal low C : Chl a ratios lead to a seasonal maximum of MPB growth rate and to an MPB spring bloom.
-P. ulvae grazing has a marked and synoptic impact on the MPB biomass accumulated during the spring bloom.
-In late spring-summer, grazing is moderate but more sustained. Both grazing and thermo-inhibition, which is limiting for MPB growth 40 % of the time in summer, contribute to maintain relatively low levels of MPB biomass.
-The model is sensitive to MPB temperature parameters (temperature optimum and maximum for photosynthesis), to the MPB light saturation parameter, and to a lesser extent to grazing parameters (the optimal temperature for grazing and the shape parameter of the temperature-related grazing function).
The seasonal MPB dynamics simulated by the model compare to time-coincident times series of remotely sensed NDVI data, hence providing a qualitative assessment of the model predictive ability. A next step would be to extend such a model-satellite data comparison to a more quantitative assessment to validate the simulated levels of MPB Chl a concentration and PP. The recent advance of multispectral and hyperspectral remote sensing allows for the development of new algorithms to retrieve products of ecological interest for MPB. Brito et al. (2013) developed local empirical relationships by relating synchronised NDVI data to in situ Chl a concentrations to retrieve from space estimates of Chl a concentration on a Portuguese intertidal mudflat. Efforts are also focused on using remote sensing reflectance from airborne hyperspectral data to assess MPB PP rates (Méléder et al., 2018). Recently, and in light of the work of Brito et al. (2013), Daggers et al. (2018 combined biomass derived from NDVI data with simulated photosynthetic capacity from environmental conditions (irradiance and air temperature) to map MPB PP on intertidal mudflats in the Netherlands. Other promising methods in the estimation of PP on intertidal mudflats at the ecosystem scale are non-invasive atmospheric and aquatic eddy covariance (EC) techniques. Atmospheric EC provides continuous and direct CO 2 flux measurements at the air-water and air-sediment interfaces during high and low tides, respectively, across different timescales from hours to years (Baldocchi et al., 1988;Aubinet et al., 1999;Zemmelink et al., 2009;Polsenaere et al., 2012). Similarly, aquatic EC measures benthic O 2 fluxes at the sedimentwater interface (Berg et al., 2003). Quantifying the MPB PP and biomass on intertidal mudflats is a prerequisite for further estimating the flux of biogenic carbon from the benthos to the pelagos. During the immersion period, MPB can be resuspended (9.7 mg C per high tide, i.e 3 % of the mean simulated production during low tides; Dupuy et al., 2014) and highly disturb the functioning of the benthic-pelagic ecosystem (Saint-Béat et al., 2014). The study of air-water and sediment-water exchanges through simultaneous atmospheric and aquatic EC measurements could allow for the quantification of metabolic fluxes during immersion and exposure periods but also the coupled processes between the benthic and pelagic compartments, such as MPB resuspension. Microphytobenthic community resuspension can significantly contribute to planktonic gross PP and, in turn, explain lower CO 2 fluxes from the water column to the atmosphere at high tide during the day than at night (Guarini et al., 2008;Polsenaere et al., 2012). To date, the modelling effort put toward the physically driven (tides and waves) resuspension processes of MPB is still limited (see Mariotti and Fagherazzi, 2012). Accounting in models for sediment bottom shear stress mediated by hydrological forcings (current and waves) along with bioturbation processes could lead to more realistic predictions of interannual MPB dynamics. Such a representation of biologically and physically driven benthic-pelagic interactions would be fully captured by the coupling of biological MPB models to high-resolution ocean models. Such an approach would open the door to an accurate assessment of the vertical and horizontal export flux of biogenic matter at the land-ocean interface and, more generally, of the contribution of productive biofilms in mudflats to the carbon cycle of the global coastal ocean. The original version of the mud temperature model of Guarini et al. (1997) is simplified by only resolving the mud surface temperature T M (z 0 , t) (K), which is governed by the following equation during emersion periods: where f (T M (z 0 , t)) is the heat energy balance (HEB, W m −2 ) at the sediment surface z 0 (m) at time t (s). This sediment surface layer is 1 cm deep. The temperature (K) is assumed to be homogeneous within the layer and is governed by the HEB (Harrison and Phizacklea, 1987;Piccolo et al., 1993). ρ M is the volumetric mass of mud (kg m −3 ). It is the sum of the water fraction and of the dry sediment fraction (ρ M = ρ W ξ + ρ S (1 − ξ ) where ρ W and ξ are the water volumetric mass (kg m −3 ) and the porosity (%), respectively. C P M is the specific heat capacity of mud at constant pressure (J kg −1 K −1 ): where η is the heat conductivity (W m −1 K −1 ) and µ the thermal diffusivity (m 2 s −1 ). Heat exchange fluxes at the sediment interface are different according to the emersionimmersion cycle. During low tide, the HEB is governed by downward fluxes of radiation from the sun (R S , W m −2 ) and from the atmosphere (R Atm , W m −2 ), by upward fluxes of radiation from the receiving surface (R M , W m −2 ), by sensible heat fluxes, by conduction due to mud-air temperature differences (S Mud→Air , W m −2 ), and by flux of evaporation (V M , W m −2 ): where ξ is the mud porosity (ξ ∈ [0, 1], %) and V W the evaporative heat flux of seawater (W m −2 ). Details about formulas and constant computation of each flux during emersion are given in Tables A1 and A2. During immersion, Guarini et al. (1997) and Harrison and Phizacklea (1987) suggested a rapid equilibrium between mud surface temperature and the temperature of the overlying water layer. The simulated mud surface temperature is therefore set to water temperature during immersion periods: The simulated seawater temperature of the whole water column (T W ) results from the mixing of the surface layer (z top ) influenced by the atmospheric forcings (i.e. equivalent to the mixed layer depth) and the bottom layer (z bot ), which remains at the seawater temperature computed at the previous time step of the model run.
The seawater temperature in the top layer of the water column is governed by the HEB at the water surface: where ρ W is the volumetric mass of water (kg m −3 ). C P W is the specific heat capacity of seawater at constant pressure (J kg −1 K −1 ). T W (z top , t) is the water temperature (K) in the surface mixed layer. The term S Air→Water is the sensible heat flux (W m −2 ) mediated by thermal conduction due to waterair temperature differences. R W (W m −2 ) is the seawater upward radiation.
The upper fraction of the water column influenced by atmospheric forcings is defined by the coefficient α top : where U is the wind speed (m s −1 ). Consequently, the simulated seawater temperature of the whole water column (T W ) results from the mixing between the fraction α top and the remaining fraction of the water column (1 − α top ): T W (K) is initialised by the following equation: T W (t) = 18.5 + 5 cos 2π day − 230 year length + 273.15, where "day" is the day of the year and the year length is in days. Details on parameters and constants are given in Tables A1 and A2.
air volumetric mass C P A : specific heat of air at constant pressure C h M→A : bulk transfer coefficient for conduction between mud and air U : wind speed measured at 10 m S Air→Water = ρ A C P A C h A→W (1 + U ) (T W (t) − T A ) C h A→W : bulk transfer coefficient for conduction between air and water L V : latent heat evaporation C V : bulk transfer coefficient for evaporation q S : specific humidity of saturated air at water temperature q A = q S H r q A : absolute air humidity L V = [2500.84 − 2.35 (T E − 273.15)] × 10 3 H r : relative air humidity formulated by Van Bavel and Hillel (1976) T E : temperature of interstitial water (in equilibrium with mud temperature)  A system of three partial differential equations describes the temporal dynamics of the MPB biomass within the surface biofilm (S), MPB biomass within the first centimetre of sediment (F ), and biomass of MPB grazer P. ulvae (Z). The system drives the MPB migration scheme according to the diurnal and tidal cycles that constrain the biological-physical coupled model (Table 1). During the daytime emersion period, where τ (h) corresponds to the potential duration of the biofilm or the potential duration of the production period. It is computed at the end of each night-time emersion and immersion period for the next daytime emersion period (Eqs. B4 and B5). When τ > 0, the MPB cells migrate upward in the sediment from F to S compartment at a transfer rate of r F (h −1 ). MPB stop migration when S reaches saturation at S max (mg Chl a m −2 ). Primary production within the S compartment regulated by the mass-specific photosynthetic rate P b (µg C (µg Chl a) −1 h −1 ) is set to zero when S = S max according to the term 1 − S S max , which represents the MPB space limitation in the S compartment. The MPB biomass produced is therefore transferred from S to F according to the term P b S S S max in the F time derivative. In order to take into account all the MPB biofilm biomass plus the biomass produced in the biofilm (S * ), the S * time derivative was computed as follows: When τ ≤ 0, the MPB cells migrate downward in the sediment from the S to F compartment at a transfer rate of r S (h −1 ). The terms m S and m F are loss rates (h −1 ) representing MPB senescence and grazing by surface deposit feeders (on S) and subsurface deposit feeders (on F ). m Z is a loss rate (h −1 ) representing P. ulvae senescence (see Appendix B3).
During the night exposure period, the MPB cells migrate downward into the sediment from S to F . P. ulvae grazes on MPB cells remaining in the biofilm (S): According to Guarini et al. (2006Guarini et al. ( , 2008, τ depends on the MPB biomass in the F compartment relative to S max and on the average time spent at the surface by a unit of biomass equal to S max (τ s , 1 h). It that suggests the higher the biomass in F , the longer S will remain at saturation S max .
During the immersion period, MPB cells remaining in the biofilm finish their downward migration and P. ulvae no longer exerts any grazing pressure: In the model, we assumed a constant rate of MPB cells resuspended during immersion periods . During immersion periods, the generic loss term (ν F , 0.003 h −1 ) includes the chronic resuspension, MPB senescence processes, and the grazing by subsurface deposit feeders. During emersion periods, the loss term is lower (m F , 0.001 h −1 ) as it only represents MPB senescence and grazing by subsurface deposit feeders. Parameter values are given in Table B1.

B2 MPB primary production
The mass-specific photosynthetic rate P b (µg C (µg Chl a) −1 h −1 ) is regulated by temperature (T , • C) and by photosynthetically active radiation (E, W m −2 ), which corresponds to 44 % of downward shortwave radiation (Britton and Dodd, 1976). The model of Platt and Jassby (1976) is used to determine the production rate as a function of light: where P b max is the photosynthetic capacity (µg C (µg Chl a) −1 h −1 ) and E k is the light saturation parameter (W m −2 ). P b max depends on the mud surface temperature T according to the relationship of Blanchard et al. (1996): Mean bulk density of sediment 520 g l −1 Present study where T max ( • C) and T opt ( • C) are the maximum and optimal temperature for photosynthesis, respectively. β is a curvature coefficient that shapes the temperature-photosynthesis relationship. P b MAX is the maximum value that takes P b max at T opt . The mass-specific photosynthetic rate P b is expressed in µg C (µg Chl a) −1 h −1 . We multiplied it to a variable Chl a : C ratio (µg Chl a µg C −1 ) on the finding of de Jonge et al. (2012) to obtain a gross growth rate in h −1 . The ratio is computed according to the formulation of Cloern et al. (1995) adapted for coastal pelagic diatoms (Sibert et al., 2010(Sibert et al., , 2011Le Fouest et al., 2013): where Chl a C min is the minimum Chl a : C ratio (g Chl a g C −1 ) and K E is the half-saturation constant for light use (in Ein m −2 d −1 ). The MPB primary production (µg C m −2 h −1 ) corresponds to the sum of the space-dependant production at the sur-face of the biofilm (i.e. the P b S 1 − S S max term) and of the biomass produced and directly transferred from S to F (i.e. the P b S S S max term). Consequently, it can be simplified by The constant values are given in Table B1.
B3 Grazer P. ulvae S is explicitly grazed by the mud snail Peringia ulvae (Z, mg C m −2 ). The grazing rate is regulated by the individual ingestion rate of snails (IR; ng Chl a ind −1 h −1 ) and Z expressed in terms of density (ind m −2 ). Density is computed as the ratio of Z (mg C m −2 ) over the mean individual weight (W mean Z , mg C) linearly interpolated on a simulation timescale (6 min; Table B2). The grazing is limited through a heaviside function (H ) including a feeding threshold (S mini , mg Chl a m −2 ). Only a fraction (γ , %) of the MPB biomass grazed by Z is assimilated into new Z biomass.
The individual ingestion rate (ng Chl a ind −1 h −1 ) by P. ulvae is calculated using a sigmoid mathematical function accounting for the effect of mud temperature T ( • C): where T opt Z ( • C) is the optimal temperature for grazing. IR max is the maximal observed individual ingestion rate. α Z (no unit) is a curvature parameter. The maximal individual ingestion rate IR max (ng Chl a ind −1 h −1 ) is calculated according to the formulation of Haubois et al. (2005) for adult snails. IR max depends on the total MPB biomass: IR max = 0.015 × (F + S) 1.72 .
The Chl a uptake rate is converted into carbon units according to the C : Chl a ratio described previously. The term (F + S) is expressed in µg Chl a g dry sed −1 . The biomass expressed in mg Chl a m −2 is converted into µg Chl a g dry sed −1 as follows: [Chl a](µg Chl a g dry sed −1 ) = [Chl a] 1.2605 (mg Chl a m −2 ) ρ S × thickness sed , where ρ S is the sediment bulk density in units of g l −1 and thickness sed is the sediment thickness, i.e. 1 cm. The Chl a concentration is scaled by the exponent 1.2605 in order to reach a maximal observed ingestion rate of 385 ng Chl a ind −1 h −1 (Coelho et al., 2011) when the Chl a concentration converges towards a maximal observed value (300 mg Chl a m −2 ; Guarini, 1998). Finally, the mortality rate of Z is a quadratic densitydependant mortality rate: where m min Z is the minimum mortality rate (h −1 ). The constant values are given in Table B1.
Author contributions. VLF, RS, PP and CD conceived and designed the study. The biological model was developed by RS and VLF. The mud temperature model was developed by KG. RS analysed the simulated data. AL and LB analysed the remotely sensed images. PB, AP and CD collected field data. All authors contributed to writing, reviewing and editing of the paper.
Competing interests. The authors declare that they have no conflict of interest.