Ecological and Environmental Effects of Estuarine Wetland Loss Using Keyhole and Landsat Data in Liao River Delta, China

: An estuarine wetland is an area of high ecological productivity and biodiversity, and it is also an anthropic activity hotspot area, which is of concern. The wetlands in estuarine areas have suffered declines, which have had remarkable ecological impacts. The land use changes, especially wetland loss, were studied based on Keyhole and Landsat images in the Liao River delta from 1962 to 2016. The dynamics of the ecosystem service values (ESVs), suitable habitat for birds, and soil heavy metal potential ecological risk were chosen to estimate the ecological effects with the beneﬁt transfer method, synthetic overlaying method, and potential ecological risk index (RI) method, respectively. The driving factors of land use change and ecological effects were analyzed with redundancy analysis (RDA). The results showed that the built-up area increased from 95.98 km 2 in 1962 to 591.49 km 2 in 2016, and this large change was followed by changes in paddy ﬁelds (1351.30 to 1522.39 km 2 ) and dry farmland (189.5 to 294.14 km 2 ). The area of wetlands declined from 1823.16 km 2 in 1962 to 1153.52 km 2 in 2016, and this change was followed by a decrease in the water area (546.2 to 428.96 km 2 ). The land use change was characterized by increasing built-up (516.25%), paddy ﬁelds (12.66%) and dry farmland (55.22%) areas and a decline in the wetland (36.73%) and water areas (21.47%) from 1962–2016. Wetlands decreased by 669.64 km 2 . The ESV values declined from 6.24 billion US$ to 4.46 billion US$ from 1962 to 2016, which means the ESVs were reduced by 19.26% due to wetlands being cultivated and the urbanization process. The area of suitable habitat for birds decreased by 1449.49 km 2 , or 61.42% of the total area available in 1962. Cd was the primary soil heavy metal pollutant based on its concentration, accumulation, and potential ecological risk contribution. The RDA showed that the driving factors of comprehensive ecological effects include wetland area, Cd and Cr concentration, river and oil well distributions. This study provides a comprehensive approach for estuarine wetland cultivation and scientiﬁc support for wetland conservation.


Introduction
Wetland loss and degradation has increased as the global demand for land [1]. Anthropic activities have globally modified wetlands [2], and more than 50% of wetlands have been altered since mid-19th century on earth [3]. Wetlands provide hydrological, biogeochemical, and ecological services (climate regulation, biodiversity protection, food production) and economic benefits (tourism, recreation) for human welfare [4][5][6]. Particularly, coastal wetland is one of the most vulnerable and economical ecosystems [7]. About 25-50% of coastal wetlands had been lost during the 20th century alone by converting into other land uses [8,9]. With the implementation of opening-up and economic policies

Data Collection and Process
The land use maps in 1988, 1998, 2008, and 2016 were interpreted based on Landsat TM/ETM images (https://earthexplorer.usgs.gov/). The land use maps in 1962 and 1972 were interpreted with Keyhole (KH5) (3.6 m × 3.6 m) and Keyhole (KH9) (0.3 m × 0.3 m) images (http://www.gscloud.cn/), including paddy field, dry farmland, built-up area, forestland, wetland, and water area, as seen in Table 1. A digital elevation model (DEM) (1:50,000) was collected from Liaoning surveying and Mapping Bureau. The images were pre-processed by registration, splicing, inlaid, and cutting for each year. The manual interpretation was used to get the classification results of land use. The annual normalized difference vegetation index (NDVI) from 1988 to 2016 was calculated using the maximum synthesis method (MVC) based on Google Earth Engine (https://earthengine.google.com/ for revision of the equivalent factors of ESVs). The spatial resolution of the data ranges from 0.3 m to 30 m, and all images were aggregated to 30 m spatial resolution by resampling. Taking the overall accuracy assessment to evaluate the accuracy of interpretations, 500 sample points were randomly selected in study area. The overall accuracy of the land use maps for 1962, 1972, 1988, 1998, 2008, and 2016 were 78.42%, 80.65%, 85.21%, 88.02%, 88.57%, and 89.35%, respectively. Manual interpretation and accuracy evaluation in 1962, 1972, and 1988 were mainly carried out by field surveys and historical records due to the smaller number of available data, which were affected by clouds, time, and other factors. The dataset of roads and oil wells was extracted from the above images. The social economic dataset, including population, GDP, agricultural and fishery output data, were collected based on the Panjin annual statistics book from 1962 to 2016, which were obtain from Liaoning statistical Bureau and field surveys. Table 1. Characteristics of multi-year images and soil sample data used in the study.

Data Date of Data Description Source
A total of 184 soil samples (surface (0-5 cm) were selected in October 2014. The geographic positions of sampling sites were recorded using portable GPS units, which is distributed as shown in Figure 1. Three duplicates were extracted and mixed at each sampling site. Six soil heavy metal pollutants were tested, including cadmium (Cd), copper (Cu), chromium (Cr), nickel (Ni), zinc (Zn), and lead (Pb). The soil standard reference material (GBW07401, GSS-1) provided by the Center of National Standard Reference Material of China was used to control the test quality. Soil samples were treated with a microwave digestion instrument (CEM Inc., Matthews, NC, USA), and then tested for the concentrations of the heavy metals. The heavy metals were measured via inductively coupled plasma mass spectrometry (ICP-MS, PerkinElmer, Waltham, MA, USA) [43]. A total of 184 soil samples (surface (0-5 cm) were selected in October 2014. The geographic positions of sampling sites were recorded using portable GPS units, which is distributed as shown in Figure 1. Three duplicates were extracted and mixed at each sampling site. Six soil heavy metal pollutants were tested, including cadmium (Cd), copper (Cu), chromium (Cr), nickel (Ni), zinc (Zn), and lead (Pb). The soil standard reference material (GBW07401, GSS-1) provided by the Center of National Standard Reference Material of China was used to control the test quality. Soil samples were treated with a microwave digestion instrument (CEM Inc., Matthews, NC, USA), and then tested for the concentrations of the heavy metals. The heavy metals were measured via inductively coupled plasma mass spectrometry (ICP-MS, PerkinElmer, Waltham, MA, USA) [43].

Land Use Dynamic Change Analysis Method
To quantitatively reflect the rate of land use change, the single dynamic index was used [44]. The land use TUPU analysis method was adopted to the land use change. The land use maps were used to build transfer matrix for information visualization of land use conversion. It was conducted to integrate the spatial information of the coded land use change value in ArcGIS. Specifically, the codes of the adjacent two phase maps were selected for the algebraic operation to obtain the value [20]. (1) where K represents the single dynamic index of land use, Ua and Ub represent the area of a certain land type at the beginning and the end of the period, respectively, and T represents the research period.

Comprehensive Evaluation Method
The comprehensive environmental effects were evaluated in this paper from three aspects, the change of ESVs, habitat suitability for water birds, and pollution risk of soil heavy mental. The ESVs and habitat suitability were assessed based on the value equivalent factor method and synthetic overlaying analysis of factors, respectively. The pollution risk was assessed by kriging interpolation and index evaluation. Then, the driving factors of comprehensive environmental effects were studied through redundancy analysis (RDA) from natural, socio-economic, and spatial geographic conditions, as shown in Figure 2.

Land Use Dynamic Change Analysis Method
To quantitatively reflect the rate of land use change, the single dynamic index was used [44]. The land use TUPU analysis method was adopted to the land use change. The land use maps were used to build transfer matrix for information visualization of land use conversion. It was conducted to integrate the spatial information of the coded land use change value in ArcGIS. Specifically, the codes of the adjacent two phase maps were selected for the algebraic operation to obtain the value [20].
where K represents the single dynamic index of land use, U a and U b represent the area of a certain land type at the beginning and the end of the period, respectively, and T represents the research period.

Comprehensive Evaluation Method
The comprehensive environmental effects were evaluated in this paper from three aspects, the change of ESVs, habitat suitability for water birds, and pollution risk of soil heavy mental. The ESVs and habitat suitability were assessed based on the value equivalent factor method and synthetic overlaying analysis of factors, respectively. The pollution risk was assessed by kriging interpolation and index evaluation. Then, the driving factors of comprehensive environmental effects were studied through redundancy analysis (RDA) from natural, socio-economic, and spatial geographic conditions, as shown in Figure 2.

Evaluation of the Ecosystem Service
The assessment method proposed by [29] has been widely calculated ESVs. The modification of the Chinese ecosystem coefficients was based on expert knowledge through over 700 ecologists [30,45], generating the table of the value per unit area of the ecosystem. The annual corrected equivalent factor was estimated at an equivalent ESVs of 410.80 US$/(hm 2 a) (based on the 2016 USD value; $1 could be converted to 6.67 RMB) based on the grain yield per unit area and the average grain price in 2016. The average equivalent ESVs coefficients of the wetlands were revised using the equivalent factors by the regional NDVI values, which may produce accurate results, especially in regions with large areas of reed wetland. Besides, built-up areas were not considered to provide ecosystem services [46]. Additionally, cultivated land was mentioned in this study, which includes paddy fields and dry farmland, for the ESVs calculation as following: where , max, and min indicate the annual average , the maximum, and minimum annual average for the local grid, respectively. VCI indicates the vegetation condition index calculated from time-series NDVI. VCI is the annual average VCI in study area, k is the average in the k county. ESV and corrected indicate the ecosystem service value before and after correction.
ij indicates the ℎ category of the coefficient for the ℎ land use type, i indicates the ℎ land use type area, n is the number of land use types, and m is the number of sub-categories of ESVs.
Sensitivity analysis was calculated the coefficient of sensitivity of the value coefficient, to test whether the results of the ESV coefficient are credible. This method was used to modify the value coefficient of each land use type by (±) 50%, and then the changes of corresponding ESVs were calculated. The reliability of the ESVs results can be verified by sensitivity analysis [24], and the index can be calculated as: where CS represents the coefficient of sensitivity. When CS < 1, it indicates that ESV lacks flexibility for VC, and the resulting reliability is high, When CS > 1, just the opposite.

Evaluation of the Ecosystem Service
The assessment method proposed by [29] has been widely calculated ESVs. The modification of the Chinese ecosystem coefficients was based on expert knowledge through over 700 ecologists [30,45], generating the table of the value per unit area of the ecosystem. The annual corrected equivalent factor was estimated at an equivalent ESVs of 410.80 US$/(hm 2 a) (based on the 2016 USD value; $1 could be converted to 6.67 RMB) based on the grain yield per unit area and the average grain price in 2016. The average equivalent ESVs coefficients of the wetlands were revised using the equivalent factors by the regional NDVI values, which may produce accurate results, especially in regions with large areas of reed wetland. Besides, built-up areas were not considered to provide ecosystem services [46]. Additionally, cultivated land was mentioned in this study, which includes paddy fields and dry farmland, for the ESVs calculation as following: where NDV I, NDV Imax, and NDV Imin indicate the annual average NDV I, the maximum, and minimum annual average NDV I for the local grid, respectively. VCI indicates the vegetation condition index calculated from time-series NDVI. VCI is the annual average VCI in study area, VCIk is the average VCI in the k county. ESV and ESVcorrected indicate the ecosystem service value before and after correction. VCij indicates the jth category of the coefficient for the ith land use type, LU Ai indicates the ith land use type area, n is the number of land use types, and m is the number of sub-categories of ESVs.
Sensitivity analysis was calculated the coefficient of sensitivity of the value coefficient, to test whether the results of the ESV coefficient are credible. This method was used to modify the value coefficient of each land use type by (±) 50%, and then the changes of corresponding ESVs were calculated. The reliability of the ESVs results can be verified by sensitivity analysis [24], and the index can be calculated as: where CS represents the coefficient of sensitivity. When CS < 1, it indicates that ESV lacks flexibility for VC, and the resulting reliability is high, When CS > 1, just the opposite. ESVi and ESV j indicate the ESVs before and after the adjustment, VCi and VCj indicate the ESV coefficients before and after the adjustment, and k represents the land use type.

Evaluation of Suitable Habitat for Birds
Estuarine wetlands in the Liao River delta are crucial bird migration stopover sites. A combination of single factor extraction and the merging method was used to establish the suitable habitat for birds [33]. In order to obtain scientific evaluation results, field surveys combined with consulting experts were used to select and determine six factors that affected the bird habitat, namely disturbance degree, food richness, water condition, shelter condition, roads, and oil wells. The classification of disturbance degree, food richness, and water condition were based on the land use classification, and NDVI analysis for shelter condition was also required. Based on previous research results and the field surveys, the characteristics of suitable habitat for birds was determined [47]. Taking the red-crowned crane, which is highly sensitive to disturbance, the nearest distance to a road needs to be 410 m, and to an oil well is 500 m [48]. These six types of habitat factor maps were carried out by synthetic overlaying analysis, which is used to combine and overlay among factor maps in spatial analysis, and generated suitable habitat maps to reflect the bird habitat changes.

Pollution Risk Assessment
The geo-accumulation index was used to assess the soil heavy metal contamination level [49]. The potential ecological risk index (RI) was evaluated the level of heavy metal pollution in topsoil, indicating the toxicity and environmental response of the soil heavy metals [50]. The kriging interpolation method was adopted to calculate the spatial distribution of the pollution. The equations are as follows: where Igeo is the geo-accumulation index of a sampling point, C i represents heavy metal i concentration in topsoil, B i is the background value of heavy metal i [51], which is obtained from Chinese soil element background value, and 1.5 is the background matrix correction factor. E i r is the single potential ecological risk factor, T i r is the toxic-response factor for a given substance, which accounts for the toxic requirement and the sensitivity requirement, and RI represents the sum of all risk factors, as shown in Table 2. C i j represents the contamination factor, C i 0 and C i n are the concentration and a reference value for soil heavy metal, respectively. In this paper, the soil background values of Liaoning province were considered as references to evaluate the present pollution conditions and the potential for ecological impacts. Table 2. Description of the Igeo, risk factor (E i r ), and risk index (RI) ranks as suggested by [43].

Driving Factors Analysis
The redundancy analysis (RDA) was applied to analyze the relationship between the wetland loss and the environmental effects. The determination of factors was chosen based on related studies [17,52]. A total of 13 factors from 3 aspects were chosen. Four variables were selected to represent socio-economic conditions, including population, GDP, agricultural, and fishery output. The five variables reflected the natural conditions, including the NDVI, soil texture, wetland area, agriculture area, and soil heavy metal contamination. Four variables were selected to represent the spatial geographic conditions, including the distance to a river, distance to a road, distance to a city, and distance to an oil well. The ESVs, suitable habitat area, and RI were selected to represent ecological services, suitable habitat, and soil heavy metal potential risk. A detrended correspondence analysis (DCA) was proceed to choose the most valid ordination method. The RDA method was chosen due to the gradient length was less than 3 (0.3). According to the values of contribution and significance, the major factors that affected the dependent variables were selected by the interactive forward-selection method. The RDA method was performed by the selected factors, which was undertaken based on the software CANOCO 4.5 [53].

Spatial and Temporal Changes in Land Use
The area of paddy fields, dry farmland, wetland, forestland, built-up areas, and water area were 1351. 30 Figure 3) showed that the built-up area had a continuous increase (0.13%, 5.92%, 7.82%, 3.59%, and 1.86%), and the change rate of built-up areas peaked from 1988 to 1998. Conversely, a continuous declining tendency (2.85%, 0.10%, 0.09%, 0.26%, and 1.18%) was observed in wetlands, and the wetland loss was the fastest in the first and the last periods. Forestland rarely appeared since the 1980s. The water area was gradually decreasing before the 1990s, and then, the water area remarkably increased to a proportion of 4.72% during 2008-2016.
based on related studies [17,52]. A total of 13 factors from 3 aspects were chosen. Four variables were selected to represent socio-economic conditions, including population, GDP, agricultural, and fishery output. The five variables reflected the natural conditions, including the NDVI, soil texture, wetland area, agriculture area, and soil heavy metal contamination. Four variables were selected to represent the spatial geographic conditions, including the distance to a river, distance to a road, distance to a city, and distance to an oil well. The ESVs, suitable habitat area, and RI were selected to represent ecological services, suitable habitat, and soil heavy metal potential risk. A detrended correspondence analysis (DCA) was proceed to choose the most valid ordination method. The RDA method was chosen due to the gradient length was less than 3 (0.3). According to the values of contribution and significance, the major factors that affected the dependent variables were selected by the interactive forward-selection method. The RDA method was performed by the selected factors, which was undertaken based on the software CANOCO 4.5 [53].

Spatial and Temporal Changes in Land Use
The area of paddy fields, dry farmland, wetland, forestland, built-up areas, and water area were 1351. 30 Figure 3) showed that the built-up area had a continuous increase (0.13%, 5.92%, 7.82%, 3.59%, and 1.86%), and the change rate of built-up areas peaked from 1988 to 1998. Conversely, a continuous declining tendency (2.85%, 0.10%, 0.09%, 0.26%, and 1.18%) was observed in wetlands, and the wetland loss was the fastest in the first and the last periods. Forestland rarely appeared since the 1980s. The water area was gradually decreasing before the 1990s, and then, the water area remarkably increased to a proportion of 4.72% during 2008-2016.  pansion, and harbor construction. The water area was converted to built-up areas (Figure 4 code 56) of 88.06 km 2 , which were distributed along the Liao River and coastal zone. The paddy fields, dry farmland, and water area were converted to built-up areas with proportions of 9.13%, 2.00%, and 4.73%, respectively. The "land use change" was slight at the core of the delta since 1988, which is mostly distributed in natural reserves, indicating that wetland resources were effectively protected in the Liao River delta through establishment of the reserve.
during 1962-2016 ( Figure 4). The wetlands were converted to paddy fields (25.30%), dry farmland (5.96%), built-up areas (13.15%), and water areas (16.43%). The conversion of "wetland-paddy field" (Figure 4 code 51) had the largest change in area of 471.79 km 2 , which indicated that the wetland was being cultivated. Meanwhile, the conversion of "wetland-water area" (Figure 4 code 56) for 306.28 km 2 and "wetland-built-up area" (Figure 4 code 53) for 245.19 km 2 took the second and third change area order, which reflected that the wetland was mainly being occupied by aquaculture activities, urban expansion, and harbor construction. The water area was converted to built-up areas (Figure 4 code 56) of 88.06 km 2 , which were distributed along the Liao River and coastal zone. The paddy fields, dry farmland, and water area were converted to built-up areas with proportions of 9.13%, 2.00%, and 4.73%, respectively. The "land use change" was slight at the core of the delta since 1988, which is mostly distributed in natural reserves, indicating that wetland resources were effectively protected in the Liao River delta through establishment of the reserve.

Ecosystem Services Changes
The ESVs results were 6.24, 5.01, 4.87, 4.67, 4.60, and 4.46 billion US$ in 1962,1972,1988,1998,2008, and 2016, respectively. The ESVs had lost about 1.78 billion US$, accounting for 28.58% since 1962 ( Table 3). The values of waste regulation and hydrological regulation were the highest ones, followed by climate regulation, biodiversity protection, soil formation and retention, recreation and culture, gas regulation, and raw material production ( Figure 5A). The ESVs values of ecological service functions showed a declining change, however the ESVs of raw material had an increasing trend (0.33%) from 1962 to 2016. The declining ESVs were mainly caused by the loss of wetland, which has the highest ecological services. Additionally, forestland and water areas were converted to cultivated land, which led to an increase in raw materials. The values of regulating services decreased (22.50%, 33.01%, 29.43%, and 31.58%), and the values of supporting services decreased (20.47% and 26.00%) 1962-2016. Regulating services and supporting services were contributed to the ESVs with the large proportion that were observed. The ESVs of wetland and water area decreased significantly, with the proportions of 36.73% and 21.47%, and the ESVs of cultivated land increased 17.89% (Table 4). The wetland contributes the most to the ESVs, followed by cultivated land and water area. Compared with other land use types, the ESVs of wetland always showed negative changes during 1962-

Ecosystem Services Changes
The ESVs results were 6.24, 5.01, 4.87, 4.67, 4.60, and 4.46 billion US$ in 1962,1972,1988,1998,2008, and 2016, respectively. The ESVs had lost about 1.78 billion US$, accounting for 28.58% since 1962 ( Table 3). The values of waste regulation and hydrological regulation were the highest ones, followed by climate regulation, biodiversity protection, soil formation and retention, recreation and culture, gas regulation, and raw material production ( Figure 5A). The ESVs values of ecological service functions showed a declining change, however the ESVs of raw material had an increasing trend (0.33%) from 1962 to 2016. The declining ESVs were mainly caused by the loss of wetland, which has the highest ecological services. Additionally, forestland and water areas were converted to cultivated land, which led to an increase in raw materials. The values of regulating services decreased (22.50%, 33.01%, 29.43%, and 31.58%), and the values of supporting services decreased (20.47% and 26.00%) 1962-2016. Regulating services and supporting services were contributed to the ESVs with the large proportion that were observed. The ESVs of wetland and water area decreased significantly, with the proportions of 36.73% and 21.47%, and the ESVs of cultivated land increased 17.89% (Table 4). The wetland contributes the most to the ESVs, followed by cultivated land and water area. Compared with other land use types, the ESVs of wetland always showed negative changes during 1962-2016. The area with high ESVs has gradually decreased, and the area with low ESVs has increased 1962-2016.
The results of ESVs were not flexible relative to the value coefficient due to the CS value was less than 1. The CS values in descending order were wetland > cultivated land > water area > forestland ( Figure 5B). The CS value of wetland in 1962 was the highest, up to 0.8, which shows a higher decreasing trend than the other land use types. Moreover, the CS values of cultivated land and water area had relatively large effects on the total ESVs. The results of ESVs were not flexible to the VC value, which was observed by sensitivity analysis, and even though the VC value was uncertain, the results were still stable.

Birds Suitable Habitat Changes
The suitable habitat for birds was divided into the suitable and unsuitable habitat. The results showed that the suitable habitat was constantly decreasing with a degradation trend ( Figure 6). The area of bird suitable habitat declined from 2359.80 km 2 in 1962 to 910.31 km 2 in 2016 due to the wetlands being cultivated and the urbanization process, which experienced a continuous decrease of 1449.49 km 2 , with a proportion of 61.42%. The wetland was dominated by natural wetland, such as reed swamp with rare built-up areas, obviously without oil wells. The bird habitat area accounted for 59% of the study area in 1962. Suitable habitat for birds was mainly distributed in the estuarine wetlands in 2016, which reflects that the wetland and tidal-flats in the Liao River delta are the main habitats of birds. Contrarily, the unsuitable habitat area was gradually increasing, mainly distributed around built-up areas and oil wells.

Birds Suitable Habitat Changes
The suitable habitat for birds was divided into the suitable and unsuitable habitat. The results showed that the suitable habitat was constantly decreasing with a degradation trend ( Figure 6). The area of bird suitable habitat declined from 2359.80 km 2 in 1962 to 910.31 km 2 in 2016 due to the wetlands being cultivated and the urbanization process, which experienced a continuous decrease of 1449.49 km 2 , with a proportion of 61.42%. The wetland was dominated by natural wetland, such as reed swamp with rare built-up areas, obviously without oil wells. The bird habitat area accounted for 59% of the study area in 1962. Suitable habitat for birds was mainly distributed in the estuarine wetlands in 2016, which reflects that the wetland and tidal-flats in the Liao River delta are the main habitats of birds. Contrarily, the unsuitable habitat area was gradually increasing, mainly distributed around built-up areas and oil wells.
The suitable habitat for birds was divided into the suitable and unsuitable habitat. The results showed that the suitable habitat was constantly decreasing with a degradation trend ( Figure 6). The area of bird suitable habitat declined from 2359.80 km 2 in 1962 to 910.31 km 2 in 2016 due to the wetlands being cultivated and the urbanization process, which experienced a continuous decrease of 1449.49 km 2 , with a proportion of 61.42%. The wetland was dominated by natural wetland, such as reed swamp with rare built-up areas, obviously without oil wells. The bird habitat area accounted for 59% of the study area in 1962. Suitable habitat for birds was mainly distributed in the estuarine wetlands in 2016, which reflects that the wetland and tidal-flats in the Liao River delta are the main habitats of birds. Contrarily, the unsuitable habitat area was gradually increasing, mainly distributed around built-up areas and oil wells.

Potential Soil Pollution Ecological Risk
The statistical results indicated that the mean concentrations of soil heavy metals were measured, Cr (154.95 ± 52.16), Zn (102.21 ± 19.44), Ni (40.71 ± 9.13), Cu (17.23 ± 6.05), Pb (11.97 ± 3.66), and Cd (1.41 ± 0.82) ( Table 5). The coefficient of variation of heavy metals was within the range of 12.65-58.23%. Particularly, the concentrations of Cd, Cr, Cu, and Pb had a high coefficient of variation (>30%), which indicated that soil heavy metals have similar sources. The Ni and Zn had a low coefficient of variation (<30%), which indicated

Potential Soil Pollution Ecological Risk
The statistical results indicated that the mean concentrations of soil heavy metals were measured, Cr (154.95 ± 52.16), Zn (102.21 ± 19.44), Ni (40.71 ± 9.13), Cu (17.23 ± 6.05), Pb (11.97 ± 3.66), and Cd (1.41 ± 0.82) ( Table 5). The coefficient of variation of heavy metals was within the range of 12.65-58.23%. Particularly, the concentrations of Cd, Cr, Cu, and Pb had a high coefficient of variation (>30%), which indicated that soil heavy metals have similar sources. The Ni and Zn had a low coefficient of variation (<30%), which indicated they were more uniformly distributed in the study area. Among them, the mean concentrations of Cd, Cr, Ni, and Zn were all above the local soil background values. The I geo values and ranks of pollutants at each sampling point are presented in Figure 7. The I geo values of Cd ranged from 0.70-3.26, with the highest value in paddy fields and the lowest value in built-up areas. Cd reached a high level. The I geo values of Cr reached light to moderate levels. Cr with highest I geo value was around the core of the Liao River estuary. The pollution hotspots of Zn were closer to the paddy fields, wetlands, and built-up areas, which reached a moderate level. Some Ni pollution hot spots were found near oil wells in the wetlands. The I geo values of Cu and Pb were almost 0, revealing that their level of pollution is relatively low. The I geo values were ranked as Cd > Cr > Ni > Cu > Zn > Pb. fields and the lowest value in built-up areas. Cd reached a high level. The I value Cr reached light to moderate levels. Cr with highest I value was around the cor the Liao River estuary. The pollution hotspots of Zn were closer to the paddy fields, w lands, and built-up areas, which reached a moderate level. Some Ni pollution hot s were found near oil wells in the wetlands. The I values of Cu and Pb were almo revealing that their level of pollution is relatively low. The I values were ranked as > Cr > Ni > Cu > Zn > Pb. The Cd posed a considerable potential ecological risk at 184 sampling points. high ecological risks of Cd are the consequences of fertilizers and pesticides used in a culture activities. Ni posed the second largest ecological risk. For Cu, Cr, Zn, and Pb risk indexes were low levels.
indicated that the degree of pollution of the soil he metals decreased as follows: Cd > Ni > Cu > Cr > Pb > Zn. The mean value of RI was 343 which indicated that all sample sites arrived at a moderate potential ecological risk l in Figure 8A. A total of 81 of 184 sites exhibited considerable or high ecological risk. The Cd posed a considerable potential ecological risk at 184 sampling points. The high ecological risks of Cd are the consequences of fertilizers and pesticides used in agriculture activities. Ni posed the second largest ecological risk. For Cu, Cr, Zn, and Pb, the risk indexes were low levels. E i r indicated that the degree of pollution of the soil heavy metals decreased as follows: Cd > Ni > Cu > Cr > Pb > Zn. The mean value of RI was 343.98, which indicated that all sample sites arrived at a moderate potential ecological risk level in Figure 8A. A total of 81 of 184 sites exhibited considerable or high ecological risk. The RI values were significantly correlated with the presence of oil wells. The RI of estuarine wetlands in the reserve were at low and moderate levels as shown in Figure 8B. The distributions of Cd and the RI were similar, indicating that the RI was impacted by Cd. RI values were significantly correlated with the presence of oil wells. The RI of estuarine wetlands in the reserve were at low and moderate levels as shown in Figure 8B. The distributions of Cd and the RI were similar, indicating that the RI was impacted by Cd.

Relationships among the Changes of Ecological Effects and Factors
The results showed that the wetland area, Cd concentrations, Cr concentrations, distance to oil wells, distance to a river, and agricultural and fishery output have power in explaining the ESV, RI, and suitable habitat (p < 0.05). Approximately 74.70% of the variation were explained by the selected variables. The first two ordination axes explained over 78.0% of the correlation with both the driving factors of wetland loss and ecological effects (p = 0.01). The results indicated that the percentage of wetland area had a greater impact on environmental effects, with contributions of 44.5%, and the Cd concentrations also had a greater impact, with contributions of 37.3% in Table 6. As shown in Figure 9,

Relationships among the Changes of Ecological Effects and Factors
The results showed that the wetland area, Cd concentrations, Cr concentrations, distance to oil wells, distance to a river, and agricultural and fishery output have power in explaining the ESV, RI, and suitable habitat (p < 0.05). Approximately 74.70% of the variation were explained by the selected variables. The first two ordination axes explained over 78.0% of the correlation with both the driving factors of wetland loss and ecological effects (p = 0.01). The results indicated that the percentage of wetland area had a greater impact on environmental effects, with contributions of 44.5%, and the Cd concentrations also had a greater impact, with contributions of 37.3% in Table 6. As shown in Figure 9, RI was positively related to Cd, Cr, agriculture and fishery output, and the distance to oil wells and negatively related to the wetland area. Wetland area has a negative relationship with the potential ecological risk. A positive correlation between wetland area and ESVs was found, which demonstrated that the change of wetland is likely to lead to ESV loss. Suitable habitat area was positively related to the distance to rivers, and negatively related to the distance to an oil well, which indicated that the water condition is an important factor for bird habitats and oil exploration is likely to lead to suitable habitat loss. The wetland loss and driving factors, such as agriculture, oil wells, and industrial exploitation have resulted in cumulative environmental effects.

Influence Process of Wetland Loss on Ecological Effects
The wetland was the primary land use type in the Liao River delta during 1962-2016. However, the wetland area has decreased dramatically, mainly being converted into cultivated land and built-up areas. The wetland has been developed as farmland, built-up areas, and aquaculture in the past few decades [54]. Some surveys were carried out on the changes of wetlands, such as the Yellow River delta, Yangtze River delta, Sanjiang plain, and Jiangsu coastal zone. Comparing the wetlands loss in the Liao River delta to those four areas showed that wetland loss was common. For example, the coastal wetland in Jiangsu has lost about 323.46 km 2 over the past 10 years [55], and wetland loss in Sanjiang plain has reached 73.3% (about 27,700 km 2 ) since 1954 [56]. The ESVs lost in the Yellow River delta has decreased by about 10.89% in 2009-2015 [57], and the ESVs of Yangtze

Influence Process of Wetland Loss on Ecological Effects
The wetland was the primary land use type in the Liao River delta during 1962-2016. However, the wetland area has decreased dramatically, mainly being converted into cultivated land and built-up areas. The wetland has been developed as farmland, built-up areas, and aquaculture in the past few decades [54]. Some surveys were carried out on the changes of wetlands, such as the Yellow River delta, Yangtze River delta, Sanjiang plain, and Jiangsu coastal zone. Comparing the wetlands loss in the Liao River delta to those four areas showed that wetland loss was common. For example, the coastal wetland in Jiangsu has lost about 323.46 km 2 over the past 10 years [55], and wetland loss in Sanjiang plain has reached 73.3% (about 27,700 km 2 ) since 1954 [56]. The ESVs lost in the Yellow River delta has decreased by about 10.89% in 2009-2015 [57], and the ESVs of Yangtze River delta has lost about 2.12 billion US$ 2009-2016 [58]. A continuous declining tendency was observed in ESVs of wetlands due to wetland being converted to cultivated land, since the ESVs of cultivated land increase significantly [56]. Conversion of wetland to agricultural land increases the material production values, however the conversion led to the ESVs loss of wetlands for pollution purification, biodiversity conservation, leisure and entertainment. In particular, this is critical to migratory birds since in the Liao River delta, suitable habitat for birds has lost from 135.77 to 74.24 km 2 along the Yangtze River 1986-2015 [59]. The suitable habitat for birds is constantly decreasing with a degradation trend in the study area, due to dam construction raised the terrain and silted up the sediment, resulting in the lack of water. Aquaculture has led to the degradation or even disappearance of wetland plant communities, which means the birds have no shelter. Potential ecological risk of heavy metal in the Yellow River delta and the Yangtze River delta are at a high level for Cd [60,61]. The comparison indicated that Cd was the primary pollutant at a high pollution level in the coastal and estuarine areas. Scholars have revealed that Pb pollution is largely related to the exploitation of mines [43]. Some Pb pollution hotspots are more evenly distributed around roads in the study area. The atmospheric Pb caused by the burning of gasoline may be the reason for the distribution of its pollution hot spots. Analysis of all these factors involved in eco-environmental changes have rarely been integrated in the past, but they are comprehensively analyzed and studied in this paper.

Policy Implications
Ecological effects always affect the human well-being and sustainable developments. The conservation of ecosystem and economic development are emphasized in national decision-making reports. Historically, the delta was a swamp area with a small human population. Agriculture and aquaculture developed rapidly after reclamation departments were established in 1960s. After oil and gas were discovered, oil production activities began in 1970. The region has experienced rapid economic growth based on the overexploitation of resources and agriculture activities, which has caused a series of environmental problems. The establishment of a national reserve and the implementation of wetland restoration projects [62], such as "returning cultivated land to wetland", has restored tidal flats and marshes, and the facilities that threaten ecosystems have been reduced. In order to achieve the coordinated development of ecological, economic, and social benefits, the "three lines" (three lines-Permanent basic farmland red line, Urban development boundary, and Ecological red line) should be promoted in China [20,63]. The "Ecological Red Line" policy formed the boundary lines of the key ecological protection area and has important strategic value for ecological security and sustainable development in the Liao River delta. Wetland protection is the prerequisite for development and utilization, so where is the red line of development? The ecological red line of wetland ecosystem is supported in the implementation of controls, which guides the economical and intensive use of wetland resources for sustainable development. The combination of the ecological red line and land use transition to build a regional ecological security pattern can better strengthen regional ecological protection and curb the degradation of the ecosystem. The results of this paper can be regarded as a reference in policy decision-making.

Validity and Limitations
The wetland changes and its ecological effects were studied in this study. The wetland ecosystem is an essential part of the Liao River delta. Which factors have the greatest impact on wetland changes? The dominant factor is still controversial due to many factors driving wetland changes [13]. Many studies have indicated that natural and socio-economic factors mainly lead to the wetland changes [7,52]. A total of 13 representative factors were chosen for analysis, and we also considered spatial geographic conditions in this study. Besides, the ESVs were calculated with modified coefficients. The price of production was decided by the mean values in Panjin city in 2016 in this study. The grain yield per unit area was selected as the revised index by related literatures [30], and the coverage coefficient was selected according to the correlation between the NDVI and the ESVs. Due to the obvious seasonal differences, NDVI data should be chosen to match the time of the study. Moreover, it is revised from the scale of nation to region, in order to ensure the comparability of the equivalents under different scales and to show the heterogeneity of the ESVs in spatial distribution. The NDVI data of 1962 and 1972 were not obtained, and the grain yield per unit area has not been revised by the coefficient. A sensitivity index was used to verify the results of the ESVs, and proved to be reliable. Applying comprehensive and explanatory evaluations such as those employed here will be key to future conservation of wetland resources in estuarine areas.

Conclusions
This paper analyzed the main land use type changes in the Liao River delta during 1962-2016, and estimated the changes of ESVs, the suitable habitat, and risk assessment of soil heavy metals. We have reached the following conclusion. The area of wetland has decreased 669.64 km 2 , with a proportion of 36.73%. The "wetland-paddy field" is the major change. The ESVs has lost about 1.78 billion US$ in 1962-2016. Regulating services and supporting services mainly contributed to the total ESVs. The suitable habitat area was decreased by 1449.49 km 2 , with a proportion of 61.42%. Cd pollution is the highest soil pollution ecological risk in the Liao River delta, followed by Ni, Cu, Cr, Pb, and Zn. Wetland area, Cd, Cr, distance to a river, distance to an oil well, and agricultural and fishery output also had impacts on ecological environmental effects. This paper put forward a set of methods to evaluate the environmental changes of estuarine wetland ecosystem. The implications of this paper provide some references for wetland conservation and the sustainable development of the region, as well as other estuaries area.