Spatiotemporal dynamics and geo-environmental factors influencing mangrove gross primary productivity during 2000-2020 in Gaoqiao Mangrove Reserve, China

Dmi Zhao, Yinghui Zhang, Junji Wang, Jianing Zhn, Zhn Shn,Kunlun Xiang, Haoli Xiang, Yongquan Wang, Guofng Wu,＊＊

a MNR Key Laboratory for Geo-Environmental Monitoring of Great Bay Area & Guangdong Key Laboratory of Urban Informatics & Shenzhen Key Laboratory of Spatial Smart Sensing and Services, Shenzhen University, Shenzhen, 518060, China

b School of Architecture and Urban Planning, Shenzhen University, Shenzhen, 518060, China

c College of Life Sciences and Oceanography, Shenzhen University, Shenzhen, 518060, China

d Key Laboratory of Wetland Ecology and Environment, Northeast Institute of Geography and Agroecology, Chinese Academy of Sciences, Changchun, 130102, China

e Guangdong Ecological Meteorology Center, Guangzhou, 510275, China

Keywords:Mangrove GPP LUE model Geodetector MGWR Spatial heterogeneity

ABSTRACT

1.Introduction

Global warming caused by carbon emissions has become a serious threat to human well-being and socioeconomic development(Raza et al.,2021).To actively respond to the trend of global warming, several countries have adopted carbon neutral policies based on “carbon reduction” (Chen et al., 2022).For example, China announced climate targets to achieve carbon peaking by 2030 and carbon neutrality by 2060(Zeng et al.,2022).

Mangrove ecosystems provide abundant ecosystem service functions(e.g.carbon sequestration, water purification and biodiversity maintenance)(Li and Lee,1997).Mangroves are one of the major“blue carbon”ecosystems in the coastal tropical region (Lee et al., 2014), with high gross primary productivity (GPP) and low soil organic carbon decomposition rate(Lovelock et al.,2015),and they store approximately 24 Tg C·y-1organic carbon around the world (Alongi, 2014).Moreover, the carbon cycle of mangrove ecosystems is sensitive to the rising concentration of CO2in the atmosphere (Liu et al., 2012; Lovelock and Reef,2020).Therefore,effectively assessing the carbon sequestration capacity of mangrove ecosystem can provide a theoretical basis for the policy formulation of carbon emission reduction and carbon neutrality in coastal areas(Sazeides et al.,2021).

GPP is defined as the total carbon assimilated through photosynthesis,which can reflect the capacity of carbon capture and migration in plants(Richardson et al.,2010).The accurate estimation of GPP plays an important role in monitoring and assessing plant growth,carbon balance and carbon conversion (Zhu et al., 2021a).In addition, GPP products have far-reaching value for monitoring the terrestrial carbon cycle and constraining the intensity of carbon sources/sinks in the Earth system(Zhu and Yan, 2022).Therefore, monitoring the spatiotemporal dynamics of GPP can capture regional changes in the carbon budget and provide a scientific basis for developing policies to address climate adaptation and carbon neutrality goals.

Currently, the main methods of estimating forest GPP include direct measurements, flux tower methods and remote sensing models (Zhang et al., 2014; Lovelock et al., 2015).The direct measurement method measures leaf photosynthesis in the laboratory; however, it cannot quantify forest GPP in large areas.The flux tower method measures the CO2exchange between ecosystems and atmosphere using eddy correlation (EC) techniques, and it requires continuous observation data; however,the current number of active flux tower sites is limited with sparse spatial distribution across the globe(Hilker et al.,2014).Estimating GPP with remote sensing data relies on the light-use efficiency (LUE) of different environmental conditions and ecosystem structures (Mu et al.,2007; Zhang et al., 2023), and global or regional GPP estimations can also be obtained from site-level GPP measurements by spatial extrapolation using remote sensing models (Huang et al.,2019).

LUE connects GPP with the photosynthetically active radiation(PAR)(Myneni et al., 1995), fraction of absorbed PAR (fAPAR) (Myneni et al.,2002) and conversion efficiency ε (Turner et al., 2003b).Several LUE models have been widely adopted to estimate global terrestrial forest GPP in the last two decades, such as the Carnegie-Ames-Stanford-Approach(CASA) model, Vegetation Photosynthesis Model (VPM), C-Fix model and MOD17 model (Pei et al., 2020; Yang et al., 2014).Although there exist several models for estimating GPP in vegetated forest ecosystems,these models have limited applicability to mangrove ecosystems with distinct tidal characteristics and high salinity levels (Barik et al., 2018).Therefore,to achieve accurate estimation of mangrove GPP,it is imperative to consider additional limiting factors such as salinity and seawater temperature (Zheng and Takeuchi, 2022).Moreover, the estimation of mangrove GPP on a large scale has been limited, mainly owing to the challenges associated with acquiring flux tower measurements and comprehensively understanding the intricate carbon exchange processes in mangrove ecosystems.

In the LUE models,GPP is affected by photosynthesis,solar radiation(Schloss et al., 1999), temperature (Zhang et al., 2009), precipitation(Zhang et al.,2009),soil properties(salinity,pH,clay and Eh)(Matin and Jalali, 2017) and forest vegetation structure, such as leaf area index(LAI), canopy height, stand basal area and crown diameter (Braghiere et al.,2019).LAI is a critical biophysical factor and is sensitive to climate change (Kovacs et al., 2005), affecting the overall GPP changes and carbon assimilation rate in mangrove (Ishtiaque et al., 2016).Although many factors(e.g.climatic factors)have been identified to influence the spatial variability of GPP, few studies have quantitatively explored the impact of geo-environmental factors (e.g.species distribution, offshore distance, elevation, slope) on time-series estimations of mangrove GPP.In recent years, Geodetector (Javi et al., 2014) and multi-scale geographically weighted regressions (MGWR) (Brunsdon et al., 1998)have been used to quantitatively investigate the spatiotemporal relationship between independent and dependent variables.However, very few studies have combined the two methods to investigate the effects of different geo-environmental factors on spatiotemporal estimations of mangrove GPP.

This study aimed to examine the spatiotemporal dynamics of mangrove forest GPP and its influencing factors in the Zhanjiang Mangrove National Nature Reserve from 2000 to 2020.We analyzed the spatiotemporal evolution of mangrove forests during 2000-2020 and investigated the spatiotemporal dynamics of GPP using the LAI-based LUE model and Google Earth Engine (GEE) cloud computation.Moreover,GPP variation patterns were analyzed using both the Mann-Kendall(MK) test and Theil-Sen median trend, and further explored the individual and interactive contributions of multiple geo-environmental factors (species distribution, offshore distance, elevation, slope, planar curvature and profile curvature) to GPP variation through both Geodetector and MGWR.

2.Methods

2.1.Study area

The study area (Gaoqiao Mangrove Reserve) is located in the core region of Zhanjiang Mangrove National Nature Reserve (21°31′-21°35′N, 109°44′-109°48′E), Guangdong Province, China (Fig.1).The main mangrove species in this region are Bruguiera gymnorrhiza(L.)Lam(BG),

Aegiceras corniculatum (L.).Blanco (AC), Rhizophora stylosa Griff (RS),Kandelia candel (L.).Druce (KC) and Avicennia marina (Forssk.) Vierh.(AM), the mangrove forests are mainly natural secondary forests and artificial forests,and only a few areas have primary forests.

The study area has a southern subtropical monsoonal maritime climate with mean annual temperature, precipitation and photosynthetically active radiation (PAR) values of approximately 22.8°C,1,700-1,800 mm and 260 μmol photons·m-2·s-1,respectively(Yu et al.,2008).The whole mangrove buffer area (yellow region in Fig.1) is dominated by mangrove forest, terrestrial forest, aquaculture zone,farming zone,building zone and water body.

2.2.Data collection

2.2.1.Field survey

The field survey was carried out on August 6-12,2020.A total of 64 sampling plots(15 m×15 m)containing a single species were randomly selected with a distance of at least 30 m between any two plots.At the center of each plot, LAI was obtained using the Hemiview digital plant canopy analysis system(Delta-T Devices,Cambridge),which uses a highdefinition digital camera (Coolpix 950, Nikon Corporation) and a 180°fisheye lens (FC-E8, Nikon Corporation) to obtain bottom-up hemispherical images(Zhao et al.,2022).

2.2.2.Carbon flux measurements

Fig.1.The true-color images of Landsat-8 images(date:10/23/2020)showing the location of study area with a 2-km buffer area.The mangrove area was divided into four zones(A,B,C and D)based on the river flow direction.(For interpretation of the references to color in this figure legend,the reader is referred to the Web version of this article.)

The carbon flux(109°45′22.32′′E,21°34′3.05′′N)and microclimate measurements of the reserve from 12/31/2012 to 12/31/2013 were obtained from the China Flux network(http://www.chinaflux.org).The ecosystem CO2flux data measured by the EC technique were analyzed following the processing steps provided by Cui et al.(2018)and Lu et al.(2019), and mainly included the GPP, air temperature (TT,°C), vapor pressure defciit(VPD,hPa)and photosynthetically active radiation(PAR,mol·m-2·d-1).

2.2.3.Remote sensing data

GEE holds freely available geospatial data and remote sensing with high-performance parallel computing services (Gorelick et al., 2017).In this study, to efficiently obtain batching geospatial data and minimize the technical difficulty of cloud detection with different Landsat sensors(e.g.Landsat-5 TM,Landsat-7 ETM+and Landsat-8 OLI/TIRS),GEE was used for data loading and processing.To avoid the effects of seasonal phenology and cloud cover on result analyses, we chose six optimal Landsat images with low cloud cover in the post-monsoon periods of 2000, 2004, 2008, 2012, 2016 and 2020.Moreover, the GPP product datasets comprise MODIS GPP (https://modis.ornl.gov/MCD15A2H), GLASS GPP (http://www.glass.umd.edu/Download.h tml) and GEOSIF GPP (https://globalecology.unh.edu/data/GOSIF-GPP.html).

2.2.4.Hydroclimatic data

Hydroclimatic data have an important impact on regional GPP (Yan et al., 2021).The temperature (TT) and precipitation (PPT) during 2000-2020 were obtained from the National Meteorological Science Data Center (http://data.cma.cn/).To better reflect the unique micrometeorology of mangroves, we fnially selected data from a site located close to the study area in the town of Gaoqiao, Lianjiang (area code:58658).The SRAD data (M2T1NXRAD) during 2000-2020 were extracted from NASA Earthdata Search (https://search.earthdata.nasa.gov/).In addition, ocean surface temperature (SST) data sources were from MODIS Aqua data (https://oceancolor.gsfc.nasa.gov/l3/); salinity data were from the hybrid coordinate ocean model(HYCOM).

2.3.Land use classification and landscape pattern analysis

Three commonly-used vegetation indices(VIs)derived from Landsat-5 TM,7 ETM+and 8 OLI/TIRS,including Normalized Difference Vegetation Index(NDVI),Modified Normalized Difference Water Index(MNDWI)and Normalized Difference Built-up Index (NDBI) (Table S1),were employed for land use classification, because they are sensitive to changes in vegetation,water and buildings(Arjasakusuma et al.,2020;Jain et al.,2020).Random forest(RF)classification using the three VIs was then employed to extract the spatial distribution of six land use types from the study area with a buffer zone of 2 km (Fig.1): mangrove forest, terrestrial forest,aquaculture zone,farming zone,building zone and water body.

We collected 1,070 samples from Sentinel-2 images with 10-m resolution using GEE, including mangrove forest (160 samples), terrestrial forest(210 samples),aquaculture zone(190 samples),farming zone(220 samples),building zone(210 samples)and water body(80 samples).We only investigated the spatiotemporal changes every four years(2000-2004, 2004-2008, 2008-2012, 2012-2020 and 2016-2020) to outline the changes in mangrove area.

The landscape-level metrics aggregate all patches in the landscape and indicate the fragmentation processes of different characteristics of spatial patterns (Koomen et al., 2012; Masnavi et al., 2021).Hence, we chose four commonly-used landscape metrics (Table S2 in the Supplementary Materials) at the landscape level, including aggregation index(AI), largest patch index (LPI), patch density (PD) and Shannon's diversity index (SHDI), to measure landscape properties, shape, aggregation extent, diversity and patch contagion.Moreover, we described the details of the above methods in the Supplementary Materials.

2.4.GPP estimation

2.4.1.LUE model

In this study, GPP was calculated based on the commonly-used LUE model by Xiao et al.(2004).Moreover,to enhance the predictive power and robustness of the LUE model, two limiting factors (surface temperature(SST)and salinity)were taken into account(Zheng and Takeuchi,2022):

where PAR is photosynthetically active radiation,fAPAR is the fraction of absorbed PAR (W·m-2) by leaf chlorophyll.εmax(g C·MJ-2) is the maximum LUE, f（Tair）, f（TSST） and f（w） are the environmental stress indices of temperature and water,respectively.

where SWRadis incident shortwave radiation;LAI(leaf area index)is the basic parameter for canopy photosynthesis and evapotranspiration(Myneni et al., 1997; Poulter et al., 2009); and the value of k (light extinction coefficient)is 0.5.

where T,Tmax,Tmin,and Toptare daily mean,maximum,minimum, and optimal air temperature for the mangrove photosynthetic activities,respectively.

Moreover, the photosynthetic response of mangroves to air temperature is similar to sea surface temperature (SST).The SST, SSTmax,SSTmin,and SSToptare daily mean,maximum,minimum,and optimal SST for the mangrove photosynthetic activities,respectively.

where LSWImaxis the maximum LSWI(Land Surface Water Index).

where msaldenotes the rate at which salinity decreases as it increases,and the value is set to 0.0047.

2.4.2.Parameter calibration

The main input parameters(εmax,Tmin,Tmax,Topt,SSTmin,SSTmaxand SSTopt) of the LUE models are generally subject to some uncertainties across different vegetation species, ages, and seasons.Optimizing these parameters is key to improving the accuracy of the model in quantifying vegetation GPP (Zhao et al., 2005; You et al., 2019).In this study, we used the Levenberg-Marquardt (LM) algorithm as a nonlinear least squares method for parameter optimization.This method locates the minimum of a multivariate function by utilizing both the gradient and the Gauss-Newton methods(Ma and Wang,2009).Previous research has demonstrated that LM has great advantages and scalability when applied to the correction and optimization of parameters for LUE (Fan and Pan,2006; Ma and Wang, 2009).Based on the LM model and flux data, we calibrated the key parameters:εmax,Tmin,Tmax,Topt,SSTmin,SSTmaxand SSTopt,which were 3.34 g C·MJ-1,10.99°C,33.19°C,27.28°C,9.40°C,31.48°C and 26.07°C,respectively.

2.4.3.GPP trend analysis and validation

According to the results of Narula and Wellington(2007),we used the multiple linear regression (MLR) model with the three VIs (Table S1,NDVI,EVI(Enhanced Vegetation Index)and GNDVI(Green Normalized Difference Vegetation Index)) derived from Landsat-8 images acquired on August 23, 2020 to estimate mangrove LAI.A total of 64 sampling plots (15 m × 15 m) were randomly selected, and a HemiView digital plant canopy analysis system (Delta-T Devices, Cambridge, UK) was employed to measure their LAI values.Six outliers were eliminated in the modeling process considering the error value, and the remaining 58 samples were randomly divided into training (60% samples) and validation subsets(40%samples).The training subset was used to train the MLR model,while the validation subset was used to validate the model.The coefficient of determination(R2)and root mean square error(RMSE)were used to evaluate the performance of the MLR model.The MLR model for LAI estimation based on the Landsat-8 image in 2020 was employed to map the spatial distribution of LAI from 2000 to 2020 with Landsat-5 TM,Landsat-7 ETM+and Landsat-8 OLI/TIRS images.

The GPP was estimated and mapped for each period using the LUE model.The trends and interannual variability of GPP were assessed using the Mann-Kendall (MK) test and Theil-Sen median trend during 2000-2020.The Theil-Sen median is a robust nonparametric statistical trend calculation method,and it is suitable for trend analysis of long time series data because it is computationally insensitive and efficient to measurement errors and niche data (Akritas et al., 1995).The Mann-Kendall test is a nonparametric time series trend test, and it is independent of missing values and outliers(Hamed and Rao,1998),which is combined with the Theil-Sen slope estimation of GPP trend lines for long time series data.

To verify the reliability of the estimated mangrove GPP,we compared the differences between the estimated GPP from the optimized LUE model and the measured GPP from the flux tower for 2012 and 2013.Furthermore,we also compared the estimated GPP within a 500-m range of the flux tower (109°45′22.32′′E, 21°34′3.05′′N) and GPP products(MODIS GPP, GLASS GPP, and GEOSIF GPP) with eight-day intervals from 2000 to 2020.

2.5.The relationship and effects of geo-environmental factors on GPP variation

Previous studies found that geo-environmental factors, such as temperature and precipitation(Zhang et al.,2009),soil water(Gebremichael and Barros,2006),photosynthetic active radiation(Myneni et al.,1995),vegetation index(NDVI,EVI)(Rodda et al.,2022)and land cover changes(Sun et al., 2018), controlled the variations in CO2fluxes.In this study,we investigated the relationship between two VIs (NDVI and EVI) and four landscape metrics(AI,LPI,PD and SHDI)against GPP.Furthermore,six geo-environmental factors (species distribution, offshore distance,elevation,slope,planar curvature and profile curvature)were employed to quantitatively explore their contributions to mangrove GPP variation by Geodetector and MGWR.In addition, the detailed description of Geodetector and MGWR methods can be found in the Supplementary Materials.

3.Results

3.1.Spatiotemporal evolution of mangrove forests during 2000-2020

3.1.1.Area changes

The RF classification model used for land use classification reported an overall classification accuracy of 85.34% (Kappa coefficient: 0.810),85.15% (Kappa coefficient: 0.792), 92.34% (Kappa coefficient: 0.882),85.46% (Kappa coefficient: 0.817), 91.70% (Kappa coefficient: 0.886)and 91.29%(Kappa coefficient:0.900)in 2000,2004,2008,2012,2016 and 2020,respectively(Fig.2).The mangrove forest experienced an area loss from 766.26 ha in 2000 to 718.29 ha in 2020 (Fig.S1 in the Supplementary Materials).The maximum reduction rate of mangrove forest was observed during 2016-2020 (13.62%), followed by the periods of 2008-2012(7.25%)and 2012-2016(6.36%),and the loss of mangrove forest was mainly converted to water bodies and aquaculture zones(Fig.S2).A significant increase in mangrove forest was found during 2000-2008 with a growth rate of 24.93%, mainly converted from aquaculture zone and farming zone with an area of 110.88 and 34.08 ha,respectively(Fig.S2).The area of buildings increased greatly from 8.91 ha in 2000 to 344.79 ha in 2020 with a growth rate of 3,769.70%,mainly converted from farming zones, terrestrial forests and aquaculture zones with areas of 244, 74.08 and 48.64 ha, respectively.The area of aquaculture zone raised from 1,660.86 to 1,840.41 ha with a growth rate of 10.81%.However, the area of the terrestrial forest and farming zones decreased at rates of 26.32% and 11.44% during 2000-2020,respectively.

3.1.2.Landscape pattern changes

The spatial distributions of the four landscape metrics(PD,LPI,AI and SHDI)at the landscape level were mapped based on 1,917 grids of 200 m×200 m(Figs.S3.1 and S3.2).Considering the whole region,the mean AI and LPI values increased from 84.37%to 73.40%in 2000 to 88.35%and 76.17%in 2020,respectively,and the mean PD and SHDI values declined from 100.21 to 0.57 in 2000 to 85.82 and 0.51 in 2020, respectively.Considering mangrove forests,the mean AI and LPI values dropped from 90.27%to 80.83%in 2000 to 89.57%and 75.60%in 2020,respectively,and the mean PD and SHDI values increased from 75.05 to 0.41 in 2000 to 84.26 and 0.52 in 2020,respectively.

The highest values of AI(Fig.S3.1)and LPI(Fig.S3.2)in each period were mainly observed in the water body and aquaculture zone,while the lowest values were found in the building zone.However, the highest values of PD (Fig.S3.3) and SHDI (Fig.S3.4) were observed in the building zone near the buffer zone of mangrove forest,while the lowest values were found in the water body and aquaculture zone.Moreover,the four landscape metrics had large differences in the four mangrove zones (A, B, C and D in Fig.1).The highest mean values of AI and LPI were found in Zone C,followed by in Zone D,Zone B and Zone A,and the highest mean values of PD and SHDI were found in Zone A,followed by Zone B,Zone D and Zone C.

Fig.2.Spatial distribution of six land use types within the 2-km buffer zone of the reserve from 2000 to 2020.

3.2.Mangrove GPP estimation and validation

3.2.1.Mangrove LAI mapping

Based on the Landsat-8 and field-measured LAI in 2020, a multiple linear regression model with NDVI, GNDVI and EVI (LAI = 1.489 +6.83NDVI --3.323GNDVI - 0.728EVI) was developed to estimate LAI,and it was applied to map LAI from 2000 to 2020(Fig.3).

The mangrove forest reported increasing LAI values with a mean LAI of 2.47, 2.99, 3.04, 3.06, 3.57 and 3.98 m2·m-2in 2000, 2004, 2008,2012,2016 and 2020,respectively.The mangrove forest in Zone A and C held higher LAI values than that in Zone B and D during 2000-2016;however,the mangrove forest in Zone A and B had higher LAI values than that in Zone C and D in 2020.The mangrove forest in Zone B showed the highest increase in mean LAI(134.88%),followed by Zone D(127.48%),Zone C(76.57%)and Zone A (50.91%)from 2000 to 2020.

3.2.2.Mangrove GPP mapping

The total GPP of mangrove forests in the post-monsoon period(September-November) was 45,611.41 (mean value = 6.35 g C·m-2·d-1),62,359.15 (mean value = 6.74 g C·m-2·d-1), 63,088.05(mean value =6.59 g C·m-2·d-1), 59,491.98 (mean value = 6.70 g C·m-2·d-1),67,768.06 (mean value = 8.15 g C·m-2·d-1) and 59,833.56 kg C·m-2·d-1(mean value=8.33 g C·m-2·d-1)in 2000,2004,2008,2012,2016 and 2020,respectively(Fig.4).When considering the period from 2000 to 2008, Zone A had a faster growth rate (10.12%) and a higher mean GPP compared to the other three zones,while Zone B had a slower growth rate(0.66%)and a lower mean GPP.However,when considering the period from 2012 to 2020, Zone A exhibited a slower growth rate(14.54%)despite having a higher mean GPP,whereas Zone D had a faster growth rate(32.10%)but a lower mean GPP than the other three zones.

When comparing GPP trends during 2000-2020, we observed that only 1.27%of mangrove forests showed a decreasing trend(Trend ＜0),while the majority of other areas exhibited an increasing trend(Fig.5a).Specifically, 75.52% and 23.21% of areas showed slightly increasing(0-0.50) and highly increasing (0.50-1.17) GPP trends, respectively.Among the four zones, Zone B had the most areas (44.97%, mainly located at river edges)with trend values greater than 0.50,followed by Zone D(24.07%),Zone A(14.42%),and Zone C(6.58%).The analysis of variance of GPP values across the six periods demonstrated that 47.35%of the entire reserve exhibited significant changes in GPP (p ＜0.05),predominantly found at the river edges of Zones B and D as well as in the outer edges of the mangrove(Fig.5b).

3.2.3.Validation of estimated GPP with eddy covariance GPP

Fig.3.Distribution map of the leaf area index (LAI, unit: m2·m-2) in the mangrove forest from 2000, 2004, 2008, 2012, 2016 and 2020.

Fig.4.The distribution map of mangrove gross primary productivity (GPP, unit: g C·m-2·d-1) in 2000, 2004, 2008, 2012, 2016 and 2020.

Fig.5.Spatial distribution (a) and significant changes (b) of gross primary productivity trends in the mangrove forest from 2000 to 2020.

Fig.6.The scatter plot of estimated GPP against eddy covariance (EC) GPP in (a) 2012 and (b) 2013.

Fig.7.The daily-scale spatial data at 8-day intervals of estimated mangrove GPP and GPP products(MODIS,GLASS and GEOSIF)in(a)2000,(b)2004,(c)2008,(d)2012, (e) 2016 and (f) 2020.

The accuracy validation results (Fig.6) between the estimated GPP and eddy covariance GPP demonstrated that the scatter plot was generally distributed close to the 1:1 line with the mean R2and RMSE value of 0.68 and 1.30 g C·m-2·d-1, respectively.Furthermore, the slope of the best-fit line was close to 1 in 2013 (R2= 0.74, RMSE = 1.06 g C·m-2·d-1), which indicates relatively consistent and accurate agreement between the estimated GPP and observed GPP from the eddy covariance data.

3.2.4.Comparison of estimated GPP with GPP products

The annual mean estimated mangrove GPP within a 500 m range of the flux tower(109°45′22.32′′E,21°34′3.05′′N)was 1041.84,1187.12,1241.84,1308.64,1533.28 and 1663.36 g C·m-2·y-1for the years 2000,2004, 2008, 2012, 2016 and 2020, respectively (Fig.7).However, the annual mean value of the three GPP products, GEOSIF GPP (1348.84,1313.12, 1485.44, 1662.64, 1789.36, 1814.93 g C·m-2·y-1) had the highest value, followed by GLASS GPP (1073.12, 1045.61, 1125.76,1433.84, 1571.92, 1771.68 g C·m-2·y-1) and MODIS GPP (1127.61,1073.44,1146.08,1300.81,1340.96,1368.72 g C·m-2·y-1)for the same years.

Moreover,the estimated ranges of GPP was 0.32-4.67(mean value=2.82 g C·m-2·d-1), 0.61-5.20 (mean value = 3.22 g C·m-2·d-1),0.42-5.25 (mean value = 3.37 g C·m-2·d-1), 0.42-5.66 (mean value =3.56 g C·m-2·d-1), 0.45-5.85 (mean value = 4.17 g C·m-2·d-1) and 0.81-6.15 g C·m-2·d-1(mean value=4.52 g C·m-2·d-1)for 2000,2004,2008,2012,2016 and 2020,respectively,using daily scale spatial data at 8-day intervals(Fig.7).When compared to the estimated mangrove GPP,the GEOSIF GPP products were found to have the highest overestimation with mean values ranging from 3.99 to 4.93 g C·m-2·d-1, followed by GLASS GPP (mean value = 2.85-4.81 g C·m-2·d-1) during 2000-2020,while MODIS GPP (mean value = 3.11-3.71 g C·m-2·d-1) had an underestimation.

Fig.8 illustrates that the fitted accuracy between the estimated GPP and MODIS GPP(mean R2=0.58)was the highest,followed by GEOSFI GPP (mean R2= 0.43) and GLASS GPP (mean R2= 0.31).During the 2000-2020, the estimation accuracy was highest in 2004 (mean R2=0.63),followed by 2012(mean R2=0.54),2008(mean R2=0.45),2016(mean R2=0.41),2020(mean R2=0.35)and 2000(mean R2=0.27).

Fig.8.The scatter plot of estimated mangrove GPP against MODIS GPP,GEOSIF GPP and GLASS GPP in(a)2000,(b)2004,(c)2008,(d)2012,(e)2016 and(f)2020.

3.3.Relationships of VIs and landscape metrics against GPP

The four landscape metrics(AI,LPI,PD and SHDI,Figs.S4-S7 in the Supplementary Materials)had much weaker correlations with GPP than the two VIs(NDVI and EVI)in all periods(Figs.9.1 and 9.2).GPP had a slightly stronger relationship with NDVI(r=0.75-0.97,mean r=0.87,p＜0.05)than EVI(r=0.58-0.94,mean r=0.82,p ＜0.05)from 2000 to 2020.Moreover,SHDI,PD,LPI and AI showed weak but significant(p ＜0.05) correlations with GPP (mean = 0.26, 0.22, 0.21 and 0.19,respectively).

Among the six periods during 2000-2020(2000, 2004, 2008,2012,2016,2020),the strongest and weakest correlations of NDVI against GPP were observed in 2008(r=0.97,p ＜0.05)and 2012(r=0.75,p ＜0.05),and the strongest and weakest correlations of EVI against GPP were found in 2008(r=0.94,p ＜0.05)and 2016(r=0.58,p ＜0.05),respectively.Moreover,there was an obvious increasing trend in the correlation of the four landscape metrics against GPP from 2000 to 2004 with r of 0.24-0.28, 0.34-0.39, 0.34-0.39, 0.39-0.43 for AI, LPI, PD and SHDI,respectively.However, there was a relatively stable and low correlation between the four landscape metrics and GPP from 2008 to 2020.Apparently, due to the urbanization and land conversion with an increasingly complex environment around the mangrove forest, the correlation between the six geo-environmental factors and GPP gradually became weaker(mean|r| =0.51-0.34)from 2000 to 2020.

Among the four mangrove zones(Table 1),Zone C(mean|r|=0.42)had the highest correlation of NDVI, EVI, AI, LPI, PD and SHDI against GPP,followed by Zone A(mean|r|=0.41),Zone D(mean|r|=0.39)and Zone B(mean|r|=0.37).The strongest correlation of the two VIs against GPP was found in 2008 with mean r of 0.97,0.96,0.97 and 0.97 for Zone A,B,C and D,respectively,and the weakest correlation was observed in 2020 with mean r of 0.71, 0.79, 0.64 and 0.60 for Zone A, B, C and D,respectively.The strongest correlation of four landscape metrics against GPP was found in 2000 with mean |r| of 0.37, 0.20, 0.32 and 0.32 for Zone A, B, C and D, respectively, and the weakest correlation was observed in 2020 with mean|r|of 0.06,0.04,0.22 and 0.02 for Zone A,B, C and D,respectively.

3.4.Significant geo-environmental factors related to GPP variation using Geodetector

With the factor detector analysis of Geodetector, we found that the mean q values (Table 2) of species distribution, offshore distance,elevation, slope, planar curvature and profile curvature from 2000 to 2020 ranged from 0.031 to 0.228,0.004 to 0.043,0.010 to 0.193,0.001 to 0.008,0.004 to 0.010,and 0.001 to 0.008,respectively.Based on the mean q values over six periods,species distribution(0.112)and elevation(0.102) showed the highest relation with GPP variation, followed by offshore distance (0.016), planar curvature (0.007), profile curvature(0.004) and slope (0.003).Compared with other geo-environmental factors, profile curvature and slope showed no significant (p ＞0.05)relation with GPP variation during 2000-2020.

Compared with the factor detector results (mean q value = 0.041)from 2000 to 2020,the impacts of the interaction detector(mean q value= 0.077) between any two factors (i.e., species distribution, offshore distance and elevation) were enhanced nonlinearly by 23.39%.Among all the combination of any two factors, the interaction of species distribution and elevation had the highest contribution to GPP variation(Fig.10)with q values of 0.130,0.185,0.170,0.372,0.257 and 0.159 in 2000,2004,2008,2012,2016 and 2020,respectively.The interaction of species distribution(mean q value=0.146)or elevation(mean q value=0.134) with the other five factors also reported relatively high explanatory power, followed by the interaction of offshore distance (mean q value = 0.065) and planar curvature (mean q value = 0.052) with the other five factors.Low q values were observed for the interaction among profile curvature and slope.

Fig.9.1.Relationships between the estimated GPP and NDVI of mangrove forests in (a) 2000, (b) 2004, (c) 2008, (d) 2012, (e) 2016 and (f) 2020.

Fig.9.2.Relationships between the estimated GPP and EVI of mangrove forests in (a) 2000, (b) 2004, (c) 2008, (d) 2012, (e) 2016 and (f) 2020.

3.5.Spatial difference of influencing factors on GPP variation using MGWR

Taking into account the results of both factor detector analysis(Table 2)and interaction detector analysis of GPP variation(Fig.10),we selected three geo-environmental factors (species distribution, offshore distance and elevation) that showed significant relations with GPP variation for further MGWR analysis to explore their strength and direction in the spatial-temporal evolution of GPP.The spatial autocorrelation analysis results indicated significant spatial autocorrelation in the distribution of GPP, with global Moran's I values of 0.653,0.699,0.673,0.836,0.407,and 0.372(p ＜0.001)for 2000,2004,2008,2012,2016 and 2020,respectively.

The standardized regression coefficients of the three explanatory variables during 2000-2020 are presented in Table S3.The standard deviations (STD) of the coefficients for offshore distance (mean STD =0.655) and elevation (mean STD = 0.493) were larger, while the coefficient for species distribution(mean STD=0.193)had a smaller degree67.01%) of areas, respectively.Species distribution showed weaker contribution to GPP variation, with minor differences in regression coefficients than offshore distance and elevation.Among the four zones,Zone B (mean value = 0.14, 0.54 and 0.36) had the highest absolute regression coefficient of species distribution, offshore distance and elevation against GPP,followed by Zone D(mean value=0.17,0.48 and 0.36), Zone A (mean value = 0.13, 0.45 and 0.41) and Zone C (mean value=0.10,0.46 and 0.34) during 2000-2020.

Table 1The correlation coefficient (r) of estimated GPP against NDVI (Normalized Difference Vegetation Index), EVI (Enhanced Vegetation Index), AI (aggregation index), LPI (largest patch index), PD (patch density), SHDI (Shannon's diversity index) of mangrove forest in Zone A,B, C and D during 2000-2020.

The contributions of species distribution, offshore distance and elevation to GPP variation gradually strengthened during 2000-2008,with mean values ranging from 0.37 to 0.39,and then weakened during 2012-2020, with mean values ranging from 0.37 to 0.28.The mean absolute values of MGWR regression coefficients for species distribution,offshore distance and elevation during 2000-2020 ranged from 0.11 to 0.17(mean value=0.15),0.36 to 0.58(mean value=0.49)and 0.31 to 0.44(mean value=0.38),respectively.In 2008,the absolute regression coefficients for offshore distance (mean value = 0.58) and elevation(mean value = 0.44) were higher than those in other years, and the factors in most areas (60.21%-74.42%) had positive contributions to GPP variation.In 2020, the absolute regression coefficients for offshore distance (mean value = 0.58) and elevation (mean value = 0.36) were lower than those in other years.In 2016, the absolute regression coefficients for species distribution(mean value=0.17)were higher than those in other years,and the factors in most areas(72.13%)had positive contributions to GPP variation.of variation.Over the six periods, elevation (mean = 0.174, median =0.171) was the most influential factor on GPP variation, followed by offshore distance (mean = 0.129; median = 0.118) and species distribution (mean = 0.076, median = 0.054).These results suggest that elevation and offshore distance had varying effects on GPP at different points in time, while species distribution had a stronger spatial heterogeneity on GPP variation to a certain extent.

All MGWR models using the three factors passed the multicollinearity diagnosis with adjusted R2values of 0.732, 0.768, 0.745, 0.842, 0.632 and 0.607 in 2000,2004,2008,2012,2016 and 2020,respectively.The spatial differences in the three factors varied for each period, with offshore distance (-3.92 to 7.36, Fig.11.1) reporting stronger spatial heterogeneity and greater fluctuations in the regression coefficient during 2000-2020, followed by elevation (-4.39 to 3.97, Fig.11.2), and species distribution (-0.40 to 1.71, Fig.11.3).During 2000-2020, species distribution, offshore distance and elevation had significant effects(p ＜0.05)on GPP variations,accounting for 18.76%-32.65%,60.61%-72.25%and 49.60%-60.27%of the total area,respectively.

During 2000-2020, positive contributions (regression coefficient ＞0) of species distribution, offshore distance and elevation to the spatial variation of GPP were observed in 56.57%-72.11% (mean = 65.32%),53.24%-68.01% (mean = 61.21%) and 55.81%-75.68% (mean =

4.Discussion

4.1.Comparing GPP estimation across different mangrove forests

During 2000-2020, the estimation of mangrove GPP for 8-day intervals (mean value = 2.82-4.52 g C·m-2·d-1) was compared with GEOSIF GPP products, which had the highest overestimation (mean value=3.99-4.93 g C·m-2·d-1),followed by GLASS GPP(mean value=2.85-4.81 g C·m-2·d-1) (Fig.8).It should be noted that the above products were primarily estimated for terrestrial vegetation, with GEOSIF GPP being derived from linear relationships between SIF and GPP(Li and Xiao,2019).Due to a wide range of stress factors that mangroves are subjected to,such as seawater erosion,salt and low oxygen content,their SIF signal is relatively lower than that of terrestrial vegetation (Li and Xiao,2019).

The GLASS GPP model parameters are estimated based on a limited range of observations and are subject to the effects of climate change,land use change,etc.,which can lead to bias in the estimates(Hu et al.,2018).However,the MODIS GPP was much lower than our results(Fig.8),which was consistent with the findings of Turner et al.(2003a) and Liu et al.(2015),who also observed an apparent underestimation in the MODIS GPP product.This could be attributed to significant errors in referenced land cover classification data,as well as the maximum light energy use values in MOD17A2H being smaller than actual values(Wang et al.,2017).

Table 2The results of factor detector analysis of GPP variation with six geo-environmental factors(species distribution,offshore distance,elevation,slope,planar curvature and profile curvature) during 2000-2020.

Fig.10.The q values of interaction detector analysis of GPP variation based on six factors in(a)2000,(b)2004,(c)2008,(d)2012,(e)2016 and(f)2020.Interaction detector analysis is one method of Geodetector.

Fig.11.1.Spatial difference in the contribution of species distribution to mangrove GPP variation in 2000, 2004, 2008, 2012, 2016 and 2020.

Based on hydroclimatic data and the LUE model, we found that our study area held an increasing tendency of GPP during 2000-2020(Fig.5)with a highly increasing trend of 23.21%(trend value ＞0.50).This result was consistent with the findings of Rasquinha and Mishra (2021), who conducted their work in India from 2000 to 2020.The mean GPP value for our study area was 6.35-8.33 g C·m-2·d-1(Fig.4), similar to the GPP values (mean value = 8.22-10.96 g C·m-2·d-1) observed in tropical forests(Rhizophoraceae family)in Asia and America(Kanniah et al.,2021).In contrast,Lele et al.(2021),who employed the LUE model to estimate GPP along the east and west coasts of India dominated by Rhizophora mucronate and Sonneratia apetala, reported a much lower mean mangrove GPP(1.20-7.70 g C·m-2·d-1) than that in our study.Additionally, Zheng and Takeuchi (2022), who estimated GPP by a spatially explicit blue carbon(BC)model,demonstrated that the mean value of mangrove GPP(4.32±2.45 g C·m-2·d-1from 2000 to 2019) in the continental United States(CONUS)was also much lower than that in our study area.

Based on the EC measurements, greatly different annual GPP values were also reported in various regions.For example,in China,Liu and Lai(2019), Zhu et al.(2021b) and Zhang et al.(2021) reported an annual mangrove GPP of 2741-2827 g C·m-2·y-1in Hongkang, 2197 g C·m-2·y-1in Yunxiao and 1335-1466 g C·m-2·y-1in Yingluo Bay.In Quintana, Mexico, Alvarado-Barrientos et al.(2021) estimated the annual mangrove GPP to be 2473 g C·m-2·y-1; while in Sundarbans,India, Rodda et al.(2016) calculated a mangrove GPP value of 1271 g C·m-2·y-1.These differences in GPP values, as estimated by the LUE model and EC measurements, might be attributed to variations in the growth environment,leading to the climate differences and differences in mangrove species, maximum LUE value, phenology and GPP models(Lele et al.,2021).

4.2.LAI-based LUE model for mangrove GPP estimation

To some extent,GPP estimations with data-driven phenology models(e.g.InTEC and SECRETS) have some uncertainties, because they are largely dependent on the quality of input model parameters (e.g.temperature, precipitation, irradiance and carbon allocation) (Zheng et al.,2018; Xie et al., 2021).However, due to the lack of high spatial and temporal resolution data in the phenology model, the GPP estimations generally differ greatly from the measured value.To overcome such limitation, several studies take LAI derived from remote sensing images into account to better estimate GPP, which can provide functional information on plant traits at an appropriate spatial and temporal resolution(Tagliabue et al.,2019).

Some studies claimed that the GPP trend was nonlinear over time(Patnaik and Biswal, 2020).Moreover, Zhang et al.(2019) argued that the growth of GPP was lower than that of vegetation quantity, and Alvarado-Barrientos et al.(2021) also showed that the GPP trend was much weaker globally than LAI when climatic factors were added to the model to estimate GPP dynamics.Although the mean LAI (Fig.3)continuously increased from 2000 to 2020, we found that the corresponding mean GPP (Fig.4) did not follow the similar trend and even decreased from 2004 to 2008.The results suggested that LAI might not be a direct proxy for mangrove GPP, because climatic (e.g.temperature,humidity and rainfall) and anthropogenic factors (e.g.land use change,deforestation and biological invasion) may also greatly affect the photosynthetic capacity of vegetation(Ferreira et al.,2007;Barik et al.,2018;Gnanamoorthy et al.,2020).

Fig.11.2.Spatial difference in the contribution of offshore distance to mangrove GPP variation in 2000, 2004, 2008, 2012, 2016 and 2020.

4.3.Spatiotemporal pattern of mangrove forest area

We found that about approximately 24.93% of mangrove forests in Gaoqiao Mangrove Reserve experienced expansions during 2000-2008,while 13.62%of areas showed losses and degradation during 2016-2020(Figs.S1-S2).Jia et al.(2018)noted that the decrease from 48,801 to 18,702 ha in mangrove area in China during 1973-2000 was attributed to economic development and agricultural reclamation.In recent years,the establishment of protected areas and the implementation of mangrove protection regulations by the government have led to a recovery of mangroves in China to 22,419 ha in 2015 (Jia et al., 2018), while the decline in mangrove area during 2016-2020 could be related to the construction of coastal dams, urban expansion and pollution of aquaculture.Ahmed et al.(2017)and Valiela et al.(2001)demonstrated that aquaculture and urbanization were the major drivers of mangrove loss.

Moreover,several studies have reported that 14%and 38%of global mangrove loss is associated with fish culture and shrimp culture,respectively(Valiela et al.,2001;Mukherjee et al.,2014),and Alam et al.(2022) even showed that mangrove defoliation could enhance the production of commercial shrimp aquaculture.Arifanti et al.(2019) and Griscom et al.(2017)reported that converting mangroves to aquaculture resulted in the loss of approximately 51%-54% of the total ecosystem carbon stocks of the existing mangroves during 1994-2015.Overall,the major drivers of mangrove loss and degradation around the world are coastal development, agricultural activities, urban land conversion,aquaculture, drainage, animal husbandry, pollution and climate change(Méndez-Alonzo et al.,2008;Valderrama et al.,2014;Lewis et al.,2016).

The quality of mangrove ecosystems is highly dependent on the landscape characteristics in surrounding habitats (Drew and Eggleston,2008; Su et al., 2022).We investigated the relationships between mangrove dynamics and surrounding environmental landscape metrics at the spatial and temporal scales(Fig.S3).During 2000-2020,we found that the mean AI and LPI of mangrove forests gradually decreased, and the mean PD and SHDI slightly increased, indicating that mangrove forests tended to be more fragmented with stronger resistance to artificial disturbance and better connectivity between patches.Such result is consistent with the results of Lele et al.(2021) who declared that the overall landscape changes of mangrove forests in Zhanjiang gradually improved over the past 40 years(1978-2018).

5.Conclusion

This study examined the spatiotemporal dynamics of mangrove gross primary productivity (GPP) by utilizing Geodetector and MGWR to quantitatively analyze the relationship between GPP variation and geoenvironmental factors.The findings revealed a reduction in the mangrove forest area within the reserve, primarily due to aquaculture and urbanization.However, despite this loss, the mean value of mangrove GPP exhibited a continuous increase during 2000-2020.Additionally, approximately a quarter of the areas analyzed displayed a significant upward trend in mangrove GPP between 2000 and 2020.Further analysis using factor detector and interaction detector methods identified species distribution, offshore distance, and elevation as significant factors influencing GPP variation.Accurately assessing mangrove GPP is crucial for effective ecological restoration efforts,hence, this study provides valuable insights for mangrove restoration efforts under the Global Framework Convention on Biodiversity and enhances our understanding of the mechanisms driving spatiotemporal dynamics in mangrove ecosystems.However,the lack of mangrove GPP products with medium spatial resolution (10-30 m) has hindered the validation of estimated mangrove GPP values from 2000 to 2020 at the same spatial resolution.Existing GPP estimation models predominantly utilize vegetation indices as proxies for light capacity,resulting in certain errors in estimating photosynthetically active radiation.Consequently,there is a need for new methods that directly incorporate plant photosynthesis, such as chlorophyll fluorescence (SIF), to accurately estimate mangrove GPP in future studies.

Fig.11.3.Spatial difference in the contribution of elevation to mangrove GPP variation in 2000, 2004, 2008, 2012, 2016 and 2020.

Funding

This work was supported by Guangdong Basic and Applied Basic Research Foundation (2019A1515010741 and 2021A1515110910),Guangdong Regional Joint Fund-Youth Fund (2020A1515111142) and Shenzhen Science and Technology Program (JCYJ20210324 093210029).

Data availability

Data are available on request from the authors.

Authors’contribution

The authors confirm contribution to the paper as follows:Guofeng Wu and Junjie Wang, study conception and design; Yinghui Zhang and Jianing Zhen, analysis and interpretation of results; Zhen Shen, Kunlun Xiang, Yongquan Wang and Haoli Xiang, field investigation and data collection; Demei Zhao and Junjie Wang, draft manuscript preparation.All authors reviewed the results and approved the final version of the manuscript.

Declaration of competing interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

Acknowledgements

We thank the two anonymous reviewers for their constructive comments and suggestions.

Appendix A.Supplementary data

Supplementary data to this article can be found online at https://doi.i.org/10.1016/j.fecs.2023.100137.