Remote Sensing-Based Mapping of Senescent Leaf C:N Ratio in the Sundarbans Reserved Forest Using Machine Learning Techniques

Carbon to nitrogen ratio (C:N) of senescent leaf is a crucial functional trait and indicator of litter quality that affects belowground carbon and nitrogen cycles, especially soil decomposition. Although mapping the C:N ratio of fresh mature canopies has been attempted, few studies have attempted to map the C:N ratio of senescent leaves, particularly in mangroves. In this study, four machine learning models (Stochastic Gradient Boosting, SGB; Random Forest, RF; Support Vector Machine, SVM; and Partial Least Square Regression, PLSR) were compared for testing the predictability of using the Landsat TM 5 (LTM5) and Landsat 8 to map spatial and temporal distribution of C:N ratio of senescent leaves in Sundarbans Reserved Forest (SRF), Bangladesh. Surface reflectance of bands, texture metrics of bands and vegetation indices of LTM5 and Landsat 8 yearly composite images were extracted using Google Earth Engine for 2009–2010 and 2019. We found SGB, RF and SVM were significant different from PLSR based on MAE, RMSE, and R2 (p < 0.05). Our results indicate that remote sensing data, such as Landsat TM data, can be used to map the C:N ratio of senescent leaves in mangroves with reasonable accuracy. We also found that the mangroves had a high spatial variation of C:N ratio and the C:N ratio map developed in the current study can be used for improving the biogeochemical and ecosystem models in the mangroves.


Introduction
Mangroves grow in the intertidal zone of tropical and subtropical countries, are highly productive and have huge carbon sequestration capacity, especially into sediments [1,2]. Thus, mangroves play a tremendous role in climate regulation. For sustaining productivity and carbon sequestration, the nutrient cycling needs to be maintained [3,4]. Through litter decomposition, a keystone ecosystem process, the nutrient supply is maintained in mangrove forests, which is essential for plant growth [3,[5][6][7][8].
Leaf litter quality such as carbon to nitrogen (C:N) ratio is often regarded as the key determinant of litter decomposition [3,[5][6][7][8]. Litter quality also influences the activities of soil microorganisms that play the direct role in litter decomposition [7,[9][10][11]. For example, poor litter quality (higher C:N ratio) with lower amount of labile C with higher amount of tannin content causes microbial immobilization [8,9,12]. In a forest ecosystem, litter mostly comes from the falling of senescent leaves [7,13]. The variability of C:N ratio in senescent leaves largely depends on species at community and ecosystem scale and Remote sensing-based modeling of functional traits in forest ecosystems has been followed by both parametric and non-parametric advance machine learning regression approaches [23,34,37,38,45]. The advantage of non-parametric machine learning over parametric modeling approaches such as multiple linear modeling is that it can handle both multicollinearity between predictors and noisy data with efficiency in selecting the key predictor from a large amount of predictors [46][47][48][49]. Moreno-Martínez et al. [23] reported the outperformance of Random Forest (RF) modeling technique over multiple linear model, neural networks, Kernel Ridge Regression, and Gaussian Process Regression in mapping leaf functional trait in terrestrial forests at the global scale. Similarly, Gara et al. [38] and Chumura et al. [37] mapped leaf N and C concentration accurately with RF. On the other hand, the application of two other machine learning modeling techniques, support vector machines (SVM) and artificial neural networks, in modeling N and C concentrations over intact and fragmented forest ecosystems have been found to have equal performance [45]. Another ensemble of learning regression method is Stochastic Gradient Boosting (SGB), which can deal with plenty of data and combine important predictors along with average predictors for strong perdition of dependent variable. It found an outperformance of RF in other forestry field (biomass mapping) [50]. However, application of SGB in functional trait such as C:N modeling and mapping is not well known.
In this study, to our current knowledge, for the first time in the mangrove ecosystem, we tested four machine learning models to map the spatial and temporal distribution of the community weighted mean of C:N of senescent in the Sundarbans Reserved Forest in Bangladesh using Landsat TM 5 and Landsat 8. Advanced remote sensing-based modeling approach (machine learning) can contribute in assessing the essential biodiversity variables for biodiversity target achievement (e.g., National Biodiversity Strategies and Action Plans related to Aichi Biodiversity Target 2020) [51]. Particularly, in remote sensing-based forest biochemical properties (C:N as an indicator of decomposition), monitoring and information systems can improve the process of measurement, reporting and verification for achieving global targets and commitments [52]. Therefore, the specific objectives of current study were (i) to select the best machine learning model to predict C:N ratio of senescence leaf using spectral and textural properties of bands and vegetation indices of Landsat TM 5, (ii) to develop a spatial map for observing the spatial variation of C:N ratio of senescent leaf as litter quality trait that controls ecosystem function (decomposition and nitrogen using strategy) in the Sundarbans Reserved Forest, and (iii) to assess the temporal change of C:N ratio of senescent leaf in Sundarbans Reserved Forest (2009-2010 to 2019) using Landsat TM5 and Landsat 8.

Description of Study Area
The Sundarbans is the world's largest contiguous mangrove ecosystem (10,017 km 2 ) shared between Bangladesh (62%) and India (38%) and was declared a Reserve Forest in 1878 [53,54]. The Bangladesh part of Sundarbans is called Sundarbans Reserved Forest (SRF). The SRF is situated in the southwestern part of Bangladesh and is one of the most diverse mangrove forests in the world [16,55] (Figure 1). The total area of SRF is 6017 km 2 , in which mangrove forests occupy about 69%, and the rest is water bodies such as rivers, small streams and canals. The mangrove extent in SRF accounts for 5% of global mangrove forests (83,495 km 2 ) [56]. It accounts for 4.07% of the total area of Bangladesh and 40% of the total forest coverage managed by the Forest Department [57]. Because of its multiple ecosystem services and biodiversity value, it is likely to influence regional and global environments under the changing climate [16,18]. This forest was declared as Ramsar Wetlands Sites in 1992 and the UNESCO (United Nations Educational, Scientific and Cultural Organization) World Heritage Site in 1997 [53].
Remote Sens. 2020, 12, 1375 4 of 21 1135 recorded species, including a large community of endangered species like the Royal Bengal Tiger, the Ganges and Irrawaddy dolphins and saltwater crocodiles [58]. The forest supports and protects the livelihoods of the local communities [59]. This forest is heterogamous in terms of species, hydrological system, and salinity gradient and also faces anthropogenic and natural disturbances similar to other mangroves, which affects the litter quality, growth, productivity and nutrient cycling [7,[16][17][18][19].  A total of 528 species of vascular plants where 24 species are true mangroves and 70 species are mangrove associates have been reported so far [54]. This forest is also rich in faunal biodiversity with 1135 recorded species, including a large community of endangered species like the Royal Bengal Tiger, the Ganges and Irrawaddy dolphins and saltwater crocodiles [58]. The forest supports and protects the livelihoods of the local communities [59]. This forest is heterogamous in terms of species, hydrological system, and salinity gradient and also faces anthropogenic and natural disturbances similar to other mangroves, which affects the litter quality, growth, productivity and nutrient cycling [7,[16][17][18][19].

Field Data
In this study, Sundarbans Reserved Forest inventory data were used, which were collected in 2009-2010 from 150 plots spanning over the whole Sundarbans Reserved Forest for calculating community weighted mean of C:N. The field survey at the 150 plots was conducted from 5 December 2009 to 30 April 2010 with two field inventory groups, each led by an Assistant Forest Conservator. Each group consisted of one Assistant Forest Conservator, one Forest Ranger/Deputy Ranger, two foresters, two students, two laborers and two armed guards. Each group was assigned a small engine boat with boatman. Each trip was seven to ten days long (an interval between two surveys) depending on stored food and availability of fresh water. Before starting each journey, both groups discussed together with detailed maps and GPS (Global Positioning System) units to plan for the next plots. Local knowledge of laborers, guards, Forest Department district staff and even fishermen aided the crews' efforts to find suitable routes to plots and minimize walking distance and time. Generally, each group completed one plot per day, but often this pre-planning activity helped the groups to complete more than one plot a day. This was due to the developed skill of our crew members and easy access to the plot. Somedays, we even completed a third plot by working both groups together. All the documentation had been reserved. Doing forest inventory in mangrove forest is always challenging. Sometimes, we had to walk more than 600 m across the dense forest and make a pool to cross a small canal (if it was unavoidable). One day, we even had to swim a canal with a width of around 10 m because after finishing the survey, when we returned to the canal, it was full of water due to the high tide.
The plots were 4 min latitude and 2 longitude away from each other and were composed of five 10-m radii nested circular subplots [16]. We also took four pictures at four directions at the center of each subplot. All the documentation had been reserved. For detailed field inventory methods see Rahman et al. [16]. This study used 93 of those 150 plots, because at least one subplot (out of five subplots) of the remaining 57 plots contained more than 10% water bodies or was underwater completely. Living tree species that had DBH (diameter at breast height) ≥ 10 cm (18 species including one unknown species) were used in this study, and species wise basal area (cross sectional area at 1.30 m height) per plot was calculated for use as the weighted value in estimating community weighted mean of C:N (Equation (1)) [25].
This study used the C:N data of Chanda et al. [13]. C:N data were averaged because they analyzed the data in three different seasons (monsoon, pre-monsoon and post-monsoon). They plucked 30 senescent leaves for each species and made a composite sample for estimating C and N. Then, they oven dried each composite sample at 75 • C for approximately 7 days until a constant weight was achieved [13]. After that, they grounded the dried samples to a fine consistency using a mortar and pestle and then analyzed for carbon and nitrogen content (Perkin Elmer 2400 Series II CHNS/O Elemental Analyzer) [13]. We used the initial C:N ratio of litter. The average C:N value (29.22; Table 1) of all the species estimated by Chanda et al. [13] was used for the species in which C:N was missing. The community weighted mean was calculated using the FD package in R [60]. Plot wise community weighted mean was estimated as the mean trait value of each weighted by the relative basal area of the species in a given plot by following the Equation (1) [61,62]. The community weighted mean of C:N varied from 18.40 to 39.99 with a mean value of 35.16 across the whole Sundarbans Reserved Forest.
where CWM is the community weighted mean of C:N ratio, S is the species richness, P i is the proportional of relative basal area of species i, and C:N i is the C:N ratio of species i.

Remote Sensing Data
The Landsat surface reflectance data from Thematic Mapper 5 (Landsat TM 5) were used in this study, which have been corrected for atmospheric, reflectance and topographic effects [63]. Both Landsat TM 5 and Landsat 8 have a 16-day repeat cycle, and their band specification is shown in Table 2 [64]. Cloud and cloud shadow pixels in all images were masked based on the pixel quality attributes provided in the products. Composite images of Sundarbans Reserved Forest for each sensor, Landsat TM 5 (2009-2010) and Landsat 8 (2019), used in this study, were aggregated from 89 and 75 images, respectively [63]. The median value of seven multispectral bands was used for modeling purposes because it was robust against extreme values and produced a good quality image which was representative of the time period [63]. One hundred and fifty-seven predictors were extracted from seven bands of Landsat TM 5 that comprised surface reflectance and 18 Gray Level Co-occurrence Metrics (GLCM) for each of the seven bands and Vegetation indices (VI). Of those 150 predictors, 48 (r ≥ 0.20 with C:N; Table 2; Supplementary Materials) were used in modeling C:N. All the image processing was done using the Google Earth Engine Platform and R environment (Version 3.5.3). Table 2. List of bands (Landsat TM 5 and Landsat 8), Grey Level Co-occurrence Metrics (GLCM) texture and vegetation indices with code used in this study.

Stochastic Gradient Boosting (SGB)
Stochastic Gradient Boosting (SGB) is a hybrid approach combining both boosting and bagging techniques [77,78]. In SGB, a small classification or regression trees are sequentially built from the residuals of the preceding tree(s). Instead of focusing on the full training set, the SGB carries out a boosting by selecting (without replacement) at each step a random sample of the data, leading to a gradual improvement of the model. Thus, increasing the chance that incorrectly classified observations will be correctly classified in the next tree. The weighted majority of classification determines the final classification of each observation across the sequence of trees [47]. Ridgeway [77] provided a comprehensive background and mathematical functions of SGB.

Random Forest (RF)
Random forest (RF) is an ensemble classifier or regression technique [79]. It is based on the well-known classification and regression tree approach (CART) [79], which consists of many decisions or regression trees, in which each tree depends on the values of a random vector sampled independently and with the same distribution for all trees in the data [80]. Random forest with a large number of trees is robust against overfitting, noise, as well as non-informative and correlated features [79]. RF is a widely used modeling technique in remote sensing-based studies in forestry such as biomass [50] and functional traits mapping [23,38,45].

Support Vector Machine (SVM)
Support vector machine (SVM) belongs to the regression model family and is a kernel-based learning method from statistical learning theory [81]. SVMs mainly involve a projection of the data into a high-dimensional feature space using a valid kernel function and then application of a simple linear regression within this enhanced space [82]. The resulting linear regression function in the high-dimensional feature space corresponds to a non-linear regression in the original input space [83].

Partial Least Square Regression (PLSR)
One of the most common methods to build a multi-independent model is Partial Least Square Regression (PLSR), and its advantage is that PLSR can handle data with strong collinearity and noise, as well as situations in which the number of variables considerably exceeds the number of available samples [84]. In modeling of functional traits, this modeling technique has also been used [85].

Model Development
In this study, we constructed repeated cross-validation (fold = 10) with ten repetitions for model fitting and model validation with the R "caret" package. In the10-fold cross validation, at each repetition, the sample data were split into 10 subsets, where nine subsets of samples were used to train models, and one was treated as validation samples. We set a seed value of seven for each modeling method, providing the same dataset for each modeling method. To select the simple and optimal model as suggested by Breiman et al. [79], this study used the "one standard error rule". By using the cross validated models, we predicted C:N ratio for 93 plots for each modeling method. Then we also fitted a simple linear regression between field C:N ratio and predicted C:N ratio for testing models fitness in terms of slope and intercept parameters against the 1:1 line. We did not perform any collinearity test or feature selection, since all the modeling methods are non-parametric, and they have their own resistance against collinearity.

Model Performance
In this study, three common performance indices, namely Mean Absolute Error (MAE), root mean squared error (RMSE) and coefficient of determination (R 2 ) were used. MAE % and RMSE % were also used as model performances which were calculated by dividing MAE and RMSE by mean value of field C:N ratio. MAE summarizes and assesses the quality of a machine learning model. The MAE Remote Sens. 2020, 12, 1375 9 of 21 of a model refers to the mean of the absolute values of each prediction error in all the test dataset. Therefore, it refers to the results of measuring the difference between all predictions and observations. RMSE (an absolute measure of fit) is a measure of the average deviation of the estimates from the observed values. It is not scaled to any particular values. RMSE shows how much the predictions deviate, on average, from the actual values. Lower values of RMSE indicate better fit. RMSE is a good measure of how accurately the model predicts the response, and it is the most important criterion for fit if the main purpose of the model is prediction. R 2 (relative measure of fit) varies between 0 and 1. R 2 is the proportion of the variance in the dependent variable that is predictable from the independent variable(s). Therefore, the value near to 1 means that the two variables are perfectly correlated and vice versa.
where P i , O i , O avg , P avg , n and p are the predicted and observed values, average of observed and predicted soil properties values at the ith point, number of data and the total number of explanatory variables in the model, respectively. Moreover, the relative improvement (RI) was calculated to compare the improvement of models and the accuracy of the prediction for C:N. The RI are defined as Equation (5): Model fitting parameter x − model fitting parameter ref Model fitting parameter ref (5) where model fitting parameters are MAE, RMSE and R 2 ; subscript x is the model whose model fitting parameter is higher (R 2 ) and lower (MAE and RMSE) than that of the reference model (least performing model).
We also performed a one-way ANOVA (Analysis of Variance) for testing the significance of the differences among different modeling techniques. For multiple comparison between modeling techniques, Least Significant Difference (LSD) was used. Furthermore, we employed simple linear regression model among different modeling techniques to see the pixel to pixel association of prediction of C:N ratio by them. After selecting the best modeling method based on Landsat TM 5 composite image for 2009-2010, we applied it to Landsat 8 composite image for 2019 for observing the temporal change of C:N.

Models Comparison
The mean MAE, RMSE and R 2 were varied from 2.06 ± 0.62 to 2.43 ± 0.68, 2.88 ± 0.85 to 24 ± 1.15, and 0.37 ± 0.23 b to 0.51 ± 0.22 a, respectively ( Table 2). A significant difference in model fitting parameters was found among different modeling techniques (p < 0.05). However, the post hoc test (LSD) revealed that stochastic gradient boosting (SGB), random forest (RF) and support vector machine (SVM) modeling techniques were not significantly different (p > 0.05; Table 3). However, the model fitting parameters of these three modeling techniques were significantly different from the least performing partial least square regression (PLSR; p < 0.05; Table 3). The relative improvement in terms of MAE, RMSE and R 2 for the SGB, RF and SVM were estimated against PLSR (the lowest performance model) is shown in Table 2. SGB modeling technique showed best performance in terms of relative improvement of RMSE and R 2 while SVM showed best performance in of relative improvement of MAE (Table 3). A similar trend in performance to RMSE and R 2 (relative improvement) was also found for MAE and RMSE percentage of field-based mean C:N (Table 3). Table 3. Comparison of model fitting parameters among four modeling methods. Relative improvement of model accuracy, precision and fitness was calculated against partial against partial least squared regression. Columns with the same letter (a or b; group of modeling methods) are not significantly different based on the LSD test (p < 0.05).

Models
Model As the three model-evaluating parameters were not varied significantly among SGB, RF and SVM, we also plotted the predicted C:N ratio based on these three models against the field-based C:N ratio for testing its fitness in terms of slope and intercept parameters against the 1:1 line (Figure 2). Observed versus predicted scatter plots showed a strong predictability of the SGB model as the slope was deviated only 0.26 from one, and points were concentrated near the 1:1 line (Figure 2a; Table 4). While the SVM had least predictability (Figure 2c; Table 4).
was found among different modeling techniques (P <0.05). However, the post hoc test (LSD) revealed that stochastic gradient boosting (SGB), random forest (RF) and support vector machine (SVM) modeling techniques were not significantly different (P > 0.05; Table 3). However, the model fitting parameters of these three modeling techniques were significantly different from the least performing partial least square regression (PLSR; P < 0.05; Table 3). The relative improvement in terms of MAE, RMSE and R 2 for the SGB, RF and SVM were estimated against PLSR (the lowest performance model) is shown in Table 2. SGB modeling technique showed best performance in terms of relative improvement of RMSE and R 2 while SVM showed best performance in of relative improvement of MAE (Table 3). A similar trend in performance to RMSE and R 2 (relative improvement) was also found for MAE and RMSE percentage of field-based mean C:N (Table 3). Columns with the same letter are not significantly different based on the LSD test (P<0.05). As the three model-evaluating parameters were not varied significantly among SGB, RF and SVM, we also plotted the predicted C:N ratio based on these three models against the field-based C:N ratio for testing its fitness in terms of slope and intercept parameters against the 1:1 line (Figure 2). Observed versus predicted scatter plots showed a strong predictability of the SGB model as the slope was deviated only 0.26 from one, and points were concentrated near the 1:1 line (Figure 2a; Table 4). While the SVM had least predictability (Figure 2c; Table 4).

Importance of Predictors
In different modeling approaches, the importance of the top ten predictors and the type of predictors were different (Figure 3). For instance, in SGB, RF and SVM models top ten predictors were a mixture of different bands, GLCM texture and vegetation indices while in PLSR model top ten predictors consist of GLCM texture of different bands (Figure 3a-d). Moisture and temperature related bands, GLCM texture and vegetation indices were enlisted in the top 10 important predictors in all the models except PLSR where GLCM texture of different bands were in the top ten predictors (Figure 3a-d).

Importance of Predictors
In different modeling approaches, the importance of the top ten predictors and the type of predictors were different (Figure 3). For instance, in SGB, RF and SVM models top ten predictors were a mixture of different bands, GLCM texture and vegetation indices while in PLSR model top ten predictors consist of GLCM texture of different bands (Figure 3a-d). Moisture and temperature related bands, GLCM texture and vegetation indices were enlisted in the top 10 important predictors in all the models except PLSR where GLCM texture of different bands were in the top ten predictors (Figure 3a-d).

Spatial and Temporal Distribution of C:N
The spatial distribution of C:N based on the stochastic gradient boosting (SGB), random forest (RF) and support vector machine (SVM) modeling techniques and its prediction difference between the three-modeling techniques is shown in Figure 4. The difference map and statistics of pixel to pixel association of prediction of C:N ratio by different modeling techniques revealed that SGB and RF predictions were slightly different ( Figure 4; Table 4). However, a large difference was observed in production of C:N ratio between SGB and SVM and RF and SVM (Table 4; Figure 4). Although, the prediction difference of C:N ratio was quite low between SGB and RF modeling techniques, the SGB modeling technique captured variation of C:N ratio that was closer to field distribution compared to RF (Table 5).
Higher C:N ratio was associated in the eastern part of Sundarbans Reserved Forest with a clear homogenous patch compared to the western part and the landward side of the central part ( Figure   Figure 3. Importance value of predictors in each of the four modeling techniques based on 10-fold cross-validation with 10 times repetition. The codes after band (e.g., B6_savg) and other predictors were explained in Table 2. (a) Stochastic gradient booting; (b) random forest; (c) support vector machine; (d) partial least square regression.

Spatial and Temporal Distribution of C:N
The spatial distribution of C:N based on the stochastic gradient boosting (SGB), random forest (RF) and support vector machine (SVM) modeling techniques and its prediction difference between the three-modeling techniques is shown in Figure 4. The difference map and statistics of pixel to pixel association of prediction of C:N ratio by different modeling techniques revealed that SGB and RF predictions were slightly different ( Figure 4; Table 4). However, a large difference was observed in production of C:N ratio between SGB and SVM and RF and SVM (Table 4; Figure 4). Although, the prediction difference of C:N ratio was quite low between SGB and RF modeling techniques, the SGB modeling technique captured variation of C:N ratio that was closer to field distribution compared to RF ( Table 5).
As the SGB modeling technique covered much spatial variation of C:N and were similar to field C:N, we applied the SGB based on Landsat 8 composite image for 2019 to observing the temporal change of C:N (Table 6; Figure 5). The difference between 2009-10 and 2019 C:N showed that C:N ratio increased mainly at the central and landward parts of Sundarbans Reserved Forest ( Figure 5). In the seaward and eastern side, it showed a decreasing trend in C:N ratio ( Figure 5).    Higher C:N ratio was associated in the eastern part of Sundarbans Reserved Forest with a clear homogenous patch compared to the western part and the landward side of the central part ( Figure 4). The streamline zone (near to river and canal bank) had lower C:N ratio compared to stable stands (forest proper zone; 100 to 200 m away from river and canal bank) ( Figure 4).
As the SGB modeling technique covered much spatial variation of C:N and were similar to field C:N, we applied the SGB based on Landsat 8 composite image for 2019 to observing the temporal change of C:N (Table 6; Figure 5). The difference between 2009-10 and 2019 C:N showed that C:N ratio increased mainly at the central and landward parts of Sundarbans Reserved Forest ( Figure 5). In the seaward and eastern side, it showed a decreasing trend in C:N ratio ( Figure 5).

Discussion
To our current knowledge, this study assessed, for the first time, the spatial distribution and temporal change of the CWM of senescent leaf C:N ratio in Sundarbans Reserved Forest, Bangladesh as a model in mangrove ecosystem with remote sensing. Four machine learning modeling approaches were compared in order to select the best modeling technique for mapping C:N ratio using spectral and textural properties of seven bands and vegetation indices of Landsat TM 5. The modeling approach that combined predictors from the above-mentioned Landsat TM 5 image properties that

Discussion
To our current knowledge, this study assessed, for the first time, the spatial distribution and temporal change of the CWM of senescent leaf C:N ratio in Sundarbans Reserved Forest, Bangladesh as a model in mangrove ecosystem with remote sensing. Four machine learning modeling approaches were compared in order to select the best modeling technique for mapping C:N ratio using spectral and textural properties of seven bands and vegetation indices of Landsat TM 5. The modeling approach that combined predictors from the above-mentioned Landsat TM 5 image properties that also related to moisture and temperate, showed better performance compared to the modeling approach that used only textural properties. The SGB modeling technique captured wide variation of C:N ratio and used a lower number of predictors compared to the other three modeling methods. Thus, this modeling technique was applied to the Landsat 8 image in 2019 for assessing temporal change of C:N ratio.

Comparison of Model Performance
The nonlinear based modeling technique showed better performance in terms of lower MAE, RMSE and higher R 2 (Table 3). This difference in performance among different modeling approaches could be due to several reasons [47,87]. First, how a specific modeling approach arranges and handles the highly correlated predictors in the top rank of importance and had high correlation with C:N. Although, all the modeling approaches used in this study have the capacity to handle the multicollinearity issue, the SGB showed strong ability to handle multicollinearity over other models [47,87]. In a study of biomass modeling using RapidEye satellite with SGB and RF in Africa terrestrials forest, Dube [50] reported the outperformance of the SGB modeling approach. PLSR showed a weak performance. A similar performance of PLSR modeling technique in comparison with other machine learning modeling approaches has also been reported in biomass modeling studies [88][89][90]. This finding could reveal that C:N ratio and remote sensing-based predictors are nonlinearly related in the Sundarbans Reserved Forest. This could not be handled well by PLSR modeling technique which is linear based [46].
The second reason is the combination of predictors from different dimensions (e.g., either band spectral properties or GLCM textural properties or vegetation indices or all of them) in modeling. The modeling approaches that aggregate both spectral and textural properties have higher predictability compared to the methods that aggregate only a single type of predictor [23,91]. The top three modeling approaches (SGB, RF and SVM) combined the most important predictors from band spectral and GLCM textural properties and vegetation indices compared to PLSR modeling technique, which mainly aggregated the GLCM texture properties of different bands (Figure 3). This incorporation of homogeneous predictors in modeling resulted in the poor performance of PLSR technique.
Finally, the type of predictors in the top rank of important value in different modeling methods are related to vegetation, canopy moisture, salinity and water content. Among the top rank predictors, Band 6 is related to moisture and temperature while the Disease Water Stress Index and Modified Normalized Difference Water Index are related to water stress and canopy moisture, respectively [92][93][94][95]. Previous studies also confirmed that in mangrove, nutrient status and dynamics vary along with water, temperature and salinity gradient [3,7]. Therefore, the modeling technique that had a combination of those predictors in the top rank of importance value had lower errors in prediction and higher capacity of explaining the variance of C:N (Table 2; Figure 3).

Spatial Distribution and Temporal Change of C:N ratio in SRF
Among the top three modeling approaches, the predicited C:N map stochastic gradient boosting (SGB) technique was closer to the field-based C:N distribution ( Figure 4; Table 6). Although, the mean and maximum values of C:N in the predicted map of C:N by the top three models was closer to field value, the SGB technique predicted the lower value of C:N better than RF and SVM ( Table 6). The spatial distribution of C:N ratio in Sundarbans Reserved Forest (SRF) was closely related to vegetation types and salinity gradient. For example, the higher C:N value (>36) in the eastern part of SRF is mainly dominated by Heritiera fomes (Figure 4) [16,22,96]. The higher C:N ratio in these areas indicates that Heretiera fomes dominated stand has higher photosynthetic nitrogen use efficiency or nutrient resorption capacity [3,16,97]. However, most of the western parts of SRF were under low C:N ratio coverage, which indicates species in those areas are in stress conditions with low photosynthetic nitrogen use or have lower nutrient resorption capacity and thereby there may be low productivity or carbon sequestration capacity [16]. These areas are mainly dominated by Ceriops decandra and Excoecaria agallocha with scattered Xylocarpus granatum, X. mekongensis [16,96]. A higher level of salinity may cause enzymatic metabolism disorder in the chloroplast; decline in K + ion, the key photosynthetic element and other philological activities, thus decreasing the rate of photosynthesis [97]. For this reason, photosynthetic nitrogen use efficiency could be lower in the high saline zone and streamlined areas in the eastern part of SRF, which could be the cause of low C:N ratio in those zones. The lower range of C:N ratio was also found in the 2010 map in the landward (the northern side of SRF), although this area is associated with freshwater or moderate saline water coverage. The reason of this lower C:N value could be the dominance of Avicenna officinalis species, which has a lower C:N ratio (24.47; Table 1) [13].
The temporal change map (difference between 2009-2010 and 2019) indicates that a gradual increasing trend was observed from seaward to landward (central and western part of SRF; Figure 5). Similar trend also observed at streamline to forest proper line ( Figure 5). This could be due to the change of species composition over the 10-years period where species with higher C:N may be shifted. In most cases it was H. fomes because other species (A. corniculatum, R. mucronata) are mainly distributed in high saline area (south-western part of SRF) and some also confined in waterlogged areas (B. sexangular) in SRF. However, temporal change map also shows that across the seaward and eastern part of SRF the C:N ratio was declined ( Figure 5). Two reasons may be caused this declination of C:N: (i) a super cyclone (Bulbul) passed over SRF on 9 November 2019 which severely damaged forest canopy and (ii) may be the change of species composition where species with low C:N ratio may come to canopy layer.

Implication for Forest Health Monitoring and Policy
In the current climate change regime, monitoring of health conditions is a key forest management task [22,33,98]. Litter quality such as C:N ratio as an essential biodiversity variable or functional trait can be used as a potential indicator of decomposition rate, plant nitrogen retention capacity, primary productivity and soil organic carbon storage [6,9,10,99,100].
We developed advance machine learning models for mapping senescence leaf C:N in SRF using Landsat TM 5 imagery and applied the best model on Landsat 8 imagery for 2019. The C:N ratio map developed in the current study can be used for improving the biogeochemical and ecosystem models in the mangroves. In this way, our developed method will contribute to evaluable ecosystem service policy and forest management at the local and regional levels [53,101]. For example, it can be used to track the endangered Red list species in Mangrove that were enlisted by the International Union for Conservation of Nature. The higher value of C:N ratio in the eastern parts of SRF is dominated by H. fomes, the endangered species in SRF. Therefore, by tracking C:N ratio as an indicator of Essential Biodiversity Variable, we can assess the progress of Biodiversity targets that related to IUCN Red List species. For example, two of 20 specific targets (targets 9 and 12) of Aichi Biodiversity Targets-2020 are related to endangered species. The aim of these biodiversity targets is to conserve the threatened species and update their conservation status as to whether they are being improved and sustained by 2020 [102,103]. Furthermore, our developed method directly contributes to monitoring the progress of the biodiversity targets (Aichi Targets 2, 3 and 15), which are related to nutrient dynamic, productivity and ecosystem health [102,104,105]. The C:N ratio is directly associated with key ecosystem functions (decomposition) that have been regarded as biodiversity indicators such as net primary productivity, nutrient retention, disturbance regime [102,104,105], helping in this way to achieve National Biodiversity Target 2 of Bangladesh [106]. While the usage of Landsat imagery in our modeling technique is related to National biodiversity Target 19 of Bangladesh which is aimed at developing a sustainable monitoring technique by means of remote sensing [106]. Ultimately, our finding will be contributed towards achieving sustainable development goal 15 and Aichi Biodiversity target 2, 11 and 15 [51,[106][107][108] for Bangladesh.

Conclusions
The current study investigated the potentiality of Landsat imagery for spatiotemporal mapping of CWM of C:N ratio of senescent leaf by comparing four machine learning techniques in SRF. SGB technique showed the best performance by incorporating key predictors from the spectral and textural properties and vegetation indices group, which had higher relative improvement against the least performing method (PLSR) compared to RF and SVM. The spatial pattern of C:N ratio map based on SGB method captured a wide variation of C:N, which was closer to field C:N ratio. The map of C:N ratio reveals that the spatial distribution of C:N mainly depends on the species type and salinity gradient. Stands close to the streamline are in a more stressed condition than forest proper or inland side. The eastern side of SRF, which is mostly dominated by H. fomes and receives more freshwater from the upstream river, had a higher C:N ratio, and thereby is in less stressed condition with higher nitrogen use or retention efficiency compared to the western and landward sides of SRF. The temporal change of C:N ratio could also be influenced by species composition change and cyclone damage. Our developed advance machine learning technique thus can be applied for future prediction of litter quality in SRF and for improving the biogeochemical and ecosystem models related to nutrient dynamic in the mangroves.
Supplementary Materials: The following are available online at http://www.mdpi.com/2072-4292/12/9/1375/s1, Pearson correlation between forty-eight predictors and community weighted mean of C:N ratio of senescent leaf used in modeling. The first column represents the type of different predictors (e.g., Surface reflectance). Other column and row represent the different variables (Dependent: CWM C:N, Predictors: all other variables) with a value of either no sign before (positive relation) or negative sign (negative relation between variables). The full names of the variables were given below at the end of Pearson correlation with type (yellow highlight).

Conflicts of Interest:
The authors declare no conflicts of interest.