Improving Remote Sensing-based Estimation of Mangrove Forest Gross Primary Production by Quantifying Environmental Stressors: Sea Surface Temperature, Salinity, and Photosynthetic Active Radiation


 Mangrove ecosystems play an important role in global carbon budget, however, the quantitative relationships between environmental drivers and productivity in these forests remain poorly understood. This study presented a remote sensing (RS)-based productivity model to estimate the light use efficiency (LUE) and gross primary production (GPP) of mangrove forests in China. Firstly, LUE model considered the effects of tidal inundation and therefore involved sea surface temperature (SST) and salinity as environmental scalars. Secondly, the downscaling effect of photosynthetic active radiation (PAR) on the mangrove LUE was quantified according to different PAR values. Thirdly, the maximum LUE varied with temperature and was therefore determined based on the response of daytime net ecosystem exchange and PAR at different temperatures. Lastly, GPP was estimated by combining the LUE model with the fraction of absorbed photosynthetically active radiation from Sentinel-2 images. The results showed that the LUE model developed for mangrove forests has higher overall accuracy (RMSE = 0.0051, R2 = 0.64) than the terrestrial model (RMSE = 0.0220, R2 = 0.24). The main environmental stressor for the photosynthesis of mangrove forests in China was PAR. The estimated GPP was, in general, in agreement with the in-situ measurement from the two carbon flux towers. Compared to the MODIS GPP product, the derived GPP had higher accuracy, with RMSE improving from 39.09 to 19.05 g C/m2/8 days in 2012, and from 33.76 to 19.51 g C/m2/8 days in 2015. The spatiotemporal distributions of the mangrove GPP revealed that GPP was most strongly controlled by environmental conditions, especially temperature and PAR, as well as the distribution of mangroves. These results demonstrate the potential of the RS-based productivity model for scaling up GPP in mangrove forests, a key to explore the carbon cycle of mangrove ecosystems at national and global scales.


31
Mangrove forest is one of the most carbon-rich ecosystems whose carbon sequestration is considerably 32 higher than terrestrial forests 1 . The estimate of the gross primary production (GPP) is important to 33 understand the carbon cycle in mangrove ecosystems. Carbon flux data measured with eddy 34 covariance (EC) techniques provide invaluable information on ecosystem productivities and can be 35 used to establish productivity models 2 . However, these models were limited to a 0.1 to 2 km spatial 36 footprint around the towers, and therefore, applying them at other sites remains challenging due to the 37 variation of GPP across species, structural features, and latitudinal locations. 38 Remote sensing (RS) provides the opportunity to characterize the ecosystem structures and 39 environmental conditions and therefore, estimate the productivity of the ecosystems 3 . The light use 40 efficiency (LUE) model was widely adopted to estimate GPP 4,5 . Currently, GPP models for terrestrial 41 forest are applicable on a global scale (e.g., C-fix, MOD17, and GLO-PEM) 6-9 , however, production 42 models have not been evaluated and employed in mangrove forests in a large scale, mainly due to the 43 lack of understanding of carbon exchange in mangrove forests and measurements from flux tower. 44 Compared to terrestrial ecosystems, mangrove ecosystems are periodically inundated by the tides 45 which contribute to the waterlogged and high salinity soil environment. Although mangroves have 46 developed special structures or tissues to adapt to such demanding surroundings such as the aerial root, 47 thick canopy, and salt-tolerance tissues, the environmental stresses remain critical to mangrove 48 productivity. In addition to being affected by air temperature (Tair) and vapor pressure deficit (VPD) 49 as terrestrial forests 10,11 , mangrove forests are also influenced by the sea surface temperature (SST), 50 salinity, and photosynthetic active radiation (PAR). Firstly, SST affects the roots and aboveground 51 metabolism of mangroves 12 . The high SST would increase the respiration rate. To minimize the water 52 loss and energy consumption, the stomatal conductance within the mangrove would be reduced, which 53 could lower the mangrove light saturation point (LSP) and constrain the photosynthesis. Secondly, the 54 salinity of surface water and porewater represents a significant control on the mangrove LUE which 55 is strongly related to the sea surface salinity, rainfall, and river discharge 13,14 . The high salinity leads 56 to the negative osmotic pressure in the environment of roots which limits the water supply, and 57 therefore inhibits the photosynthesis net photosynthetic rate 15,16 . Thirdly, the relatively low LSP makes 58 the mangrove easy to reach light-saturated status 17 . Hence, the high PAR condition would bring excess 59 light absorption and heat to the canopy that reduces the LUE of the canopy. These typical 60 environmental stressors are not well understood and quantified. Studies have typically focused on the 61 seasonal dynamics and interannual variation of carbon fluxes through modeling GPP 17-20 or tidal 62 effects on CO2 exchange based on in-situ measurements 21-24 . Barr Currently, there are no RS-based productivity products for mangrove forests globally. Therefore, 68 scaling up carbon fluxes from flux tower to national and global scales considering the coastal 69 environment is of great importance and challenge. Modeling the GPP of mangrove forests provides 70 the first step in using RS to discover the role of mangrove ecosystems in global carbon budgets. 71 Therefore, the objectives of this study are 1) to improve the LUE model for mangroves 72 considering environmental stresses in coastal zone (SST, salinity, and PAR), 2) to estimate the GPP of 73 mangroves in the whole coastal zone of China combining flux tower-based measurements and RS, 3) 74 to map the spatiotemporal distributions of mangrove productivity and analyze the possible affecting 75 factors. 76

77
Effects of environmental stressors on mangrove LUE.

78
LUEmax: The GPPmax and Re for two mangrove forests are summarized in Table 1 based on the flux 79 tower data. The NEE-PAR fitted curves are displayed in Figure 1. The initial slope of each rectangular 80 hyperbola was calculated as the LUEmax and listed in figure, we can clearly see that RS-based salinity differed significantly from the measured salinity. Barr,106 et al. 29 ) found that the surface water salinities above 28 ppt result in reduced NEE and LUE of 107 mangroves. As the surface water salinity is usually lower than the sea water salinity due to the rainfall 108 and river discharge, we assumed that surface water salinity is lower than 28 ppt if the RS-based sea 109 water salinity is below 28 ppt. Therefore, a constraint was finally added to the calculation of Salscalar, 110 which is equal to 1 when the salinity is below 28 ppt. 111 PARscalar: LUE declined with the increasing PAR as shown in Figure 2. However, the decreasing rates 112 changed with the increase of PAR. When PAR values were less than 1 mmol/m 2 /s, the decreasing rate 113 was high. While when PAR exceeds 1 mmol/m 2 /s, the decreasing rate became lower. Therefore, we 114 set the threshold of PAR to get the decreasing rate of PARscalar. The mpar for PAR <= 1 mmol/m 2 /s and 115 PAR > 1 mmol/m 2 /s were 0.5171 mmol/PAR and 0.3080 mmol/PAR.

LUE validation.
119 Figure 3 shows the validation results of the terrestrial LUE model and mangrove LUE model 120 considering coastal environments, with hourly data (Figure 3a-b) and daily data (Figure 3c- Figure 3. Validation of LUE estimated from terrestrial (b&d) and mangrove model (a&c) with daily-131 (c&d) and hourly-scale (a&b) data. 132 GPP validation. 133 The validation results for the GPP model are shown in Figure 4. The results reveal that the GPP has 134 relatively high accuracies with Pearson's r around 0.5 and RMSE less than 7 μmol/m 2 /s. The GPP 135 estimated was generally lower than the measured value. Figure

166
The improved performance of the mangrove LUE model considering coastal environments in this 167 study was mainly attributed to the determination of environmental scalars. Parameters determining 168 environmental stressors (e.g., Topt, Tmin, Tmax, VPDmin, and VPDmax) were set based on the general 169 characteristics of mangroves worldwide. It may not be as accurate for the mangroves in our study sites, 170 but it generally reflects the response of mangroves to environmental changes and is applicable to other 171 study sites. Despite the specific characteristics of each mangrove ecosystem at different sites being 172 preferred, however, this study first offers the possibility to estimate mangrove productivity at a larger 173 scale to track GPP, thus emphasizing the role of mangrove ecosystems nationally or worldwide. 174 The validation results showed that the LUE values of the mangrove model agreed well with 175 those estimated by EC method (Figure 3)  were verified by field observations to be more similar to mangrove productivity near the flux towers. 227 The MODIS GPP and EC-estimated GPP showed that the MODIS GPP had a large fluctuation 228 and weakly reflects productivity, being overestimated in 2012 and underestimated in 2015. The 229 response of mangrove productivity to Tair has not been well-calibrated in the MODIS GPP product, 230 which may partly account for the poor correlation between the MODIS GPP and EC estimates. 231 However, the GPP model generated in our study showed similar trends to the field measurements, 232 capturing seasonal variations. Additionally, the increased difference in MODIS GPP versus EC 233 estimates may be due to the structure of the MODIS, which assumes a linear increase of GPP with 234 PAR. In our model, the response of GPP to PAR is suppressed, resulting in seasonal changes in GPP 235 that better match the observations. 236 Most studies provide EC-based estimates of GPP that are measurements from a limited footprint. 237 It is possible to extrapolate results across similar vegetation types and geographic settings, but not to 238 areas of heterogeneous vegetation. The RS-based GPP model offers spatial-scale estimates that can be 239 directly incorporated into ecosystem-type models. PAR, SST, and salinity are the key environmental 240 parameters of this RS-based mangrove GPP model. SST and salinity data were derived from the 241 satellite images, while PAR was generated from the reconstructed PAR data, given that is more 242 accurate than the available RS data and available for historical years of data. However, there are 243 already some RS-based PAR products from Hamawari-8, MERIS, and SeaWiFS, which provide an 244 opportunity to obtain large-scale PAR data using RS in the future. In addition to this, GPP of two 245 mangrove forests was assessed and validated with three-year measurements. Validation at different 246 sites and years showed similar results, which indicated the model has similar performance across 247 mangrove forests. Nonetheless, these estimates need to be corroborated with EC databases, which are 248 relatively accurate and provide many additional variables that are currently beyond the scope of higher 249 spatial-resolution RS estimates. The proposed GPP model considering coastal environments was well 250 suited to extend the study area by incorporating RS information and meteorological data. 251 The LUE model considering the effects of SST, salinity, and PAR performed well, however, the 252 GPP estimated from the LUE, fAPAR, and PAR showed discrepancies and were generally lower than 253 the measured values. Although the results are better than MODIS products, limitations exist still. 254 Firstly, the effects of salinity and SST on mangrove productivity were directly related to tide 255 activities. The soil pore water and surface water salinity could affect the osmotic pressure of 256 mangroves especially for the submerged parts which would control the stomatal conductance. In the 257 same way, SST could roughly represent the temperature of mangrove root systems and soil sediments 258 which has impacts on mangrove roots' respiration and transpiration. Therefore, the tide duration, tide 259 height, and tide cycle would determine the effect of salinity and SST on the mangrove LUE and GPP. 260 However, quantifying the influence from the tidal cycle remains a challenging task, which could 261 influence the performance of salinityscalar and SSTscalar. 262 Secondly, mangroves of different species and ages exhibit diverse structural and physical 263 conditions, resulting in different LUEmax, and optimal growing conditions such as Tair opt and VPDmin. 264 However, we have not specified the variables for different mangrove species and ages which could be 265 improved in the future. 266 Thirdly, the relatively low spatial and temporal resolution of the environmental data from RS 267 would influence the accuracy of the model. The datasets have a relatively coarse resolution (usually 268 500 m -1 km) and are thereby less suitable for smaller nature reserves, especially in the narrow patches 269 of mangrove areas that are rapidly being exploited in coastal China. Besides, VPD is on a monthly 270 scale, which is less consistent with in-situ measurements. Porewater salinity is controlled by sea 271 surface salinity, precipitation, and river discharge. However, currently, pore water salinity was 272 expressed in terms of sea surface salinity, which may lead to an underestimation of Salscalar. Finer 273 resolution and meteorological data are highly needed to improve the model performance more 274

significantly. 275
Lastly, the RS-derived fAPAR only considers the absorptions by living green vegetation 276 elements, whereas the ground measured fAPAR refers to the contributions from all absorbing 277 components 49 . The lower fAPAR-S2 values in mangrove forests may be due to the exposed-to-air root 278 systems which absorb the radiation. Moreover, the spatial distribution of PAR was determined by Co-279 Kriging interpolation. The elevation was taken as the covariate to estimate spatial PAR. There are 280 many other variables affecting the incoming PAR (e.g., slope and clearness) 50 . A more comprehensive 281 set of variables needs to be included in the Co-kriging interpolation to improve the PAR estimation. 282 The spatial and seasonal variation of the mangrove GPP was related to environmental changes 283 along the shoreline. The low summer GPP was explained by the lower fAPAR in summer compared 284 with other seasons, which was principally due to the underestimation of fAPAR in summer.

285
Furthermore, PARscalar took a mean value of LSP as 1 mmol/m 2 /s, however, LSP varied with different 286 species and environmental conditions. In summer, mangroves are more likely to obtain light saturation, 287 and thus PARscalar may lead to an underestimation of LUE and thus GPP. On the contrary, PAR values 288 in winter were relatively low but increased slightly with decreasing latitude. Thus, the inhibitory effect 289 of PAR on LUE was not significant, and GPP increased with decreasing latitude. Salinity and VPD 290 were more stable across years and locations and had no noticeable effect on the mangrove LUE and 291 GPP. The seasonal latitudinal patterns and effects on mangrove productivity were similar for Tair and 292 SST. Tair and SST were lower in winter, especially at high latitudes where mangroves were more 293 sensitive to cold weather. Therefore, the GPP of mangroves at high latitudes in winter was the lowest 294 throughout the year. However, hot weather in summer also limited the photosynthesis in mangroves, 295 especially at low latitudes, where Tair and SST were higher. Nevertheless, there were some correlations 296 among these environmental constants. For example, the Tair affects the vapor pressure and SST. There 297 was a positive correlation between PAR and Tair. The multicollinearity among these variables and the 298 various conditions of mangroves may affect the performance of the model and show variations along 299 the coastline, which would be improved in future studies. 300 Additionally, the GPP of mangroves increased from 2007 to 2018, which was mainly due to the 301 expansion of mangrove forests in the coastal areas. As mangroves grow, canopy size and tree density 302 increase, which may lead to higher LUE and less underestimation of fAPAR, thus contributing to high 303 productivity. However, Zhejiang province (27°02′N-31°11′N) experienced extremely cold weather in 304 January 2016 caused by the East Asia cold wave 51,52 , and large areas of mangrove forests died or 305 became sick, leading to a decline in the mangrove GPP at high latitudes in 2018. 306

307
In conclusion, we presented a RS-based productivity model to estimate the GPP of mangrove forests 308 in China. The model considered the environmental stressors induced by tidal inundation, therefore, 309 involving SST, sea surface salinity, and PAR as environmental scalars to develop the LUE model. SST 310 was first-ever included in the mangrove LUE model and parameterized by a similar model for Tair. In 311 addition, it was the first to indicate the downscaling effects of PAR on the mangrove LUE and 312 determine the LUEmax for mangroves under different temperatures. Consequently, the mangrove GPP 313 was estimated based on the mangrove LUE model, fAPAR generated from Sentinel-2 images and 314 reconstructed PAR from meteorological stations. The results revealed that PAR, Tair, VPD, SST, and 315 salinity are clearly drivers of diurnal and seasonal variations in the mangrove LUE and CO2 fluxes. 316 Among them, PAR, SST, and salinity are unique to mangrove ecosystems. The LUE model developed 317 for mangrove forests had higher overall accuracy (RMSE = 0.0051, R 2 = 0.64) than the previous LUE 318 model for terrestrial forests (RMSE = 0.0220, R 2 = 0.24). GPP estimated in this study generally agreed 319 with in-situ measurements from two carbon flux towers. Although there are still limitations, the 320 modeled GPP maintained higher accuracies compared with MODIS GPP products. These results 321 demonstrated the potential of RS-driven productivity models for the large-scale mangrove GPP 322 estimation and provided fundamental data and scientific methodological support for future mangrove 323 blue carbon potential assessment and restoration policy development. 324

325
Study area. 326 Two carbon flux towers have been established in Zhangjiang Estuary Mangrove National Nature 327 Reserve (117°24'53.02"E, 23°55'26.63"N) in Fujian and Zhanjiang Mangrove National Nature 328 Reserve (110°09'44.67"E, 20°56'24.08"N) in Guangdong (Fig. S4), China. The mangroves in 329 Zhangjiang and Zhanjiang are mainly composed of Kandelia obovate and Sonneratia apetala. The 330 forest structures and microclimate in these two sites were various and listed in Table S3. 331 Site-specific data from carbon flux towers 332 Half-hourly carbon fluxes between the canopy and atmosphere were obtained from flux towers and 333 processed by the EC method. Meteorological and tidal information was measured with multiple 334 instruments near the flux tower. All data were provided by the ChinaFlux network 335 ( http://www.chinaflux.org ). Table S4 summarizes the data availability. More details of the EC system 336 structure and data processing can be referred to in published papers 17,30 . 337 Remote sensing data. 338 Historical climate data were derived from different satellites or based on the reanalysis data. The 339 summary of the dataset can be found in Table 2. Climate data were obtained from the Google Earth 340 Engine platform. Sentinel-2 L1C images were adopted for computing the fraction of absorbed 341 photosynthetic active radiation (fAPAR). PAR was derived from the reconstructed PAR dataset 53 342 which was derived from the meteorological data, MODIS AOD data and NASA/GSFC O3 data. Zheng  for Tair and SST. The low temperature may freeze the water-conducting xylem vessels of the mangrove, 380 and the high temperature would reduce the stomatal conductance 58 . So, we assumed that the mangrove 381 LUE increases with the increase of temperature, however, it will start to decrease after a certain value. 382 Therefore, Tscalar was defined as equation (4) where Tmin, Tmax, and Topt are minimum, maximum, and optimal temperatures for the mangrove 385 photosynthetic activities, respectively. If Tair is below Tmin, Tscalar is set to be zero. The daily mean 386 temperature (Tmean) and daily maximum temperature (Tmax) were employed to calculate the daytime 387 Tair following equation (5)  where VPDmax and VPDmin are the maximum and minimum daytime VPD. If VPD is less than VPDmin, 397 VPDscalar is set to be 1. If VPD is larger than VPDmax, VPDscalar will be set as 0. We determined VPDmin 398 and VPDmax by summarizing from previous studies and verified with in-situ data. 399 Salinityscalar: The decline in LUE with increasing salinity was quantified by equation (7)  Gross primary production (GPP) modeling. GPP was calculated as equation (9) 4,64 and the overall 408 flowchart can be summarized as Figure 8: Firstly, the LUE of mangroves in these two mangrove reserves was calculated based on the LUE model 410 proposed in section 2.4.1. Resampling of RS data was carried out to keep the spatial resolution at 500 411 m and the temporal resolution at daily. Since SST and salinity were derived from the sea surface, the 412 nearest sea surface pixel to the mangrove was adopted to represent the effects of SST and salinity on 413 mangrove LUE. Then, fAPAR was computed using the biophysical processor in SNAP software 65 . The 417 processed Sentinel-2 fAPAR products (fAPAR-S2) represent the daily integrated fAPAR values, 418 following the assumption that the instantaneous fAPAR value at 10:00 (or 14:00) solar time is close 419 to the daily integrated value under clear sky conditions 66 . Besides, outlier pixels in fAPAR-S2 were 420 eliminated and only pixels with "QA = 0 0 0" were adopted. 421 After that, reconstructed PAR data were obtained from 724 meteorological stations provided by 422 Tang, et al. 53 ). We further interpolated the PAR data from meteorological stations to the whole coastal 423 zone by the Co-Kriging interpolation method 67 taking surface elevation as covariate 68,69 . Finally, GPP 424 in these two mangrove reserves was estimated based on the derived LUE, fAPAR, and PAR. 425 Model validation and application 426 The LUE results modeled with hourly and daily environmental data were validated with the LUE 427 values from the carbon fluxes tower in the Zhangjiang mangrove reserve. Moreover, LUE was 428 estimated following the LUE model for terrestrial forests considering only the effects of Tair and VPD 429 as shown in Eq.
(1). The experimental results were also compared with in-situ LUE to evaluate the 430 performance. GPP estimated using the proposed model was validated with the flux tower 431 measurements. The GPP derived considering coastal environments was in turn converted to a 432 cumulative 8-day composite and compared with MODIS GPP product at the same resolution 70 . 433 After validation, the GPP model was applied to estimate the GPP of the mangrove forests in the 434 whole coastal zone of China for the years 2007, 2010, and 2018. Besides, seasonal variations were 435 displayed to reflect the different productivity of mangroves under various environmental conditions. 436 Data availability 437 The datasets generated during and/or analyzed during the current study are available from the 438 corresponding author on reasonable request. 439