Impact of Land Use and Land Cover Changes on Urban Ecosystem Service Value in Dhaka, Bangladesh

: Urban ecosystem services provide a wide range of services to sustain life, social relation, health, etc., and address most of the challenges, including climate change and environmental pollution. While it is recognized that the urban ecosystem substantially contributes to human well-being in cities, there is less attention to consider the value of urban ecosystem service in urban planning and policymaking. This study analyzed the land use and land cover (LULC) dynamics of city of Dhaka over the past three decades (1990–2020) to evaluate the impact of LULC on ecosystem services value (ESV). The estimation of ESV in relation to LULC has been done using the globally used benefits transfer method (BTM). Findings of the study show that built-up area has increased by 188.35% from 1990 to 2020, with an average annual growth rate is about 6.28%. The analysis of ESV shows that it has decreased by 59.55% (85 million USD) from 142.72 million USD in 1990 to 57.72 million USD in 2020 due to the development of the built-up area through conversion of agricultural land, waterbodies, and forest and vegetation land. This study also identified that waterbodies are the greatest contributor to ESV. The result on the elasticity of ESV in relation to LULC implies that about 1% transition in LULC would result in about 0.33% change in total ESV during the study period. We believe that the findings of this study would serve as a reference for the policy maker and urban planner to devise appropriate land use decision to ensure sustainable urban development of Dhaka.


Introduction
Ecosystem services (ES) can be defined as the conditions and processes through which natural ecosystems, and the species that make them up, sustain and enrich human life [1]. Ecosystem services can be considered as the benefits human populations derive, directly or indirectly, from ecosystem functions [2,3]. Ecosystem services are categorized into four major groups: (a) provisioning services (e.g., food production); (b) regulatory services (e.g., climate regulation); (c) supporting services (e.g., pollination); and (d) cultural services (e.g., recreation). These ecosystem services provide many functions to human life [4]. Costanza et al. [3] identified 17 functions that are derived from ecosystem services necessary for the survival and well-being of human life. Urban ecosystem services play an essential role in connecting cities with the biosphere and to reduce the ecological footprint and ecological debt of cities and at the same time enhance urban public health, resilience, and quality of urban life [5]. Urban ecosystem services contribute directly or indirectly to human well-being providing many functions, including food supply, water supply, waste treatment, regulation of the urban heat island effect, clean air, water filtration, noise reduction, pollination, and climate regulation, etc. [5][6][7][8]. A healthy markable LULC changes [35][36][37]. It is argued that unplanned changes in LULC will dismantle the urban ecosystem, leading to environmental degradation, air and water pollution, the urban heat island effect, urban disturbance, social and economic burden, all hampering human well-being and public health [38,39]. For example, in past decades, Dhaka has developed a large amount of built-up area at the expense of wetland. Due to this transformation, the natural landscape was disturbed, and as a result, Dhaka faces waterlogging and rainfall-induced urban flooding, which cause enormous problems to urban people. Similarly, the present scenario of air pollution and water pollution in Dhaka is the consequence of an endangered urban ecosystem [40][41][42].
In response to the above, some researchers studied the nexus between LULC change and ecosystem services in Bangladesh. For example, Hoque et al. [43] studied the future impact of land use/land cover changes on ecosystem services in the lower Meghna river estuary of Bangladesh. Their findings revealed that rural and urban settlement in the lower Meghna estuary had increased at the expense of agricultural and mangrove forest lands, leading to a total ESV decrease of US$105.34 million during 1988-2018. Akber et al. [44] assessed the impact of land use change on ecosystem services (ES) of southwest coastal Bangladesh. Their results revealed that agricultural land decreased by 253,928 ha, resulting in a loss of ESV US$ 1.41 billion during the period 1980-2016. The study of Huq et al. [45] analyzed the interactions between freshwater ecosystem services and land cover changes in southern Bangladesh. Their results showed that during the period 1973 to 2014, substantial land cover transformations occurred, e.g., an increase of agricultural land, rural vegetation with settlement (RVS) in exchange for wetlands resulting in USD 1.06 billion worth of freshwater ecosystem services loss. The above studies mainly focused on coastal areas of Bangladesh but not on the cities. Zinia and McShane [46] evaluated urban ecosystem services in Dhaka on the basis of household surveys, field-level observation, interviews, and discussion with experts and local people. They mainly identified and mapped available urban ecosystems and conducted inventory analysis. Their work did not consider the impact of LULC change on ESV.
According to the literature, there is insufficient work evaluating the impact of LULC changes on urban ecosystem services, especially the use of satellite image-derived LULC change for urban ecosystem service assessment is very limited locally and globally. Although there are some studies conducted in urban areas or cities globally, the LULC characteristics of those urban areas or cities are different from that of Bangladesh. For example, the study of Lin et al. [30] examined the influence of LULC changes on ecosystem service in Chengdu city, where the agricultural land and built-up area occupied about 51.86% and 18.16% of the total land respectively in 2018. On the other hand, Dhaka occupied 8.08% and 51.85% of Bangladesh's agricultural land and built-up area respectively in 2020. Although there are many separate studies on LULC change and ecosystem services valuation in Bangladesh, there is limited research that focuses on the influence of LULC change on urban ecosystem services. According to the literature, there is no previous study in Dhaka analyzing the impact of LULC change on ecosystem services. Thus, this study is the first of its kind to assess the impact of LULC changes on ecosystem services using satellite images. In the view of above, this study aims (1) to analyze the LULC change dynamics of Dhaka during the period from 1990 to 2020; (2) to estimate and map the total ESV over the study period; (3) to assess the value of individual ecosystem service functions; and (4) to assess the impact of LULC change on ESV through determining the elasticity of ESV in response to LULC change.

Study Area
The study was conducted in Dhaka, Bangladesh. The city is centrally located in Bangladesh and lies between 23.67 0 and 23.90 0 N in latitude and 90.33 0 and 90.51 0 E in longitude, and has an area of 305.82 Km 2 . The study area covers the Dhaka metropolitan police (DMP) area ( Figure 1). The topography of the city is predominantly flat, having a mean surface elevation of 15 m. The city is surrounded by three major rivers, namely the Turag to the northwest, Buriganga-Dhaleshwari to the south, and Shitallakhya-Balu to the east. Dhaka is the capital of Bangladesh and the center of the national government, trade, and culture. It is the largest city in Bangladesh and in Eastern South Asia. It is also considered as one of the major cities in South Asia due to its urban character and socioeconomic and political diversity. Dhaka is one of the largest and densely populated megacities of the world [47]. The city is the home of more than 11 million people within a total land area of about 305.82 km 2 . According to Demographia World Urban Areas, as of May 2021, Dhaka is the firstranked high-density city in the world, having 3 6,941 people per square kilometer. The population growth rate in Dhaka is about 3.5% in 2021, which is one of the highest growth rates amongst the most Asian cities [48]. This high level of population growth rate is the reflection of rural-urban migration in Dhaka. The vibrant culture, business opportunities, income opportunities, health, and education facilities, etc., have contributed to the increasing migration and population growth in Dhaka.
Dhaka experiences a tropical wet, dry, and monsoon climate, having an average temperature of 25 °C (77 °F) and monthly means varying between 18 °C (64 °F) in January and 29 °C (84 °F) in August. The annual rainfall of Dhaka ranges between 1169 mm to 3028 mm, and the yearly average rainfall is about 2076 mm. About 80% of the average annual rainfall of 1854 mm occurs through the monsoon period, spanning from May to the end of September [49]. Due to the ever-increasing population in Dhaka, the air and water are becoming polluted beyond the acceptable level, wetland and green space are being occupied by multistory buildings and real estate developments to fulfill the needs of the increasing population, which, in turn, threatened the overall environment, biodiversity, and urban ecosystem of Dhaka.

Datasets
This study focuses on the valuation of urban ecosystem services based on land use land cover (LULC) change. Thus, LULC is the core data for this study, and this LULC data has been derived from satellite images. This study uses time-series datasets from Landsat-5, Landsat-7, and Landsat-8. Landsat satellite images for the years 1990, 1995, 2000, 2005, 2010, 2015, and 2020 were downloaded from US Geological Survey (USGS) official website (https://earthexplorer.usgs.gov/). Table 1 presents the detailed description of Landsat images used for this study. The whole study area lies at the intersection of Landsat path 137 and row 44. The downloaded images had a built-in projection system of the Universal Transverse Mercator (UTM) projection within Zone 46 North based on the World Geodetic System (WGS)-1984 datum. The spatial resolution of the images used for this study is 30 x 30 m per pixel. Since the temporal variation might affect the spectral reflectance of the Earth's surface, similar dated images were downloaded to avoid such variation. We selected Landsat data for the month of November in each year except 1990 because the cloud-free images were available in the month of November. In addition, there was no rain during November, this ensures the actual presence of waterbodies. If images were taken during rainy season, some low land could have been identified as wetland/waterbodies although those were not actual waterbodies. Therefore, we had selected data for the month of November. In 1990, there were no cloud-free images available for the month of November, so in that case, we took image for the month of December.

Landsat Image Processing
Landsat data product has a different level of processing. In this study, we have downloaded Level-1 precision and terrain corrected (L1TP) and cloud-free multispectral images. L1TP images are radiometrically calibrated and orthorectified using ground control points (GCPs) and digital elevation model (DEM) to correct for relief displacement. The highest-quality level-1 products are suitable for pixel-level time series analysis. L1TP is geometrically corrected; hence no further processing for geometric correction is required for Landsat images [50]. The raw images contain a digital number (DN) value for each pixel. However, these DN values do not reflect the actual earth surface characteristics due to atmospheric effects and solar angle. The reflected or emitted electromagnetic radiation (EMR) from features on the Earth's surface area is influenced by distortions due to the sensor, solar, atmospheric, and topographic effects. Thus, the satellite sensor does not receive the actual reflected signal from the earth's surface [51]. Therefore, preprocessing is required to minimize these effects, depending on remote sensing applications. Most preprocessing focuses on radiometric and atmospheric correction to improve the accuracy of surface reflectance, emittance, and back scattering measurements. Many atmospheric models are available to correct the radiometric error, such as 6S (second simulation of the satellite signal in the solar spectrum) and MODTRAN (moderate resolution atmospheric radiance and transmittance), atmospheric correction for flat terrain (ATCOR2), fast lineof-sight atmospheric analysis of spectral hypercubes (FLAASH), simple dark object subtraction (SDOS), dark object subtraction (DOS), and cosine estimation of atmospheric transmittance (COST) [52]]. In this study, we have applied COST model for radiometric correction of images because the COST model is entirely based on image characteristics and was developed to account for multiplicative effects of atmospheric scattering and absorption [51]. The COST model also performs well compared to other models [53] and is available in many image processing software. Three-step procedures were performed to convert DN values to surface reflectance through radiometric correction. These steps are a) conversion of DN values to spectral radiance at the sensor, b) conversion of spectral radiance to reflectance at the sensor, and c) atmospheric correction and conversion of sensor reflectance to surface reflectance. We a applied COST model for radiometric correction using the radCor() function of RStoolbox in R software [54].

Classification of Land Use and Land Cover (LULC)
The surface reflectance value acquired after preprocessing through a radiometric and atmospheric correction in the earlier step was used to derive LULC maps. We used an unsupervised classification method to classify Landsat images into five broad LULC classes. These LULC classes are a) built-up area; b) vegetation; c) bare land; d) waterbody; e) agricultural land. The definition of these LULC classes are presented in Table 2. These LULC categories have been selected considering biome types used to calculate ecosystem value. We used the unsupervised classification method because this method has some advantages over the supervised classification method. This method is comparatively quick and easy to run and does not require extensive prior knowledge about the area. In the unsupervised method, classes are created completely based on spectral information of the images and, thus, are not subject to visual interpretation [55]. In this method, the interpreter has to define the number of classes, and the software automatically groups the pixels of similar spectral properties into the number of classes previously defined by the user [56]. One of the disadvantages of this method is that spectral reflectance does not always correspond to the informational category. This mostly happens to the pixels having values like two adjacent classes. This problem can be solved by determining the optimal number of classes through an iterative process [57].
Different clustering algorithms are available to classify land cover using an unsupervised method. Popular clustering algorithms include interactive self-organization data analysis (ISODATA), clustering large applications (CLARA), mean-shift clustering, agglomerative hierarchical clustering (AHC), partitioning around medoids (PAM), random forest, and k-means clustering, etc. [58][59][60][61]. Amongst, K-means clustering is one of the most commonly used unsupervised clustering algorithms used to classify LULC. Considering its popularity, we have used a k-means clustering algorithm using R software [54] to classify processed Landsat images into five LULC classes. In K-means clustering, data are partitioned into a set of k groups or clusters where k represents the number of clusters predefined by the analyst. K-means clustering classifies data in such a way that data within the same cluster are as similar as possible, whereas data from the different clusters are as dissimilar as possible [62].
As mentioned earlier, determining the optimal number of clusters is a fundamental issue in k-means clustering. Because the performance and accuracy are associated with determining the number of clusters. An optimal number of clusters prevents misclassification of the dataset. There are several methods to determine the optimal number of classes. Among them, three methods are popular; these are elbow, silhouette, and gap statistic methods [63]. We used the elbow method to determine the optimal number of clusters and identified that eight is the optimal number of classes for our dataset, although our required number of LULC classes is five. So, we first applied k-means clustering with k=8, and then we regrouped the class into five LULC classes.
'Salt-and-pepper' noise was observed in the classified images. Some misclassifications were identified in the classified images. Since Landsat images have low resolution, such misclassification and 'salt-and-pepper' noises are common. Hence, some degree of postprocessing is required to eliminate 'salt-and-pepper' and misclassification. To remove the 'salt-and-pepper' noise, we applied a 3 x 3 majority filter in ArcGIS 10.8 software. The classified raster images were then vectorized to manually rectify the localized error due to misclassification. This output has been used to create the final version of LULC maps for the different time periods. After removing the localized error, the vectorized images were converted to raster images again for accuracy assessment.
Accuracy assessment was done using the Google Earth platform. Classified images were compared to corresponding land cover in Google Earth for accuracy assessment. It is recommended that about 50 random points for each land use class are required to assess the classification accuracy if the study area is less than 1 million acres and the land cover category is less than 12 classes [64]. In our case, the size of the study area was 75569.77 acres, and the number of LULC classes was five. The minimum sampling point is 250 for the assessment of classification accuracy. We generated 300 random sampling points using the stratified sampling technique in ArcGIS 10. In addition to quantifying the LULC type, we also analyzed the LULC change over the study period since LULC changes influence ecosystem services and their function. We used the land-use dynamic index to understand the LULC changes over time [65]. Landuse dynamic index was calculated using the following Equation (1).
where K indicates the land-use dynamic index for a single LULC type, and are the initial and final area of specific LULC type, respectively, and T is the study period.

Assignment of Ecosystem Services Values (ESV)
Ecosystem service valuation has received much attention, and it has been used globally to understand the multiple benefits offered by ecosystem services. A lot of studies have been conducted to estimate the value of ecosystem services in the past few decades, including the economic valuation of tropical forests, endangered species management, and protected areas. However, the ecosystem service valuation model presented by Costanza et al. [3] is considered as the robust method to estimate the economic value of ecosystem services. They identified 17 ecosystem services from 16 biomes in order to calculate ESV. Using the benefits transfer method (BTM), they estimated global economic values for those ecosystem services based on existing literature and original calculations. Since then, this method has been used by many studies globally to assess ecosystem services. Although the method presented by Costanza et al. [3] received some critiques due to low estimate of farmland unit value and uncertainties, the model presents the most comprehensive set of initial approximations for ecosystem service valuation. Therefore, we used this method in our study for the valuation of ecosystem services for Dhaka. To obtain the values of ecosystem services for each of the five LULC categories in Dhaka, we compared the five categories of LULC with 16 biomes identified in Costanza et al. [3]. The most relevant biome for each category was assigned as the proxy for that LULC type. For example, 'cropland' biome for 'agricultural land', 'wetland' for 'water body', 'tropical forest' for 'forest and vegetation', 'urban' for 'built-up area', 'desert' for 'bare land' were assigned for the proxy to calculate the value of ecosystem service for each land use type. In our case, most of the waterbodies are wetlands, and a negligible portion belongs to shallow lakes/river. Therefore, we assigned waterbodies as 'wetlands' instead of lakes/rivers. Although the LULC categories and equivalent biomes do not perfectly match, their use proved feasible in many other studies [17,[66][67][68]. LULC types, equivalent biome, and corresponding value coefficient are presented in Table 3.

Calculation of Ecosystem Services Values
Based on the value coefficient in Table 3, the total value of ESV was calculated by multiplying land area with value coefficient according to the following Equation (2) [69,70].
Where ESV is the total estimated value of ecosystem service, is the area in hectare and is the value coefficient (US $ ℎ ) for each LULC category 'k'. In order to determine the change in ecosystem service value, we calculated the differences between the estimated ESV for each LULC category for the year 1990, 1995, 2000, 2005, 2010, 2015, and 2020. We estimated the ESV change rate using the following Equation where refers to the yearly ESV change rate for a single LULC type, and are the initial and final ecosystem service value, and T is the study period. In addition to calculating the ESV for each LULC category and estimating the effect of LULC change on total ESV, we also estimated the value of 17 ecosystem functions provided by each biome using the following Equation (4) Where, is the estimated value of ecosystem service function f, is the area (ha) and is the value coefficient of function f (US $ ℎ ) for each LULC category 'k'. The value coefficient for each individual function was obtained from Table 2

Calculating Elasticity of ESV due to LULC Changes
Elasticity describes how well one variable responds in relation to the changes in another variable. To understand how total ESV changes due to the change of LULC, we calculated the elasticity of ESV. The elasticity of ESV in relation to LULC changes measures the percentage change of ESV in response to the percentage change in LULC. The elasticity of ESV was calculated using the following Equations (5) and (6) In the above equations, EEL is the elasticity of total ESV in relation to LULC changes, and are the initial and final ecosystem service value, and T is the study period. LTP is the percentage of land transition which implies the degree of LULC changes, △ is the changes of LULC category during the study period and is the area of lULC category of .

Sensitivity Analysis
In the early step, we made a similarity check of our five LULC categories with 16 biomes identified by Costanza et al. [3], and the most similar biome was assigned to the corresponding LULC category as the proxy, but the biomes used as the proxies are not perfect matches of each LULC category; thus, there exist uncertainties in the coefficient values used to calculate ESV. Hence, sensitivity analysis is necessary to determine the level of dependence of total ESV due to a change in the coefficient value of a specific LULC type. Coefficient of sensitivity (CS) is a measure to verify the effective representation of biomes and the accuracy of the corresponding value coefficient in various LULC categories. We used the standard economic concept of elasticity to calculate CS as per the following Equation (7) where ESV and VC refer to estimated total ecosystem service value and value coefficient, respectively; i and j represent the initial and adjusted value, and k represents the specific LULC type. The value coefficient of each LULC category was adjusted by ±50% to calculate the adjusted ESV and CS. By calculating the value of CS, we can understand the dependence level of change in the total ESV due to the change in the value coefficient of a specific LULC type. If CS is less than 1, then the estimated ESV is inelastic with respect to value coefficient, and if CS is greater than 1, the calculated ESV is considered to be elastic. If the estimated ESV with respect to value coefficient is inelastic, then the representation of proxies and associated estimated ESV are reliable. In many previous studies, sensitivity analysis was widely used to understand the change in ESV due to change in value coefficient [71][72][73][74].

Land Use Land Cover Change
Spatial distribution and trend of LULC changes of Dhaka for the last three decades have been presented in Figures 2 and 3. Figure 2 illustrates that, in 1990, the built-up area was spatially distributed along the southern part of Dhaka, which extends towards the northwest direction over the past 30 years. Initially, agricultural land was distributed over the northern part of the city, but in the course of time some portion of the forest and vegetation area in the southeast part was converted into agricultural land, and later most agricultural land was converted to the built-up area. In 1990, forest and vegetation was distributed in the central part of the city, but most parts of it have been replaced by built-up area in the past 30 years. Figure 2 illustrates that most parts of the different LULC types were converted to the built-up area during the period 1990 to 2020. This transition is the result of the high urbanization rate in Dhaka.    We also calculated the land use change dynamic index (k) according to equation 1 to understand the annual rate of changes of different LULC types in the different study periods. Table 5 shows the value of k in which the "-" sign indicates the decrease and others indicate an increase in LULC changes. From Table 5, it can be seen that agricultural land, waterbodies, forest and vegetation annually decreased by 2.06%, 1.98%, and 2.05% respectively from the period 1990 to 2020. On the other hand, built-up area and bare land annually increased by 6.28% and 4.61%, respectively, during the period 1990 to 2020. It is also observed that the annual decreasing rate of agricultural land use change (4.66%) was highest during the period 1990-1995 and that of waterbodies (6.49%) and forest and vegetation (5.65%) during the period 2015-2020. The annual increase of built-up area change (11.01%) was highest during the period 1990-1995 and that of bare land (6.71%) during the period in 2015-2020. Figures 2,3, and Table 4 show the spatial distribution and trend of LULC changes but do not provide the idea of how a particular LULC type was redistributed over the period. To understand this redistribution mechanism, we created a transition matrix using ArcGIS 10.8 software, which is presented in Table 6. This transition matrix shows how the different LULC types of the study area were redistributed during the period from 1990 to 2020.  Table 6 shows that built-up area occupied 9.18% (28.07 km 2 ), 5.46% (16.70 km 2 ), 15.58% (47.65 km 2 ), and 5.39% (16.48 km 2 ) area from agricultural land, waterbodies, forest and vegetation, and bare land during the period from 1990 to 2020, whereas it lost only 1.71% from its own. It can also be noted that, except waterbodies, the major portion of different LULC types was converted to the built-up area; a major portion of waterbodies (11.52%) was converted to bare land. Table 6 also implies that the forest and vegetation and waterbodies were the highest contributors (15.58% and 11.52%, respectively) of builtup area and bare land increase.

Change in Ecosystem Service Value
In our study, we estimated the ESV of Dhaka for each LULC type based on Equation (2), ecosystem service value coefficient in Table 7, and area of different LULC types in Figure 3. Please note that we converted the measurement unit of land area from km 2 to hectare since ecosystem service value coefficients were calculated based on the later unit of measurement [3]. The result of the estimated ESV has been presented in Table 7, Table  8, Figure 4, and Figure 5. According to result, as presented in Table 7 and Figure 4 Table 7 clearly shows that waterbodies were the highest contributor to total ESV, followed by forest and vegetation, and agricultural land. On the other hand, built-up areas and bare land have no contribution to ESV.   According to Figure 4, the ESV trend of Dhaka shows a declining change process. A gradual decrease in ESV was observed for forest and vegetation, and agricultural land over the study period, and a sharp decrease in ESV was observed during the period 1990-2000 and 2010-2020. Since waterbodies were the highest contributor to total ESV, the characteristics of the total ESV change process are dominated by the characteristic change process of ESV derived from waterbodies.
To understand the yearly changes of ESV, we calculated the rate of annual changes for ESVs according to Equation (3). Table 8 shows the annual rate of changes of ESV for different LULC types during the period 1990 to 2020. The negative sign in the values represents the decrease in ESV. Table 8 shows the annual rate of ESV for agricultural land, waterbodies, forest and vegetation only; ESV changes for the built-up area and bare land are not present because they derive no ESV. It can be seen that ESV derived from agricultural land, forest and vegetation, and waterbodies annually decreased by 2.06%, 2.05%, and 1.98%, respectively, from the period 1990 to 2020. It can also be observed that annual change in ESV from agricultural land use (4.66%) was highest during the period 1990-1995 and that of waterbodies (6.49%) and forest and vegetation (5.65%) during the period in 2015-2020. ESV derived from agricultural land decreased by 0.13 million USD during the period from 1990 to 1995. From 2015 to 2020, ESV derived from waterbodies and forest and vegetation reduced by 24.89 and 2.3 million USD, respectively. Table 8 shows that the annual rate of changes in total ESV of Dhaka were −2.56%, −1.65%, −0.89%, −1.37%, and −3.

Estimated Value of Individual Ecosystem Service Function
In this study, we estimated the value of 17 individual ecosystem functions based on Equation (4), the area of different LULC types in Figure 3 , and the value coefficient of individual ecosystem functions provided in Table 2 in Costanza et al. [3]. In addition, we calculated the impact of LULC changes on the change of individual ecosystem function during the period 1990 to 2020. The values of 17 individual ecosystem functions are presented in Table 9. These 17 ecosystem services are grouped into four broad categories: provisioning, regulating, supporting, and cultural. The value of individual ecosystem function helps to understand the contribution of each ecosystem function to overall ESV. As shown in Figure 6, the regulating services are the major contributor to total ESV in the study area, followed by the provisioning, cultural and supporting services. The regulating services contributed about 56.43%, 56.34%, 55.91%, 56.42%, 56.85%, 56.80% and 56.60% of the total ESV in the year 1990, 1995, 2000, 2005, 2010, 2015, and 2020. In the case of provisioning services, they showed the highest contribution (27.42%) in 2020. Compared to regulating and provisioning services, the contribution of supporting services and cultural service to the total ESV is low: on average 6.79%, 9.37%, respectively.  According to Table 9, water supply, waste treatment, and disturbance regulation are three major contributors and share 22.90%, 25.60%, and 27.30% of total ESV, respectively, in 1990. Due to the large quantity of waterbodies and higher value of the coefficient, these three ecosystem functions generate higher ESV and contribute to the major portion of total ESV. Change in the value of individual ecosystem function were calculated during the period 1990-1995, 1995-2000, 2000-2005, 2005-2010, 2010-2015, 2015-2020 and are presented in Table 10.  Tables 9 and 10, although the percentage of contribution to total ESV of major three ecosystem functions remained similar during the period 1990 to 2020, the quantity of the major ecosystem functions decreased by 59.31% (19.38 million USD), 59.34% (21.67 million USD), and 59.31% (23.13 million USD) respectively in 2020. The highest declines in the value of ESV for water supply (32.43%), waste treatment (32.37%), and disturbance regulation (32.43%) were observed during the period 2015-2020. In the case of supporting and cultural groups of ESV, nutrient cycling function and cultural function are the major contributors to the total ESV. They contributed about 4.89% (6.98 million USD) and 5.31% (7.58 million USD) to the total ESV in 1990. However, their contributions declined by 61.57% (4.30 million USD) and 59.31% (4.49 million USD) during the period 1990 to 2020.

Impact of LULC Change on ESV
To understand the impact of LULC changes on the total ESV and their relationship, we have calculated the elasticity of ESV in response to LULC changes. The elasticity of ESV in response to LULC changes was calculated using Equation (5) and is presented in Figure 7. The figure shows that the elasticity of ESV in response to LULC changes were 0.29, 0.40, 0.15, 0.27, 0.65, 0.79, and 0.33 during the period, 1990-1995, 1995-2000, 2000-2005, 2005-2010, 2010-2015, 2015-2020, 1990-2020 respectively. The elasticity implies that 1% transition of LULC would result respectively in on average 0.29%, 0.40%, 0.15%, 0.27%, 0.65%, 0.79%, and 0.33% change of total ESV during the study period ( Figure 7). As illustrated in Figure 7, the elasticity of ESV during the different time periods over 1990 to 2020 shows a fluctuating trend. The elasticity of ESV increased from 0.29 to 0.40 during the period 1990-2000; however, it declined during the period 2000-2005. After that, an increasing trend in the elasticity of ESV was observed during the period 2005 to 2020. The increasing trend of elasticity implies that ESV is becoming more sensitive to LULC changes with time. The highest elasticity of ESV was observed during the period 2015 to 2020. The elasticity of ESV implies the sensitivity of ESV in relation to LULC change. The higher the value of elasticity, the higher the sensitivity of ESV to LULC change. The higher value of elasticity of ESV during 2015-2020 implies that the greatest change in the ESV occurred during that period. The impact and relationship between ESV and LULC can also be understood in Figure 8. Figure 8 shows that waterbodies, forest and vegetation, and agriculture had a positive impact upon ESV. On the other hand built-up areas and bare land negatively affect ESV. As shown in Figure 8, the trend between ESV and waterbodies is more similar compared to others which implies that waterbodies have a higher positive impact on ESV change compared to other LULC types. The higher impact of the water body is also justified due to its higher coefficient value and large area. Built-up areas and bare land have a negative impact upon ESV, causing the value of ESV to decrease. Although both of them have a negative impact, the built-up area is more responsible for reducing total ESV due to its larger area compared to bare land.

Sensitivity Analysis
In order to evaluate the reliability of the results, we performed a sensitivity analysis. We calculated the value of the coefficient of sensitivity (CS) based on Equation (7), considering ±50% adjustment in the original value coefficient proposed by Costanza et al. [3]. For the estimated ESV to be reliable, the value CS must be less than 1. If the value of CS is consistently lower than 1, then it is confirmed that the estimation of ESV was relatively inelastic in relation to the value of the coefficient, and the assignment of ecosystem service value coefficient was perfect. In our study, the value of CS was less than 1 (Table  11), which confirms that estimation of total ESV was inelastic in relation to the value coefficient (Table 3, proposed by Costanza et al. [3]. The CS for waterbodies was higher, about 0.9%, compared to agricultural land and forest and vegetation due to its high value of coefficient and large area. The lower value of CS of agricultural land and forest and vegetation implies that the adjustment in their value coefficient had very little impact on the estimation of ESV. Less than ±0.2% change was observed in total ESV due to ±50% change in the value coefficient for agricultural land, and about ±5% change was observed due to the ±50% change in the value coefficient for forest and vegetation.  On the other hand, the comparatively higher CS value of waterbodies implies that change in the value coefficient has a greater impact upon the estimation of ESV, ±50% change in the value coefficient resulted in about ±44% change in the estimated ESV. However, the sensitivity analysis in our study showed that the estimated ESV using the value coefficient in Table 3 is robust despite the uncertainties in the value coefficient.

Land Use Land Cover Change Dynamics
This present study attempts to quantify the LULC dynamics and their impacts on urban ecosystem services in Dhaka during the past three decades from 1990 to 2020. Our findings regarding LULC change dynamics are in close conformity with the findings of other studies: a) agricultural land, waterbodies, and forest and vegetation show a gradual decrease; and built-up area and bare land show a gradual increase over the study period [75]; b) initially, the proportion of waterbodies/wetland was higher, and it decreases over time [35]; c) built-up area continues to increase over time and dominates all other LULC types [75]; and d) a considerable portion of agricultural land, waterbodies, and forest and vegetation have been replaced by built-up area and bare land during the study period [76,77].
LULC changes are dynamic and nonlinear, which means that transition of one land use from another does not follow similar pattern over the year due to many natural and manmade factors. During the study period, the built-up area increased by 188.35%. This increase was only possible at the expense of valuable agricultural land, waterbodies, and forest and vegetation. During the past three decades, the fringe area of Dhaka was rapidly developed through the conversion of wetland/ and agricultural land into built-up area, including housing, commercial activities, service facilities, industries, and associated urban infrastructures [36]. This conversion of wetland/water body, agricultural land was done through land filling. This conversion was the obvious consequence of urbanization and migration in Dhaka [78].
LULC change dynamics in Dhaka was influenced by many factors. Population migration, urbanization, and industrialization are the main factors contributing to LULC change in Dhaka. Being the capital city of Bangladesh, Dhaka is characterized by many urban facilities compared to other cities in Bangladesh. Most importantly, Dhaka is the main hub of income generation in Bangladesh due to the presence of many industries, business opportunities, commercial activities, and private job opportunities. In addition to income generation facilities, Dhaka has many pull factors, including education and health facilities. Due to Dhaka's pull factors and push factors of other cities and rural areas, a large number of people migrate towards Dhaka every year [79]. The increase in population requires many residential, commercial, industrial activities and many urban facilities, including education, health, civic services, which directly affect the dynamics of the LULC change. For example, a high rate of urbanization led to the conversion of wetland to the built-up area of Dhaka, causing rainfall-induced waterlogging and urban flooding in many parts of the city [80].

Assessment of Urban Ecosystem Service Values
Based on the LULC change, we have quantified the ESV of Dhaka during the period from 1990 to 2020. Our findings identified that LULC change is one of the important drivers of ecosystem services. However, due to the differences in value coefficient and extent of land area, the impacts upon ESV of different LULC types have a different magnitude . Our findings suggest that the ESV trend of Dhaka shows a declining change process due to the transition of LULC during the study period. While the decreases in forest and vegetation, agricultural land, and waterbodies were responsible for reduction of ESV, waterbody reduction was the main contributor among them due to its high-value coefficient and large area compared to the other two. During the period from 1990 to 2020, the total estimated ESV of Dhaka decreased by 59.55%, in which 52.74% ESV decreased due to the reduction of waterbodies. This is in line with the findings of other studies in the wellestablished literature [17,[81][82][83]. For example, the study of Wang et al. [84] shows a declining trend of ESV in Sanjiang Plain, northeast China, during the period from 1980 to 2000. Their study identified that the substantial decline in ESV was essentially due to the 53.4% loss of wetlands. During the period from 1980 to 2000, a considerable amount of wetland was lost through the conversion to agricultural land in Sanjiang Plain, northeast China. Our study also identified that, from 1990 to 2020, the built-up area increased from 17.98% to 51.85% in Dhaka at the expense of agricultural land, waterbodies, and forest and vegetation land. Since the built-up area derives no ecosystem service value, the total ESV declined due to the development of the built-up area during the past three decades. A similar scenario is also evident in some other studies showing the development of the built-up area through the conversion of agricultural land [85], forest, and vegetation [86] and waterbodies [87].
In the case of individual ecosystem services functions, water supply, waste treatment, and disturbance regulation are three major contributors over the study period. Their importance is also supported by other studies. For example, according to the study by Hasan et al. [18], water supply, waste treatment, soil formation and retention, and biodiversity protection were the most valuable ecosystem services, affecting the total ESV Guangdong,  Figure 6). The reason behind the highest contribution of regulating services is the presence of a large amount of waterbodies and their high-value coefficient for individual ecosystem services functions. For example, two individual functions of waterbodies have a higher value coefficient: waste treatment and disturbance regulation with a high-value coefficient (Table 3). However, the value of individual ecosystem functions decreased over the study period. For example, the individual value of water supply, waste treatment, and disturbance regulation functions decreased by 59.31% (19.38 million USD), 59.34% (21.67 million USD), and 59.31% (23.13 million USD) respectively in 2020. Some other studies also identified a similar trend. For example, the study of Hasan et al. [18] identified the highest decline in the water supply value (−22.30 billion CNY, −14.72%) between 1986 and 2017, followed by waste treatment (−20.77 billion CNY, −14.63%) in Guangdong, Hong Kong, and Macao in South China.
Our findings on the elasticity of ESV in response to LULC change show that change of ESV to LULC change is highly elastic. The elasticity result implies that about 1% transition in LULC would result in about 0.33% change in total ESV during the study period. Apart from a fluctuation during 2000-2005, the elasticity seems to be increasing ( Figure  7), which implies that ESV is becoming more dependent upon LULC changes over time. A similar result was also identified in other studies. For example, the elasticity of the ecosystem service value change responding to LULC changes during the periods of 1992-2001, 2001-2009, and 2009-2018 were 0.28, 0.34, and 0.50, respectively, indicating an increasing trend over time [30].
The sensitivity analysis was done to test the reliability of our estimation of ESV, and it confirmed that our estimation of ESV is robust since the value of the coefficient of sensitivity (CS) is lower than 1 in all the cases, which is much similar to other studies by Tolessa et al. [17], Gashaw et al. [74], Hasan et al. [18], and Wang et al. [88]. However, we understand and believe that Costanza et al. [3] developed the value coefficient for different land cover types for global estimation and might underestimate real services of land cover types in the local context. In addition, ESV for some LULC types was estimated using equivalent biomes (Table 3), which, in turn, may not be able to derive the actual ESV from that particular LULC type. In addition, ecosystem service value coefficients were assumed to be spatially homogeneous, but, in reality, value coefficients may be different in different places based on local context [17]. Hence, the development of ecosystem service value coefficients based on local context, physical and socioeconomic characteristics is very important. Owing to this reality, different researchers have tried to modify the ecosystem service value coefficients originally developed by Costanza et al. [3] to fit the local context. For example, Xie et al. [89], on the basis of a questionnaire survey of 700 ecological experts, developed an equivalent weighting factor per hectare ESV to be applicable for Chinese terrestrial ecosystems. However, in Bangladesh, there has been no such study that worked on adjusting the ecosystem service value coefficients for urban areas.

Limitations and Future Scope of the Study
In this study, ESV was calculated based on Landsat satellite images derived from land cover and the benefit transfer method proposed by Costanza et al. [3]. However, there are several uncertainties and limitations with the assessment of ESV using this method. First, we classified different land cover types using time-series Landsat images, which are not high-resolution. Since the estimation of ESV in this study is completely based on different land cover types, the accuracy of the estimated ESV entirely depends upon the accuracy of land cover classification. However, due to the medium spatial resolution of Landsat images, the accuracy of land cover classification is also low compared to some other high-resolution images, such as quick bird and world view images, but most high-resolution images are not freely available. Therefore, Sentinel images can be used to derive LULC because Sentinel images have comparatively higher resolution than Landsat images. However, Sentinel images are not available before 2015. In this case, image fusion technique (for Landsat and Sentinel images) can be tested to derive high-resolution LULC before 2015. Second, benefit transfer methods assume the spatial homogeneity of the ecosystem service value coefficient, but human-environmental systems are much complex and heterogeneous; the same LULC may derive different ecosystem service values in different countries since the social, physical, and economic context is inevitably different. So the generalization of the value coefficient would inevitably generate erroneous ESV [90]. Hence, we recommend developing a value coefficient considering the local context. For Bangladesh, there is no specific value coefficient of different LULC types. So, future study can focus on developing value coefficient considering the local contest of Bangladesh. Third, since LULC categories are not a perfect match to equivalent biomes, several studies calculated the coefficient of sensitivity (CS) to examine the reliability and robustness of the estimation. Most of the studies adjusted value coefficients by ±50% to calculate CS. However, testing the robustness and reliability of the result was criticized by many researchers because of the tendency for CS value to be often less than 1 even when the value coefficient was adjusted by ±25% [91]. This indicates that the CS is robust but erroneous. Fourth, we have assumed that the value coefficient of ESV remains constant during the study period, but, in reality, it is unlikely that the value coefficient of ESV would remain constant over time [92], so the ESV change rate over time would not reflect the actual ESV change rate. Therefore, we recommend adjusting the value coefficient over time in future studies.

Conclusions
This study shows how ESV changes over time due to the change of LULC. The study identified that built-up area is the most dominant LULC type during the whole period except in 1990. The built-up area increased by 188.35% from 1990 to 2020, with an average annual growth rate is about 6.28%. The built-up area increased at the expense of agricultural land, waterbodies, and forest and vegetation land due to the high rate of urbanization, population migration, industrialization, and socioeconomic development in Dhaka. The analysis of ESV shows that it has decreased by 59.55% (85 million USD) from 142.72 million USD in 1990 to 57.72 million USD in 2020 due to the development of the built-up area through conversion of agricultural land, waterbodies, and forest and vegetation land. The elasticity result implies that 1% transition in LULC would result in about 0.33% change in total ESV during the study period. However, the highest decline in total ESV (32.20%, 27.32 million USD) occurred during the period 1990 to 1995. The study identified that waterbodies have a higher positive impact on ESV change compared to other LULC types. It has also been identified that water supply, waste treatment, and disturbance regulation are three major contributors to ESV in our study area.
This study is the first of its kind to understand the impact of LULC change on urban ecosystem services and their functions in Dhaka. The valuation of urban ecosystem services has many significances in sustainable land use planning because it has a great impact on the human well-being of city people. For example, urban vegetation controls the urban heat island effect. Loss of urban vegetation can require additional energy demand for cooling and heating buildings, increase the number of people with respiratory diseases, and many more impacts. All of these are related to an economic cost, so it is important to explore how ecosystem services respond in relation to LULC changes to facilitate decision-making towards urban sustainability. We strongly believe that this study would help policymakers and urban planners to devise appropriate land use decisions to ensure sustainable urban development of Dhaka considering natural resources and the trade-off between urban development and reduction of ecosystem services.