Mapping the Intertidal Microphytobenthos Gross Primary Production Part I: Coupling Multispectral Remote Sensing and Physical Modeling

Introduction

Tidal mudflats are soft sediment coastal marine ecosystems that undergo regular tidal inundation. Their bare surface is colonized by biofilms constituted of photosynthetic microorganisms (microalgae and cyanobacteria), collectively known as microphytobenthos (MPB) (e.g., MacIntyre et al., 1996; Paterson and Hagerthey, 2001). Under temperate latitudes, MPB assemblages are mainly dominated by diatoms that form a brown dense biofilm at the sediment surface during daytime low tides (MacIntyre et al., 1996; Underwood and Kromkamp, 1999). Intertidal flats provide important ecosystem services such as biodiversity depositories, storm protection, and shoreline stabilization (Murray et al., 2018; Legge et al., 2020). They also provide essential food resource for higher trophic levels, from benthic fauna to birds (Herman et al., 2000; Kang et al., 2006; Jardine et al., 2015), and for pelagic organisms when MPB is resuspended in the water column (Perissinotto et al., 2003; Krumme et al., 2008). As such, they contribute to the so-called blue carbon (Otani and Endo, 2019; Legge et al., 2020). With sand and rocky flats, tidal mudflats are one of the most extensive coastal ecosystems, with a recent global area estimation of at least 127,921 km² (Murray et al., 2018). With a global annual gross primary production (GPP) estimated to be in the order of 0.5 Gt C y^–1 (Cahoon, 1999), MPB are of key importance for local and global coastal ecosystem functions, including carbon budget. However, their actual contribution remains unknown, due to the lack of estimation at ecosystem scale (i.e., the entire mudflat). A more comprehensive mapping of these intertidal mudflats is therefore needed to improve the accuracy of coastal carbon budgets (Legge et al., 2020).

The spatial heterogeneity of intertidal mudflats, as well as the high degree of temporal variability in process rates, adds to the challenges of accurately quantifying carbon stocks and flows in coastal areas (Legge et al., 2020). MPB spatiotemporal distribution is highly variable, as it is driven by highly variable physical [photosynthetically active radiation (PAR), mud surface temperature (MST), tides, and waves] and biological (grazing, biostabilization, and bioturbation) factors (e.g., Pinckney and Zingmark, 1991; Kingston, 2002; Cohn et al., 2003; Consalvey et al., 2004; Spilmont et al., 2007; Coelho et al., 2009, 2011; Serôdio et al., 2012; Savelli et al., 2018). Such a variability impedes on an accurate and robust assessment of MPB contribution to the coastal and global marine carbon cycle. To assess MPB biomass and/or GPP, tidal mudflat measurements are usually limited to single-point sampling (e.g., Vieira et al., 2013; Orvain et al., 2014; Cartaxana et al., 2015; Pniewski et al., 2015). Station-based sampling makes it possible to locally assess MPB temporal dynamics, but it fails to describe MPB spatial and temporal variations at the scale of the entire mudflat (Forster and Kromkamp, 2006). Whereas only few studies have resolved MPB variability using sampling campaigns requiring time and important logistical resources (i.e., Guarini et al., 1998; Ubertini et al., 2012), satellite remote sensing appears to be the most efficient upscaling tool. Since the end of the last century, starting with Jobson et al. (1980) initiative, airborne and spaceborne remote sensing methods have been developing increasingly and are now more widely used for MPB studies (Méléder et al., 2003b; van der Wal et al., 2010; Brito et al., 2013; Benyoucef et al., 2014; Echappé et al., 2018). Remote sensing data can cover large spatial scales (from one meter to several kilometers), and vegetation indices such as the normalized difference vegetation index (NDVI) were successfully applied to multispectral broadband satellite sensors to map MPB biomass at the scale of a whole mudflat (e.g., Méléder et al., 2003b; van der Wal et al., 2010; Brito et al., 2013; Benyoucef et al., 2014; Echappé et al., 2018). However, although it is common use to estimate terrestrial vegetation GPP (e.g., Goetz et al., 1999; Huemmrich et al., 2010) and oceanic GPP (e.g., Babin et al., 2015) from space, remote sensing has never been used to map MPB GPP over an entire mudflat before the study by Daggers et al. (2018). While this recent study represents a major contribution to the field, their GPP model does not take into account the seasonal variability of photophysiology (Kromkamp and Forster, 2006) and also strongly depends on the relationship between NDVI and sediment MPB chlorophyll-a (Chl a) concentration, which is known to be highly sensitive to Chl a sampling depth, MPB vertical distribution, and MPB small-scale horizontal variability (Jesus et al., 2006). In order to resolve the NDVI versus Chl a issue, we propose here an original method where GPP is directly calibrated to NDVI. First, a GPP algorithm was obtained using laboratory measurements of NDVI and carbon fluxes (mg C h^–1 m^–2) fitted on Production–Irradiance (P–E) curves. Second, the seasonal variability of the photophysiological parameters was taken into account in a series of laboratory experiments performed during winter, spring, and summer. Third, the NDVI-calibrated GPP algorithm was applied to high-resolution satellite images acquired during the three seasons and coupled to emersion time, mud surface temperature (MST), and light intensity (PAR) obtained from the physical model for Applications at Regional Scale (MARS-3D). Finally, we compared the remotely sensed GPP with field observations, and we discussed the ability of our algorithm to map MPB GPP at mudflat scale and to provide new insights on the role of MPB in the coastal carbon cycle.

Materials and Methods

Study Site: Brouage Mudflat

The study site was located in the Pertuis Charentais Sea, a shallow semi-enclosed sea located on the French Atlantic coast (Figure 1). The tidal regime is semi-diurnal and macrotidal. The tidal range reaches up to ∼6 m during spring tides. This study focused on the Brouage mudflat which extends over 42 km² in the southeastern part of the area (Figure 1). The mudflat sediment is composed of very fine and cohesive grains (median grain size is 17 μm and 85% of grains have a diameter <63 μm; Bocher et al., 2007) distributed on a gentle slope (∼1/1000; Le Hir et al., 2000).

Periods of Investigation

The investigated periods for MPB primary production estimation by multispectral remote sensing and modeling were selected in accordance to the seasonal cycle of MPB biomass (Figure 2). This cycle was extracted from a time series from 2000 to 2015, obtained by the processing of 582 low-tide images from the Moderate Resolution Imaging Spectroradiometer (MODIS) onboard the Terra satellite. MODIS Terra data [namely the Surface Reflectance Daily L2G Global 250 m SIN Grid product (MOD09GQ)] were used because the morning pass [10 to 11 h Universal Time (UT)] was better suited than the MODIS Aqua satellite afternoon pass to observe the emerged mudflat during spring low tides at our study site. The MOD09GQ product provides surface reflectance data in a red band (R_red, from 620 to 670 nm) and a near-infrared (NIR) band (R_NIR, from 841 to 876 nm) at a spatial resolution of 250 m and with a revisit time of 1 to 2 days. This medium spatial resolution has been previously demonstrated to be valid for the study of MPB seasonal dynamics in large intertidal areas (van der Wal et al., 2010). For the present study, 582 MODIS low-tide images were initially downloaded, from which 343 cloud-free scenes were eventually selected. An NDVI (Eq. 1) was computed as a proxy of MPB biomass. The NDVI quantifies the changes in the reflectance’s spectral shape due to Chl a absorption in the red band and to the absence of absorption by pigments in the NIR plateau (Méléder et al., 2003b; Brito et al., 2013; Benyoucef et al., 2014; Echappé et al., 2018). The NDVI is a widely used vegetation index, relatively robust to the variability in sediment backgrounds (Barillé et al., 2011) and mainly driven by changes in algal biomass (i.e., the higher vegetation biomass, the higher the NDVI). Several methods have been developed to discriminate the MPB from other Chl-bearing intertidal vegetation (e.g., seagrass and macroalgae), from simple reflectance thresholds (Méléder et al., 2003b) to more complex clustering methods (Hossain et al., 2015). As previous field studies did not report significant areas colonized by macrophytes in the MPB-dominated Brouage mudflat (Lebreton, personal communication), the NDVI-derived biomass was non-ambiguously assigned here to MPB biofilms.

From this NDVI time series, three periods were selected for field campaigns and for laboratory P–E curve calibration (Figure 2): March, when MPB biomass reaches its highest level; July, corresponding to the lowest level; and May, corresponding to an intermediary level. For the three periods, the weather was also expected to be contrasted: low temperature and light intensity in March, higher temperature and light intensity in July, and intermediate temperature and light in May. Field campaigns occurred the 5th and 6th of May, 2015; the 2nd and 3rd of July, 2015; and the 5th of March, 2018. They included in situ measurements and sediment sampling for further laboratory experiments.

Gross Primary Production (GPP) Measurements

Field Campaigns

The sampling station was located at the “Merignac” site, in the Brouage mudflat (Figure 1, red dot 45°53′11.20^″N; 1°7′538^″W). CO₂ fluxes were measured at the air/sediment interface (enclosed sediment area of 165 cm² down to 5-cm depth) using the closed-chamber method described in Migné et al. (2002). Air CO₂ concentration (ppm) changes were monitored in the benthic chamber (0.8 L) continuously over an incubation period of 20 to 30 min, using an infrared gas analyzer (IRGA Li-840A, LI-COR, Lincoln, NE, United States) connected to a datalogger (Li-1400, LI-COR) with a 30-s frequency. CO₂ flux was calculated as the slope of the linear regression of CO₂ concentration (μmol mol^–1) against time (min) and expressed in mg C m^–2 h^–1.

Transparent chambers were used to estimate the net benthic community production [NCP; balance between the community GPP and the community respiration (CR)]. Dark chambers were used to estimate the CR. Light and dark incubations were performed successively. Due to the tidal cycle duration, a maximum of four incubations were done per day (two transparents and two darks). The GPP expressed in mg C m^–2 h^–1 was computed following Migné et al. (2002) (Eq. 2):

The biomass-specific productivity (P^b) was then computed from GPP. Here, the NDVI was used as a proxy of MPB biomass, and P^b was directly expressed in C m^–2 h^–1 ndvi^–1 (Eq. 3):

The NDVI was computed from MPB reflectance spectra acquired using a JAZ (Ocean Optics, Largo, FL, United States) spectroradiometer (200–1100 nm; sampling: 0.3 nm; spectral resolution: 0.3–10 nm FWHM) pointed at an area of the MPB biofilm close to the benthic chambers. Because the spectral resolution of this detector was higher than a satellite multispectral detector, NDVI calculation was adapted (Eq. 1): the red reflectance, R_red, was considered as the averaged value of data at 675 nm ± 3 nm and the NIR reflectance, R_NIR, as the average at 750 nm ± 3 nm (Méléder et al., 2003a, 2010). NDVI was measured during the same period of each light incubations (equal to NCP estimation), and a time-averaged NDVI was used to standardize GPP in Eq. 3.

Synchronously, the incident photosynthetically active radiation (PAR, from 400 to 700 nm, in μmol photons m^–2 s^–1) and mud surface temperature (MST, in °C) were continuously (every 30 s) measured (LS-C sensor plugged to a ULM-500 data-logger, Walz, Effeltrich, Germany) at the sediment surface, near to the chambers and to the area of reflectance acquisition and measurement of MPB Chl a biomass content. The latter was measured in the first 250 μm of sediment continuously sampled (every 5–10 min) by the “crème brulée” technique (Laviale et al., 2015). This technique, derived from the contact-core technique (Ford and Honeywill, 2002), consists of freezing by contact the top surface of sediment (here 250 μm) with a metal surface (1.5 cm²) previously immersed in liquid nitrogen. The obtained sediment disks were stored in liquid nitrogen during field campaign and were kept at −80°C in the laboratory until pigment analysis. After freeze-drying of the sediment disks, pigments were extracted in a cold mixture (4°C) of 90% methanol/0.2 M ammonium acetate (90/10 vol/vol) and 10% ethyl acetate. Injection, HPLC device (Hitachi Lachrom Elite, Tokyo, Japan), pigment identification, and quantification were detailed before (Roy et al., 2011; Barnett et al., 2015). The Chl a amount in sediment was standardized to the sampled surface (1.5 cm²) in order to be expressed in mg Chl a m^–2.

In an aim to assess MPB photophysiological status at tide time scale and check if changes occurred during incubation time, several photophysiological parameters (F_v’/F_m’, rETR, α, and E_k) were measured continuously (every 5 to 10 min) using a Water-PAM fluorimeter (Fiber version, Walz, Effeltrich, Germany). For details, see Supplementary Appendix A. All in situ data (GPP, PAR, MST, biomass, and photophysiological parameters) were used as ground-truthing for laboratory P–E curve calibration and the remotely sensed GPP validation.

Laboratory Experiment for Production–Irradiance Curve Calibration

During the three field campaigns, the upper layer (approx. the top first centimeter) of sediment was collected and brought back to the laboratory. The mud was cleaned of fauna by sieving through a 500-μm mesh. The sediment was homogenized by thoroughly mixing and was spread as plane layer in plastic trays of 4-cm depth (Serôdio et al., 2012). A water layer was added for the night and sediment was left undisturbed overnight. The next morning, the water layer was manually removed by a syringe 3 h before the lowest water level timing expected in situ (i.e., at sampling site) and trays were kept in darkness at 22 ± 1°C. Experimentation started 1 h later, when the MPB biofilm started to darken in color at the sediment surface. It consisted of lighting up the trays one by one with an LED panel (LED Light SL 3500-E, Photo System Instrument, Tøeboò, Czechia) at a given PAR (=E), varied from 5 to 2,200 μmol photons m^–2 s^–1 for 30 min at 22 ± 1°C, exactly in the optimal temperature range for MPB production (according to Blanchard et al., 1997). During the lightening, NCP was estimated through transparent benthic chamber incubations, following the same method as the field measurements (Migné et al., 2002). At the beginning and end of each light incubation, NDVI, sediment Chl a content, and photophysiological parameters using PAM-fluorimetry were measured the same way as in situ (see above). After illumination, CR was measured during the 20- to 30-min dark benthic chamber incubations. The same process was repeated as long as MPB was present at the surface of sediment (i.e., with similar measured NDVI) and for 8 to 11 Es. For each E, GPP was calculated as the sum between NCP and CR as done in situ (Eq. 2, Migné et al., 2002). GPP was then standardized by NDVI in order to directly express P^b in mg C m^–2 h^–1 ndvi^–1 (Eq. 3), so that the GPP algorithm could be consistently applied to satellite-derived NDVI maps (see section “Coupling Remote Sensing Data and Modeling: the GPP-Algorithm”).

For the three campaigns (March, May, and July), a season-dependent P^b was obtained. For each season, several P–E models widely used in the literature [namely Platt et al. (1980), Eilers and Peeters (1988), Steele (1962), Platt and Jassby (1976), and the modified version of Platt and Jassby (1976)] were fitted to the P–E curves to select the best model to be integrated into the GPP-algorithm. The selection of the most appropriate model was done using the determination coefficient (r²) and residual standard deviation (RSD) calculated using the in R-software. For details, see Supplementary Appendix B.

Data Analysis for Field Campaign and Laboratory Experiments

All data are available for download: 10.5281/zenodo0.3862068. Changes in MPB biomass (NDVI and Chl a content), photophysiological parameters measured by PAM-fluorimetry (F_v’/F_m’, α, rETRm and E_k) and GPP and P^b were detected at a tidal time scale (i.e., during incubation) and at a monthly scale (March, May, and July). Knowing the potential high variability of environmental conditions (i.e., PAR, MST, and light dose calculated from PAR and emersion time), the objective was to check if biomass and photophysiological status changed drastically or not during the light incubation estimating GPP as well in situ and at the laboratory. In this aim, Spearman or Pearson correlations, and ANOVA or Kruskal–Wallis (KW) test were performed on all data using R software, after a Shapiro test, to test data normality. For details, see the Supplementary Appendix A.

Coupling Remote Sensing Data and Modeling: the GPP-Algorithm

The successive steps in the remotely sensed GPP procedure are described below and synthesized in Figure 3.

Remote Sensing Data

High-resolution satellite remote sensing was used to upscale the GPP algorithm to the whole mudflat, using the season-dependent, NDVI-specific P^b parameter. Multispectral images were selected from the SPOT, Landsat, and Pléiades archives following three acquisition criteria: (1) acquisition day as close as possible to the field campaign day, (2) acquisition time as close as possible to low tide timing (i.e., the lowest water height, when mudflat was the most exposed), and (3) cloud-free observation (cloud cover < 10%) with an almost zenithal sun. These criteria allowed us to select three images from the SPOT6 and Pléiades archive (one per field campaign, Table 1). The top-of-atmosphere data were converted into surface reflectance using the fast line of sight atmospheric analysis of spectral hypercubes (FLAASH, Matthew et al., 2000) atmospheric correction using the ENVI software. For a consistent analysis of the SPOT and Pléiades images, the same FLAASH parameters were applied to each image: United States atmospheric model with a visibility of 40 km, and a maritime aerosol model. The spectral images were registered in the WGS 84 UTM 30N coordinate system. Finally, the NDVI was calculated from the surface reflectance following Eq. (1) to map MPB biomass (Méléder et al., 2003b). As the SPOT6 and Pléiades data have similar spectral characteristics, the NDVI was not recalibrated between the two sensors (Echappé et al., 2018). The satellite-derived NDVI maps were used as inputs for the GPP-algorithm, in complement to other data (Figure 3).

Tidal Height and PAR Modeling

The tidal height and PAR were simulated over the Brouage mudflat by the 3D hydrodynamical MARS-3D (Figure 3). The bathymetry was extracted from the model numerical grid. The Navier–Stokes primitive equations were solved under assumptions of Boussinesq approximation, hydrostatic equilibrium, and incompressibility (Blumberg and Mellor, 1987; Lazure and Dumas, 2008). The numerical domain of the Pertuis Charentais Sea consisted of 100- × 100-m grid cells discretized over 20 sigma-levels. For details, see Lazure and Dumas (2008). The meteorological forces (i.e., 10 m wind speed, 2 m air temperature and relative humidity, atmospheric pressure at sea level, nebulosity, and solar fluxes) used to constrain MARS-3D were extracted from the Meteo France AROME model¹. The tidal model CST FRANCE, developed by the SHOM (Simon and Gonella, 2007), forced MARS-3D at the domain boundaries. The tidal model solved the amplitude and phase of 115 harmonic constituents. Initial and boundary conditions of seawater temperature, salinity, current velocity, and sea surface height were extracted from the MANGAE 2500 Ifremer model of 2.5-km lateral resolution (Lazure et al., 2009). Simulated tidal height associated to bathymetry allowed to estimate the beginning, the end and thus the duration of emersion period of each grid cells. During emersion period, PAR intensity varied every 10 min at 100- × 100-m spatial resolution, but values were interpolated on the horizontal grid of satellite data (2 or 6 m, see Table 1) and used as inputs in the GPP-algorithm (Figure 3).

Mud Surface Temperature Modeling

The simulated mud surface temperature (MST) was obtained from the coupling of the MST model of Savelli et al. (2018) with MARS-3D. The simulated heat fluxes in a 1-cm-deep sediment layer were solved by thermodynamic equations detailed in Savelli et al. (2018). The horizontal fluxes of heat were neglected. During exposure periods, the simulated MST resulted from heat exchanges between the sun, the atmosphere, and the sediment surface; from the conduction between mud and air; and from evaporation. The simulated MST of immersed mud was set to the temperature of the overlying seawater simulated by MARS-3D since MST was not used in this study. The MST differential equation was solved by the MARS-3D numerical scheme. For more details, see Savelli et al. (2018). During emersion period, MST varied at the same time resolution than PAR intensity (10 min and 100 × 100 m). But, as PAR, the MST values were interpolated on the horizontal grid of satellite data (2 or 6 m, see Table 1) and was used as input in the GPP algorithm (Figure 3).

GPP Algorithm

Finally, the GPP algorithm (Figure 3) coupled NDVI maps from SPOT and Pléiades scenes with forces by hydrodynamical MARS-3D. Whereas the NDVI estimates the horizontal distribution of MPB biomass, the MARS-3D model simulates the emersion time over the whole mudflat using the bathymetry and the tidal height. Coupled with the MPB biomass, the emersion time determined the photosynthetically active biomass at the mud surface. The MPB biomass detected by satellite was assumed to correspond to the fully established biofilm during the daytime low tide (total photosynthetically active biomass). The migration behavior of MPB was introduced in the model through a progressive establishment of the total photosynthetically active biomass at the sediment surface, known to take place during 20 min at Brouage mudflat (Herlory et al., 2004). The MPB started to migrate at the sediment surface just after water removal to reach 50% of the total amount of the photosynthetically active biomass (=NDVI/2) after 10 min of emersion of respective grid cells. After 20 min, the MPB biofilm at the sediment surface was fully formed. 20 min before the immersion, downward migration started and only the half of the total photosynthetically active biomass was still at the surface 10 min before immersion. No biomass was at the surface when water overlaid the sediment. PAR and MST simulated by MARS-3D were used to constrain the algorithm: the selected P–E model and its respective parameters values fitted on the laboratory measurements were used to compute the P^b according to the simulated light conditions. The effect of MST was simulated using the Blanchard et al. (1996) model to compute the P^b variations according to the simulated temperature. Finally, combined with the horizontal distribution of the MPB biomass of the NDVI maps, the P^b (mg C m^–2 h^–1 ndvi^–1) was further used to map the remotely sensed GPP (mg C m^–2 h^–1). The time resolution was 10 min (following PAR and MST variations), whereas the spatial resolution was the one of the SPOT or Pléiades image (2 or 6 m, see Table 1).

The current configuration of the GPP-algorithm was consistent only for intertidal mudflats composed by very fine cohesive grains over the Brouage mudflat (Poirier et al., 2010). Consequently, no-muddy areas were excluded of the GPP-algorithm and MST model used was developed and validated for the mud only (Savelli et al., 2018). Therefore, laboratory and in situ measurements were conducted only on muddy sediment. The GPP maps thus corresponded to MPB assemblages known to be dominated by epipelic diatoms at the study site (Haubois et al., 2005; Du et al., 2017).

GPP Maps Validation

The validation of the GPP model was performed in three steps for each seasonal experiment, using the field data acquired at the study site (see section “Gross Primary Production (GPP) measurements”). First, the simulated MST and PAR were compared to in situ measurements using a Mann Whitney (MW) test to assess the accuracy of the physical model. Second, the satellite derived NDVI was validated against field measurements. The GPP was then computed from satellite-derived NDVI maps coupled with the modeled PAR and MST data. The remotely sensed GPP maps were averaged hourly (mg C m^–2 h^–1) over the emersion period, as well as daily-integrated period (mg C m^–2 d^–1). Finally, the hourly remotely sensed GPP data was compared to in situ GPP measurement using a MW test.

To ensure that the delay between measurements for P–E calibration and image acquisition (2, 13, and 14 days in March, May, and July, respectively) was not an issue, simulated PAR were averaged during daytime emersion periods 2 weeks before each sampling/measurement session and image acquisition to be compared (MW tests). The main limitation was a change in physiology and metabolic acclimation status of MPB cells due to a change in light and temperature conditions during this delay, preventing the use of P–E photophysiological parameters retrieved from laboratory experiments to calibrate images acquired several days later.

Results

Microphytobenthos GPP Variation Between Months

Field Campaigns

During field campaigns, PAR intensity and MST were significantly different between months (Figures 4a,b; KW test, p ≤ 0.001). Minimum values were mainly measured in March for PAR and MST with 624 ± 10 μmol photons m^–2 s^–1 (mean ± SE) and 15.7 ± 0.1°C, respectively, whereas maximum values were reached in May for PAR (1,509 ± 27 μmol photons m^–2 s^–1) and July for MST (30.7 ± 0.2°C) (Figures 4a,b). According to its seasonal cycle (Figure 2), MPB biomass was higher in March than in May and July (KW test, p≤0.001) with averaged NDVI values and Chl a sediment content of 0.61 ± 0.01 and 95.2 ± 3.1 mg m^–2, respectively (Figures 4c,d). The GPP did not vary significantly among campaigns, although the minimum value was measured in March, and the highest was measured in July (Figure 4e). When GPP was standardized by NDVI to be expressed into biomass-specific productivity (P^b), this difference was more visible (Figure 4f) even if it was not significant (KW test, p = 0.1). Regarding photophysiological parameters measured by PAM-fluorimetry (F_v’/F_m’, α, rETRm, and E_k), see Supplementary Appendix A. At the tidal scale, biomass (NDVI and Chl a content), but also PAM photophysiological parameters (F_v’/F_m’, α, rETRm, and E_k) changed with PAR, light dose, and MST, illustrating the rapid responses of MPB (behavioral migration and/or physiology) to environmental conditions.

Laboratory Measurements for P–E Curve Calibration

During laboratory experiments, whereas temperature and PAR were controlled, the MPB biomass, expressed in NDVI or Chl a sediment content, changed with the sampling campaign date; as for in situ, it was the highest in March, with an averaged value of 0.57 ± 0.01 for NDVI and 97.7 ± 5.5 mg Chl a m^–2 (three-way ANOVA, p ≤ 0.001; Figure 5 and Supplementary Table A3). The biomass at the sediment surface of the plastic trays globally did not change between PAR tested and during incubation time for a given month (see Supplementary Table A3 for details). Consequently, GPP measured from benthic chamber incubations for each sediment tray were considered to correspond to a same MPB biomass. Biomass specific productivity, P^b was then obtained by dividing GPP by the averaged NDVI value over the sediment trays for each month: 0.57 ± 0.01 in March, 0.08 ± 0.01 in May, and 0.20 ± 0.01 in July (Supplementary Figure A3a). The shape of the relationship between P^b and the irradiance provided by artificial lighting (i.e., P–E curves) varied with seasons (Figure 5), as well the fitted photophysiological parameters from the five P–E models tested (Supplementary Table B1). The high r² and low RSD values demonstrated that all P-E models fitted well with the laboratory measurements (Supplementary Table B2). Because it exhibited the best fit to P–E laboratory measurements, the Eilers and Peeters (1988) model was selected to estimate and map remotely-sensed GPP using the GPP-algorithm.

Mapping Microphytobenthos GPP From NDVI Maps

Overall, SPOT and Pléiades provided consistent spatial distribution of MPB biomass, with a higher NDVI in the middle and lower areas of the mudflat than in the upper shore, especially in March and May (Figure 6). The Pléiades-derived NDVI varied from 0 to 0.4, with an averaged value of 0.2 ± 0.09 in March and of 0.14 ± 0.05 in July (Figures 6a,c), whereas the SPOT6-derived NDVI varied from 0 to 0.3, with an averaged value of 0.14 ± 0.05 in May (Figure 6b). Besides the seasonal variability, the difference in the sensors’ spatial resolution could also partly explain the differences between SPOT6 and Pléiades. Compared to Pléiades (2 m), the lower spatial resolution of SPOT6 (6 m) smoothed out the fine-scale distribution of MPB biofilms, thus resulting in lower NDVI averages.

Spatially resolved GPP rates were then modeled from the satellite-derived NDVI maps for each season (Figure 7). The short-scale temporal variability of the production factors was taken into account using the hourly PAR and MST MARS-3D simulations and the laboratory P–E curves (Figure 5 and Supplementary Table B1). For each date of satellite acquisition, remotely sensed GPP was averaged over the daytime emersion period (Figures 7a–c) and eventually integrated to yield daily GPP maps (Figures 7d–f). GPP was at its maximum in March and its minimum in July (Figures 7a,c). The most productive areas of the mudflat were the middle and lower shores, with values up to 14.4 mg C m^–2 h^–1 in March and 10.8 mg C m^–2 h^–1 in May (Figures 7a,b). The upper shore was less productive with an hourly GPP of ∼7.5 mg C m^–2 h^–1 in March and ∼6.0 mg C m^–2 h^–1 in May (Figures 7a,b). The hourly GPP exhibited no spatial pattern in July, and GPP was low over the entire mudflat (∼1.8 mg C m^–2 h^–1, Figure 7c). The mean daily-integrated GPP, reaching 50.2 ± 30.1 mg C m^–2 d^–1 in March, 40.9 ± 26.8 mg C m^–2 d^–1 in May, and 22.3 ± 20.5 mg C m^–2 d^–1 in July, allowed to integrate GPP over the mudflat, reaching, respectively, 2.06 (for an emerged surface of 41 km²), 1.42 (for 34.5 km²), and 0.80 tC d^–1 (for 36 km²) in March, May, and July (Figures 7d–f).

Validations

A Priori Physical Environment Validation

The MST and PAR conditions simulated by the MARS-3D model satisfactorily compared to the in situ conditions, even if there were some significant differences (Figures 8a,b). The delay between field campaigns and image acquisition was not an issue (Figures 8c,d): for the three periods, the averaged PAR and the MST 2 weeks before measurements and before images acquisition were not significantly different (Figures 8c,d). These observations mean that MPB cells were expected to be in a similar acclimation status during laboratory experiments and image acquisition. After these first validations, simulated MST and PAR, satellite images and calibrate P–E model were used in the GPP algorithm to predict and map GPP.

Measured Versus Estimated NDVI and GPP

The NDVI measured in situ was always higher than the remotely sensed NDVI at the study site (Figure 9a). In March, the NDVI measured in situ (0.61 ± 0.03) was almost 5-fold higher than the remotely sensed NDVI at the pixel corresponding to the study site (0.14). The NDVI measured in situ in May (0.14 ± 0.02) was 1.5-fold higher than the remotely sensed NDVI at the pixel of the study site (0.09). In July, the NDVI measured in situ (0.18 ± 0.09) was almost 1.8-fold higher than the remotely sensed NDVI at the pixel of the study site (0.1). Regarding the GPP, the measured value in situ varied from 4.8 ± 2.1 mg C m^–2 h^–1 in March to 6.3 ± 0.3 mg C m^–2 h^–1 in July (Figures 4e, 9b), whereas the GPP remotely-sensed at the respective grid cell averaged during daytime emersion varied from 2.2 ± 1.4 mg C m^–2 h^–1 in July to 7.8 ± 3.1 mg C m^–2 h^–1 in March (Figure 9b). However, there was not significant difference between in situ and remotely sensed GPP (MW tests, Figure 9b).

Discussion

NDVI Versus Chl a GPP Standardization

In order to be able to upscale MPB biomass over an entire mudflat, NDVI has been used as a proxy of MPB abundance (e.g., Méléder et al., 2003b; van der Wal et al., 2010; Brito et al., 2013; Benyoucef et al., 2014; Echappé et al., 2018). Here, the NDVI 2000–2015 time series demonstrated a seasonal cycle characterized by a maximum of MPB biomass in winter–spring and a minimum in summer. This MPB seasonality observed in the Brouage mudflat is consistent with the seasonal pattern reported before for the same mudflat (Cariou-Le Gall and Blanchard, 1995; Savelli et al., 2018) and for other Northern European mudflats (van der Wal et al., 2010; Echappé et al., 2018).

Usually, on a local scale, GPP is standardized to the Chl a sediment content (e.g., Cahoon, 2006). Ideally, in order to convert our NDVI maps into GPP, a conversion of NDVI in Chl a sediment content should be performed. However, a direct relationship between NDVI and Chl a is a real issue for three main reasons. First, the relationship is known to be non-linear and to saturate for high values of Chl a (Méléder et al., 2003a; Serôdio et al., 2009; Barillé et al., 2011; Daggers et al., 2018). Second, there is an ongoing debate on the sediment depth to be sampled for Chl a content estimation in order to be in accordance with the sediment depth detected by sensors for NDVI calculation. Jesus et al. (2006) suggested the sampling of the first 150 μm, however, all depths tested along the first 2 mm were highly correlated, leading to similar NDVI for different Chl a contents (i.e., different depths). Moreover, the optical depth varies with the MPB biomass at the sediment surface, the sediment texture, the organic and water content, and the incident light wavelengths (Kühl et al., 1994; Jesus et al., 2006; Barillé et al., 2011; Kazemipour et al., 2011; Fisher et al., 2018). In the current study, the length of the path of the reflected light (i.e., the NDVI) is assumed to correspond to the photic zone, where the biomass is photosynthetically active, which rarely exceeds 500 μm for muddy sediments (Cartaxana et al., 2011). Third, the dilution of the reflectance signal from the sediment surface to the satellite sensor is variable. MPB-specific NDVI data obtained in this study from MODIS-Terra, SPOT6, and Pléiades satellites and measured in situ with a handheld field spectroradiometer reach, respectively, maximal values of 0.25, 0.30, and 0.40 for satellite and 0.60 in situ, meaning that maximum NDVI value decreases with spatial resolution. Even if such values are in the range of the MPB-specific NDVI derived from satellite data over temperate mudflats (Méléder et al., 2003b; van der Wal et al., 2010; Brito et al., 2013; Benyoucef et al., 2014; Echappé et al., 2018), their variability illustrates the different sensitivity between devices, but also the dilution of the reflectance signal for larger surfaces, mainly due to the patchiness distribution of the MPB biofilm (Saburova et al., 1995; Jesus et al., 2006; Spilmont et al., 2011). Méléder et al. (2010) and Launeau et al. (2018) demonstrated previously that MPB patchiness is responsible for a reflectance signal, due to the non-linear mixing of individual MPB patches and apparent mud. This generates one Chl a content corresponding to diverse NDVI values and conversely, one NDVI value could correspond to diverse Chl a content. This non-linear mixing increases the difficulty to use the Chl a sediment content measured at local scale (a few millimeters squared) to map NDVI over several meters squared or kilometers squared. In the present study, the direct calibration of the GPP-algorithm in NDVI allowed to decrease bias due to uncertainties of the NDVI versus Chl a relationship. Nevertheless, and in spite of this major improvement, some issues remain, such as the NDVI saturation for high MPB biomass leading to a potential non-linear GPP-NDVI relation, which will need further investigation.

GPP Algorithm Physical Setting

The MST and PAR conditions simulated by the MARS-3D model compared well to the conditions measured in the field. The model-versus-in situ data comparison suggested that the 3D model can resolve with confidence the physical environment experienced by MPB at the sediment surface. In regards to the frequency of the atmospheric AROME model (1 h), the simulated PAR conditions varied less than the in situ observations. The 3D model could not reproduce the observed synoptic variations of irradiance at the sub-hourly scale that can induce a substantial variability in the remotely sensed GPP over an emersion period. In addition, the horizontal resolution of the 3D model (100 × 100 m) may also translate into model–data discrepancies. In the GPP-algorithm, the bathymetric level and simulated water height originating from 100- × 100-m grid cells delay the emersion timing by ∼30 min. However, considering this delay, the comparison of in situ and simulated physical conditions and GPP were made on the corresponding low tides. Most importantly, the preservation of the horizontal resolution of satellite data in order to capture the MPB patchiness suggests that the GPP algorithm can resolve with confidence the overall dynamics of MPB GPP at the tidal scale.

Ability of the GPP Algorithm to Map the Current Productive State of the Mudflat

Our study is not the first to assess MPB primary production coupling remotely sensing and physical–biological modeling. Daggers et al. (2018) proposed a first approach using (i) remotely sensed information on MPB biomass and on sediment mud content, (ii) surface irradiance and ambient temperature, (iii) directly measured photophysiological parameters by PAM fluorimetry, and (iv) a tidal model. The current GPP algorithm improves the (Daggers et al., 2018) approach by (i) the use of NDVI rather than the conversion of NDVI into Chl a, which introduces uncertainties (see section “Discussion” above); (ii) the use of MST rather than ambient temperature, which was identified as a weakness of their model; and (iii) photophysiological parameters derived directly from benthic chamber CO₂ exchange measurements on natural MPB communities in sediment under controlled conditions, rather than an averaged electron requirement (EE) rate derived from PAM fluorimetry. The EE used by Daggers et al. (2018) corresponds to the ETR efficiency for C fixation to translate ETR (μmol electrons m^–2 s^–1) into C (mg C m^–2 h^–1). However, it is known to vary with season, species, and site (Barranguet and Kromkamp, 2000), and the relationship between ETR and C-fixation can be non-linear, especially at irradiances exceeding E_k. The difficulty to use photophysiological parameters derived from PAM fluorimetry to predict C fixation was confirmed by the current study. PAM photophysiological parameters could vary rapidly even during the time of benthic chamber incubation (30-min duration) rendering the use of an averaged EE rate weakly representative of a given season, a given day and even a given emersion. To overcome this issue, we suggest to directly calibrate the P–E model with C-fixation standardized to NDVI values.

Additionally, the environmental conditions 2 weeks before measurements for the GPP algorithm calibration and 2 weeks before the acquisition of satellite images used to apply GPP algorithm were checked in the aim to support the representativeness of the photophysiological parameters. The conditions were similar, assuming the same photosynthetic capabilities of MPB during experiments for the calibration and during image acquisition. Otherwise, it would have been hazardous to apply the GPP algorithm on MPB that would have been differently acclimated between lab experiments and satellite image acquisitions.

The GPP algorithm uses a vertical migration scheme of MPB biomass within the upper layer of sediment, which is represented through the modulation of the total photosynthetically active biomass detected from the remotely sensed NDVI. Such a migration scheme in the GPP algorithm was set according to the observation of the progressive superficial sediment covering by MPB during emersion at our study site (Herlory et al., 2004). However, the migration speed can be faster [a few minutes; see Méléder et al. (2003b)] or slower [one hour; see Paterson et al. (1998)], and it is mainly controlled by the tidal cycle and the light climate and spectral quality (Pinckney and Zingmark, 1991; Spilmont et al., 2007; Coelho et al., 2011; Barnett et al., 2020; Prins et al., 2020), but also by temperature (Cohn et al., 2003), nutrient availability in the sub-surface of sediment (Kingston, 2002), and desiccation (Coelho et al., 2009). Currently, the GPP-algorithm does not include the short-term variations of MPB photosynthetically active biomass at the sediment surface (i.e., “micro-migrations”), as it has been also observed some days/timings during our field campaigns. As suggested previously by Daggers et al. (2018), further research on the mechanisms and triggers that determine the vertical phototaxis of diatom cells in sediment is required to better predict and model intertidal MPB migration patterns and therefore changes in MPB photosynthetically active biomass at the sediment surface.

When constrained by the PAR and MST conditions simulated by the MARS-3D model, the remotely sensed GPP predicted by the GPP-algorithm reached rates up to 14 mg C m^–2 h^–1 or 100 mg C m^–2 d^–1, similar to rates reported for other European mudflats (Barranguet et al., 1998; Underwood and Kromkamp, 1999; Cahoon, 2006; Hubas et al., 2006; Daggers et al., 2018; Frankenbach et al., 2020). These values further supports the key role of MPB in supporting the high productivity of temperate intertidal bare mudflats (e.g., Cahoon, 1999; Underwood and Kromkamp, 1999; Barranguet and Kromkamp, 2000) and its paramount support to local socio-economics (Lebreton et al., 2019). This growing recognition also supports the necessity for the worldwide consideration of mudflats as key ecosystems in the marine global carbon budget (Ciais et al., 2014; Legge et al., 2020).

The predicted and measured GPP reasonably compared. Both GPP show low seasonal variability, whereas the MPB biomass displays a pronounced seasonal cycle. This leads to lower P^b values during the period with the highest MPB biomass (i.e., March). It has been shown before that high MPB biomass does not always generate high production (Barranguet et al., 1998). During winter, the MPB biomass standing stock increases due to lower grazing activity (Thompson et al., 2000). As a consequence, MPB photosynthetically active biomass is mainly concentrated in an extremely thin layer at the sediment surface. However, such high concentration induces strong competition for light and nutrients, limiting the biomass specific productivity (P^b) (Barranguet et al., 1998; Stal, 2010; Vieira et al., 2016). This competition, compensated by high biomass, suggests a bottom-up regulation of the GPP in winter (by light, temperature, and nutrients), whereas a top-down control (by grazing activity) occurs in spring and summer.

Conclusion and Perspectives: The Use of Hyperspectral Remote Sensing

The GPP algorithm developed in this study combines satellite remote sensing, laboratory measurements, and a 3D physical model. The algorithm is constrained by realistic simulated 2D fields of tidal heights, MST, and PAR. It is standardized by photophysiological parameters estimated from laboratory measurements on natural MPB communities in sediment and expressed in C fixation rates. In addition, the direct calibration of the algorithm in NDVI is a step forward to limiting outstanding bias due to NDVI to Chl a conversion. Moreover, the calibration of the GPP algorithm to NDVI allowed us to consistently apply the algorithm to satellite images. This study shows that:

• The NDVI data retrieved from the SPOT6 and Pléiades sensors were consistent with the seasonality of the MPB biomass previously reported for the study site, and their range was comparable to NDVI data from other European mudflats;

• The GPP-algorithm yields MPB GPP rates in the range of in situ GPP measurements, including the seasonal variability of GPP;

• The GPP algorithm was well-adapted to intertidal mudflats mostly composed by fine cohesive sediments dominated by motile epipelic diatoms, and it could be applied for similar habitats.

However, this study highlights several challenging issues that need to be tackled to better estimate MPB production on regional and global scales from in situ information:

• Photosynthetic ability changes over a range of time scales, from the emersion to the season via the tidal fortnight cycle. To overcome this issue, photophysiological parameters, and not only MPB biomass, have to be measured at the ecosystem level (i.e., entire mudflat scale). Currently, only hyperspectral remote sensing is able to capture such a detailed information based on fine pigment absorption features, as recently proposed for terrestrial vegetation (DuBois et al., 2018; Lees et al., 2018). This is an approach we have started to develop successfully on MPB (Méléder et al., 2018).

• Up-scaling remains the main issue. It could be overcome also using hyperspectral remote sensing as demonstrated recently by Launeau et al. (2018). The use of the optical properties retrieved from hyperspectral images to predict MPB biomass allows to remove the patchiness effect (=non-linear mixing) which build a linear relationship between MPB biomass and optical properties independent of the size of the analyzed surface. This approach could be applied for GPP mapping to improve the up-scaling bias.

Data Availability Statement

All the data, inlcuding maps from GPP-algotithm, are avaible on a data repository (Zenodo): doi: 10.5281/zenodo.3862068.

Author Contributions

VM, RS, VL, and JL developed the theoretical formalism. VM, RS, AB, PP, AL, and JL carried out the experiments. PG, AL, and ALB processed the remote sensing data. VM, RS, and VF designed the model and the computational framework. VM, RS, AB, and JL analyzed the data. VM and RS wrote the manuscript with inputs from all authors. All authors contributed to the final version of the manuscript.

Funding

This work was supported by (i) the DYCOFEL project, funded through the 2015 Fondation de France call “Quels littoraux pour demain?”; (ii) the MIMOSA project, funded through the 2018 CNRS EC2CO-LEFE call; (iii) the HYPEDDY project, funded through the 2018–2020 Tosca-CNES call; (iv) the BIO-Tide project, funded through the 2015–2016 BiodivERsA COFUND call for research proposals, with the national funders BelSPO, FWO, ANR, and SNSF; (v) the public funds received in the framework of GEOSUD, a project (ANR-10-EQPX-20) of the program “Investissements d’Avenir,” managed by the French National Research Agency; and (vi) the projects Littoral 1 and ECONAT funded by the Contrat de Plan Etat-Région (CPER) and the CNRS and the European Regional Development Fund. Pléiades and SPOT images were acquired by CNES’s ISIS program, facilitating scientific access to imagery. Pléiades CNES 2015, 2018, Distribution Airbus DS, all rights reserved. Commercial uses forbidden. This research was part of fulfillment of the requirements for a Ph.D. degree (RS) at the Université de La Rochelle, France. RS was supported by a Ph.D. fellowship from the French Ministry of Higher Education, Research and Innovation.

Conflict of Interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Acknowledgments

We thank the CNRS, Ifremer, the Université de Nantes, and the Université de La Rochelle for their support and V. Ouisse and N. Spilmont for their help with setting up the benthic chambers and with analyzing the CO₂ data.

Supplementary Material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fmars.2020.00520/full#supplementary-material

Footnotes

1. ^ https://donneespubliques.meteofrance.fr

References