The Integrated Valuation of Ecosystem Services and Trade (InVEST) model was jointly designed and developed by Stanford University, the World Wildlife Fund, and the Nature Conservancy to quantitatively assess ecosystem services. The model adopts the production function and empirical modeling methods to quantify the impact of land use and water body changes on the output value of ecosystem services. The evaluation results are visualized in the form of spatial distribution maps, which helps researchers quickly and accurately understand the changing characteristics of ecosystem services. Based on meteorological data, land use data, and other data, this study selected the indicators—soil conservation, water production services, carbon storage, and habitat quality—from the InVEST model to evaluate the ecosystem service value of four wetland nature reserves in the Tianjin Ancient Coast area, Beidagang, Tuanbo, and Dahuangbao.

(1) Soil Conservation Service Module

The soil conservation function of ecosystems in the InVEST model is represented by soil conservation, which is defined in two parts: erosion reduction and sediment retention. Erosion reduction refers to the potential erosion reduction of different land use methods and is expressed as the difference between the potential erosion and actual erosion. Sediment interception capacity is the retention of sediment by a certain land use method and is expressed as the product of sediment volume and sediment retention efficiency. Based on the Soil Loss Equation (USLE), this module considers the ability of the plot to intercept upstream sediments, which leads to more scientific calculation results (Gong et al., 2019). The formulas are as follows (Equations 1–4): S E D R E T x = R K L S x − U S L E x + S E D R x (1) R K L S x = R x × K x × L S x (2) U S L E x = R x × K x × L S x × C x × P x (3) S E D R x = S E x ∑ y = 1 x − 1 U S L E y ∏ z = y + 1 x − 1 1 − S E x (4) where SEDRETx is the soil conservation amount of the x grid; RKLSx is the potential soil erosion amount of the x grid (t); USLEx and USLEy represent the actual erosion amount (t) after considering management and engineering measures for the x grid and its uphill y grid, respectively; SEDRx is the sediment retention capacity (t) after considering management and engineering measures for the x-grid; Rx is the precipitation erosion factor of the x grid; Kx is the soil erodibility factor of the x grid; LSx is the slope length factor of the x grid; Cx is the vegetation cover management factor for the x grid; Px is the soil and water conservation engineering measurement factor for the x grid; and SEx is the sediment retention efficiency of the x-grid.

The precipitation erosion factor (R) is the main external force of hydraulic erosion and reflects the potential ability of rainfall to cause soil erosion. This can be obtained using various calculation methods based on precipitation data. This article selects the empirical formula for precipitation erosivity proposed by Wischmeier (Wischmeier and Smith, 1965), and the specific calculation is as follows (Equation 5): R = ∑ i = 1 12 1.735 × 10 1.5 g p i 2 p − 0.8188 (5) where P is the annual average precipitation (mm), and Pi is the average monthly precipitation (mm). The unit of R is 100 ft t·in/(ac·h·a), which needs to be multiplied by a coefficient of 17.02 to convert the unit to the international unit MJ·mm/(ha·h·hr).

The soil erodibility factor (K) characterizes the difficulty of hydraulic separation and transport of soil particles and can reflect the sensitivity of the soil itself to erosion. Commonly used calculation methods include the Nomo corrosion resistance equation proposed by Wischmeier and the EPIC model equation established by Williams. This study used Williams' model formula (Singh and Jain, 2025), and the specific formula is as follows (Equations 6, 7): K = 0.2 + 0.3 ⁡ exp − 0.0256 S A N 1 − S I L 100 × S I L S I L + C L A 0.3 × 1 − 0.25 C C + exp 3.72 − 2.95 c × 1 − 0.7 S N 1 S N 1 + exp 22.9 S N 1 − 5.51 × 0.1317 (6) S N = 1 − S A N 100 (7) where SAN is the soil sand content (%); SIL is the soil silt content (%); CLA is the soil clay content (%); and C is the soil organic carbon content (%). The unit K reflects t·ha·h/(ha·MJ·mm) from the American system, and it needs to be divided by a coefficient of 7.59 to convert it to the international unit t·ha·hr/(MJ·ha·mm).

(2) Water Production Service Module

The water production calculation module in the InVEST model can not only study the impact of land use/vegetation cover changes on surface water production but also evaluate the relative contribution of various watersheds in the study area to water production. It primarily utilizes the Budyko water heat coupling equilibrium assumption (1974) and annual precipitation data. Based on the principle of water balance, the formula for calculating water production ((Y _xj )) within each grid unit in the study area is as follows (Equation 8): Y x j = 1 − A E T x j P x × P x (8) where AET _xj is the actual evapotranspiration of land use type j in grid unit x (mm), and P _x is the precipitation amount (mm) of grid unit x. A E T x P x = 1 + P E T x P i − 1 + P E T i P i ω 1 ω (9) ω x = Z A W C x P x + 1.25 (10) A W C x = min m a x l a y e r _ d e p t h x , R o o t _ d e p t h x × P A W C (11) The variables involved in (Equations 11–13) are defined as follows: PET _x is the potential evapotranspiration of grid unit x (mm); ω is a non-physical parameter of natural climate soil characteristics, dimensionless; Z is the seasonal precipitation distribution and precipitation depth parameter, which is an empirical constant that characterizes precipitation characteristics. In regions dominated by precipitation in winter, Z approaches 10, whereas in regions with a uniform precipitation distribution or dominated by precipitation in summer, Z approaches 0. In addition, AWC _x is the effective soil moisture content of the land use type on grid unit x, which is the amount of water provided by the soil for plant growth (mm); maxlayer_depth is the maximum burial depth of soil roots (mm); Root_depth is the depth of the plant root system (mm); PAWC is the plant available water content. P E T x = K c l x × E T 0 x (12) In Equation 12, K _c is the plant evapotranspiration coefficient of land use type l _x on grid unit x, and ET _0x is the reference crop evapotranspiration (mm) on grid cell x. P A W C = 54.509 − 0.132 × S s a n − 0.003 × S s a n 2 − 0.055 × S s i l − 0.006 × S s i l 2 − 0.738 × S c l a + 0.007 × S c l a 2 − 2.688 × S o r g + 0.501 × S o r g 2 (13) In Equation 13, S _san represents soil sand content (%); S _sil is soil silt content (%); S _cla is soil clay content (%); and S _org is the soil organic carbon content (%).

(3) Carbon Storage Module

Carbon storage refers to the total amount of carbon stored by terrestrial ecosystems, as well as the natural processes by which ecosystems (including soil and vegetation) fix carbon. The InVEST carbon storage module divides carbon storage into four parts: aboveground biochar (carbon in living plants on the ground), underground biochar (carbon in the roots of living plants), soil carbon (organic carbon in organic and mineral soils), and dead organic carbon (carbon in the litter, standing, or inverted wilted plants). Based on the land use data, the area of each land use type was multiplied by the average carbon density and summed to obtain the total carbon storage capacity of the study area. The calculation formulas are as follows (Equations 14, 15): C = C a b o v e + C b e l o w + C s o i l + C d e a d (14) C t o t a l = ∑ k = 1 n A k × C k k = 1 , 2 , . . . , n (15) where C _above represents the carbon density of above ground organisms, C _below is the underground biological carbon density, C _soil is the soil carbon density, C _dead is the carbon density of dead organic matter, C _total is the total carbon storage, A _k is the area of land use type k, and n is the total number of land use types.

(4) Habitat Quality Module

Habitat quality is an important indicator for evaluating whether an ecosystem can provide a suitable environment and sustainable development for biological reproduction and population survival; it can reflect the biodiversity status of the region to some extent. This study used the habitat quality module in the InVEST model to calculate and evaluate the habitat quality of the study area and implement a visual analysis. This module includes four functions: the relative impact of each threat factor, the relative sensitivity of each type of land use to the threat factor, the distance between the threat source and land use, and the degree to which land use is legally protected. The module operates based on raster data and first calculates the degree of habitat degradation. The formulas are as follows (Equations 16–18): i r x y = 1 − d x y d r ⁡ max (16) i r x y = exp − 2.99 d r ⁡ max d x y linear regression (17) D x j = ∑ r = 1 R ∑ y = 1 Y r ω r ∑ r = 1 R ω r r y i r x y β x S j r index recession (18) where i _rxy represents the impact of the threat source r’s habitat in grid cell x on grid y; d _xy is the linear distance between grid cells x and y; d _{r max} is the maximum impact distance of threat r; D _xj is the habitat degradation index; R is the number of threat factors; ωr is the weight of the threat factor r; Y _r is the grid number of threat factors; r _y is the value of the threat factor on the grid; and β _x is the degree of legal protection (0 for areas protected by law and one for other areas). The larger the calculation result of the habitat degradation degree, the higher the degree of threat posed by threat factors to the habitat and the more susceptible the habitat is to the influence of threat factors causing degradation.

Based on the degree of habitat degradation D _xj , the formula for evaluating habitat quality is as follows (Equation 19): Q x j = H j 1 − D x j z D x j z + K z (19) where Q _xj is the habitat quality index of grid unit x in land use type j, H _j is the habitat suitability of land use type j, K is a semi-saturation constant (usually 1/2 of the maximum value of habitat degradation), and Z is the scaling factor (usually set to 2.5) (Feng et al., 2018). Q _xj ranges from 0 to 1, where one indicates the highest habitat quality, which is beneficial for maintaining biodiversity, and 0 indicates the lowest habitat quality.

The ecosystem services trade-off degree (ESTD) method is based on linear data fitting that reflects the interrelationships between ecosystem services. The formulas are as follows (Equations 20, 21): E S C I i = E S i a − E S i b E S i b (20) E S T D i j = E S C I i / E S C I j + E S C I j / E S C I i 2 (21) where ES _ia and ES _ib are the values of the ith ecosystem service at times a and b, respectively; ESCI _i is the ith ecosystem services change index (ESCD); ESCI _j is the change index of the jth ecosystem service; and ESTD _ij represents the degree of trade-off and synergy between the ith and jth ecosystem services. A negative ESTD indicates a tradeoff relationship between the ith and jth ecosystem services, whereas a positive ESTD indicates a synergistic relationship between the two. The absolute ESTD value reflects the level of trade-off/synergy.

Soil conservation services assess the ability of ecosystems to control soil erosion by calculating and quantifying the amount of soil retained due to natural vegetation, land management, and land use types. Figure 2 shows the spatial distribution of soil retention in the four reserves in 2000, 2010, and 2020. I represents Tianjin Dahuangbao Wetland Nature Reserve, II represents Tianjin Ancient Coast and Wetland National Nature Reserve, III represents Tianjin Tuanbo Bird Nature Reserve, and IV represents Tianjin Beidagang Wetland Nature Reserve.

As shown in Figure 2, soil retention in the four nature reserves changed to varying degrees from 2000 to 2020. The highest soil retention exhibited significant interannual variation, showing an initial increase followed by a decline. Specifically, the highest soil retention increased from 43.4 t/hm² in 2000 to 64.80 t/hm² in 2010, representing a 49.31% increase. However, it then decreased to 56.91 t/hm² in 2020, reflecting a 12.17% decline compared to 2010, but still 31.08% higher than the level in 2000. Adding these percentage changes provides a clearer understanding of the magnitude of variation over time, making the temporal trends more intuitive for readers. From the perspective of spatial distribution, the soil retention of the four nature reserves was generally dominated by low values and low average values, and the high-value areas were unevenly distributed in the form of dense patches. The retention capacities of the waterside areas of the Beidagang Wetland Reserve and Tuanbo Bird Nature Reserve were the most significant, and the retention capacities of the land areas were similar. Among them, soil retention in the Dahuangbao Wetland Nature Reserve and Tuanbo Bird Nature Reserve was relatively high. From the perspective of the different functional zones (Table 2), the average soil retention in the core area was the highest, and it increased steadily from 2000 to 2015 and decreased slightly in 2020. Soil retention in the buffer and experimental zones was lower than that in the core zone, which was related to the land use type and control method of the region.

The water production service calculates the amount of available water generated in each area, thus reflecting the amount of water resources produced by each grid unit. The spatial distribution is shown in Figure 3. In general, the water production distribution of each reserve showed obvious spatial heterogeneity, and the differences in water production gradually expanded over time. From the perspective of spatial distribution, the spatial distribution of water production in the Tianjin Beidagang Wetland Nature Reserve and Tianjin Tuanbo Bird Nature Reserve was relatively stable, with a small change range, and the overall water production characteristics were relatively uniform. In contrast, the Dahuangbao Wetland and Ancient Coast and Wetland National Nature Reserve showed large spatial changes, especially in 2010, when the high water production areas of these two reserves increased significantly and the spatial heterogeneity was enhanced. This shows that the water production capacity of different wetlands was affected by the geographical location, wetland characteristics, and climate change. From the experimental results, the average water production of the four reserves reached the highest value in 2010, which was 123.99 m3. From the perspective of functional zoning (Table 3), unlike soil conservation services, water production in the buffer and experimental zones was generally greater than that in the core zone and fluctuated and increased between 2000 and 2020. This is because the water area accounts for a large proportion of the land use structure in the core zone, and the model does not calculate water production in the water area.

Carbon storage services calculate the amount of carbon stored in ecosystems, especially carbon captured and stored by natural processes, such as plants, soils, and water bodies. As shown in Figure 4, the spatial distribution patterns of carbon storage in the four protected areas in Tianjin showed significant spatial heterogeneity. The highest carbon storage density in the four protected areas was 20 t/hm2, and the lowest was 0.75 t/hm2. Specifically, the carbon storage density of the Tianjin Dahuangbao Wetland Nature Reserve was concentrated in the lower range in 2000 and 2010. By 2020, the carbon storage density of the reserve had increased significantly, particularly in the core and buffer zones (Table 4). These changes were particularly evident. Carbon storage in the Tianjin Tuanbo Bird Nature Reserve and Tianjin Beidagang Wetland Nature Reserve showed a downward trend, and low-value areas expanded. The Tianjin Tuanbo Bird Nature Reserve was completely covered by low-value areas in 2020 with an average of 2.1 t/hm2. From the perspective of different functional zones, the carbon storage in the core area is generally higher than that in the buffer zone of the nuclear experimental area. However, mean carbon storage in the core area exhibited a fluctuating downward trend and was relatively unstable. From 2015 to 2020, carbon storage in the core and buffer zones increased, whereas that in the experimental area continued to decline, indicating that human control had a positive effect on the carbon storage function of the ecosystem.

Habitat quality services assess the quality and suitability of ecosystems and habitats, with a focus on whether a habitat can support the survival, reproduction, and migration of species. According to the data distribution characteristics of the habitat quality index in the study area, the results were divided into four levels using the equal interval method: low-quality area (0, 0.2), medium-quality area (0.4, 0.6), relatively high-quality area (0.6, 0.8), and high-quality area (0.8, 1). The habitat quality index in this study area does not involve the (0.2, 0.4) interval. This classification is illustrated in Figure 5.

According to the spatial distribution characteristics, the habitat quality of each nature reserve shows significant spatial differences and agglomeration. In general, the high-quality and relatively high-quality areas increased significantly from 2000 to 2020, whereas the low-quality areas gradually decreased, indicating that the quality of wetland habitats has shown a trend of gradual improvement over the past 20 years. Among them, the habitat quality of the Tianjin Dahuangbao Wetland Nature Reserve showed an upward trend by 2010 but declined from 2010 to 2020. The highest growth in high-quality areas occurred in the Tianjin Tuanbo Bird Nature Reserve and the Tianjin Beidagang Wetland Nature Reserve. In particular, the high-quality area of the Tianjin Tuanbo Bird Nature Reserve covered the entire reserve area, which is related to the different main protection objectives of the different reserves. From the perspective of functional zoning (Table 5), the areas of medium quality and above in the buffer and core zones accounted for the largest proportion, showing a fluctuating upward trend, among which the increase from 2010 to 2020 was relatively large. As of 2020, the area of medium quality and above in the core zone was 35.99%.

Based on the evaluation results of the above four key ecosystem services (i.e., soil conservation, water production, carbon storage, and habitat quality), maximum and minimum normalization was performed. The same weights were used to construct a comprehensive ecosystem service index to quantitatively evaluate the ecosystem service capacity of different regions under a unified framework. The calculation results are shown in Figure 6.

According to the spatial distribution, the four reserves showed significant spatial differentiation, and the positions of the regional high- and low-value areas were relatively stable. Among them, the Tianjin Dahuangbao Wetland Nature Reserve and Tianjin Ancient Coast National Wetland Nature Reserve showed a significant upward trend during the study period; Tianjin Beidagang Wetland Nature Reserve showed a downward trend, as shown by the reduction of high-value areas and decrease in the highest value; and Tianjin Tuanbo Bird Nature Reserve showed a relatively small change trend. According to functional zoning (Table 6), the mean values of the three functional zoning areas were relatively small, the interannual changes were small, and the regional mean was stable.