Léa Tardieu
email: tardieu@centre-cired.fr

Laetitia Tuffery
email: laetitia.tuffery@inra.fr

From supply to demand factors: what are the determinants of attractiveness for outdoor recreation?

Keywords: Forest recreational attractiveness, Supply and demand factors, Travel cost method, Function benefit transfer, Regional park strategy Q26, Q51, R12, R58

The concept of ecosystem services (ES) remains underused in supporting practical decisions in conservation/development plans and programs. One of the most important factors explaining this non-consideration is the lack of spatial information describing the naturesociety relationship in environmental and economic analyses. In this paper, we developed a novel method to predict, in spatially explicit terms, the recreation attractiveness potential combining supply and demand factors. Our method is based on the combination and transfer of a Lancasterian function of biophysical aspects and a travel cost model based on agents' sociodemographic characteristics. We further validate the usefulness of the proposed recreation model by using it in the evaluation of a regional park charter pursuing two main objectives: recreational attractiveness and habitat quality (modeled with InVEST). The results demonstrate first that the biophysical context plays a large role in the recreational trip choice and thus should not be ignored in travel cost studies. Second, from a policy guidance perspective, we show that providing spatial information appears particularly critical for ES to be a useful lever for action in day-to-day decision-making.

INTRODUCTION

Despite the growing political attention paid to social dependencies of human beings on natural systems, these dependencies remain under-considered in practical decisions with regard to conservation and management plans [START_REF] Laurans | Use of ecosystem services economic valuation for decision making: Questioning a literature blindspot[END_REF]. One of the most important factors explaining this non-consideration is the lack of spatial information describing the nature-society relationship in environmental and economic analyses. Ecosystem services (ES) modeling and mapping is now recognized as a unifying approach for different research communities to support the consideration of this link for decisionmaking (Bateman et al. 2014, Burkhard and[START_REF] Burkhard | Mapping Ecosystem Services[END_REF]. Refined techniques for mapping the supply side (representing ES flows produced by the ecological functioning of ecosystems), the demand side (flows to beneficiaries) or different values associated with ES have been developed. By providing spatially explicit representations of ES, land managers have the opportunity to broaden the goals that can be pursued when evaluating land management scenarios [START_REF] Geneletti | Handbook on biodiversity and ecosystem services in impact assessment[END_REF]. For instance, winners and losers of potential land management decisions affecting ES can be more clearly identified and located by using these types of analysis.

Nevertheless, shortcomings remain in this field because of a lack of combined consideration of both ecological and social aspects of ES. This leads spatial models to be based predominantly on biophysical attributes, ignoring individuals' preferences for these attributes, or to a lesser extent, based on socioeconomic drivers while ignoring ecological flows of ES [START_REF] Tardieu | The need for integrated spatial assessments in ecosystem service mapping[END_REF]. A good example for the latter case is the recreational service for which, in contrast to other ES in the literature, the modeling of the demand side is much more developed than that of the supply side [START_REF] Termansen | Modelling and mapping spatial heterogeneity in forest recreation services[END_REF]. Since Hotelling's letter in 1947, recreation demand has been historically assessed by modeling visits with the travel cost method. The method enables a recreation demand function to be derived for predicting a participation decision and the number of trips taken to a site according to an implicit price (the travel cost) and a set of other variables (time available, income, other available sites as substitutes, visitors' socioeconomic characteristics, etc.) (e.g., [START_REF] Shresta | Value of recreational fishing in the Brazilian pantanal: a travel cost analysis using count data model[END_REF]Parsons, 2003;[START_REF] Martínez-Espiñeira | Recreation demand analysis under truncation, overdispersion, and endogenous stratification: An application to Gros Morne National Park[END_REF][START_REF] Bujosa Bestard | Estimating the aggregate value of forest recreation in a regional context[END_REF][START_REF] Jones | Estimating Arrival Numbers for Informal Recreation: A Geographical Approach and Case Study of British Woodlands[END_REF][START_REF] Binner | Valuing the social and environmental contribution of woodlands and trees in England, Scotland and Wales[END_REF]. The advantage of the method is that it relies on real individuals' behaviors (revealed preference) and allows the assessment of the actual benefits derived from existing policies or management types (from an ex-post perspective). However, the produced results are often (1) aggregated numbers of visits that do not allow a spatial differentiation of recreation trips within the destination site and (2) focused on individuals, while underestimating the importance of the biophysical characteristics of destination sites.

Other approaches enable the recreational service to be modeled and mapped. [START_REF] Colson | Map and determinants of woodlands visiting in Wallonia[END_REF] relied on a survey in which forest managers reported the visitation frequency to forests in Wallonia, enabling the attributes explaining forest attractiveness to be determined. Other new methods based on volunteered geographic information such as self-geotagged photographs on web applications directly uploaded by visitors (e.g., Flickr©) may allow the development of spatially sensitive recreation maps. Finally, [START_REF] Zulian | ESTIMAP: A GIS-based model to map ecosystem services in the European Union[END_REF] developed the ESTIMAP model to map the recreation potential across Europe. However, these models do not rely on individuals' recreation preferences. Instead, recreation potential often assumes implicitly that visitors' preferences are the same regardless of the sociodemographic characteristics of potential users, with the main driver of recreation being the distance/accessibility, and, more importantly, that every individual (located at the right distance from the recreational site) is a potential user. According to the wide literature in economics on this subject, this assumption is far from being verified. Moreover, no demand evolution, consequent to a change in sociodemographic context, can be predicted from these models. [START_REF] Termansen | Modelling and mapping spatial heterogeneity in forest recreation services[END_REF], [START_REF] Johnston | Introduction to Benefit Transfer Methods[END_REF] and [START_REF] Czajkowski | Spatial Heterogeneity of Willingness to Pay for Forest Management[END_REF] mapped recreational values by transferring a function derived from a discrete choice experiment (DCE) study applied from national surveys. De [START_REF] De Valck | Outdoor recreation in various landscapes: Which site characteristics really matter?[END_REF] also mapped site quality scores from a distance-based DCE. However, although these studies consider visitors' characteristics and preferences, they do not characterize the actual demand in terms of visits, provided that the DCE is based on hypothetical situations of policy initiatives or management alternatives (from an ex-ante perspective). In one exception, [START_REF] Sen | Economic Assessment of the Recreational Value of Ecosystems: Methodological Development and National and Local Application[END_REF] 3 took into account the spatial variability of recreational services, based on actual visits, and combine supply and demand factors in their study. The authors developed a model to predict visits, obtained from a national recreation survey, according to destination site characteristics, the attributes of outset locations and the distance traveled by visitors to reach the site. They further applied the estimated model to the national scale in the United Kingdom. Nevertheless, whether these types of models are applicable to more local territories, requiring more precise output data to develop accurate strategies at the local scale remains an open question.

In this paper, we rely on existing and observed data to derive recreation functions, according to supply and demand factors, in spatially explicit terms at the local level to inform decisionmaking. We develop a methodology to predict, in a spatially sensitive way, the recreation supply and demand of forest ecosystems. The developed method can be applied to any other ecosystem type. The predictions are based on a Lancasterian function of biophysical supply and a travel cost method to inform visitation preferences. We further upscale the supply and demand functions to the studied area. In this way, this study constitutes one of the first attempts to develop a function benefit transfer based on a travel cost method enabling transferring a recreation function adjusted for on-site attributes and individual preferences. We apply the method to a real case study of a regional park in France to test its applicability. Then, we validate the usefulness of the created output by using it to evaluate the regional park charter. Because the park charter is developed around two objectives, biodiversity conservation and recreational attractiveness, we add a habitat quality index in the evaluation, which is computed with the InVEST 4 model. We show that site attractiveness is driven both by the biophysical attributes of forests and the socioeconomic characteristics of individuals, the former having a stronger effect. This finding calls for a better accounting of these attributes in travel cost methods that tend to ignore key factors of recreational attractiveness or to better account for individual preferences in studies focusing on supply factors. Second, we show that the created indexes, expressed as geographic information, can provide policy guidance for territorial management planning by allowing overlapping layers to inform multiple objectives (i.e., biodiversity conservation and region attractiveness).

MATERIALS AND METHODS

Snapshot of the methodological framework

The methodological framework can be divided into three steps (Figure 1). The first step is a recreation model, described in section 2.3. This model is based on an initial survey conducted in the Lorraine region (section 2.3.1), further combined with spatial data. A Lancasterian function of attractiveness based on supply factors and a travel cost model are estimated to derive a combined attractiveness index. Simultaneously, a habitat quality model is computed using the InVEST software (section 2.4). The results of the two models are detailed and analyzed with regard to the Ballons des Vosges Regional Park (BVRP) orientation strategy to assess whether it is accurate and efficient (section 3).

Case study: The Ballons des Vosges Regional Park, France

The Ballons des Vosges Regional Park (BVRP) is a 2700 km 2 territory situated in two regions of eastern France: Grand-Est and Bourgogne Franche-Comté. Forest covers approximately 61% of the total area of the BVRP, and the remaining territory is covered by wetlands, lakes, agricultural areas and meadows (Figure 2). The mountainous relief supports a high diversity of natural habitats governed by different bioclimatic stages. The High Vosges are obstacles to oceanic disturbances that come from the west, producing heavy rains in the western part of the BVRP (more than 2 m/year) to the crest line and, conversely, approximately 50 cm/year in the eastern part, making this side of the park one of the driest regions in France. One-third of BVRP territory has been placed under "remarkable natural site" status and 22% under Natura 2000 conservation commitments5 . A few remarkable species are present in the site, such as Western capercaillie (Tetrao urogallus) and the lynx (Lynx lynx), which are subject to particular conservation programs. In terms of recreation, the BVRP is a popular site for skiing and hiking activities.

Natural Regional Park is a protection status in France 6 . Such parks are governed by a charter approved by the states and municipalities composing the area. The charter sets the development strategy in terms of the maintenance/improvement of environmental and cultural heritage and the means to implement it over a 15-year period. According to the charter, the environmental strategy must be focused on two principal activities: (1) recreational attractiveness and (2) habitat and biodiversity conservation. To that end, land managers have developed a management strategy, specified in the 2024 horizon charter, for three different territories (Figure 2). First, in the High Vosges territory, the objective is to ensure habitat conservation while maintaining good recreational attractiveness. More specifically, managers rely on the European Charter for Sustainable Tourism in Protected Areas7 and try to define tourist attendance strategies. Second, the key challenge in the Valleys and Piedmont is to control urbanization and habitat fragmentation. Finally, the 1000 Ponds Plateau benefits from exceptional habitat richness; however, declines resulting from industrial and farmland activities are weakening the attractiveness of this area. Thus, the objective in this area is to sustain its vitality.

Recreation model

Original survey and scales of analysis

We rely on an online survey initially conducted in the former Lorraine region (now the Grand-Est region), covering one-third of the BVRP territory (Figure 3). This survey has been used in different cases to study local recreation by Abildtrup et al. (2015a and2015b). The survey (a) (b) was carried out between July and August 2010 by email in the former Lorraine region. A total of 1144 respondents completed the survey, and 526 had visited a forest and provided information about which specific forest they had visited in the 12 months before the survey (i.e., July/August 2009 -July/August 2010). All of the survey details are reported in Abildtrup et al. (2015a and2015b).

In the initial survey, Lorraine's forests were divided into "forest units" representing relevant recreational units equal to or greater than 5 ha. A total of 5568 forest units were delineated in the initial survey. From this total sample, we selected the forest units included in the BVRP, which are in the Vosges department in the Lorraine region, resulting in 256 forest units. In total, 1236 visits were recorded in 94 of the 256 forest units and distributed throughout the entire area of the BVRP, surrounded by a 10 km buffer zone8 (Figure 3). Our recreation model is developed based on those 94 forest units. However, because we had no clear rules to define forest units in the rest of the BVRP area, we defined recreational units through a raster distributed homogeneously according to a kilometric mesh (1*1 km). Based on the French forest database (BD Forêt®, IGN), we assigned the surface of forest areas to each mesh. Only meshes including at least 50% of the closed canopy forest were considered as forest recreational units in this work. Meshes cover the entire area of the park, expanded through a 10 km buffer zone around its perimeter to take into account possible visitors who live just outside the park and, therefore, to avoid border effects. Ultimately, we obtained 3774 forest recreational units over the entire BVRP (Figure 3). These recreational units were used to upscale recreational functions.

Attractiveness model based on supply factors

In a typical travel cost model, visited sites are directly revealed by visitors during the survey. However, in this paper, such locations are unknown in a large part of the BVRP and have to be predicted. The attractiveness model based on supply factors characterizes the supply of recreational forests in the area with a Lancasterian function of biophysical characteristics and forecasts the most interesting characteristics. The prediction model achieves this forecast by taking into account information on the number of visits to destination sites from the initial survey and relating this information to (i) the forest structure characteristics, (ii) the amenities available at the site, and (iii) the distribution of the population around the site (Table 1). The expectations concerning these variables are that elevation; coniferous, broadleaved and mixed forests; and water surfaces are more attractive, as are sites strongly endowed with amenities. The population living within 2 km of the forest represents the forest recreational pressure. We expect that when the pressure is substantial, the forest is less attractive. These variables were not present in the initial survey and have been reconstituted for the purposes of the paper. The attractiveness based on supply factors of each recreational unit is predicted with a count data model. These models are particularly accurate when the dependent variable is an integer that takes few values, such as visitors' trips to a destination site (Shaw, 1988;[START_REF] Englin | Estimating social welfare using count data models: an application to long run recreation demand under conditions of endogenous stratification and truncation[END_REF][START_REF] Baerenklau | Spatial allocation of forest recreation value[END_REF][START_REF] Roussel | Recreation demand analysis of sensitive natural areas from an onsite survey[END_REF].

Many pixels are not visited. We addressed this issue by using a zero-inflated count data model, where the probability of participation is estimated simultaneously with the visit function [START_REF] Cameron | Econometrics Models based on count data: comparisons and applications of some estimators and tests[END_REF][START_REF] Gurmu | Excess Zeros in Count Models for Recreational Trips[END_REF]. Zero-inflated count data models are more general than typical count models in that they relax the restriction that an identical process generates both zeros and positive integers. Moreover, as [START_REF] Englin | Estimating Forest Recreation Demand Using Count Data Models[END_REF] argue and as it is the case in our dataset, zero trips can be generated by both a binomial process (for people not in the market, i.e., not visiting any forests) and a Poisson process (for people in the market, i.e., surveyed visitors who did not take any trips in the specific forest in question).

In contrast to simple-hurdle models, zero-inflated models account for both types of zeros by following a two-step procedure. First, a participation model determines whether the zero observation belongs to the group where the dependent variable is always null (i.e., nonvisitors) or whether it belongs to the group in which the dependent variable can be positive or null (i.e., visitors). Furthermore, a Poisson or a negative binomial regression estimates the visit function. Thus, the density function with a Poisson distribution for each forest is as follows [START_REF] Englin | Estimating social welfare using count data models: an application to long run recreation demand under conditions of endogenous stratification and truncation[END_REF]:

= | = + 1 - e = 0 1 - e λ ! > 0 1
where is the probability that a visitor coming from an outset location will take no trip to forest and 1 -is the probability that follows a Poisson law with a λ parameter (see [START_REF] Roussel | Recreation demand analysis of sensitive natural areas from an onsite survey[END_REF] for more explanations on the standard Poisson model).

For the λ parameter representing the average number of visits, is dependent on a vector of explanatory variables that can be different from , the vector of the explanatory variables for λ . Participation can be either assessed with a logit or a probit model, depending on the distribution we assume.

Attractiveness model based on demand factors

Attractiveness based on demand is derived from the initial survey by applying a travel cost model. This model aims to predict the trips taken to each forest recreational unit from any given outset area. To allow the trip function to be transferable, we limited the number of variables to those that can be reconstituted on non-surveyed individuals (e.g., income, distance from the outset location to forest). Moreover, as [START_REF] Bateman | Making Benefit Transfers Work: Deriving and Testing Principles for Value Transfers for Similar and Dissimilar Sites Using a Case Study of the Non-Market Benefits of Water Quality Improvements Across Europe[END_REF] suggest, this type of parameterization excluding ad hoc variables is recommended because of the multiplicative role of coefficients resulting in major transfer errors in over-parameterized models. We thus derive a recreation demand function from the travel cost model determined by an implicit price (the trip cost) and a set of independent variables (generic predictors of individual demand based on economic theory such as the availability of substitutes and income).

The trip cost is a combined cost between trip costs for individuals using motorized means of travel and an opportunity cost of time. The trip cost is computed as follows:

"# = 2 × & ' ( × # ) + *#"+ 2
where "# is the trip cost and ' ( is the distance between the outset location and the site9 in km. # ) is a cost per km based on the vehicle fiscal power published annually by the fiscal administration. Because we did not have this information, we made the assumption that vehicles had a mean fiscal power of 4 fiscal horsepower (the highest proportion of registered vehicles in France since 2000 have between 4 and 5 fiscal horsepower10 ), corresponding to a cost of 0.493/km11 . This cost takes into account vehicle depreciation, maintenance costs, tire expenses, fuel consumption and insurance premiums. is the number of individuals in the group. The costs are multiplied by two to consider the entire round trip. We did not consider foot and bicycle travel, even though some studies considered material depreciation in terms of cycle or trekking shoes (e.g., Bertram et al., 2017). *#" is the opportunity cost of time. An individual who visits a recreation site has an opportunity to use his time differently (working, for example) and is therefore subject to an opportunity cost. This cost relies on the assumption of an individual's trade-off between labor and leisure. However, we chose not to consider an *#" for two reasons: (i) because we rely on very local recreation, implying short trips to forests, and (ii) because this approach assumes that individuals have flexible jobs and are able to substitute work for leisure time at the margin, and this assumption is rarely verified 12 .

As in the attractiveness model based on supply factors, visits are predicted with a count data model where the dependent variable is the number of observed visits in each forest unit .

The observations come from the initial survey. The dependent variables are recomputed to serve our trip function transfer and are described in Table 2. The expectations regarding these variables are that the trip cost between the respondent's outset location and the visited site lowers visit demand as well as the availability of substitutes around the outset location and that the income increases recreation demand.

Benefit transfer and construction of the combined attractiveness index

Two transfer functions are completed to map the combined attractiveness index (CAI) of BVRP forests:

13 IRIS (Ilots Regroupés pour l'Information Statistique) is the smallest level of French national statistics, and it includes approximately 2000 persons. 14 The 5 km buffer was chosen following [START_REF] Sen | Economic Assessment of the Recreational Value of Ecosystems: Methodological Development and National and Local Application[END_REF], who tested different radiuses (1 km, 2.5 km, 5 km, and 10 km). The 5 km radius had the best fit according to the Akaike Information Criterion (AIC).

(1) A transfer of the function developed in the attractiveness model based on supply factors that leads to the development of an attractiveness supply index (ASI); and (2) A trip transfer function predicting the visits to BVRP territory, allowing the development of the attractiveness index based on demand factors (ADI).

In the case of the attractiveness model based on supply factors, the significant coefficients are applied to each 1 km 2 mesh of recreational units in the park. The result, after normalization, gives the ASI based on the characteristics sought by visitors. The ASI for each cell varies between [0; 1] and is calculated as follows:

/01 = /0 -/0 2 3 /0 245 -/0 2 3 3

where

/0 = exp 9 , ∀ 4
/0 is the attractiveness score based on the supply of recreational cell , with ∈ =1, … , 1? before normalization. is a vector of the supply independent variables (the same variables as those described in Table 1) describing recreational forest unit in the entire park, and 9 is a vector of the estimated parameters associated with the vector in the attractiveness model based on supply factors (specified as a zero-inflated Poisson (ZIP) model, which is why the functional form of the equation is semi-logarithmic).

Furthermore, the demand function is applied to construct the ADI applied to the variables described in Table 2. /' is the attractiveness based on demand factors of recreational cell from a series of outset locations present in the BVRP, with ∈ =1, … , 1? before normalization. As in the common travel cost, it also denotes the number of visits to cell . /' is the demand for cell from a given visitor coming from an outset location , with ∈ =1, … , P?. F is the vector of demand independent variables describing visitors and the outset area characteristics, i.e., the percentage of various land use types within a set radius of the outset location. "# is the travel cost between recreational cell and outset . All the variables are described in Table 2.

E and D are the vector of parameters associated with the trip cost and other variables estimated in the attractiveness model based on demand factors.

Finally, the CAI is attributed to each recreational unit by computing the geometric mean of the two indexes as follows:

#/1 = Q/01 × /'1 7

Using a geometric mean makes it possible to consider non-compensation between the two indexes. That is, the ASI and the ADI cannot completely compensate for each other in the CAI if the ASI is strong and the ADI is weak or vice versa. Therefore, a strong combined index reveals strong attractiveness of supply factors combined with a strong attractiveness of demand factors (which is symmetrical for a low CAI).

Habitat quality index

Several tools and models to spatially assess natural habitat quality are already available; they are usually based on ecological indicators [START_REF] Maes | Mapping ecosystem services for policy support and decision making in the European Union[END_REF] or on ecological modeling, such as Globio and Marxan [START_REF] Alkemade | GLOBIO3: A Framework to Investigate Options for Reducing Global Terrestrial Biodiversity Loss[END_REF][START_REF] Chan | Ecosystem Services in Conservation Planning: Targeted Benefits vs. Co-Benefits or Costs?[END_REF], or on land use land cover change modeling, such as InVEST [START_REF] Kareiva | Natural capital: Theory and practice of mapping ecosystem services[END_REF][START_REF] Terrado | Model development for the assessment of terrestrial and aquatic habitat quality in conservation planning[END_REF][START_REF] Sallustio | Assessing habitat quality in relation to the spatial distribution of protected areas in Italy[END_REF][START_REF] Salata | Mapping Habitat Quality in the Lombardy Region, Italy[END_REF]. We use the InVEST model of habitat quality as a proxy to represent biodiversity richness. The aim is to estimate the quality of habitats in spatially explicit terms by measuring the appropriate conditions for the survival and reproduction of species generally based on the potential threat and degradation of natural habitats (Tallis et al., 2011).

Land cover is composed of seven categories extracted from the Corine Land Cover database (CLC, 2012): artificial areas (CLC1), agricultural lands (CLC2), coniferous forests (CLC311), broadleaved forests (CLC312), mixed forests (CLC313), semi-natural areas (CLC32), and water bodies (CLC4-CLC5). Thus, only forests are considered and more precisely identified (i.e., three categories of forestland cover), as our analysis is based on only the habitat quality of forest areas.

Habitat degradation and quality are functions of threats from different sources, e.g., artificial areas, agricultural lands, main roads, secondary roads, and trails. Artificial areas and agricultural lands are included as sources of landscape fragmentation and habitat degradation [START_REF] Girvetz | Integration of landscape fragmentation analysis into regional planning: A statewide multi-scale case study from California[END_REF]. In addition, a few studies have revealed the important impact of forest roads on biodiversity and forest habitat (e.g., [START_REF] Marcantonio | Biodiversity, roads, & landscape fragmentation: Two Mediterranean cases[END_REF]. Consequently, analysis of road networks according to their categorization (primary, secondary and trails, derived from the French Topographic database, BDTOPO®) is critical to understanding the role of linear infrastructure on habitat loss, fragmentation, and degradation [START_REF] Underhill | Effects of roads on wildlife in an intensively modified landscape[END_REF][START_REF] Der Lippe | Do cities export biodiversity? Traffic as dispersal vector across urban-rural gradients[END_REF].

The habitat quality index depends on four factors: ST, the impact of the threats on habitat (i.e., relative impact of each threat to other threats on habitat quality, ranging from 0 to 1); UVW. ', the distance between habitat and threats (i.e., maximum distance, in km, under which each threat affects habitat quality) 15 ; X , the habitat suitability (classification, from 0 to 1, of land-covers as habitat); and 0 T, the sensitivity of habitats to threats (ranging from 0 to 1).

Based on these different factors, the degradation level ' is defined as follows: 15 The distance decay is defined by the InVEST model as exponential [START_REF] Tallis | InVEST 2.2.2 User's Guide, The Natural Capital Project[END_REF].

' = J J Y ST ∑ [ \ ] \LM ^T_ `a _LM ] \LM W \_ \ 8
where c is a mesh of the threat raster T, dT represents all the meshes of the raster T, ST is the relative weight of the impact of threat T, W \_ is the impact of threat T that originates in grid cell c on the habitat in grid cell (function of the distance), and 0 \ is the sensitivity of habitat to threat T.

After computing the degradation level, the habitat quality e is defined as follows:

e = X f1 -Y ' ' + g ^h 9

X ranges from 0 to 1, where 1 indicates LULC classes with the highest suitability for species and g is the half-saturation constant, i.e., 50% of the maximum level of degradation that our territory can support. The degradation (Equation 8) is then translated into habitat quality level according to the half-saturation constant (Equation 9), with UVW' = 0.8; g = 0.4 in this case study. Thus, the higher the degradation index ' and the half-saturation constant g are, the lower the quality of habitats e .

To define the parameters of the model, we used a combined approach based on a literature review [START_REF] Polasky | The Impact of Land-Use Change on Ecosystem Services, Biodiversity and Returns to Landowners: A Case Study in the State of Minnesota[END_REF][START_REF] Terrado | Model development for the assessment of terrestrial and aquatic habitat quality in conservation planning[END_REF][START_REF] Salata | Mapping Habitat Quality in the Lombardy Region, Italy[END_REF][START_REF] Sallustio | Assessing habitat quality in relation to the spatial distribution of protected areas in Italy[END_REF] and participatory meetings with local stakeholders and managers. More specifically, we first defined our scores based on their absolute and, more importantly, their relative values according to the literature. Then, our scores were further discussed and validated by our panel of experts (a public forest manager, a biodiversity officer in the BVRP, and a forest officer of the National Geographic Institute). These experts are specialists in biodiversity protection issues, ecological corridors (e.g., forest continuum, network of wetlands) and forest management. Inputs are presented in Tables 3 and4.

RESULTS AND DISCUSSION

Recreation model

Attractiveness model based on supply factors

We estimated four count data models: a Poisson model, a Zero Inflated Poisson (ZIP) model, a negative binomial model and a zero-inflated negative binomial model. The likelihood ratio test on alpha, representing the dispersion parameter, showed that our dataset was not overdispersed, justifying the use of a Poisson model over a negative binomial model.

Vuong's test confirmed that the ZIP model was preferable to a standard Poisson model (z>2) [START_REF] Vuong | Likelihood ratio tests for model selection and non-nested hypotheses[END_REF].

We tested various model specifications, including the following:

-Hiking paths as the number of paths instead of meters; -Waterways and forest surfaces instead of coverage share variables; -Public versus private forests; -Distance to the nearest forest; and -Different variables explaining non-visited forests. We present the best fitting model according to the Akaike information criterion (AIC) and Bayesian information criterion (BIC) model selection criteria in Table 5. We can observe that the most powerful predictor is elevation, followed by the presence of water and the presence of natural and cultural points of interest at the site. This result is also observed in previous studies based on stated preferences (e.g., [START_REF] Nielsen | An economic valuation of the recreational benefits associated with nature-based forest management practices[END_REF][START_REF] Colson | Map and determinants of woodlands visiting in Wallonia[END_REF][START_REF] Abildtrup | The effect of forest land use on the cost of drinking water supply: A spatial econometric analysis[END_REF][START_REF] Termansen | Modelling and mapping spatial heterogeneity in forest recreation services[END_REF][START_REF] De Valck | Benefits of clearing forest plantations to restore nature? Evidence from a discrete choice experiment in Flanders, Belgium[END_REF][START_REF] Zulian | ESTIMAP: A GIS-based model to map ecosystem services in the European Union[END_REF][START_REF] De Valck | Outdoor recreation in various landscapes: Which site characteristics really matter?[END_REF]. Forest types are integrated as logarithms and can thus be interpreted as elasticities. Mixed forests provide the greatest attraction compared to coniferous and broadleaved forests. These types of forest characteristics are, to the best of our knowledge, rarely integrated in models explaining recreational visits, making comparisons difficult. However, some studies have integrated some of these characteristics in their regressions. For instance, [START_REF] Sen | Economic Assessment of the Recreational Value of Ecosystems: Methodological Development and National and Local Application[END_REF] found a positive relationship between visits and the presence of woodland (they did not separate forest types). [START_REF] Termansen | Modelling and mapping spatial heterogeneity in forest recreation services[END_REF] and [START_REF] Colson | Map and determinants of woodlands visiting in Wallonia[END_REF] found a preference for broadleaved forest, and [START_REF] Schägner | Mapping recreational visits and values of European National Parks by combining statistical modelling and unit value transfer[END_REF] found that the forest type had no effect. The presence of a population at 2 km around the site has no impact on-site attractiveness. This finding invalidates the potential presence of a crowding effect in the area [START_REF] Pröbstl | Management of Recreation and Nature Based Tourism in European Forests[END_REF], that is, a conflict between users and a congestion effect influencing attractiveness. The same result applies for the variable hiking path. This result contradicts our first intuitions and the findings in the literature (e.g., De Valck et al., 2017) but may be because the case study includes many hiking paths so that the variable does not discriminate between different supply units. Furthermore, elevation has a negative effect on attractiveness. Nevertheless, the square of elevation has a positive impact, giving an overall non-linear effect of elevation. That is, elevation decreases attractiveness to a certain point (823 m, in our case study), above which it increases attractiveness at higher elevations. This non-linearity is also observed in [START_REF] Colson | Map and determinants of woodlands visiting in Wallonia[END_REF] for the slope effect 16 . This result makes sense in a mountainous environment, such as the BVRP; indeed, it shows that individuals are more attracted by lowland and highland forests and have lower preferences for forests located in medium-altitude mountains.

The participation function is explained by only hiking path presence. Sites with fewer hiking paths are more likely to be visited, but the predictor is weak.

Attractiveness model based on demand factors

We estimated four count data models. However, for the same reasons, we ultimately relied on a ZIP model to explain attractiveness based on demand factors. The results are presented in Table 6. Other variables have been included to explain participation decision; however, none of them was significant. 16 [START_REF] Colson | Map and determinants of woodlands visiting in Wallonia[END_REF] modeled two slope variables as a proportion of woodland on slopes less than 10° and greater than 30°. They showed that woodlands located on slopes attract more visitors. However, overly steep slopes, that is, slopes greater than 30°, have the opposite effect, as is the case for many of our forests located below 823 m in our case study. As expected, distance traveled (and thus trip cost) negatively affects trip demand, as shown in Figure 4 representing the distance decay effect on trip demand. The figure shows that the negative trip cost effect plays a major role, greatly decreasing the probability of a visit to the first kilometers that individuals must travel. The income variable, positively but weakly, affects trip demand. This result is unusual in travel cost models, as income is commonly not significant in explaining visitors' trips [START_REF] Shresta | Value of recreational fishing in the Brazilian pantanal: a travel cost analysis using count data model[END_REF]Parsons, 2003;[START_REF] Martínez-Espiñeira | Recreation demand analysis under truncation, overdispersion, and endogenous stratification: An application to Gros Morne National Park[END_REF]Garcia andJacob, 2010, Roussel et al., 2016). The availability of substitutes around a potential outset location also influences visit numbers, which is common in recreation models [START_REF] Sen | Economic Assessment of the Recreational Value of Ecosystems: Methodological Development and National and Local Application[END_REF][START_REF] Schägner | Mapping recreational visits and values of European National Parks by combining statistical modelling and unit value transfer[END_REF]. The presence of urban areas, agricultural areas, water bodies and broadleaved forests negatively influences visits at more distant sites, in line with previous studies 17 . More surprisingly, the presence of mixed forests and wetlands around outset locations tends to positively influence 17 We tested for a potential multicollinearity bias between substitute types. It appeared that they are weakly correlated, except for coniferous forest and agricultural lands. However, coniferous forest has been included as a reference in the model, and thus is not included as a variable, reducing strongly the risk of a multicollinearity bias. visits to forests located at more distant sites. This means that almost all land cover types are more attractive than coniferous forest cover (i.e., the reference in the model) except for mixed forest and wetlands.

Combined attractiveness index (CAI)

The CAI is finally computed and mapped in the BVRP, as shown on the left side of Figure 5. Graphically, we can observe that attractiveness is principally concentrated in the center of the BVRP, which is also the most elevated part of the area. A hotspot analysis (right side of Figure 5) statistically confirms this graphical insight. Hotspot areas are characterized by highdensity clusters of the CAI and surrounded by low-density clusters of CAI, referred to as coldspots. Hotspot analysis is performed by applying a combination of statistical analysis and spatial procedures using high/low clustering18 . Thus, the results allow us to presume a high influence of elevation, which is confirmed by the Pearson correlations presented in Table 7. The correlations show that the CAI is primarily determined by the supply index ASI, which, in turn, is highly influenced by elevation and, to a lesser extent, by the ADI even though the correlation is also strong. This means that biophysical attributes are particularly critical to take into account in outdoor recreation models. Moreover, we can observe that demand and supply are significantly correlated but not very strongly. This result confirms the fact that models driven by individuals may have a truncated view of recreation, as the analysis of both sides of this ES flow reveals different and potentially complementary information.

Habitat quality

The InVEST output on the habitat quality index is presented in Figure 6 and is based on two indexes produced: degradation and suitability. Habitat quality can be interpreted as a relative level of habitat suitability to the current landscape. A high score in grid cells means that the habitat quality in the cell is high relative to that in other cells (Tallis et al., 2013). Figure 6 shows spatial heterogeneity in terms of quality across the study area. The BVRP seems separated into two negatively correlated parts. The western part shows the highest quality index of habitats; inversely, the lowest index appears in the eastern part of the site. The two threats that appear to be the most important determinants of habitat quality are trails and agricultural land. Indeed, the blue part of Figure 6 corresponds to the densest part of hiking paths and follows the large agricultural belt located at the eastern part of the BVRP that can also be observed Figure 2. These two factors correspond to the lowest weights in terms of impact and sensitivity (see Tables 3 and4), but their effect on habitat quality is relative to their high concentration in the northeast of the territory. The other threat factors (urban environments and roads) appear to be less important for the relative habitat quality because of the low concentration of these elements in the territory.

Relationship between habitat quality and recreational attractiveness

The relationship between habitat quality and the recreation attractiveness (CAI) indexes can be analyzed with Spearman correlations (Table 8). We find a negative correlation between the attractiveness index and habitat quality, showing statistically that strong recreational attractiveness is opposed to a high quality of habitats. [START_REF] Pröbstl | Management of Recreation and Nature Based Tourism in European Forests[END_REF] highlighted this type of conflict between outdoor recreation versus environmental quality in the context of a multifunctional forest. This result is not surprising because the habitat degradation characterized by InVEST is a function of anthropogenic pressure proximity (such as the presence of roads and thus accessibility) to natural habitat, and conversely, habitat quality is related to habitat preservation (a proxy of low accessibility). However, a synergistic effect can also be expected between environmental quality and recreation attractiveness as has been shown by various papers (e.g., [START_REF] Dehez | Une approche multi-attributs de la demande de loisirs sur les espaces naturels : l'exemple de la forêt publique[END_REF][START_REF] Abildtrup | The effect of forest land use on the cost of drinking water supply: A spatial econometric analysis[END_REF]. Our correlation does not show this relationship, and this may be directly related to the way the habitat quality index is assessed by linking distance to threats and quality. To verify this hypothesis, we measured the correlation between another type of habitat quality index and our CAI. The second index is the level of environmental protection status in every pixel, obtained by overlapping the different layers of protected zones19 . Correlations show a positive association between the protection statuses and recreation attractiveness (CAI), highlighting that a higher level of protection increases the attractiveness and no relationship with the habitat quality index computed by InVEST. The latter result can be explained by the fact that the two indicators do not necessarily represent the same type of biodiversity. The protection statuses are related primarily to remarkable biodiversity, and the habitat quality computed by InVEST includes suitable habitat for both remarkable and ordinary biodiversity.

Using the indexes to evaluate the regional park charter

Regarding BVRP management strategy, the Spearman statistics presented in Table 9 tend to validate the choices made in terms of recreation; however, the results are more ambiguous with regard to habitat quality.

In the High Vosges, the objective is to combine biodiversity conservation with recreational attractiveness. Our results on attractiveness show that this territory is the most attractive among the three territories. However, the correlation with habitat quality is negative; even though habitat quality seems high in the western part of the High Vosges and it is also the territory with the higher level of protection status. The efforts to conserve and protect habitats in the High Vosges seem thus to be insufficient.

Regarding the 1000 Ponds Plateau, the charter objective is to preserve habitat quality and to limit the decline from agricultural and industrial activities. We can observe that the territory effectively benefits from high habitat quality. This area is the sole territory showing a positive correlation with the index. The goal seems relevant here, even though it does not give the impression of being ambitious, as it aims to preserve the strength of the area and not to ameliorate a weakness. Limiting the territory's economic decline could be achieved by developing a recreation policy around agricultural products and by implementing a better communication plan around the habitat richness of the area. Indeed, it has been shown that promoting agritourism enhances and promotes the development of rural areas [START_REF] Greffe | Is rural tourism a lever for economic and social development?[END_REF][START_REF] Barbieri | Perceived Benefits of Agricultural Lands Offering Agritourism[END_REF]. However its success highly depends on the support local policies grant to farmers [START_REF] Greffe | Is rural tourism a lever for economic and social development?[END_REF][START_REF] Sgroi | Agritourism and local development: A methodology for assessing the role of public contributions in the creation of competitive advantage[END_REF].

Finally, the aim of BVRP managers in the Valleys and Piedmont area is to control urbanization and habitat fragmentation. The habitat degradation index, which is also an indicator of habitat fragmentation, demonstrates that the objective does not tend to be reached in this territory (the correlation being positive). Efforts to control urbanization should continue. Therefore, PNRBV managers and all stakeholders have to pursue their current efforts. For the 1000 Ponds Plateau area, this means reviving local attractiveness and developing a greater economic autonomy thanks to the valorization of local resources (e.g. energy, food and forage resources). In the Valleys and Piedmont area, limiting natural habitats fragmentation must remain a priority. The balance between the attractiveness and the rich natural habitats of the High Vosges has to be maintained. This can be achieved by implementing integral biodiversity reserves (prohibiting public access). In this sense, our work helped PNRBV managers to layout quiet areas for the Western Capercaillie (Tetrao urogallus), as part of the Park's new forestry charter.

Limitations

Three principal limitations and refinements of the study are possible to provide better modeling and policy orientations.

First, the recreational attractiveness index developed is valid for the local population, that is, the population living in and just around the park, but it does not include individuals living outside regions covered by the park (Grand-Est and Bourgogne Franche-Comté). Although a large part of the visits are made by local people, 86% of the French people visiting the park come from the studied regions (ORTA, 2011), and a significant part of recreationists, approximately 35%, are tourists coming from other countries (principally Germany and Belgium). Their preferences should also be studied to complete the recreational attractiveness index. This can be done with the same type of model but based on an international online survey or in an on-site survey capturing this type of tourist visits.

Second, the proxy for testing crowding effect characterized by our variable "Population 2 km around the site", could be approached by a seemingly more appropriate proxy such as the distance to the nearest city. This may allow avoiding looking in all directions, as population density can be highly isotropic. However, the territory of our case study comprises many small cities. This implies that distances between each pixel (i.e. forests) and the nearest city would be very similar, resulting in a non-significant coefficient related to this variable. Third, we consider only motorized trips, which is typical in travel cost studies. Although 84% of visitors visit forests by using their car or their motorbike (ORTA, 2011), the remaining part of visitors who use others means of transport should be accounted for to improve our model. Including non-motorized visitors, by considering the opportunity cost of time for these visitors, may exacerbate the bias of spatial sorting effect inherent to the travel cost approach, i.e., people with a strong preference for outdoor recreation tend to live closer to forest areas. To address this issue, the combination with a stated preference approach such as a choice experiment study is possible, as in Abildtrup et al. (2015b).

CONCLUSION

In this paper, we developed an inventive method to predict, in spatially explicit terms, the recreation attractiveness potential combining supply and demand factors to inform decisionmaking. This approach is innovative in the sense that travel cost models are usually a-spatial and specified with only the sociodemographic characteristics of individuals. The method is based on the combination and transfer of a Lancasterian function of biophysical aspects and a travel cost model based on agents' sociodemographic characteristics. Functions are derived from a previous survey (comprising a part of the studied area) and are fed by GIS data. The method developed for recreational service mapping is reproducible for any extrapolation of a single-site or multiple-site travel cost model, allowing a better apprehension of the spatial distribution of the service. The results demonstrate that travel cost models should include the biophysical context of visited ecosystems, as they play a large role in the recreational trip choice. Not considering biophysical characteristics of different destination sites ignores the sole adjustment variable that management planners have for their recreation policy planning; and more importantly, it tends to ignore an important part of the trip decision choice.

We further validate the usefulness of the created recreation model by using it in the evaluation of a regional park charter. Because the park charter is developed around two objectives, i.e., biodiversity conservation and recreational attractiveness, we add a habitat quality index (computed with InVEST) to the evaluation. Recreational attractiveness and habitat quality are then statistically compared with the spatial planning strategy pursued by park managers. This comparison helps to highlight and locate the strengths and weaknesses of the established planning strategy. Moreover, because land planners principally rely on GIS technology for the definition of planning strategy, producing spatial information appears particularly critical for ES to be a real lever for action in day-to-day decision-making. The relevance of this work has already been proven by the practical use of results in the design of two policies in the park (i.e., the territorial forestry charter20 , and the special protection areas designed for Tetrao urogallus).

We conclude with directions that can be taken in future research. Here, we evaluate the state of the art of recreational service and habitat quality in a positive analysis. An interesting extension would be to adopt a normative analysis and thus to optimally allocate different services with spatial optimization modeling (e.g., by using production possibility frontiers). Doing so would allow the emphasis of areas that benefit from comparative advantages in the provision of different ES. The spatial assessment of wood production can enrich the analysis because it will allow the analysis of a bundle of ES in terms of trade-offs and synergies. The functions developed in this paper will serve as inputs in the optimization model. Furthermore, the inclusion of biophysical attributes is critical, as we have observed in this paper. However, what visitors actually perceive in regard to those attributes remains an open question. Indeed, here, we make the implicit assumption that individuals perfectly perceive the biophysical attributes of visited sites using GIS data. Further research should investigate the accuracy of this assumption, for instance, by estimating the model with perceived variables on characteristics of the forest or perceived distance traveled, to compare results with the GIS-based model.

Figure 1 :

 1 Figure 1: Snapshot of the methodological framework

Figure 2 :

 2 Figure 2: The land cover types (a) and the three main territories of the Ballons des Vosges Regional Park (b)

Figure 3 :

 3 Figure 3: Area covered by the initial survey in the BVRP, the forest units visited, and the kilometric mesh used for the benefit transfer

Figure 4 :

 4 Figure 4: Effect of travel cost on travel distance

Figure 5 :

 5 Figure 5: Spatial distribution of the combined attractiveness index (left) and hotspot analysis (right)

Figure 6 :

 6 Figure 6: Habitat quality in the BVRP

Table 1 :

 1 Description of the variables for the attractiveness model based on supply factors

				Descriptive statistics
	Variable	Description	Data source	Mean	Median	Std. err.
	Population at 2 km Population within 2 km buffer	Census data, INSEE	2165	1827	1472
		zone of the mesh
	Share of forest	Share of coniferous forests at	GIS calculation from the	0.41	0.43	0.22
	types	the site	French Forest database,
			BD Forêt®, IGN
		Share of broadleaved forest at		0.17	0.09	0.15
		the site
		Share of mixed forests at the		0.27	0.25	0.14
		site
	Hiking path	Hiking path (in m/km 2 ) including	French Topographic	11369	6998	14781
		biking paths	database, BDTOPO®,
			IGN -paths, roads
	Natural and	Number of natural and cultural	BDTOPO®, IGN -nature	0.3	0	0.86
	cultural points of	points of interest in the mesh	and culture
	interest
	Waterways, water	Share of waterways in the	BDTOPO®, IGN -water	0.09	0.06	0.14
	surfaces	mesh with a buffer zone of 200	surfaces
		m around the courses
	Elevation	Altitude (in m)	French Digital Terrain	669	630	177
			Model (DTM), BD Alti®,
			IGN

Table 2 :

 2 Variables used in the attractiveness model based on demand factors and in the transfer

	Variables	Description	Data source	Descriptive statistics
				Mean Median	Std. err.
	Visits	Visits reported by visitors in the year	Survey	15.63	6	23.5
		preceding the survey
	(Trip cost) TC	Initial survey
		Calculated from the centroid of the visitors' IRIS centroid 13 to the recreational unit	ArcGIS GIS calculation, ESRI	16.20	10.46	23.93
		Reconstituted variable in the transfer
		Calculated from each IRIS centroid to each
		recreational unit centroid with Equation 1
	Income (€)	Initial survey
		Median individual net income in different	Census data INSEE	33007 33750	14581
		income classes
		Reconstituted variable in the transfer
		Median individual net income at the IRIS level
	Availability of	Initial survey
	substitutes
		Share of other land use types around a 5 km 14 buffer around the visitor's IRIS centroid	GIS calculation (ESRI ArcGIS) from Corine Land Cover (2012)
		% Urban areas (CLC1)		0.38	0.25	0.34
		% Agricultural areas (CLC2)		0.25	0.24	0.27
		% Wetlands (CLC4)		0.001 0	0.008
		% Water bodies (CLC5)		0.008 0	0.014
		% Coniferous forests (CLC311)		0.19	0.14	0.18
		% Broadleaved forests (CLC312)		0.05	0	0.15
		% Mixed forests (CLC313)		0.04	0	0.11
		Reconstituted variable in the transfer
		Share of other land use types within a 5 km buffer around the IRIS centroid

Table 3 :

 3 The characteristics of the threats to habitat quality in the BVRP

	Type of threat	Maximum distance	Impact weight
	T Artificial areas	UVW. ' 1.5	ST 0.9
	Main roads	1.5	0.9
	Secondary roads	1	0.7
	Agricultural land	1	0.56
	Trails	0.3	0.35

Table 4 :

 4 Habitat suitability and the sensitivity to each threat in the BVRP

		Habitat suitability	Urban Agriculture	Road1	Road2	Trail
	Forest cover Broadleaved forest	X 0.9	0.8	0.5	0 T 0.8	0.6	0.4
	Coniferous forest	0.8	0.72	0.45	0.72	0.54	0.36
	Mixed forest	0.9	0.8	0.5	0.8	0.6	0.4

Table 5 :

 5 Zero-inflated Poisson model for the attractiveness model based on supply factors of destination sites

	Variables

Table 6 :

 6 Zero-inflated Poisson model for the attractiveness model based on demand factors of outset locations

	Variables	ZIP coefficients Std. Err.
	Trip demand
	Trip cost	-0.0189***	(0.00212)
	Income	4.16e-06**	(1.88e-06)
	Log (% Urban substitutes availability)	-0.0862***	(0.0226)
	Log (% Agricultural substitutes availability)	-0.0480**	(0.0186)
	Log (% Wetland substitutes availability)	0.203***	(0.0775)
	Log (% Water bodies substitutes availability)	-0.0621***	(0.0181)
	Log (% Broadleaved substitutes availability)	-0.140***	(0.0219)
	Log (% Mixed forest substitutes availability	0.0386**	(0.0171)
	Constant	2.267***	(0.109)
	Participation
	Trip cost	-38.23	(11,902)
	Constant	3.989***	(0.583)
	Log likelihood LR chi 2 (9)	-1072.682 244.97
	Prob>chi 2	0.000
	AIC	2167.36
	BIC	2206.36
	Number of zero observations	162
	Number of non-zero observations	94
	Total observations	256	256
	Standard errors in parentheses *** p<0.01, ** p<0.05, * p<0.1 Vuong's test of ZIP vs. standard Poisson: z=2.46, Pr>z=0.0069 Alpha test on ZINB=1.082131, Pr>z=0.580 Coniferous forests are set as the base case coverage share

Table 7 :

 7 Pearson correlation of the CAI, ASI and ADI

		ASI	ADI	CAI
	Supply Index -ASI	1.0000
	Demand Index -ADI	0.3926***	1.0000
	Combined Index -CAI	0.9538***	0.6224***	1.0000

Table 8 :

 8 Spearman correlations between recreational attractiveness (CAI), habitat quality index and level of protection status

		Recreation	Biodiversity indicator
		CAI	Protection status	Habitat Quality
	CAI	1.0000
	Protection status	0.1779***	1.0000
	Habitat quality	-0.2471***	n.s	1.0000
	***p<0.01 n.s: not significant

Table 9 :

 9 Spearman correlations between recreation and habitat quality indexes and territories

		Recreation		Habitat quality
		CAI	Suitability	Degradation
	High Vosges	0.6373***	-0.2393***	0.2114***
	1000 Ponds Plateau	-0.1606***	0.5203***	-0.5147***
	Valleys and Piedmont	-0.5534***	-0.1106***	0.1358***
	***p<0.01 n.s: not significant

The paper by[START_REF] Sen | Economic Assessment of the Recreational Value of Ecosystems: Methodological Development and National and Local Application[END_REF] is mostly based on various seminal papers developing the use of GIS in travel cost techniques:[START_REF] Jones | Estimating Arrival Numbers for Informal Recreation: A Geographical Approach and Case Study of British Woodlands[END_REF],[START_REF] Brainard | Using isochrone surfaces in travel-cost models[END_REF][START_REF] Brainard | Integrating geographical information systems into travel cost analysis and benefit transfer[END_REF] and[START_REF] Lovett | Improving Benefit Transfer Demand Functions: A GIS Approach[END_REF].

Integrated valuation of ecosystem services and trade-offs from the Natural Capital Project. InVEST is an open-source software available at: https://www.naturalcapitalproject.org/invest/ (webpage access date 02/14/19)

More information on the park charter is available online (webpage access 02/11/14): https://issuu.com/parcdesballons/docs/charte_2012-2024_fevrier2014

French environmental code (Article L333-1) modified after the Law on Biodiversity 2016 (L2016-1087 8 août 2016, Article 48).

http://www.europarc.org/library/europarc-events-and-programmes/european-charter-for-sustainabletourism/ (webpage access date 01/25/19)

We considered this buffer zone because forest units are sometimes juxtaposed with BVRP territory or are very close to the limits of the BVRP.

D, the distance between the outset location and the forest, is calculated by using the osrmtime command from Stata[START_REF] Huber | Calculate travel time and distance with OpenStreetMap data using the Open Source Routing Machine (OSRM)[END_REF].

https://www.statistiques.developpement-durable.gouv.fr/donnees-sur-les-immatriculations-des-vehicules (webpage access date 12/03/18)

K takes into account vehicle depreciation, maintenance, fuel and insurance costs and is published annually by the fiscal administration in France (webpage access date

12/01/16) https://www3.impots.gouv.fr/simulateur/calcul_impot/2017/pdf/baremekm.pdf This inclusion of vehicle depreciation is typical in the TCM (seeParsons, 2003). However, this inclusion has been discussed in the literature (see[START_REF] Earnhart | Do travel cost models value transportation properly?[END_REF] because individuals may not explicitly perceive these types of costs.12 For a discussion regarding the OCT, refer to[START_REF] Roussel | Recreation demand analysis of sensitive natural areas from an onsite survey[END_REF].

This is done with the command Getis-Ord General G from the ESRI ArcGIS© software package of spatial statistics.

The protection index is based on French biodiversity protection layers ranging from 0 to 8, including Natura

areas, Nature Reserves, Biological Reserves, Biosphere Reserves, Regional Nature Reserves, Biotope Protection or Natural Zone of Interest for Ecology Flora and Fauna (ZNIEFF types I and II).

The territorial forestry charter brings together all the actors of a territory to define a program of actions to enhance their forest areas. It takes into account all the forest uses in a territory: economic, environmental and social.

ACKNOWLEDGMENTS

This research has primarily been funded by the AFFORBALL project (PSDR4). This work has also been supported by grants from the French National Research Agency (ANR) as part of the "Investissements d'Avenir" program (ANR-11-LABX-0002-01, Lab of Excellence ARBRE). The authors wish to thank the following for research assistance and constructive discussions regarding the issues raised in this paper: Jens Abildtrup, Antoine Fargette, Emeline Hily, Claude Michel (BVRP director), the different members of the CIRED and BETA seminars, and reviewers/participants at the LEO conference (April, 2018), 6 th WCERE (June 2018), the Ulvön-WONV conference (June 2018), the FAERE and ESP-EU conference (October 2018) for their useful remarks. We particularly thank Jens Abildtrup for sharing the database. We also thank the editor Irene Ring and two anonymous reviewers for their valuable comments that significantly improved the manuscript. The usual disclaimer applies, and any shortcomings in the article are our responsibility.