Valuing water purification services of forests: a production function approach using panel data from China's Sichuan province

Abstract

The water purification functions of forests represent one of the most frequently invoked examples of nonmarket ecosystem services that are economically valuable. This study quantifies the monetary value of forests' water purification services in the form of the ensuing cost savings of municipal drinking water treatment, using a rich panel dataset from China's Sichuan province. Moreover, this study has undertaken a novel spatial piecewise approach to investigate the spatial patterns of such cost savings delivered by forests at different distances from the water intake point. The estimation results find that forests within a 2 km radius upstream from the water intake point have the most sizeable and statistically significant cost saving effect. For forests within a 3 km radius, this effect becomes somewhat smaller but remains statistically significant. Beyond a 4 km radius, this effect becomes notably smaller and statistically equal to zero. Our analysis facilitates the optimal spatial targeting of forest conservation.

1. Introduction

The water purification functions of forests represent one of the most frequently invoked examples of nonmarket ecosystem services that are economically valuable (Freeman et al., 2014; Vincent et al., 2016). Forests help enhance water quality by reducing soil erosion (and hence reducing silt) as well as filtering out nutrients and pollutants carried in water, which allows the municipal water supply sector to simplify or expedite many costly water treatment procedures and thereby save on operating costs (Millennium Ecosystem Assessment, 2005; Holmes et al., 2017). Despite that, there has been a paucity of statistical estimates that robustly quantify such cost savings (Ferraro et al., 2012; Price and Heberling, 2018; Ovando and Brouwer, 2019), which has precluded inclusion of this essential ecosystem service in large-scale programmes intended to mainstream ecosystem services into national accounts, such as the UK National Ecosystem Assessment (Bateman et al., 2013).

There is a recent and limited body of formal econometric literature that attempts to bridge this gap, such as the studies of Abildtrup et al. (2013, 2015), Fiquepron et al. (2013), Singh and Mishra (2014), Vincent et al. (2016, 2020), Lopes et al. (2019), Westling et al. (2020) and Piaggio and Siikamäki (2021).¹ Among these studies, only Vincent et al. (2016, 2020), Westling et al. (2020) and Piaggio and Siikamäki (2021) undertook fixed effects panel data estimation, and none of these studies adopted the further quasi-experimental approaches or supplementary analyses (such as placebo tests) recommended by Greenstone and Gayer (2009), Imbens and Wooldridge (2009) and Athey and Imbens (2017) to more convincingly identify the water treatment cost savings ‘caused’ by forests, as opposed to those that are in fact induced by other unobserved ‘confounding factors’ but are likely to be mistakenly attributed to forests.²

Moreover, the spatial patterns of such cost savings delivered by forests at different distances from the water intake point remain less understood. Existing studies typically focus on the average effect of forest cover in the entire catchment area (e.g., Singh and Mishra, 2014; Vincent et al., 2016, 2020; Piaggio and Siikamäki, 2021).³ Such estimates suffice for the valuation and accounting of the average or aggregate value of water purification services of forests in each catchment area. However, the provision of such services (via forest conservation) nearly always entails costs to society such as administrative costs of conservation activities and forgone benefits of alternative land use (Pearce, 2004; Armsworth, 2014). Economically optimal land use decision making critically depends on whether the benefits of such services (savings of drinking water treatment costs) outweigh the costs of forest conservation.

Yet, such benefits and costs tend to be spatially heterogenous (Polasky et al., 2008; Vincent et al., 2016). Therefore, the cost-effective provision of water purification services via forest conservation would necessitate a better understanding of the presumably location-specific implications of forest cover for drinking water treatment costs. For instance, intuition might suggest that upstream forest cover within a certain radius of the water intake point is likely to have a larger effect on water treatment costs, relative to forest cover in a farther upstream location, because solids and other pollutants entrained by surface runoff in farther upstream locations may have been naturally removed before they could reach the water intake point, due to sedimentation and other self-purification functions of water bodies. In that case, it might be worthwhile to undertake forest conservation actions within that radius of the water intake point, but the opposite might be true for more distant locations or the entire catchment area. This is particularly relevant if water treatment works extract water from large rivers with massive catchment areas.

This study seeks to contribute to this literature through a spatial piecewise analysis of the effects of forest cover on drinking water treatment costs. We undertook fixed effects estimation using a rich panel dataset that contains annual observations on the drinking water treatment works of 170 county-level administrative divisions of China's Sichuan province during 1992–2015. This allows us to control for many unobserved confounding factors that might bias the estimated water treatment cost savings attributable to forest cover. Further, we performed a Heckman correction and a placebo test to formally assess the implications of missing data and potential confounders for the estimates of interest. This study thus enriches the currently thin evidence base around the monetary value of water purification services of forests.

In particular, our study contributes to the geographic representativeness of this literature by providing perhaps the first formal case study from China, a country that, in the past decades, experienced significant changes in forest cover and rapid expansion of drinking water supply on a large scale. Moreover, the spatial piecewise approach⁴ enables us to identify the farthest upstream distance beyond which forest cover no longer directly affects water treatment costs. We first delineated a sequence of equal distance concentric circles or buffers (with a 1 km step length) surrounding each water intake point. We next measured the forest cover within the overlapping areas between the catchment and each buffer, giving rise to a vector of variables that represent forest cover within different radiuses of the water intake point. These forest cover variables enter the regression model individually, which reveals the threshold radius where the effect of forest cover on water treatment costs just disappears.

Such spatially explicit quantification of forests' water purification services has pronounced implications for developing countries in the tropics that are subject to continuing net losses of natural forests (Song et al., 2018) and limited access to safe drinking water (UNICEF and WHO, 2019). For instance, it was estimated that only 35 per cent of the population in the world's least developed regions was covered by safely-managed drinking water services in 2017 (UNICEF and WHO, 2019). In the upcoming decade, these regions will need heavy investments in sanitised drinking water supply in order to achieve the United Nations' Sustainable Development Goal 6: ensure safe drinking water for all by 2030. Preserving forests' water purification services may help rein in such costs, which might be particularly pertinent for developing regions under tight constraints on public spending.

However, local households in these regions tend to rely heavily on forest extraction as an essential source of livelihoods (Wunder et al., 2018). In these contexts, the within-catchment optimal targeting of forest conservation actions takes on particular importance in terms of reducing drinking water treatment costs on the one hand, and protecting local livelihoods on the other. This study, admittedly, did not collect data from the tropics and the results may not be globally generalisable. Despite that, this study demonstrates a viable empirical approach that facilitates the identification of the forests (within a catchment) that provide the most materialisable water purification services to local residents. This empirical approach constitutes the primary contribution of this study and is transferrable to other geographic contexts.

Our discussion so far has yet to mention other types of ecosystem services provided by forests because our empirical analysis primarily focuses on the water purification service. However, when conducting a cost-benefit analysis for forest conservation or restoration, it is important to bear in mind that the water purification service is only one of a wide range of ecosystem services provided by forests, as discussed by the Millennium Ecosystem Assessment (2005), Binder et al. (2017), Holmes et al. (2017), and Ovando and Brouwer (2019). For example, forests offer provisioning services such as timber and non-timber forest products. The water purification service of forests is typically regarded as a type of regulating service; other examples under this broad category include carbon sequestration, regulating streamflow, recharging groundwater, and mitigating flood risk. Cultural services such as recreational benefits represent another well-known category of ecosystem services delivered by forests.

Moreover, forests provide supporting services such as soil formation, nutrient cycling and the provision of habitats for wildlife. Many of these ecosystem services exhibit spatial heterogeneity or distance decay similar to the water purification service (Siikamäki et al., 2015). For example, provisioning and recreational services are more likely to be utilised in locations closer to human settlements or roads, due to high transport/travel costs in more distant or less accessible locations (Bateman et al., 2006; Hyde, 2012; Siikamäki et al., 2015; Busch and Ferretti-Gallon, 2017). The travel cost approach is a classic methodology that captures the distance decay of forests' recreational services by modelling the negative relationship between distance and the number of visits (Freeman et al., 2014).

Moreover, a few recent hedonic pricing studies (i.e., Conway et al., 2010; Sander et al., 2010) performed the spatial piecewise analysis mentioned above and found that the recreational services of urban forests are capitalised into the prices of houses, but only within a certain radius. Another example is that forests provide habitats for natural pollinators which contribute to the production of pollinator-dependent crops, and such benefits tend to be larger if the forests are in close proximity to farmland plots, because natural pollinators usually have a limited foraging distance. Tibesigwa et al. (2019) explored this spatial pattern using the spatial piecewise approach mentioned above. Our study represents the first empirical study on the spatial pattern of forests' water purification service, which provides a novel addition to this broader literature around the spatial heterogeneity of forests' ecosystem services.

The remainder of this paper is structured as follows. Section 2 describes the study area, data sources and measurement of variables. Section 3 performs the econometric analysis and reports the results. Section 4 discusses these findings and concludes.

2. Study area, data and variables

The geographic focus of this study is China's Sichuan province (figure A1, online appendix A). In 2018, the province had a population of over 80 million and its economy was similar in size to Switzerland in terms of GDP (National Bureau of Statistics of China, 2019). The province's forests are widely considered to have a pivotal role in the delivery of watershed services and biodiversity. Sichuan is one of the only two provinces of China that accommodate the upper courses of both the Yangtze River and the Yellow River, the two largest rivers of the country. Its forests contribute greatly to water inflows to both rivers. Moreover, they provide vital habitats for a variety of endangered species, including the giant panda.

Therefore, central and local government bodies have been heavily investing in forest conservation and restoration programmes, such as the Sloping Land Conversion Programme that retires highly erodible agricultural lands and converts them to woodlands. Public spending on forest conservation and restoration in Sichuan province in 2018 amounted to CNY9.5 b (USD1.4 b) (National Bureau of Statistics of China, 2019), which is a considerable sum, using the annual cost of the US Conservation Reserve Program (USD1.8 b) as a reference point (Hellerstein, 2017). However, less is known about the economic returns of forest conservation and restoration, which adds to the difficulty of setting such activities at the optimal level (where the benefits outweigh the costs). The postulated water purification services of forests would represent a highly perceivable nonmarket benefit of forests that continuously accrues to society.

We next outline the theoretical framework which guided the selection of variables and data collection. The water purification effect of forests is modelled as an input to the production of drinking water, following the theoretical framework of Vincent (2011) and Freeman et al. (2014). The production function is assumed to take the Cobb-Douglas functional form:⁵

(1)\begin{equation}z = \boldsymbol{x}_1^{\boldsymbol{\alpha }}\boldsymbol{x}_2^{\boldsymbol{\beta }}{\boldsymbol{e}^{\boldsymbol{\gamma }}}. \end{equation}

In this production function, x₁ represents a vector of variable inputs such as raw water, chemicals, labour and electricity; x₂ denotes a vector of fixed inputs such as machinery; e consists of a vector of environmental factors that are exogenous to production decisions, such as forest cover and rainfall. Producers (water treatment works) choose the amounts of variable and fixed inputs that minimise the long-run production costs (since our dataset shows considerable adjustments in production assets during the 24-year period it covers), given exogenous levels of output $\tilde{z}$ and environmental factors $\tilde{\boldsymbol{e}}$. Solving this cost minimisation problem gives rise to the conditional factor demands:

(2)\begin{align}\boldsymbol{x}_1^{\boldsymbol{\ast }} & = {\left[ {{{\left( {\frac{{{\boldsymbol{w}_1}\boldsymbol{\beta }}}{{{\boldsymbol{w}_2}\boldsymbol{\alpha }}}} \right)}^{ - \boldsymbol{\beta }}}{{\tilde{\boldsymbol{e}}}^{ - \boldsymbol{\gamma }}}\tilde{z}} \right]^{\frac{1}{{\alpha + \beta }}}}, \end{align}

(3)\begin{align}\boldsymbol{x}_2^{\boldsymbol{\ast }} & = {\left[ {{{\left( {\frac{{{\boldsymbol{w}_2}\boldsymbol{\alpha }}}{{{\boldsymbol{w}_1}\boldsymbol{\beta }}}} \right)}^{ - \alpha }}{{\tilde{\boldsymbol{e}}}^{ - \boldsymbol{\gamma }}}\tilde{z}} \right]^{\frac{1}{{\alpha + \beta }}}},\end{align}

where w₁ and w₂ represent prices of x₁ and x₂, respectively. The unit production cost can thus be expressed as:

(4)\begin{equation}\bar{c} = \frac{c}{{\tilde{z}}} = \frac{{{\boldsymbol{w}_1}\boldsymbol{x}_1^{\boldsymbol{\ast }} + {\boldsymbol{w}_2}\boldsymbol{x}_2^{\boldsymbol{\ast }}}}{{\tilde{z}}} = \Omega \boldsymbol{w}_1^{\frac{\boldsymbol{\alpha }}{{\boldsymbol{\alpha } + \boldsymbol{\beta }}}}\boldsymbol{w}_2^{\frac{\boldsymbol{\beta }}{{\boldsymbol{\alpha } + \boldsymbol{\beta }}}}{\tilde{\boldsymbol{e}}^{ - \frac{\boldsymbol{\gamma }}{{\boldsymbol{\alpha } + \boldsymbol{\beta }}}}}{\tilde{z}^{\left( {\frac{1}{{\boldsymbol{\alpha } + \boldsymbol{\beta }}} - 1} \right)}}, \end{equation}

where $\Omega = {\left( {\frac{\boldsymbol{\alpha }}{\boldsymbol{\beta }}} \right)^{\frac{\boldsymbol{\beta }}{{\boldsymbol{\alpha } + \boldsymbol{\beta }}}}} + {\left( {\frac{\boldsymbol{\beta }}{\boldsymbol{\alpha }}} \right)^{\frac{\boldsymbol{\alpha }}{{\boldsymbol{\alpha } + \boldsymbol{\beta }}}}}$. Taking the logarithm of both sides yields a linear form of the model:

(5)\begin{equation}ln\bar{c} = ln\Omega + \frac{\boldsymbol{\alpha }}{{\boldsymbol{\alpha } + \boldsymbol{\beta }}}ln{\boldsymbol{w}_1} + \frac{\boldsymbol{\beta }}{{\boldsymbol{\alpha } + \boldsymbol{\beta }}}ln{\boldsymbol{w}_2} - \frac{\boldsymbol{\gamma }}{{\boldsymbol{\alpha } + \boldsymbol{\beta }}}ln\tilde{\boldsymbol{e}} + \left( {\frac{1}{{\boldsymbol{\alpha } + \boldsymbol{\beta }}} - 1} \right)ln\tilde{z}, \end{equation}

where the coefficient on $ln\tilde{\boldsymbol{e}}$ represents the elasticity of the unit drinking water treatment cost with respect to environmental factors, which captures the monetary value of forests' water purification services. This theoretical framework provides a basis for the selection of variables in our data collection and empirical analysis.

A full catalogue of the geographic coordinates of 276 urban drinking water intake points in Sichuan province was obtained from the Ministry of Ecology and Environment of China. This has enabled the delineation of the upstream catchment area of each water intake point, within which the surface runoff drains into the water intake point (as illustrated in figure 1). The upstream catchment area (the polygon that encompasses the river networks in figure 1) was delineated from a Digital Elevation Model (DEM) in ArcGIS. We next drew a sequence of equal distance concentric circles (with a 1 km step length) surrounding each water intake point. (The smallest circle surrounding the water intake point has a radius of 1 km, the second smallest circle has a radius of 2 km, and so forth.) These circles were used to define the segments of the catchment at different distances to the water intake point. We then measured the area and percentages of different land use types within the overlapping areas between the catchment and each circle (the shaded areas in figure 1).

Figure 1. Schematic diagram of a water intake point and its catchment.

In contrast, forests in the remaining fractions of the circles outside the catchment (the unshaded portions of the circles in figure 1) will be used as a placebo test, as will be detailed shortly. The land use dataset was sourced from the Climate Change Initiative of the European Space Agency (https://www.esa-landcover-cci.org/?q=node/175), which consists of consecutive annual land use data from 1992–2015 with a spatial resolution of 300 m.⁶ Moreover, following Singh and Mishra (2014) and Vincent et al. (2016), we controlled for rainfall in these buffers, which is likely to correlate with both land use changes and water quality (and hence drinking water treatment costs). The rainfall dataset was obtained from the Climate Research Unit of the University of East Anglia (Harris et al., 2014; Osborn and Jones, 2014), which contains monthly rainfall data for the same period (1992–2015). This dataset has an original spatial resolution of 0.5 degrees and was resampled to a 30 m resolution. The land use and rainfall variables constitute the vector of environmental inputs ($\tilde{\boldsymbol{e}}$) in equation (5).

Annual water treatment data were extracted from the Statistical Yearbook of City Water Supply and the Statistical Yearbook of County Water Supply. These data provided the remaining quantity variables in equation (5), namely the unit water treatment cost $\bar{c}$, the level of water supply $\tilde{z}$ and the value of fixed assets ${\tilde{\boldsymbol{x}}_2}$. Each observation in these yearbooks refers to a water supply firm, which may own multiple water treatment works. This has added to the difficulty of performing the analysis at the water treatment work level. Moreover, many water supply firms were renamed (even more than once) due to reorganisation, acquisition and expansion so forth. over the 24-year-long time span of this study. Therefore, we have opted to conduct the analysis at the county level to circumvent the uncertainties around identifying the same water supply firms under different names. We thus aggregated the land use, rainfall, water supply and asset variables to the county level,⁷ and derived the county-level unit water treatment cost variable as each county's total treatment cost divided by its total water supply. In fact, these ‘aggregated’ observations still mostly contain firm-level information, since only a very small proportion (less than 7 per cent) of the counties in our dataset have multiple water supply firms in one year. Turning to the prices of variable inputs and production assets (w₁ and w₂), the wage rate (or the price of labour) was measured at the county level using data from the Sichuan Statistical Yearbook, whereas the prices of raw water, electricity, chemicals and production assets were controlled for by year fixed effects, which will be further discussed in section 3. In addition, the privatisation of China's urban water supply, as opposed to public management, has been expected to enhance the sector's performance in terms of cost-effectiveness (Jiang and Zheng, 2014; Li, 2018). We have therefore controlled for the proportion of private water supply firms in estimating the water treatment cost function. Panel 1 of table A2 (online appendix A) describes these variables, although for brevity it contains land cover and rainfall variables only for the 3 km radius as an example, since the preferred regression models did not find a statistically significant effect of forest cover beyond the 3 km radius on water treatment costs. Figure A2 in online appendix A presents the means of the land cover variables for all different radiuses involved in the regression analysis.

3. Estimation methods and results

We performed a series of econometric analyses on the spatial patterns of forests' water purification services. This section reports our estimation procedures and main findings.

3.1 Visual evidence

Figure A3 in online appendix A visualises the spatial distribution of unit water treatment cost levels and forest cover within a 3 km radius upstream of the water intake point.⁸ The two maps exhibit a reasonably distinguishable negative association between the two variables: unit water treatment costs tend to be higher (darker areas in figure A3a) where water intake points are surrounded by less forest cover (lighter areas in figure A3b), such as the dashed square areas. In the opposite direction, water treatment costs appear to be lower where water intake points have more forested buffers, such as the dashed circular areas. Moreover, figure A3c suggests a timewise negative correlation between water treatment costs and forest cover. These visually observed negative correlations can be quantitatively substantiated by the correlation coefficient between the two variables (−0.075, p-value < 0.01).

These findings are suggestive of the postulated water purification effect of forests. Yet, the strength of such evidence is rather limited, as the observed spatial and temporal heterogeneity of water treatment costs might be induced by certain unobserved factors other than forest cover. For example, less forested locations might imply higher levels of urbanisation and hence higher prices of labour, in which case the observed higher water treatment costs in these places might be caused by higher labour costs instead of lower source water quality associated with less forest cover. These patterns will be further tested via the more rigorous regression analysis reported below.

3.2 Main analysis

The departure point of our regression analysis is a county-level fixed effects model specified as per the water treatment cost function (equation (5)):

(6)\begin{align}arcsinh({{{\bar{c}}_{it}}} )& = {\boldsymbol{\lambda }_{\boldsymbol{l}}}\; arcsinh({{\boldsymbol{l}_{\boldsymbol{it}}}} )+ {\lambda _{r1}}\; arcsinh({{r_{it}}} )+ {\lambda _{r2}}\; {[{arcsinh({{r_{it}}} )} ]^2}+ {\lambda _z}\; arcsinh({{z_{it}}} )\notag\\ & \quad + {\lambda _w}\; arcsinh({{w_{it}}} )+ {\lambda _p}\; arcsinh({{p_{it}}} )+ {\boldsymbol{\theta }_{\boldsymbol{i}}} + {\boldsymbol{\theta }_{\boldsymbol{t}}} + {\varepsilon _{it}}. \end{align}

(The notation ‘arcsinh’ refers to the inverse hyperbolic sine transformation which approximates the logarithmic transformation. We will further discuss this later.) This specification explicitly contains the following explanatory variables for county i in year t: the level of annual water supply (${z_{it}}$), the wage rate (${w_{it}}$), the percentage of private water works (${p_{it}}$), the percentages of different land use types (${\boldsymbol{l}_{\boldsymbol{it}}}$), and annual rainfall (${r_{it}}$). In particular, runoff from cropland and urban areas tends to carry considerable amounts of agricultural fertilisers and household sewage, which is likely to induce higher water treatment costs.

On the other hand, these land use types are often converted from natural land cover such as forests, and therefore may have a negative correlation with the percentage of forest cover. We thus explicitly controlled for the percentages of cropland and urban areas to avoid potential confounding bias.⁹ Moreover, we attempt to account for other price variables using year fixed effects (${\boldsymbol{\theta }_{\boldsymbol{t}}}$). Prices of raw water and electricity are regulated by the government to be the same across different counties of the province in a given time period. Chemicals (such as clarifying agents and disinfectants) and production assets are likely to be purchased from a single and sufficiently competitive market. We thus assume that prices of chemicals and production assets vary over time but not across water treatment works. These year fixed effects also control for a variety of other time varying factors that are homogeneously faced by all water treatment works, such as inflation and changes in national and provincial environmental and resource policies.

In addition, this specification includes county fixed effects (${\boldsymbol{\theta }_{\boldsymbol{i}}}$) to eliminate potential confounding factors that are location specific but do not vary over time, such as various topographical and geological characteristics including the slope of forestland. Standard errors are clustered at the county level to address unobserved within-county correlation (Cameron and Miller, 2015). Lastly, the logarithmic transformation is approximated using the inverse hyperbolic sine (IHS) transformation, which allows zero values of the land use variables. All estimates are expressed as elasticities derived using the approach recommended by Bellemare and Wichman (2020).¹⁰

Figure 2a presents the estimated elasticities of the unit water treatment cost with respect to forest cover within different radiuses upstream from the water intake point, where the solid circles represent the point estimates, and the crosses and diamonds give the confidence intervals at the 5 and 10 per cent significance levels, respectively. These estimates were derived from 10 regression models that each contain a single forest cover variable for a certain radius, as exemplified by Model 1 in table 1. Figure 2a shows that forest cover within a 2 km radius upstream from the water intake point has the most sizeable effect on drinking water treatment costs (estimated elasticity = –2.48 × 10^–2, p-value < 0.05). For a 3 km radius, this effect becomes somewhat smaller but remains statistically significant (estimated elasticity = –1.84 × 10^–2, p-value < 0.10). The estimate for the 4 km radius is similar to that for the 3 km radius, although the p-value goes slightly above 10 per cent after correcting for the unbalanced panel dataset, which will be further discussed shortly. Beyond 4 km, the estimates become much weaker in terms of both the magnitude and statistical significance (the absolute value of the estimated elasticity <9.06 × 10^–3, p-value > 0.43). We therefore opted to focus on the estimate for the 3 km radius (for the moment).

Figure 2. Elasticity estimates with respect to forestland percentages within different radiuses.

Table 1. Estimated water treatment cost function using land cover percentage variables

Notes: Estimates that are statistically significant in both models are highlighted in bold italics (up to the 10% significance level). Asterisks indicate statistical significance: *p-value < 0.10, ***p-value < 0.01. Standard errors are in parentheses.

Table 1 reports the full estimates for the corresponding model (Model 1). As mentioned above, forests in farther upstream locations (e.g., beyond 4 km) are less likely to influence water treatment costs because sediments and nutrients from those locations are more likely to have settled before reaching the water intake points. That said, it is worth noting that the estimate for the closest radius (1 km) is small in size and statistically insignificant. Many ecological studies (e.g., Allan, 2004; Zhao et al., 2015) found that river water quality tends to have a higher correlation with land use configurations in a radius larger than a few hundred metres because sediments and nutrients can be transported longer (but not unlimited) distances, which could be 2–3 km in our case. Therefore, land use at a very local scale (e.g., within 1 km of the water intake point) may not be able to explain sediments and nutrients that are transported downstream from reasonably close upstream locations (e.g., within a 2 or 3 km radius). This implies that forest cover within a 1 km radius may have lower explanatory power of water quality at the water intake point and hence a smaller effect on water treatment costs. Comparing figures 2a and A5a (online appendix), we have qualitatively similar findings regardless of whether forest cover is measured in percentages or area units.

As can be seen in Model 1, the magnitude of the estimate for the variable ‘IHS of percentage of forestland 0–3 km’ implies that a 1 per cent increase in forest cover within a 3 km radius upstream from the water intake point would decrease drinking water treatment costs by almost 0.02 per cent (at the means of the two variables). (Recall that all estimates have been converted to elasticities following Bellemare and Wichman (2020).) Another way to interpret this elasticity estimate is that a 1 km² increase in forest cover would reduce water treatment costs by CNY0.006/m³ (USD0.001/m³). This is much smaller than the elasticity estimates reported by tropical case studies such as Singh and Mishra (2014) and Vincent et al. (2016, 2020), but comparable to those from temperate regions such as Abildtrup et al. (2013) and Lopes et al. (2019), as shown in table A1 in online appendix A. It appears that forests in the tropics help reduce drinking water treatment costs to a greater extent than those in temperate regions, which is consistent with ecological and hydrological evidence (e.g., Bruijnzeel, 2004; Wei et al., 2005). The aggregate cost savings derived using the total water supply in 2018 amounts to CNY63 m (USD9.5 m). This value, although lower than the province's forest investments in the same year, is provided by a small proportion of the province's forests (those within 3 km upstream of drinking water intake points), and represents only one part of forests' ecosystem services. Moreover, such benefits will continue to accrue in increasing amounts over time, in light of the ongoing rapid urbanisation of Sichuan province and the accompanying expansion of municipal water supply.

In addition, the estimated elasticity of the unit water treatment cost with respect to the scale of water supply is negative and less than one, as shown by the estimate for the variable ‘IHS of water supply’ in Model 1. This implies that the unit water treatment cost decreases less than proportionally with the scale of water supply; in other words, the water supply firms in our dataset exhibit increasing returns to scale.

3.3 Heckman correction for missing data

Nonetheless, our county-level panel dataset is still unbalanced, which might bias the estimates if certain types of counties are more likely to have missing data, in which case the available observations with complete information would become an unrepresentative sample of the study area. For instance, the complete observations in our dataset have an average urban population of 110.66 thousand, which is notably higher than that of incomplete and missing observations (71.60 thousand), and the difference is strongly significant with a p-value below 0.001. This implies that counties with a smaller urban population are more likely to have missing data; this is not surprising because water supply firms smaller than a threshold size (probably serving a smaller urban population) are not legally obliged to report to statistics authorities.

As can be seen in Model 1, the negative and statistically significant estimate on the variable ‘IHS of water supply’ suggests that counties with a higher scale of water supply (which is likely associated with a larger urban population) tend to have lower unit water treatment costs. Therefore, the complete observations in our dataset are likely to underestimate the province's unit water treatment costs, since counties with a smaller urban population (and hence potentially a smaller scale of water supply and higher unit water treatment costs) are more likely to be missing. There are similar yet more subtle implications for the estimates on forest cover. If the available observations constitute an unrepresentative sample of the province, it would be possible that the estimated water purification effect of forests deviates from the overall situation of the province.

We performed a Heckman correction following a two-stage procedure described by Wooldridge (2010) to formally assess whether the estimates in Model 1 are indeed biased due to missing data. The first stage is a probit sample selection equation (Model A1 in table A3 in online appendix A) estimated using data for the entire sample (170 counties × 24 years = 4,080 observations):

(7)\begin{equation}Pr({{y_{it}} = 1\textrm{|}{\boldsymbol{v}_{\boldsymbol{it}}}} )= \Phi ({{\delta_0} + \boldsymbol{\delta }{\boldsymbol{v}_{\boldsymbol{it}}}} ). \end{equation}

The binary dependent variable ${y_{it}}$ equals one if the data point for county i in year t is observed, and zero otherwise. $\Phi ({\cdot} )$ represents the standard normal cumulative distribution function, and ${\boldsymbol{v}_{\boldsymbol{it}}}$ consists of a vector of variables that explain the probability of a data point being observed, which are listed in panel 2 of table A2 (online appendix A). The variable ‘urban population’ captures the scale-dependent obligation of statistical reporting mentioned above. In addition, we included several variables that account for communication costs within a centralised economic system, as per Huang et al. (2017). These variables include each county's distance to the capital city of the province (Chengdu), number of phones per capita, road density and percentage of private water works. Furthermore, the first-stage model contains a binary variable that indicates whether a county is an ethnic minority autonomous county, as China's ethnic minority autonomous divisions are typically less integrated with the country's governing system (Han and Paik, 2017) and therefore may be less responsive to the centralised statistical bureaucracy. These variables were obtained from the Sichuan Statistical Yearbook and hence were available for the entire sample. Lastly, we estimated the inverse Mills ratio¹¹ from the sample selection equation and inserted it into the water treatment cost function estimated using the selected sample.

Returning to table 1, the second column of results (Model 2) presents the selection-corrected estimates of the water treatment cost function that contains forest cover within a 3 km radius upstream. The statistically significant coefficient on the inverse Mills ratio term provides corroborating evidence of the conjectured sample selection bias associated with missing data. However, the consequences of such bias for this particular model are somewhat limited, as the selection-corrected estimates on other regressors mostly closely resemble the uncorrected estimates in Model 1. In particular, there is no substantial change in the estimate for the variable ‘IHS of percentage of forestland 0–3 km’, although its statistical significance becomes slightly weaker. Figure 2b reports the selection-corrected estimates for forests at other distances. As mentioned above, the 4 km radius estimate becomes statistically insignificant, since the upper bound of its 10 per cent level confidence interval goes above zero. Aside from that, there is still no evidence that forests at farther distances have a statistically distinguishable effect on water treatment costs. Therefore, the previous findings and interpretation pertaining to the monetary value of forests' water purification services remain largely robust.

3.4 Placebo tests and robustness checks

We conducted a placebo or falsification test in an attempt to examine whether the foregoing findings stem from some unobserved time-varying factors that are systematically correlated with variation in forest cover. In this test, we replaced the forest cover variable in the water treatment cost model using forests inside the same radius but outside the catchment area (the unshaded portions of the circles in figure 1). These out-of-catchment forests are either downstream from the water intake point or belong to another drainage system, and therefore, intuitively, should have no direct implications for water quality at the intake point. If there exist certain unobserved local factors that covary with both forest cover and water treatment costs, these factors would likely be correlated with forests inside and outside the catchment alike. In that case, regressing water treatment costs against out-of-catchment forests would (falsely) pick up a significant effect. Reassuringly, the estimates of the placebo regressions, which are statistically insignificant throughout all radiuses (as shown in figure 2c), ease this concern to some extent.

Admittedly, the estimates for out-of-catchment forests are not negligible in size relative to those for within-catchment forests. In particular, the negative estimate on out-of-catchment forests in the 2 km radius suggests that higher levels of out-of-catchment forest cover correlate with lower water treatment costs, which cannot be explained by forests' water purification effect and is hence likely confounded by certain unobserved factors. In that case, the estimates on within-catchment forests in the 2 km radius (as shown in figures 2a and 2b) are likely subject to similar bias and should be taken with caution, since these estimates have likely overestimated forests' water purification effect. In contrast, the positive estimate on out-of-catchment forests in the 3 km radius implies that the estimates on within-catchment forests in the same radius have likely underestimated forests' water purification services due to omitted factors. This is another reason that we opted to focus on estimates for the 3 km radius (rather than those for the 2 km radius which have higher magnitudes and statistical significance levels), since conservative estimates help reduce Type I error (or false positive findings) and are thus preferable for statistical hypothesis testing. Despite the high correlation between forest cover inside and outside the catchment area (as shown in figure A4 in online appendix A), out-of-catchment forest cover was not found to have a statistically discernible effect on water treatment costs, which provides additional credibility for our main results.

In addition, we formally assessed the implications of potential heterogeneity in the quality of treated water. Drinking water supply firms in China are required to treat water to the same set of national standards, and the quality of treated water is subject to routine self-monitoring and government inspections. Despite that, these standards and regulations tend to be loosely enforced, and the quality of treated water is likely to vary across regions. This is because economic growth is geographically unbalanced throughout the country, and less developed regions tend to have financial difficulties treating drinking water to the national standards, since treating drinking water to higher quality usually incurs higher treatment costs (Browder et al., 2007; Jiang and Zheng, 2014; Li, 2018). If these less developed regions are also less urbanised and therefore have higher levels of forest cover, omitting the quality of treated water in the regression analysis would confound the estimate on forest cover.

Our dataset contains water supply firms' self-reported percentage of tested water samples that achieve the quality standards, which shows some degree of heterogeneity: this percentage ranges from 80.71 to 100 per cent, yet is above 98 per cent for 90 per cent of the observations that contain complete information for the regression analysis. Grieco and McDevitt (2017) developed a novel form of the production function which accounts for both the quantity and quality of output. As described in online appendix B, this model basically assumes that producing the same amount of output with higher quality requires a higher level of aggregate production input. This production function allows the quality of treated water to be explicitly controlled for in the regression models as an explanatory variable. We next re-estimated all regression models controlling for the quality of treated water, and the estimates turned out to be largely stable, as can be seen in figure A7 in online appendix B; they are almost indistinguishable from their counterpart estimates from the models that do not account for the quality of treated water (figure 2). Table A5 in online appendix B presents an example of these new models (Model B1) which has the same specification as Model 2 in table 1 (and hence focuses on forest cover in a 3 km radius), except that Model B1 controls for the quality of treated water. It can be seen that the two models give very similar estimates for all the common regressors. Admittedly, our data on the quality of treated water were self-reported by water treatment firms and hence may not fully reflect the actual variation in quality. Despite that, the placebo test mentioned above has implicitly accounted for unobserved heterogeneity in the quality of treated water: if it has a substantial confounding effect (jointly with other unobserved factors), this should have been captured by the estimates for out-of-catchment forest cover.

So far, we have been following the approach of Tibesigwa et al. (2019) and focusing on a series of models which separately look at the effects of forest cover in different radiuses one at a time (i.e., 0–1 km, 0–2 km and 0–3 km, etc.). This approach differs from another type of spatial piecewise analysis which includes in a single model all forest variables for different distances (i.e., 0–1 km, 1–2 km and 2–3 km, etc.). We did not opt for the latter approach, primarily due to concerns about multicollinearity issues among forest variables for adjacent ‘rings’ or ‘bands’ (e.g., 1–2 km and 2–3 km, or 5–6 km and 6–7 km).¹² High multicollinearity would lead to implausible coefficient estimates and inflated standard errors (Greene, 2020). That said, we estimated another version of Model 2 which still focuses on forest cover for the 3 km radius but controls for forest cover in farther rings (i.e., 3–6 km and 6–10 km). The band-widths of the farther rings are wider than 1 km, which may help reduce multicollinearity. The estimates of this model are shown by Model A2 in table A4 (online appendix A). The estimate for the 3 km radius (−2.91 × 10^–2, p-value < 0.01) remains reasonably stable compared to its counterpart in Model 2 (table 1), which leaves out forest cover in farther rings (−1.67 × 10^–2, p-value < 0.10), although the estimate in Model A2 appears stronger in terms of the magnitude and statistical significance. Regarding forest cover in farther rings (i.e., 3–6 km and 6–10 km), the estimates in Model A2 are statistically insignificant (p-value = 0.31 and 0.55, respectively). However, the magnitudes of those estimates are much larger (in absolute value) than those in figure 2b, where the models do not contain multiple forest variables, and the estimate for the 3–6 km ring in Model A2 has the opposite sign. These swings in estimates likely stem from high multicollinearity and correlation between the two forest variables for farther rings: the variance inflation factor of the two forest variables is above the conventional rule of thumb (10) and the correlation coefficient is close to one (0.83, p-value < 0.001), which suggests substantial multicollinearity issues.

Finally, we tested the robustness of our findings using the logarithms rather than the IHS transformations of the variables. The IHS transformation is becoming increasingly popular in applied econometrics due to its properties that approximate those of the logarithmic transformation whilst being straightforwardly applicable to zero-valued observations (Bellemare and Wichman, 2020). Despite that, using logarithms in our econometric analysis was better in line with the theoretical model presented in section 2. In addition, many previous econometric studies on the value of forests' water purification services used the logarithms of the variables and accommodated zero-valued observations by adding a very small value to the original variables (e.g., Vincent et al., 2016, 2020; Westling et al., 2020; Piaggio and Siikamäki, 2021).

Therefore, it is useful to investigate whether our results are sensitive to the choice between the IHS and the logarithmic transformation. We repeated all our main regression models and the placebo tests using the logarithms of the variables, where we added a small value (0.001) to the original variables that contain zero values. The estimates pertaining to forest cover are graphically reported in figure A6 in online appendix A. It can be seen that these estimates are highly similar (in terms of both the magnitude and statistical significance) to their counterparts derived from the regression models using IHS transformations (as shown in figure 2). For example, for the model that uses logarithms and corrects for sample selection, the estimate for forest cover within a 3 km radius is −1.42 × 10^–2 (p-value < 0.10), which closely resembles its counterpart from the model using IHS transformations (−1.67 × 10^–2, p-value < 0.10). In fact, the two sets of models yield highly comparable estimates for all the explanatory variables, which demonstrates the stability of our results regardless of the choice between the IHS and the logarithmic transformation. (Further details are available upon request.)

4. Discussion and conclusions

This study robustly measures the monetary value of forests' water purification services in the form of the ensuing cost savings of municipal drinking water treatment. This was enabled by a rich panel dataset from China's Sichuan province, which allowed us to adopt the fixed effects approach to control for a variety of observed and unobserved factors that might otherwise have biased the estimates of interest. This study thus adds to the currently thin formal econometric evidence base in this regard, which was recently reviewed by Price and Heberling (2018) and Ovando and Brouwer (2019), and typically represented by the study of Vincent et al. (2016). Moreover, this study has undertaken a novel spatial piecewise approach to investigate the spatial patterns of such cost savings delivered by forests in different segments of the catchment area.

We found that forests inside a 2 km radius upstream from the water intake point have the strongest effect of reducing drinking water treatment costs. For a 3 km radius, this effect becomes somewhat smaller but remains statistically significant. Forests within a 4 km radius have a similar sized effect to the 3 km radius, but the p-value is marginally above the conventional significance levels. Beyond 4 km, the estimated effect becomes notably smaller and statistically equal to zero. These findings provide suggestive practical implications for the optimal spatial targeting of forest conservation efforts. The magnitudes of our estimates are comparable to the results of previous studies from temperate regions such as Abildtrup et al. (2013) and Lopes et al. (2019), but smaller than those from the tropics such as Singh and Mishra (2014) and Vincent et al. (2016, 2020). Our results are robust to a series of alternative methods of specifying the regression model and measuring the key variables. The credibility of our estimates is further strengthened by a placebo test which found no statistically discernible effect for forests outside the catchment area.

The aggregate water treatment cost savings delivered by forests in the study area amount to CNY63 m (USD9.5 m) in 2018. This value is only moderate relative to the scale of the province's economy and forest investments. However, the primary contribution of this study is to demonstrate the spatial piecewise approach, which helps identify the heterogenous water purification services provided by forests at various distances from the water intake point. This approach also facilitates a placebo test utilising out-of-catchment forest cover in the same radius of the water intake point, which provides important insights as to the direction and magnitude of potential omitted variable bias. This methodological twist is transferrable and applicable to other regions where the optimal spatial targeting of forest conservation is particularly relevant, such as less developed tropical regions where local livelihoods tend to heavily rely on extraction of forest resources (which implies a higher opportunity cost of forest conservation).

Moreover, forests' water purification services may provide other benefits aside from savings of drinking water treatment costs. For instance, in our study area, more than two-thirds of the province's population still rely on untreated water from natural sources. These untreated water users are likely to benefit from forests' water purification services in the form of, for example, reduced exposure to waterborne diseases (Herrera et al., 2017). Yet such benefits are not captured by the production function approach in this study, which focuses on centralised drinking water treatment. Additionally, as mentioned in the introduction, forests provide a wide range of ecosystem services aside from water purification. Therefore, the monetary value of water purification does not represent the entirety of forests' benefits when assessing whether forest conservation efforts are economically worthwhile.

Lastly, the focal point of this study is the spatial piecewise analysis which seeks to identify the spatial patterns of forests' water purification services. Due to resource constraints, we were unable to fully explore other important research questions around this topic. For example, previous studies have conjectured that forests' water purification services depend on rainfall levels (Calder et al., 2007; Vincent et al., 2016), which implies that the future scale of this ecosystem service may be affected by climate change. We conducted some preliminary analysis (see online appendix C) but did not find any evidence for this hypothesis. However, this is likely associated with the nature of our annual-level dataset which has precluded us from making full use of the variation of rainfall and water treatment costs at a higher temporal resolution. Moreover, we did not investigate whether different types of forests (e.g., natural versus planted forests) may deliver different levels of water purification services. This would have been particularly informative in the context of Sichuan, because the province has been heavily investing in both forest conservation and restoration. Further research is warranted in light of the importance of these research questions.