Responses of ecosystem services to natural and anthropogenic forcings: A spatial regression based assessment in the world's largest mangrove ecosystem

• Multiple valuation approaches were used to quantify Ecosystem Services (ESs). • Six regression models were utilized for spatial modeling of ESs. • Six major driving factors were considered for spatial regression modeling. • Climate change is the most crucial component of ESs degradation. • Socio-economic and developmental factors have negligible effects on ESs. ⁎ Corresponding author. E-mail address: srikanta.arp@iitkgp.ac.in (S. Sannigrah https://doi.org/10.1016/j.scitotenv.2020.137004 0048-9697/© 2020 Elsevier B.V. All rights reserved. a b s t r a c t a r t i c l e i n f o


H I G H L I G H T S
• Multiple valuation approaches were used to quantify Ecosystem Services (ESs). • Six regression models were utilized for spatial modeling of ESs. • Six major driving factors were considered for spatial regression modeling. • Climate change is the most crucial component of ESs degradation. • Socio-economic and developmental factors have negligible effects on ESs.

G R A P H I C A L A B S T R A C T
a b s t r a c t a r t i c l e i n f o

Introduction
Ecosystem Services (ESs) are a variety of supports and benefits that humans obtain from the natural environment (Costanza et al., 1997;Braat and Rudolf, 2012;Nelson et al., 2009;Daily et al., 2009). Ecosystem Service Function (ESF) refers to a bundle of ecological processes that function at a varying environmental setups and maintain fundamental ESs such as production of raw material, soil formation, decomposition of dead materials, and nutrient recycling (Braat and Rudolf, 2012;Seppelt et al., 2011). Another concept termed ecosystem service values (ESVs) encompasses computational approaches that involve quantification of both biophysical and monetary values of natural goods and services (Schmidt et al., 2017). The ESVs of natural capitals demonstrate the monetary importance of both ecological assets and socio-ecological statuses of an ecosystem (Sannigrahi et al., 2018(Sannigrahi et al., , 2019a(Sannigrahi et al., , 2019b. Most of the Earth's ESs have experienced a decreasing trend over the last few decades (Millennium Ecosystem Assessment, 2005;Chen et al., 2019). According to Millennium Ecosystem Assessment, 2005, approximately 60% of the global ESs are either degraded or used in an unsustainable way (Xu et al., 2018). Among the key driving forces (e.g., climatic factors, land-use change, and socio-economic factors) that are responsible for the changing structure, process, and function of ESs, research has focused mostly on the factors pertaining to land degradation (Cowie et al., 2018). However, the significance and directional (both synergies and trade-offs) effect of the other driving forces on spatially varying ESs have not been substantially explored in the existing literature (Keesstra et al., 2018;Chen et al., 2019;Wu et al., 2019). Thus, this paper presents a thorough evaluation that explores the causal impacts of driving forces on the provision of ESs using as a pilot study in the Sundarbans region, India.
Among the key driving factors, Land Use Land Cover (LULC) change plays a significant role in ES provision Zambon et al., 2019;Luo et al., 2019). For instance, food production, soil formation, nutrient, climate, and gas regulation services are directly connected to the abundance and health of ecological land, including forest, wetland, and grassland (Arowolo et al., 2018). The destruction and degradation of such ecological lands is a widely spread phenomenon in the world, and it is predominantly due to urbanization and agricultural expansion to meet the demand of the ever-growing population that would eventually lead to a substantial loss in biodiversity . Chen et al. (2019) showed that the changes in pattern and configuration of the ecologically important natural capitals have severely affected the fundamental ecological processes. Robertson and Swinton, 2005 found that the expansion of agricultural lands at about 13 million ha yr −1 globally was mainly at the expense of forest land, and 40% of the global Earth's surface had already been converted to cropland in order to meet the increasing food demand.
Climate change can significantly impact the provision of ESs through altering the functions of ecological systems (Nelson et al., 2013;Boone et al., 2018;Chiabai et al., 2018;van der Geest et al., 2019;Smith et al., 2019;Luo et al., 2019). The climate change-induced impact, both positive and negative, on welfare-enhancing ESs is expected to increase rapidly around the world (Nelson et al., 2013;Schäfer et al., 2018, Dow et al., 2013Staudinger et al., 2012). The negative influences of climate change on ESs may affect and change behavioural patterns and sensitivity of different species and biotic/abiotic organisms that are beneficial for providing many regulating, supporting, and cultural services (Schäfer et al., 2018). Dow et al. (2013) found that, in South Asia, every 1°C increase above the 26°C threshold of the night-time temperature would cause a decline of crop production by 10% due to its impact on the rice pollination and flowering services. Additionally, if the night time temperature exceeds the maximum threshold limits (35°C), the current rice varieties would face hardship to adjust the increased temperature.
Socio-economic proxies, including governmental policies, demographic and economic structure of the society, also play crucial roles in the destabilization of ESs in a given ecosystem (Wang et al., 2012, Maes et al., 2012Hauck et al., 2013;Lü et al. (2012)). These factors can be intertwined with the aforementioned LULC and climate changes in ES dynamics. Lü et al. (2012) assessed four key ESs in the Loess Plateau under China's ecological rehabilitation policies and found that the policies with embedded economic incentives enhanced soil conservation and carbon sequestration but reduced water yield at certain climate conditions. Zoderer et al., 2016 investigated the link between perceived social-cultural values of ESs and socio-demographic background in South Tyrol (Italy), where the results demonstrated the importance of cultural background of people in ES valuation. Wang et al. (2018) used a dynamic-CLUE model to project different land-use scenarios for ES valuation in Wuhan city (China) and revealed that the socio-economic development under urban expansion would lead to degradation of all individual ESs, but ecological protection policies could mitigate these adverse impacts.
In the last few decades, the Sundarbans mangrove ecosystem has suffered from a variety of anthropogenic and natural adversities. The physical disturbances include sea-level rise and accelerated coastal erosion, periodic cyclonic storms and coastal flooding, increasing salinity due to the shortage of freshwater and siltation in Bidyadhari river located in the Sundarbans, destruction of embankments and degradation of mangrove, etc. (Sánchez-Triana et al., 2018;Mukherjee et al., 2013;Manna et al., 2010;Mukherjee et al., 2013). The Sundarbans mangrove ecosystem provides the necessary habitat for several commercially and ecologically significant aquatic organisms that support inland and deepsea fisheries, which is one of the primary economic activities of the Indian Sundarbans. Additionally, the mangrove ecosystem is highly effective in fixing (15-46 × 10 12 mol yr −1 ) and storing (3 × 10 14 mol yr −1 ) carbon (Alongi, 2009;Alongi, 2012), while the Sundarban mangrove ecosystem alone sequesters nearly 25 × 10 10 mol yr −1 of atmospheric CO 2 Ray et al., 2014). Apart from these benefits, the Sundarbans mangroves contribute substantially to the provision of other ESs, such as reduction of coastal erosion and flood protection, nutrient storing and recycling, atmospheric gas regulation, biomass production, and preventing soil erosion, which are all essential for maintaining the livelihood of millions of people in this region. Despite the profound ecological and economic significance, limited efforts are offered so far to explore the socioecological significance of this vibrant ecosystem adequately.
This study has made an effort to identify the key controlling factors and assess their causal effects on ESs in the Sundarbans region. Specific objectives of this study are (1) to identify the key driving factors, i.e., biophysical factors, climatic factors, land degradation factors, environmental stress factors, development and socio-economic factors that have substantial impact on ESs; (2) to analyze synergies and trade-offs among the driving factors and ESs; (3) to evaluate the effectiveness of the proposed hybrid approach for policy implication and decision-making. The outcome of this study can not only enrich the current knowledge of ES research but also provide a more inclusive and in-depth view of the complex socioecological nexus between human development and ecosystem management, in line with the 12th Sustainable Development Goal stated as "Sustainable consumption and production aims at doing more and better with less by empowering sustainable and responsible consumption and improving the status of ecosystem services." (https://www.undp.org/content/undp/en/home/ sustainable-development-goals.html).

Study area
The Sundarbans mangrove is the world's largest single tract halophytic-mangrove ecosystem that covers nearly 3% total area of the world mangroves (Mitra, 2015;Akhand et al., 2016). It is also the world's largest coastal wetland ecosystem, lying in the convergence of Rivers Ganga, Brahmaputra, and Meghna and covering approximately 10,000 km −2 geographical area, of which Bangladesh shares 62% and India shares 38%. The physiography of the deltaic Sundarbans includes sand beaches, tidal creeks and inlets, sand flats, mudflats, dunes, estuaries, salt marsh, and mangrove littoral swamps (Fig. 1). This area experiences three meteorological seasons, pre-monsoon from February to May, monsoon from June to September, and post-monsoon from October to January (Ray et al., 2014). The average yearly rainfall is 1920 mm and the average yearly humidity is around 82%. The elevation of the deltaic complex of Sundarbans varies from 3 to 8 m from the mean sea level (Ray et al., 2014). The Indian Sundarbans consists of 102 connected and dispersed islands, of which 52 are habited, 48 are uninhabited, and the remaining 2 have been eroded due to sea-level rise (Mitra and Zaman, 2015). The Sundarbans area was enlisted as a UNESCO world heritage site in 1997 for its profound socio-ecological and biodiversity importance.

Data source and processing
This study incorporates multi-source remote sensing, climatic, biophysical, socio-economic, demographic, and ancillary data products for approximating six major driving factors: biophysical, climatic, land degradation, environmental stress, socio-economic, and developmental factors that can have substantial effects on the provision of ESs in the Sundarbans (Table S1). The smoothed and filtered Enhanced Vegetation Index (EVI) and Normalized Difference Vegetation Index (NDVI) data provided by the University of Natural Resources and Life Sciences, Vienna (http://ivfl-info.boku.ac.at/) were used as biophysical variables. A biophysical scale factor was used to retrieve the actual EVI and NDVI values. The Net Primary Productivity (NPP) was derived from the Moderate Resolution Imaging Spectroradiometer (MOD17) ecosystem model The data required for the MOD17 model include absorbed photo-synthetically active radiation (APAR), fraction of photosynthetically active radiation (fPAR), photosynthetically active radiation (PAR), solar radiation, vapour pressure deficit, temperature, and water stress scalar factors. The topographic variables, including elevation and slope information, were derived from the 90 m Shuttle Radar Topography Mission (SRTM) Digital Elevation Model (DEM) data. The Land Surface Water Index (LSWI) and Soil Adjusted Vegetation Index (SAVI) variables were retrieved from the Landsat Thematic Mapper (TM) data. The Soil Moisture (SM) and other key climatic variables were extracted from the gridded climate data provided by Terra Climate. The development variables, including educational, infrastructure, transport, and power usage facility, were taken from district statistical handbook (DSH), Census of India (http://censusindia.gov.in/). Land Surface Temperature (LST), Palmer Drought Severity Index (PDSI), and Standardize Precipitation and Evapotranspiration Index (SPEI) are the key environmental stress variables used in this study. Monthly PDSI data were aggregated and converted into the annual unit to calculate the yearly drought severity of the region. Additionally, 12-year SPEI values were gathered for estimating the annual water balance and water availability of the region. The aridity index (AI) and evaporative index (EI) values were derived from the ratio of precipitation and evapotranspiration. These two indices designate the climatic in-favorability of the region. The land degradation indicator include landscape fragmentation and connectivity indices. These were calculated using the Fragstat v4.2.1.603 software. Landsat TM data was utilized for this purpose. Several demographic and economic proxies were incorporated to establish the linear and non-linear associations between the socio-economic changes and ESs. These variables were mainly derived from the census of India report and district statistical handbook provided by the Department of Planning & Statistics, Government of West Bengal (https:// www.wbpspm.gov.in/). Tables S1 and S2 provide all the aforementioned details of data processing.

Identification of key ecosystem services of the Sundarbans
A primary survey was conducted in early 2018 to identify the most relevant ESs of the Sundarbans region. The intention of the participatory public survey was to get acquainted with the perception and awareness of the local residents about ESs and how these varieties of ecosystem and landscape components shape their daily livelihoods in this complex and dynamic deltaic lobe of the Sundarbans. In the beginning, 24 provisioning, 18 supporting, and 10 regulating services were included in the initial questionnaire for the public participatory survey. A total of 14 villages, mainly located in the Patharpratima, Gosaba, and Kultali blocks in the Sundarbans, were selected for this analysis. The responses of the local residents (N = 160) in connection with their perception and awareness about key ESs were collected using a five-point scale, ranging from 0 to 4. For provisioning and cultural services, responses were collected using the following scale: 0 = never used, 1 = used once; 2 = used few/several times; 3 = used regularly; 4 = do not know. For regulating services, response modules were categorized as 0 = no relevance; 1 = low relevance; 2 = moderate relevance; 3 = high relevance; 4 = very high relevance. In this selection process, the 5th option (i.e., 4 = do not know) has not been considered in the final analysis on provisioning and cultural services. This option (4 = do not know) was included in the questionnaire only to know how many people have absolutely no idea what ESs stand for and what was the reason behind this unawareness. Additionally, service providing capacities of the main ecosystem types, including rural/urban settlement, grassland/ mixed vegetation, cropland, aquaculture, mangrove, social forestry, wetland/pond/lake, and river/stream, were evaluated through a landscape-based ESs scoring approach proposed by Burkhard, 2009;Burkhard, 2014 (Tables S3, S4, S5).

Quantification of ecosystem services
To overcome the uncertainty and bias in the valuation and calculation of ESs, three different approaches, namely biophysical, economic, and hybrid methods, were incorporated in this study. Both biophysical and monetary values of the key ESs were estimated. Furthermore, a hybrid ESs valuation method was developed by aggregating the biophysical and monetary values of ESs. The details of each valuation method are described below.

Biophysical methods
The biophysical estimates of ESs were based on the calculation of Net Primary Productivity (NPP). NPP was derived from five ecosystem models including the Carnegie-Ames-Stanford-Approach (CASA) (Potter et al., 1993;Field et al., 1995), Eddy Co-variance Light Use Efficiency (EC-LUE) (Yuan et al., 2007(Yuan et al., , 2014, Global Production and Efficiency Model (GLO-PEM) (Prince and Goward, 1995), Moderate Resolution Imaging Spectroradiometer Model (MOD17) (Zhao et al., 2006Running et al., 2004, and Vegetation Photosynthesis Model (VPM) (Xiao et al., 2004). The overall performances of the five NPP models were thoroughly evaluated to identify the best performing model. The biophysical and economic values of 10 key ESs including biological control, climate regulation, cultural, disturbance regulation, genetic, habitat, nutrient cycling, raw material provision, water supply, and waste treatment services were estimated from NPP and other biophysical inputs (e.g., precipitation, evapotranspiration, runoff, elevation, slope, water body occupancy ratio). Additionally, a spatially explicit Integrated Valuation of Ecosystem Services and Tradeoffs (InVEST) model was adopted for quantifying the biophysical values of six major ESscarbon storage, sediment retention, habitat service, water supply, and nutrient retention (Barral and Oscar, 2012;Song et al., 2015, Song andDeng, 2017).

Economic methods
The economic valuation of ESs consisted of five successive steps: (1) determining equivalent weight coefficient; (2) parameter adjustment and rectification; (3) determining standard invariant equivalent value factor; (4) dynamic correction and invariant/comparable economic valuation; and finally (5) estimating regional ESV using adjusted coefficient (Sannigrahi et al., 2019a(Sannigrahi et al., , 2019b. For determining the adjusted equivalent weight coefficient of each ES of interest, the global equivalent weights developed by Costanza et al. (1997Costanza et al. ( , 2014 were adopted. As the food production function of agricultural land is considered to be the most direct ES, the weight coefficient of food production service of cropland was used as the base for approximating the weight coefficients of other ESs. In addition, for adjusting the global weight coefficients and making it functional for local-or regional-level assessment, six dynamic biophysical and climatic variables, including NPP, NDVI, crop yield, precipitation, fractional vegetation cover (FVC), and NPP/NDVI were applied. Using the average crop production and yield statistics of the Sundarbans, the economic value of the food production services of farmland were estimated by assuming that the projected monetary value of food production service could be 1/7 of the real food production estimates (Xie et al., 2008(Xie et al., , 2017Liu et al., 2012). After calculating the per unit food production services of cropland, the economic values of the rest of the ESs were estimated. Additionally, Pearl's S-shaped Growth Curve (PGC) model, Engel coefficient, Inflation Rate (IR), and Consumer Price Index (CPI) data were used for calculating the invariable and comparable economic values of the ESs. The details about the economic valuation methods and approaches are discussed in detail in Sannigrahi et al., 2019aSannigrahi et al., , 2019b.

Hybrid methods
A hybrid method was developed in this study, which is a combination of biophysical and economic valuation methods. The biophysical and economic valuation approaches, which were used for valuation and mapping of ESs, are not exempted from inevitable biases and uncertainties. The biophysical approach mostly depends on variation and quality of key biophysical variables, including NPP, EVI, NDVI, and climatic variables like precipitation, temperature, evapotranspiration, etc. The biophysical models included Integrated Valuation of Ecosystem Services and Trade-off (InVEST), ARtificial Intelligence for Ecosystem Services (ARIES), which rely solely on current and future LULC, temperature, precipitation, evapotranspiration, soil texture, elevation, and slope information. These variables were predominantly derived from the secondary data sources (remote sensing images and satellitebased reanalysis gridded climate datasets). On the contrary, the economic valuation of ESs is connected to several direct and indirect proxy methods including Contingent Valuation (CV), Payment for Ecosystem Services (PES), Travel Cost (TC), Damage Cost (DC), Benefits Transfer Method (BTM), and statistical value transfer model. Pearson correlation coefficient and Coefficient of determination have been implemented to evaluate the consistency between biophysical and economic valuation approaches. The outcomes of biophysical and economic indices were integrated to form a hybrid method. Additionally, to assess the effects of the selected biophysical, climatic, environmental stresses, land degradation, socio-economic, and developmental factors on ESs, a spatial regression was performed for all three valuation methods (biophysical, economic, and hybrid). The details of the driving factors and spatial regression models are discussed in Sections 2.5 and 2.7.

Selection of driving factors
Six relevant driving factors, including biophysical, climatic, environmental stresses, land degradation, socio-economic, and developmental factors, were used in this study to evaluate their effects on ESs. The biophysical factors consist of 8 variables, the climatic factors include 12 variables, the development factors have 4 variables, the environmental stress factors include 9 variables, the land degradation factors are comprised of 16 variables, and a total of 11 variables was selected to represent the socio-economic status of the region. The details about the selection and sources of these driving factors are given in Table S2. The zonal mean values of raster layers were calculated using the ArcGIS zonal statistics tool.

Data dimensionality reduction and variance partitioning analysis
The data dimensionality reduction (DDR) techniques including principal component analysis (PCA), factor analysis (FA), cluster analysis enable to reduce the dimension of data and hence make the data-mining process more efficient by simplifying the classification, visualization, communication, and storage of large arrays of data. Among the widely used DDR methods, PCA is most frequently used and appears to be the most consistent DDR technique. The PCA finds the direction of principal components that explains maximum variances in the data set and accordingly rotates each data point to its coordinates using the orthogonal, oblique, and varimax rotation approaches (Holden et al., 2006). The overall functions of DDR are categorized into two main types: feature extraction and feature selection (Fu and Wang, 2003). Although the variety of DDR methods are intended to uphold the main characteristics of the original data, it is sometimes challenging to prevent data loss mainly to misjudgment of original data layers (Fu and Wang, 2003). In this study, the PCA method was utilized to adjust the dimension and magnitude of the original data layers used for approximating the six major driving factors. Subsequently, six variance partitioning methods were performed, including Partial Least Square Regression (PLSR), Principal Component Regression (PCR), Variable Inflation Factor (VIF), Variable Importance Analysis (VIP), Canonical Correspondence Analysis (CCA), and Redundancy Analysis (RDA), to remove the redundant data layers and identify the most relevant driving factors. The PLSR and PCR models consist of descriptive statistics, correlation and standardized coefficients of estimates, actual and predicted estimates, redundancy and VIP information, confidence intervals, factor loadings, factor correlation, factor scores, and data outliers. Using the PLSR approach, the key driving factors that had maximum explanatory power and thus explained the maximum model variances were identified. The PCR and PCA were used to recognize the key components, which explained the maximum model variances and thus increased the model performances. Using the multicollinearity test, the variable inflation factor, correlation, regression, and tolerance values of the key explanatory variables were calculated. Higher collinearity indicates higher dependency among the control variables and vice versa. The PCA, PCR, PLSR, VIP, and VIF tests were performed in XLSTAT Version 2016.02.28451. The associations between the six groups of factors and the ESs provided by natural capitals were examined through hybrid Redundancy Analysis (RDA) and Canonical Correspondence Analysis (CCA). These tests were performed in the CANOCO v4.5 statistical software. Both RDA and CCA evaluated the synergetic and trade-off association between the driving factors and ESs. Additionally, for CCA and RDA tests, 11 final explanatory variables, i.e., elevation, LSWI, slope, SM, ET, TRANS_FAC, POWER_FAC, SPEI, FOREST_A, AREA_NA, and CULT_W, were used. The final 11 variables were derived from the variance partitioning and data reduction analysis. These variables had the acceptable multicollinearity values (below 4) and based on the collinearity statistics of these variables, several multivariate regression and redundancy analysis were performed (Table S6). The aforementioned PCA, PLSR, PCR, VIF, VIP, RDA, and CCA tests were implemented for all three types of ESs methods (biophysical, economic, and hybrid) to identify the best method for ES evaluation (Tables S7, S8, S9, S10, S11).

Spatial regression and autocorrelation models
The spatial autoregressive and multivariate regression methods have been the focus of recent scientific studies in spatial decision models (Kelejian and Prucha, 1998;Drukker et al., 2013). During the past decades, this approach was used extensively in urban growth modeling and environmental suitability studies (Qu and Lee, 2015;Qiu and Turner, 2013;Li et al., 2018). The spatial regression models (SRM) are highly crucial to evaluate the spatial effects (Chakraborti et al., 2019;Fotheringham et al., 1998;de Lima et al., 2005) including spatial autocorrelation, spatial stationarity, and heterogeneity of spatial objects on regression models (Schröter et al., 2015;Tenerelli et al., 2016;Maes et al., 2012). Spatial autocorrelation measures the spatial dependency of the objects in a feature space. This indicates the similarity of spatial patterns over space. Conversely, the random spatial pattern exhibits no spatial autocorrelation. Six spatial regression models were performed in this study to evaluate the spatial dependency and association between the driving factors and ESs. These include Geographical Weighted Regression (GWR), Spatial Error Model (SEM), Spatial Lag Model (SLM), Spatial Error Lag model (SEM_SLM), Ordinary Least Square (OLS) model, and Spatial Autoregressive Model (SAM). Among these six regression models, the OLS model was utilized for determining the global interaction between the six driving factors and 10 key ESs. However, the OLS model does not consider any spatial autocorrelation or homogeneity in the modeling. Furthermore, to examine the individual and combined effects of the driving factors on ESs, the Geographical Detector Model (GDM) was implemented in this study. The GWR is a spatial model that exhibits a spatial non-stationarity in modeling processes (Fotheringham et al., 1998). Unlike OLS, the GWR model produces varying local attributes throughout the feature space after integrating the spatially referenced data layers (Fotheringham et al., 2003;Brunsdon et al., 2002;Lugoi et al., 2019). The GWR model is expressed as follows: where Y i is response variable (ES in this case); β o , β i , and ε are model parameters; a and b are geographical coordinates (latitude and longitude) of the jth point, X i is explanatory variable i (climate, biophysical, land degradation, environmental stress, socio-economic, and development factors). Additionally, several spatial autoregressive models, including SAR, SEM, SLM, SEM_SLM were also incorporated. The SAR model showed great potential in resolving spatial econometric problems as it consists of a define structure that often leads the model inferences more straightforward and facilitates replication (Qu and Lee, 2015;Kelejian and Prucha, 1998). The GWR model was performed using the GWR v4.0 software. The SAM, SEM, SLM, SEM_SLM models were performed using GeoDa v1.12.1.161 and GeoDaSpace v1.2 software modules. The GDM was utilized using the GeoDetector software package (Wang et al., 2010). Additionally, the OLS model was performed using ArcGIS spatial statistics toolbox.
To analyze the spatiotemporal association between explanatory and response variables (ESs), the Local Geary statistic was used in this study. The Local Geary statistic is generally used to estimate the local indicators of spatial association for assessing the attribute similarity of feature space (Anselin, 1995, Anselin, 2019. The value of Local Geary statistic ranges from 0 to 2: a value close to 0 indicates a positive spatial autocorrelation and clusters, while a value close to 2 indicates a negative spatial autocorrelation and outliers; whereas a value close to 1 indicates a random pattern of distribution.

Spatial regression between driving factors and ESs using GWR
The spatial interactions between the six driving factors and ESs were analyzed using the GWR model (Fig. 2). For all the three ES valuation methods, biophysical and climatic factors were highly correlated with the ESs, characterized by a moderate to very high local R 2 values (R 2 = 0.83-0.97). Additionally, very high local R 2 values (R 2 = 0.84-0.97) were obtained over the southern region (Gosaba, Kultali, Patharpratima, Namkhana, Sagar, Kakdwip, Mathurapur, Kulpi, Jayanagar). The same association was observed between the climate factors and ESs, where moderate to very high local R 2 values (R 2 ≥ 0.84) were accounted for the southernmost blocks. For the economic method, a reverse association was observed between climate factors and ESs. While examining the spatial association between the development factors and ESs, moderate local R 2 values (R 2 = 0.46-0.70) were estimates for Gosaba, Patharpratima, Kultali, Kakdwip, Basanti, Hingalganj blocks, while very low regression values were accounted for the northern region. The association between the climatic factors and ESs derived from the economic method was very high (R 2 = 0.84-0.86) over the northern region, and considerably lower local R 2 values (R 2 = 0.81-0.84) were counted for the southern part of the region. For environmental stress, land degradation, and socio-economic factors, most of the blocks adjacent to Bay-of-Bengal (BOB) have produced a high to a very high spatial association between the said driving factors and ESs, mainly those are derived from biophysical and hybrid methods (Fig. 2). The spatial interaction between all six factors, VIF based filtered factors with ESs were also analyzed and presented in Fig. 3. Among the three ESs valuation methods, moderate to very high regression values were accounted between the all/VIF factors and ESs, while, comparatively lower estimates were calculated between the all/VIF factors and ESs derived from the economic methods (Figs. 3,  4). The summarized results of the GWR model is reported in Table 1. From the GWR estimates, it was found that among the six driving factors, the climatic factors produced the highest association with the ESs, characterized by a maximum local R 2 value, followed by environmental stress, biophysical, land degradation, development, and socioeconomic factors.

Effects the driving factors on ESs
The individual and joint effects of the driving factors on ESs were analyzed using the GDM ( Table 2). The interaction effects (q) of biophysical and climatic factors on ESs were the greatest, followed by biophysical/environmental stress, biophysical/land degradation, and biophysical/socio-economic factors. Additionally, the biophysical factors demonstrated very weak interaction effects with the development factors. Climate factors exhibited very strong interaction effects with socio-economic factors, followed by environmental stress, development, and land degradation factors. The interaction between development factors and other factors was the highest for land degradation and environmental stress factors ( Table 2). The q statistics between the final 11 variables and ESs were estimated and presented in Table 3. For elevation, the interaction effects were the largests for ELEV/AREA_NA, followed by ELEV/SM, ELEV/LSWI, ELEV/SPEI, ELEV/ET, ELEV/SLOPE, ELEV/CULT_W, ELEV/TRANS_FAC, ELEV/POWER_FAC, and ELEV/FOREST_A, respectively. Similarly, the interaction effects between the 11 explanatory factors were estimated. Among the 66 pairs of interaction results, the highest interaction was obtained between ET and SPEI, followed by ET/LSWI, ET/TRANS_FAC, ET/SLOPE (Table 3).  The interactions between the six driving factors and ESs were evaluated by six spatial regression models (Table 4). Among the six driving factors, the climatic factors are significantly associated with the ESs, characterized by a very high coefficient of determination value (R 2 = 0.74-0.81). The other driving factors also impacted the ESs significantly. The biophysical factors produced a very high regression estimate (R 2 = 0.59-0.77), followed by environmental stress (R 2 = 0.57-0.79), land degradation (R 2 = 0.71-0.78), socio-economic (R 2 = 0.45-0.66), and development factors (R 2 = 0.27-0.44), respectively (Table 4). Table S7 shows the effects of the six explanatory factors on individual ESs derived from PCR and PLSR models. The interaction between the explanatory factors and ESs showed the highest values for the climate regulation and raw material provision services, followed by genetic service, nutrient cycling, water supply, habitat, waste treatment, disturbance regulation, cultural service, and biological control services, respectively (Table S7).
To examine the directional associations between the 11 driving variables and ESs, two redundancy tests, including hybrid redundancy analysis (hRDA) and hybrid canonical correspondence analysis (hCCA), were performed and presented in Fig. 5. Among the 11 variables, SM, LSWI, ET, AREA_NA, FOREST_A were highly correlated with the ESs derived from biophysical and hybrid methods. There was no significant association between the driving factors and ESs estimated using the economic method. The other variables, i.e., POWER_FA, SLOPE, ELEV, CULT_W, and SPEI, did not exhibit any significant association with the ESs. Therefore, the moisture factors characterized by SM, LSWI, and ET and socio-economic factors (the area not available for agriculture and forest area) were the most significant driving factors with a substantial contribution to the degradation of ESs (Fig. 5).

Discussion
While considering the interaction between each biophysical factor and ESs, the greenness factors, i.e., EVI, NDVI, were observed to be highly associated with gas regulation and climate regulation functions. Both gas regulation and climate regulation services were calculated from the results of NPP, which is the outcome of biophysical (EVI, NDVI), climatic (solar radiation, temperature, soil moisture), and bioclimatic (PAR, LUE, fPAR, APAR) variables. The NPP was found very high over the southern BOB region, which is mostly covered by dense mangrove forests. This suggests that the Sundarbans mangroves are sequestering a substantial amount of gaseous carbon; therefore, strengthening the protection and preservation of Sundarban mangroves would be an effective strategy for the reduction of carbon dioxide emissions and maintaining the carbon balance in this deltaic ecosystem (Rodda et al., 2016). Several research has reported about the carbon sink capacity of the Sundarban mangroves. Rodda et al. (2016) estimated the net carbon influx of the Sundarbans using the Eddy covariance measurement, which is 249 ± 20 gC m −2 year −1 from April 2012 to March 2013. Ganguly et al. (2008) found that the estimated carbon sink of entire Sundarbans was 206 Gg day −1 while the mean net flux of Lothian Island and Sajnekhali region was 48.3 gCO 2 m −2 day −1 . The study conducted by Chanda et al. (2013) measured the CO 2 influx of Sundarban mangroves at different locations viz. Jharkhali, Henry Island and noted that the fluxes varied from 33.69 to 114.91 gCO 2 m −2 day −1 . The findings of this study are in accordance with all these reported statistics. The biophysical factors, i.e., EVI, NDVI, NPP, SAVI, LSWI, showed positive effects on ESs, especially over the southernmost region adjacent to coastal BOB. This region is mostly covered by mangrove and truncated riverine network originating from Hooghly, Bidyadhari, Matla rivers stretched from east to west of Sundarbans. Therefore, the proxies that resemble the biophysical characteristics and richness of the region can be good estimators for evaluating the effects of ecosystem health and vigor on ESs.
Changes in the key climatic components (e.g., temperature, precipitation, evapotranspiration) can have significant impacts on ESs (Geest et al., 2018;Bangash et al., 2013;Buytaert et al., 2011;Walther, 2010). Study conducted by Nelson et al. (2013) thoroughly evaluated the plausible consequences of climate change on different ecosystem functions and suggested that climate change will modify the ability of different key ecosystem functions that support the provision of multiple ecosystem services including food production, wildfire regulation, hazard reduction and coastal flood protection, marine and inland fishery production, water supply, and nature-based tourism and recreational services. Climate change also has a significant impact on human health by increasing the intensity of urban heat islands, amplifying the risk of flooding and resulting in mortality due to the spread of infectious and other diseases (cardiovascular and respiratory) (Chiabai et al., 2018). In this study, the GWR model evaluated the spatial interaction between the aforementioned climatic factors and ESs in Sundarbans (Fig. 2). The climate change impact on ESs has a high variability with geographical space as the corresponding climatic variables produced very high local R 2 values over the BOB region while they have lower regression estimates over the northern region. This suggests that the region and islands (both human-occupied and isolated) that are closer to the coastal stretch have a higher climatic vulnerability than the landward region. The strong association between the climatic and biophysical factors with ESs is observed throughout the study region which suggests that these two factors can significantly determine the ESs provision of an ecosystem. This also implies that climate change is the most critical factor for the ESs changes in Sundarbans. Therefore, among the six relevant controlling factors integrated into this study, the climate change is considered to be the highest contributing elements, followed by biophysical, environmental stress, land degradation, socio-economic, and developmental factors, which is in agreements with the findings of Nelson et al. (2013) and Chiabai et al. (2018).
The Sundarban mangrove ecosystem provides many livelihood alternatives and ecosystem services to the coastal communities living in this region and contributes substantially to improving the overall socio-ecological status of this deltaic ecosystem. Many livelihood options of this region, such as honey collection, crab collection, fishing, aquaculture, and cultivation practices, will be directly or indirectly affected by the prolonged climate change and the associated adversities. This could be even more prominent when anthropogenic and development factors would be aligned with the physical factors such as sealevel rise (SLR), coastal erosion, increasing sea surface temperature and salinity, coastal flood and storms, etc. It is also anticipated by earlier research that the present biodiversity region of Sundarbans will be reduced from 60% to 30%, mainly due to sea-level rise and resulted in the loss of mangroves and coastal land (CEGIS, 2006;Uddin et al., 2013). Due to the loss of mangroves, it is expected that sediment discharge and nutrient load of the riverine network will be affected, salinity (both soil and water) will be increased, and the disproportion between precipitation and evaporation will expedite the formation of cyclones and storms. This would pose serious environmental and livelihood threats to millions of peoples living in this region. Additionally, according to the study of CEGIS (2006) in Bangladesh Sundarbans, the area suitable for Heritiera species (the most abundant mangrove species of Sundarbans) would be decreased by 14% in 2050 (in 32 cm SLR scenario) and by 45% in 2100 (in 88 cm SLR) from the base year of 2001; while for Excoecaria species, the suitable area will be decreased by 7% in 2100 (Uddin et al., 2013). Additionally, the monetary loss due to climate change would be much visible in the southern part of the Sundarbans, where most people are dependent on the forest-based provisioning and supporting services. Therefore, the formulation of effective policies, increasing adaptability, resilience, and promoting alternative livelihoods for the coastal communities of the Sundarbans are required to cope with the environmental problems pertinent to this ecosystem.
As the prolonged climatic extremities have been witnessed in many places of the world, it is expected that coastal storm surge/flood intensities and frequencies will increase across the scale, especially in the tropical coastal region, where climatic extremes designate severe socio-economic and livelihood threats to billions of peoples located in the coastal stretch in India, Bangladesh, China, Thailand, Malaysia, Myanmar and other coastal countries (Brouwer et al., 2007;Douglas et al., 2008;Mirza, 2011). The disturbance regulation service is one of the key regulating services in the Indian Sundarbans. The natural coastal habitats of the Sundarbans, including mangroves, inland wetlands, coastal estuary, salt marshes, sand dunes, act as the first defense system against the prevalent and frequent coastal storms and wave surges that cause the exaggerated shoreline erosion and coastal flooding due to increasing sea-level rise and lack of sediment and water inputs into river system under the construction of Farakka Dam along the Hooghly river. Strong protection measures should be incorporated in the coastal storm management policies for restoring the ecosystem of these natural habitats. This can be done by promoting the soft engineering options provided by these natural habitats. Moreover, the soft engineering supports provided by the natural habitats demand little ongoing maintenance costs, turning this option economically viable and cost-effective (Nelson et al., 2013). In this study, a strong spatial agreement was observed between the driving factors and ESs. Most of the regulating services, i.e., climate regulation, disturbance regulation, nutrient cycling, raw material provision, water supply, and waste treatment services have exhibited very strong synergic associations with climatic and biophysical variability, especially over the southern coastal region. The nexus between the development factors, land degradation factors, socio-economic factors, and ESs were also evaluated (Figs. 2,3,4,5). The development factors, represented by four major development indicators such as educational facility, transport facility, infrastructure facility, power usage facility, were considered as controlling factors for spatial regression and spatial interaction modeling. The spatial associations between development factors/land degradation factors/socioeconomic factors and ESs were not uniformly distributed, as higher local R 2 values were evident over the southern region, and lower regression estimates were accounted in the northern part of the region. These tendencies were reversed while the ESs derived from the economic valuation methods were taken as response variables in spatial regression modeling. However, the results of this study show that the development and socio-economic factors did not produce any strong (non)spatial association with the ESs. The results of the six regression models and GDM are in accordance with this finding. While the land degradation factors represent the land fragmentation and connectivity of an ecosystem, it exhibited a moderate to strong association with the ESs, especially the ESs calculated using the biophysical and hybrid valuation methods (Table 4). The socio-economic factors incorporated in this study, such as forest area, total unirrigated area, total irrigated area, pisciculture area, pisciculture production, etc. unveils the socio-ecological status of the region. The higher the socio-ecological stability, the higher the ESs we observed. Perhaps, this might have been the cause of having a higher local R 2 value over the southernmost administrative block of the Sundarbans, which is mostly covered by forest and agriculture land, while executing the spatial regression models between the socioeconomic/development factors and ESs.
The environmental limiting factors such as temperature stress, water stress have significant impacts on the provision of ESs. Both spatial regression and spatial interaction models have advocated a strong (non)spatial association between the limiting factors and ESs. This association is more prominent over the southern blocks (Gosaba, Patharpratima, Namkhana, Basanti, Kultali, Sagar), and lower local R 2 values of the same are evident in the northern region. The outcome of this test indicated that apart from the profound effects of climatic, biophysical, land degradation factors on ESs provision, the environmental limiting conditions, which refer to both climatic, hydrological, and biogeochemical stress scalars, have a substantial impact on the ecosystem service production.
There are uncertainties associated with the methods for ES quantification and valuation. Among the three methods, the biophysical and hybrid methods exhibited strong spatial accuracies (Fig. 2). The ESs quantified using the economic method did not produce correlation with those using biophysical and hybrid methods. These exceptions could be due to the structural differences between the three valuation methods. As a variety of modeling and calculation approaches were adopted for each valuation method, it is obvious that there might be some uncertainties and biases that exist in the modeling. The methodological differences and uncertainties embedded in the valuation and mapping of ESs manifest complications for the evolving sustainable ES valuation framework and inclusion of the ES concept effectively in national capitals accounting. This also flags some issues in broadening the relevant natural resource management policies for effective decision making and policy implication (Crossman et al., 2012;Boerema, 2017;Wong et al., 2015). Additionally, Crossman et al. (2012) stated that the methodological issues that exist in ES valuation also turns the commodification (ecosystem service productions such as carbon offsets and taxation, auction of ecological conservation, payment for ecosystem services, banking of natural capitals such as wetland, forest) and trading of (non)marketable ecosystem goods and services highly ambiguous as the valuation markets require a certain clarity and transparency.

Conclusion
The interaction effects and complex nexus between the driving forces and ESs supply in the Sundarbans region have revealed that climate change is the most crucial component of ES degradation. In order to maintain the ecological stability and flow of multiple ESs from this dynamic deltaic lobe of the Sundarbans, it is optimal to analyze the sensitivity and responses of natural capitals to any adversarial effects. The series of spatial regression models have explicitly discussed the spatially varying interaction among the natural and anthropogenic forcing and ESs. The results of all six regression models reveal that the socio-economic and development factors have weak to moderate negative effects on ESs. The q statistics derived from the GDM model suggest that the joint effects of the driving factors are much higher than their individual effects. Furthermore, this study has proposed a hybrid valuation method, and it has been observed throughout the study that most of the regression models produced a better estimate for the ESs derived from the hybrid methods. The findings of this research could be useful to the land administrators, environmentalists, policymakers for adopting suitable land resource conservation and management plans for strengthening and protecting the natural capitals, thereby improving the overall socio-ecological status of the region. This study included the most suitable and identical driving factors, which were found to have substantial positive or negative effects on ESs. There is a scope for future work to improve the methods and approaches adopted in this study.

Declaration of competing interest
The authors whose names are listed in this manuscript certify that they have NO affiliations with or involvement in any organization or entity with any financial interest or non-financial interest in the subject matter or materials discussed in this manuscript.