Assessing the Impact of Marine Tourism and Protection on Cultural Ecosystem Services Using Integrated Approach: A Case Study of Gili Matra Islands

Cultural ecosystem services (CES) are intangible benefits people obtain from an ecosystem through physical and cognitive interactions. Understanding CES provides vital insights into how activities impacting ecosystem services also impact people. Gili Matra Islands, a set of three small tropical islands located in West Nusa Tenggara Province, Indonesia, are an increasingly busy marine tourism destination and a marine protected area. By integrating a hedonic monetary value model with a eudaemonic non-monetary value model, this study examines the impacts of tourism and marine protected area management on cultural ecosystem services in the Gili Matra Islands. Results showed that the distance had significantly influenced property prices to coastlines, beach spots and coastlines with sunset views. In addition, the property prices of each individual island showed significant correlations with particular marine tourism and protection features. Less restricted marine protected zones and coastlines were the most significantly influencing variables to the strong eudaemonic well-being dimensions expressed by residents. The Spiritual dimension produced the highest score and was most significantly affected by several features. This study utilised higher accuracy of properties and residents’ location, enabling more accurate assessments of interaction between CES and the features. This study also discusses how these novel insights in the small island’s CES case can inform vulnerability assessments, reviews of recreation taxes, and spatial planning for marine protected areas and help optimise beach nourishments.


Introduction
Ecosystem services are defined as "benefits people obtain from ecosystems" [1]. Marine and coastal environments provide a wide range of ecosystem services, including fisheries, food, tourism, spiritual opportunity, genetic material, coastal protection, biodiversity, fish habitat, climate regulation, and so forth [2,3]. Marine and coastal environments also provide cultural ecosystem services (CES) that are largely unexplored. CES arises through physical and cognitive interactions between the natural environment and people and includes intangible benefits people obtain, such as spiritual, aesthetic, and recreational experiences [4,5]. CES is more challenging to assess than other ecosystem services, which might make CES the least explored [4]. However, CES research is emerging, with CES becoming more recognised for providing benefits to people [6], particularly when more complex conflicts occur when 1.
Use a hedonic approach to assess the monetary value of CES influenced by marine tourism and marine protection in the Gili Matra Islands.

2.
Use a eudaemonic approach to assess well-being values associated with CES in marine and coastal environments in Gili Matra Islands. 3.
Investigate the interrelations between CES-associated well-being value to marine tourism and marine protection in Gili Matra Islands.

Study Area: Gili Matra Islands
Gili Matra Islands (Figure 1), West Nusa Tenggara, Indonesia, were chosen as the case study to represent a marine protected area and marine tourism destination [23]. The name "Gili Matra" is an abbreviation of Gili Meno, Gili Trawangan, and Gili Air-three islands that provide a home for many marine tourism sites, such as coral reefs, scenic beaches, etc. To help sustain the marine and coastal resources, the government of the Republic of Indonesia designated a marine protected area (MPA) for the islands-a geographical space defined and managed through legal or other effective means to achieve the long-term conservation of nature and associated ecosystem services and cultural values [24]. Its primary function as a Marine Tourism Park is to promote sustainable naturebased tourism while conserving the natural environment so the society can benefit from the area [25][26][27].  Since tourism emerged in 1986 [28], Gili Matra Islands have experienced increasing tourism activities, predominantly landscape changes for residential and tourism accommodation [23,29,30]. Gili Matra Islands have developed unique and diverse CES like other coastal tropical islands that have undergone increasing tourism activities [10]. The number Since tourism emerged in 1986 [28], Gili Matra Islands have experienced increasing tourism activities, predominantly landscape changes for residential and tourism accommodation [23,29,30]. Gili Matra Islands have developed unique and diverse CES like other coastal tropical islands that have undergone increasing tourism activities [10]. The number of tourists and residents reportedly increased, as shown in Figure 2 [23,30,31], except during the COVID-19 pandemic. Unfortunately, these increases imposed direct (e.g., physical damage) and indirect stressors (e.g., pollution waste, as shown in Figure 3) that degraded ecosystem and seawater quality [23,29,30,[32][33][34]. This problem becomes even more critical given the increasing anthropogenic value of tourism interacting with CES. That is, negative impacts may degrade the ecosystem and the CES, such as marine recreational, aesthetic value, spiritual value and other cultural uses [4].

Questionnaire
A questionnaire was prepared to collect information about residential property prices, locations, and well-being values associated with CES. The questionnaire aimed to extract hedonic and eudaemonic primary data from the residents of Gili Matra Island through 4 general questions, 6 hedonic well-being open-questions, and 23 eudaemonic a 7-point Likert scale questions. The survey managed to collect 131 responses for hedonic Since tourism emerged in 1986 [28], Gili Matra Islands have experienced increasing tourism activities, predominantly landscape changes for residential and tourism accommodation [23,29,30]. Gili Matra Islands have developed unique and diverse CES like other coastal tropical islands that have undergone increasing tourism activities [10]. The number of tourists and residents reportedly increased, as shown in Figure 2 [23,30,31], except during the COVID-19 pandemic. Unfortunately, these increases imposed direct (e.g., physical damage) and indirect stressors (e.g., pollution waste, as shown in Figure 3) that degraded ecosystem and seawater quality [23,29,30,[32][33][34]. This problem becomes even more critical given the increasing anthropogenic value of tourism interacting with CES. That is, negative impacts may degrade the ecosystem and the CES, such as marine recreational, aesthetic value, spiritual value and other cultural uses [4].

Questionnaire
A questionnaire was prepared to collect information about residential property prices, locations, and well-being values associated with CES. The questionnaire aimed to extract hedonic and eudaemonic primary data from the residents of Gili Matra Island through 4 general questions, 6 hedonic well-being open-questions, and 23 eudaemonic a 7-point Likert scale questions. The survey managed to collect 131 responses for hedonic

Questionnaire
A questionnaire was prepared to collect information about residential property prices, locations, and well-being values associated with CES. The questionnaire aimed to extract hedonic and eudaemonic primary data from the residents of Gili Matra Island through 4 general questions, 6 hedonic well-being open-questions, and 23 eudaemonic a 7-point Likert scale questions. The survey managed to collect 131 responses for hedonic data and 132 responses for the eudaemonic instrument. The questionnaire can be found in Supplementary Materials Questionnaire S1. The hedonic pricing model examined the impact of the features of interest on people's selections, which were property prices representing the monetary value of CES [4]. The features of interest included marine tourism sites (i.e., coastline, sunset view coastline, beach spots, and dive spots), marine protected area (MPA) core zone (strict zone) and another zone (less strict), public facilities, and road networks [29,30]. To determine the impact of these features of interest, spatial-environmental attributes such as distance to and visibility were considered [4,13]. They were employed through standard multiple regression statistical analysis using SPSS Statistic 26 [39], as described in Equation (1).
where P is property price, representing the CES's monetary values; α is Intercept (constant term); β is the slope coefficient; E is attributes of the features of interest. Unlike previous research [4], this study used property price data from both residential and tourism accommodations (e.g., hotels and cottages). Consequently, separate models were employed for different property types and other variations (i.e., set of islands and The hedonic pricing model examined the impact of the features of interest on people's selections, which were property prices representing the monetary value of CES [4]. The features of interest included marine tourism sites (i.e., coastline, sunset view coastline, beach spots, and dive spots), marine protected area (MPA) core zone (strict zone) and another zone (less strict), public facilities, and road networks [29,30]. To determine the impact of these features of interest, spatial-environmental attributes such as distance to and visibility were considered [4,13]. They were employed through standard multiple regression statistical analysis using SPSS Statistic 26 [39], as described in Equation (1).
where P is property price, representing the CES's monetary values; α is Intercept (constant term); β is the slope coefficient; E is attributes of the features of interest. Unlike previous research [4], this study used property price data from both residential and tourism accommodations (e.g., hotels and cottages). Consequently, separate models were employed for different property types and other variations (i.e., set of islands and individual islands) that provided further insights. Property price and location data for the residential properties were obtained through the questionnaires along with well-being questionnaires (with assistance from a GPS device for high accuracy of location). For primary data/residential property prices, a total of 131 observations were collected. Subsequently, tourism accommodation property prices and addresses were compiled from online markets. Keywords included 'Gili Matra', 'Gili Trawangan', 'Gili Air', 'Gili Meno', and so forth. Data from 46 tourism accommodation properties were collected. Prior to statistical analysis, the attributes of each feature of interest, i.e., distance and visibility, were derived through spatial analysis using QGIS 3.18.2 (QGIS Association) and ArcGIS Pro (ESRI, Redland, CA, USA), which required location data and Digital Elevation Model (DEM). The DEM was collected from Indonesia Geospatial Portal [37] and the United States Geological Survey database [38]. The inverse distance weighted (IDW) method was then used to interpolate property price data and visualise data distribution.

Eudaemonic Approach to Assess Well-Being Values Associated with CES in Marine and Coastal Environments in Gili Matra Islands
To investigate the value of marine tourism and protection features to marine users, the well-being value questionnaire was based on an instrument previously developed by the UK National Ecosystem Assessment (NEA) [16] and other studies [40][41][42]. The instrument incorporated seven dimensions of well-being influenced by CES, including Place identity, Therapeutic value, Engagement and interaction with nature, Social bonding, Memory/transformative value, Spiritual value, and Challenge and Skill. Each dimension was associated with several indicators investigated by exploratory factor analysis [4,16]. Using the instrument, respondents assigned a value to each indicator on a 7-point Likert scale regarding their agreement with natural places nearby, reflected from 1 (strongly disagree) to 7 (strongly agree). Target respondents were the residents of the Gili Matra Islands, and a total of 132 responses were collected.
Statistical analysis of the collected data was performed through confirmatory factor analysis (CFA) using SPSS Amos 25 and SPSS Statistic 26 (IBM, Armonk, NY, USA). CFA was chosen because well-being indicators and dimensions were based on previous literature rather than through own exploration [4]. This analysis investigated relationships between these indicators and well-being dimensions associated with marine and coastal CES. The model was evaluated through several fitness tests, i.e., Chi-Square statistic and degree of freedom (CMIN/df), comparative fit index (CFI), Tucker-Lewis index (TLI), root mean squared error of approximation (RMSEA), and standardised root means squared residual (SRMR). The indicators are provided in Appendix A, Table A1. The questionnaire template for the eudaemonic and hedonic model is available in Supplementary Materials Questionnaire S1.

Interrelations between CES Associated Well-Being Value to Marine Tourism and Marine Protection in Gili Matra Islands
Correlations between CES-associated well-being scores and features of interest were performed through standard multiple regression modelling using SPSS Statistic 26 [39]. Correlations determined if there were influences from the features of interest on respondent valuation over their CES-related well-being related. The dependent variable was the factor scores for each dimension of CES well-being, and the independent variables were the distance to and visibility of features of interest. Since seven dimensions of CES well-being were adapted [4], the models were run seven times, corresponding to each factor.

Influence of Marine Tourism and Marine Protection on the Monetary Value Associated with CES
The first model (model HR1) of residential property prices used price data for the three islands as a dependent variable. Table 1 shows the statistical summary. The minimum price of residential properties was found at 75 million Indonesian Rupiah (IDR) per are (1 are equals to 100 square meters (100 m 2 )) (100 m 2 ), while the highest price was 1000 million IDR per are. The closest feature was the road with almost no distance (only 0.97 m), while the furthest feature was the core zone (MPA's strict zone) with 1620.35 m. The roads themselves are always visible to every residential property. An ANOVA/F-test showed that the model reached statistical significance (Sig. = 0.000), and the R 2 value was considerably high (0.641), implying that the independent variables collectively explained 64.1% of the variance in the dependent variable [39]. Three variables had a significant individual contribution, i.e., distance to the coastline, distance to beach spots, and visibility of the coastline with a sunset view. An increase of one meter in the distance to the coastline decreased the price by 0.697 IDR per are (100 m 2 ). Similarly, an increase of one meter to a beach spot lowered the price by 0.128 million IDR per are. Meanwhile, visibility of at least one coastline facing the sunset increased price by 143.95 million IDR per are. The regression equation and t-test for Model HR1 are shown in Supplementary Materials (Equation (S1) and Table S1).  Figure 5 reveals that high prices were generally located near the coastline. Spatial interpolation also indicated that lower values or price within the first quartile (<25%) resides in the middle of the island, far from the beach spot, coastline, or other marine tourism sites, as shown in Table 2. On the other hand, higher property prices, i.e., prices in the third (50-75%) and fourth (75-100%) quartile, were mostly located near the coastline. Each island also displayed a unique distribution pattern for properties with higher prices.          IDR 1000 million per are.
The second model, Model HR2, used residential property prices for Gili Trawangan Island. ANOVA/F-test showed that the model was statistically significant (Sig. = 0.000) [39], with the R 2 value was 74.4%. Two independent variables had a significant individual contribution, i.e., distance to coastline and distance to facilities. An increase of one meter in the distance to the coastline decreased property prices by 1.011 million IDR per area. Furthermore, a one-meter increase in the distance to public facilities decreased residential property prices by 0.464 million IDR per are. Spatial distribution is presented in Figure 6. The regression equation and t-test for Model HR2 are provided in Supplementary Materials (Equation (S2) and Table S2). IDR 1000 million per are.
The second model, Model HR2, used residential property prices for Gili Trawangan Island. ANOVA/F-test showed that the model was statistically significant (Sig. = 0.000) [39], with the R 2 value was 74.4%. Two independent variables had a significant individual contribution, i.e., distance to coastline and distance to facilities. An increase of one meter in the distance to the coastline decreased property prices by 1.011 million IDR per area. Furthermore, a one-meter increase in the distance to public facilities decreased residential property prices by 0.464 million IDR per are. Spatial distribution is presented in Figure 6. IDR 1000 million per are.
The second model, Model HR2, used residential property prices for Gili Trawangan Island. ANOVA/F-test showed that the model was statistically significant (Sig. = 0.000) [39], with the R 2 value was 74.4%. Two independent variables had a significant individual contribution, i.e., distance to coastline and distance to facilities. An increase of one meter in the distance to the coastline decreased property prices by 1.011 million IDR per area. Furthermore, a one-meter increase in the distance to public facilities decreased residential property prices by 0.464 million IDR per are. Spatial distribution is presented in Figure 6. IDR 1000 million per are.
The second model, Model HR2, used residential property prices for Gili Trawangan Island. ANOVA/F-test showed that the model was statistically significant (Sig. = 0.000) [39], with the R 2 value was 74.4%. Two independent variables had a significant individual contribution, i.e., distance to coastline and distance to facilities. An increase of one meter in the distance to the coastline decreased property prices by 1.011 million IDR per area. Furthermore, a one-meter increase in the distance to public facilities decreased residential property prices by 0.464 million IDR per are. Spatial distribution is presented in Figure 6. The regression equation and t-test for Model HR2 are provided in Supplementary Materials (Equation (S2) and Table S2).
The second model, Model HR2, used residential property prices for Gili Trawangan Island. ANOVA/F-test showed that the model was statistically significant (Sig. = 0.000) [39], with the R 2 value was 74.4%. Two independent variables had a significant individual contribution, i.e., distance to coastline and distance to facilities. An increase of one meter in the distance to the coastline decreased property prices by 1.011 million IDR per area. Furthermore, a one-meter increase in the distance to public facilities decreased residential property prices by 0.464 million IDR per are. Spatial distribution is presented in Figure 6. The regression equation and t-test for Model HR2 are provided in Supplementary Materials (Equation (S2) and Table S2).
The second model, Model HR2, used residential property prices for Gili Trawangan Island. ANOVA/F-test showed that the model was statistically significant (Sig. = 0.000) [39], with the R 2 value was 74.4%. Two independent variables had a significant individual contribution, i.e., distance to coastline and distance to facilities. An increase of one meter in the distance to the coastline decreased property prices by 1.011 million IDR per area. Furthermore, a one-meter increase in the distance to public facilities decreased residential property prices by 0.464 million IDR per are. Spatial distribution is presented in Figure 6. The regression equation and t-test for Model HR2 are provided in Supplementary Materials (Equation (S2) and Table S2).  The third model, Model HR3, analysed residential property prices for Gili Air Island. An ANOVA/F-test showed that the model was statistically significant, and the revealed R 2 value was 0.729. Four independent variables had significant individual contributions (t-test with Sig. < 0.05, [39]), i.e., distance to beach spots, distance to dive spots, distance to core zone, and distance to the road. An increase of one meter in the distance to beach spots lowered residential property prices by 0.467 million IDR per are, as did the distance to dive spots decreased prices by 1.003 million IDR per are for every one-meter increase. Increased distance to the core zone of MPA also lowered prices by 0.454 million IDR per are for every one-meter increase. The most substantial impact on property prices was observed in the distance to the road, which reduced prices by 2.004 million IDR per are for every one-meter increase. The spatial distribution of these findings is presented in Figure 7. The fourth model, Model HR4, analysed residential property prices for Gili Meno Island. An ANOVA/F-test showed that the model was statistically significant, and the R 2 value was 0.876. Two independent variables had significant individual contributions. Increased distance to the coastline decreased prices by 0.271 million IDR per i for every one-meter increase. Similarly, increased distance to the MPA core zone decreased prices by 0.140 million IDR per are for every one-meter increase. The spatial distribution of these findings is presented in Figure 8. The regression equation and full t-test for Model HR3 are available in Supplementary Materials in Equation (S3) and Table S3, whereas for Model HR4 are in Equation (S4) and Table S4.
The final model in the hedonic analyses, Model HT1, analysed property prices of tourism accommodation derived from secondary data. Model HT1 included 46 properties across all three islands. Unlike the residential property, tourism accommodation property prices were expressed per whole property rather than per are (100 m 2 ). Table 3 shows a statistical summary for Model HT1. An ANOVA/F-test showed that the model reached statistical significance (Sig. = 0.000), revealing that all independent variables collectively significantly influenced the prices. The R 2 value of Model HT1 was 0.826. However, only one independent variable, bedroom number, had a significant individual contribution (t-test Sig. < 0.05) to predicting accommodation property prices. Subsequently, due to marine tourism and protection features not significantly contributing to property prices in Model HT1, the regression model was repeated with gradual eliminations of the independent variables. Five omitted variables were selected in order from least significant to most significant; visibility of the road, distance to beach spots, visibility of beach spots, distance to dive spots, and visibility of dive spots. The new regression produced four significant independent variables. Distance to the coastline facing the sunset decreased tourism accommodation property prices by 43.686 million IDR for every one-meter increase. The visibility of at least one core MPA zone increased prices by 29,446.45 million IDR, while the visibility of at least one other MPA zone increased prices by 30,042.83 million IDR. Lastly, the addition of one-bedroom increased prices by 2445.998 million IDR. Figure 9 shows a spatial interpolation of Model HT1. The regression equation and full t-test for Model HT1 are provided in Supplementary Materials (Equation (S5) and Table S5).
The third model, Model HR3, analysed residential property prices for Gili Air Island. An ANOVA/F-test showed that the model was statistically significant, and the revealed R 2 value was 0.729. Four independent variables had significant individual contributions (t-test with Sig. < 0.05, [39]), i.e., distance to beach spots, distance to dive spots, distance to core zone, and distance to the road. An increase of one meter in the distance to beach spots lowered residential property prices by 0.467 million IDR per are, as did the distance to dive spots decreased prices by 1.003 million IDR per are for every one-meter increase. Increased distance to the core zone of MPA also lowered prices by 0.454 million IDR per are for every one-meter increase. The most substantial impact on property prices was observed in the distance to the road, which reduced prices by 2.004 million IDR per are for every one-meter increase. The spatial distribution of these findings is presented in Figure 7. The fourth model, Model HR4, analysed residential property prices for Gili Meno Island. An ANOVA/F-test showed that the model was statistically significant, and the R 2 value was 0.876. Two independent variables had significant individual contributions. Increased distance to the coastline decreased prices by 0.271 million IDR per i for every one-meter increase. Similarly, increased distance to the MPA core zone decreased prices by 0.140 million IDR per are for every one-meter increase. The spatial distribution of these findings is presented in Figure 8.   The final model in the hedonic analyses, Model HT1, analysed property prices of tourism accommodation derived from secondary data. Model HT1 included 46 properties across all three islands. Unlike the residential property, tourism accommodation property prices were expressed per whole property rather than per are (100 m 2 ). Table 3 shows a

Eudaemonic Well-Being Values Associated with CES in Marine and Coastal Environments
Kaiser-Meyer-Olkin (KMO) and Bartlett's test results showed that the KMO Measure of Sampling Adequacy value was 0.864 (threshold of >0.6) and Bartlett's Test of Sphericity Sig. value was 0.000 (threshold of <0.05), which meant that the factor analysis model (Model E1) was appropriate and statistically significant [39]. Next, the model fit measure and validity tool test by the SPSS Amos plugin [43,44] showed that the model could be improved by removing several indicators, including Eng3, Eng5, Plid 8, and Soc17 (see Appendix A Table A1). Following the removal of these indicators in Model E1, communalities values were robust (>0.3), ranged from 0.446 to 0.852 (see Supplementary Materials Table S6), showed that these indicators fitted well with other indicators and significantly explained total variances in Model E1 [39]. Several goodness of fit tests were used and are presented in Supplementary Materials Table S7. Three tests showed good results in standardised root means squared residual (SRMR), Chi-square (CMIN), and Degree of Freedom (df). The SRMR presented a value of 0.074, which was less than 0.08, denoting that the model was a relatively good fit between the observed data and the hypothesised model [45]. The CMIN/df value was 2.155, less than 3, indicating an acceptable [46] and reasonable fit (<5) [47] between the sample data and the model.
Indicator loadings values were mostly above 0.6, implying a valid indicator and its solid association with the respective well-being factor [4], as shown by the values on arrows from each dimension to corresponding indicators in the structure model in Figure 10. For instance, the loading factor from the Spiritual dimension to indicator "Spir14" was the highest, indicating that this indicator most represents the well-being dimension expressed by the local residents. The CES well-being factor scores for each response were then calculated using the factor score weight produced by the CFA analysis (see Supplementary Materials Table S8). For instance, the Skill factor score for response number One was calculated by the equation in Supplementary Materials Equation (S6). After calculating each factor and response, the average and standard deviation of the factor score were calculated for overall responses, as shown in Figure 11.  The average CES well-being dimension scores were high. The lowest factor score was 3.92, still above the midpoint (3.5). These results suggested that the residents of Gili Matra Island experienced positive impacts from CES associated with marine and coastal environments. Among all dimensions, Spiritual value received the highest score with a value of 6.51, inferring that the residents benefited most from CES through their spiritual aspect. The next highest scores were Social bond, Place identity, Memory value, Therapeutic value, Engagement and interaction with nature, and Skill and achievement. Notably, each island also showed the same patterns as shown in Figure 12.   The average CES well-being dimension scores were high. The lowest factor sco 3.92, still above the midpoint (3.5). These results suggested that the residents of Gili Island experienced positive impacts from CES associated with marine and coasta ronments. Among all dimensions, Spiritual value received the highest score with a of 6.51, inferring that the residents benefited most from CES through their spiritual a The next highest scores were Social bond, Place identity, Memory value, Thera value, Engagement and interaction with nature, and Skill and achievement. Notably island also showed the same patterns as shown in Figure 12.

Interrelations between CES Associated Well-Being and Marine Tourism and Protectio
Subsequent analyses investigated the correlation between well-being values an viously utilised marine tourism and protection features in the hedonic monetary analyses. The regression ran seven times, one for each of the CES well-being dime used. The first model of this analysis, Model R1, used the first dimension, Skill and lenge value, as the dependent variable from 132 responses. The result of the ANO test suggested that the model was not statistically significant in collectively affecti factor score of Skill values. The independent variables also could only explain 16.4 was 0.164) of the variances in the factor score. However, two independent variabl have significant individual contributions to Skill values-distance to coastline an tance to other zones. One-meter increases in the distance to the coastline increas Skill values score by 0.002, while one-meter increases in the distance to other MPA decreased the Skill values score by 0.002.
Correlations between all well-being dimensions and marine tourism and prot

Interrelations between CES Associated Well-Being and Marine Tourism and Protection
Subsequent analyses investigated the correlation between well-being values and previously utilised marine tourism and protection features in the hedonic monetary value analyses. The regression ran seven times, one for each of the CES well-being dimensions used. The first model of this analysis, Model R1, used the first dimension, Skill and challenge value, as the dependent variable from 132 responses. The result of the ANOVA/F-test suggested that the model was not statistically significant in collectively affecting the factor score of Skill values. The independent variables also could only explain 16.4% (R 2 was 0.164) of the variances in the factor score. However, two independent variables did have significant individual contributions to Skill values-distance to coastline and distance to other zones.
One-meter increases in the distance to the coastline increased the Skill values score by 0.002, while one-meter increases in the distance to other MPA zones decreased the Skill values score by 0.002.
Correlations between all well-being dimensions and marine tourism and protection features were similar. Table 4 summarises the impact of each independent variable on each well-being factor/dimension, along with their R 2 values. Even though the models did not collectively explain the influences on the factor scores, there were consistent significant individual contributions from several independent variables to the respective well-being dimension. That is, there were consistent significant individual contributions by distance to other zones and distance to the coastline. As the distance to other zones (i.e., the less strict zones) increased, the score for every CES well-being factor decreased. Conversely, the further distance to the coastline, the higher the factor score. For the core zone (no-take zone), the influence only emerged for Spiritual value, where the score dropped as the distance increased. Meanwhile, minimal impacts (i.e., a coefficient value of zero or almost zero) were shown by several independent variables, including distance to the coastline facing the sunset, distance to the dive spot, distance to the core zone, and distance to public amenities.

Influence of Marine Tourism and Marine Protection on the Monetary Value Associated with CES
The first model, Model HR1, indicated that the independent variables used were mostly comprehensive and covered various marine tourism and protection features. Additionally, the infrastructures that support tourism activities in the area include a road network, road accessibility, and amenities [30,[48][49][50][51]. Housing prices can indeed represent the benefits gained from close proximity to coastlines, like beach aesthetics, therapeutic, and others [4,52]. That is, house owners are willing to pay (WTP) a higher price to enjoy the CES provided [52][53][54]. This observation is supported by Table 2, where the highest prices were found to be located adjacent to the coastline. The strong impact of beach spots and coastline on residential property prices should also be a warning about depleting the coastal environment by mass utilisation. Following a WTP appraisal, mechanisms like zoning plans or recreation fees could be improved to hinder the depletion [53]. Visibility of the coastline facing the sunset also significantly influenced residential property prices. The sunset view might enhance property prices [55] as buyers are willing to pay more for the premium views [56,57]. However, other visibility features did not significantly influence. This may be due to the flat terrain of Gili Matra Islands [58,59] reducing visibility. Subsequently, the model for each island displayed particular characteristics where independent variables (i.e., features) significantly affected residential property prices. In Model HR2, Gili Trawangan, amenities were distributed mainly on the island's east coast (Figure 7), including the tourism commercial centre, school, etc. [30]. This finding supported the hypothesis and previous research that considered amenities to support coastal tourism infrastructure [48,49,51]. These findings provide evidence of 'hot spots' for economic activities and settlement. For instance, based on the high property prices dominating the eastern coast of Gili Trawangan, a cost-benefit analysis may recommend beach nourishment to prevent erosion in premium locations [53,60].
Model HR3 suggested that there were four significant variables in Gili Air Island: distance to dive spots, distance to road, distance to beach spots, and distance to the MPA core zone. Eleven dive points on Gili Air were spread evenly around the island, as shown in Figure 8 [29]. This is likely due to local awareness of the recreational activities, biodiversity, and other values that dive spots offer. The distribution of high residential property prices and their correlation with nearby dive spots (coral reefs) could be used as a novel input for updating vulnerability assessments. The latest vulnerability assessment on the island used parameters such as coastline changes, live coral changes, coral reef changes, and developed areas [29]. An accessible road network is a major determinant of nearby property prices [61][62][63]. Roads are considered crucial infrastructure for supporting tourism development [49,51], particularly in Gili Matra Islands, where gas-fueled vehicles are prohibited, making accessible roads essential for pedestrians and cyclists.
Model HR4 showed that distance to the coastline and core zone significantly influenced property prices in Gili Meno Island (Figure 9). The core zone is a "no-take" zone in MPA with strict protection and very minimum human activity [64][65][66], thus, is more likely to have more biodiversity [66]. Although these benefits may not be visible and could not be directly enjoyed, the residents in Gili Matra were most likely aware of these protected areas and related conservation activities [29,67], which may convince them to pay more for properties closer to marine conservation sites [4,68]. The impact of core zones observed in this study may prompt updating of MPA zoning.
According to Model HT1, it is reasonable that hotel or tourism accommodation properties that provide tourists with water views would receive higher prices [69,70]. Figure 10 shows how higher prices were generally situated closer to the coastline in Gili Trawangan. However, this pattern was not obvious in the other two islands, suggesting other variables were dominant. Lastly, the influence of bedroom number on property price was strongly evident as it is one of the primary features that any tourist seeking accommodation considers. The results from hedonic models used in this study showcased various patterns across different islands and property types. These multiple models managed to reveal the impact of distinct marine tourism and protection features on people's valuation of their properties. Each island displayed its own price pattern and what unique features influenced this pattern, which subsequently became highly crucial for island and water management.

Eudaemonic Well-Being Values Associated with CES in Marine and Coastal Environment
For all six CES well-being dimensions, the high score displayed by the results indicated that residents of Gili Matra Island have been experiencing positive impacts from CES. The residents of Gili Matra Islands experienced the most benefit from Spiritual value ( Figure 11). This result varied from previous research by Spanou et al. [4] where the highest score was Therapeutic value. The respondents in Spanou's research were local residents in Argyll and Bute, a region on Scotland's scenic west coast [4,71]. The picturesque coastal area most likely contributed to the therapeutic value experienced by the locals as they can benefit every day. Meanwhile, for locals in Gili Matra, it appeared that the scenery of the islands impacted their well-being mainly on a spiritual level. Spirituality is strong among the people of Gili Matra Islands and Lombok Island, who are well known for practising their religion and living harmoniously [72][73][74]. A coastal environment may enhance spirituality by promoting outdoor activities that enable people to enjoy a scenic natural environment, away from the hustle of cities [75].
The social bond is one of the most well-preserved values in most Indonesian communities. Furthermore, Gili Matra Islands cover only 665 hectares [76], possibly promoting closer and more frequent social interactions between residents. The residents were also found to have community social capital and resilience through iterative and communal activities [77]. This solid social network grew from groups in the islands and even affected governance institutions that consider social and ecological states [78]. Subsequently, the locals were found to highly value the economic aspect of the coastal environment and its surrounding MPA [67]. In addition, a number of marine tourism events have been held, providing avenues for locals and tourists to exchange with each other. These activities have created memories, inspirations, and experiences, both for locals and the tourists in Gili Matra Islands [79].
Despite many respondents showing strong connections with the surrounding natural environment, some expressed uncertainty in taking care of the environment and having their skills improved. These results differed to previous findings [67], suggesting that locals in Gili Matra Islands believed that they value the learning aspect the most. This may be due to different methods. Nevertheless, the eudaemonic models in this study disclosed the residents' well-being as the result of interaction between CES and surrounding tourism and management features. The discrete pattern between models or dimensions suggested that the eudaemonic models revealed the favoured CES impact on residents' well-being. In contrast, the similar pattern between islands indicated that the overall results are valid and can be used for the case of small tourism islands. To better understand well-being values in Gili Matra Islands, additional dimensions should be investigated using exploratory factor analysis. Lastly, revisiting the same approach in future research with a larger sample size could also help investigate data consistency and temporal patterns.

Interrelations between CES Associated with Well-Being and Marine Tourism and Protection
The seven regression models for CES-associated well-being dimensions showed significant correlations with two features of interest ( Table 4). The other zones in MPA (less strict zone) allow for a wide range of regulated activities, including tourism and education [64]. Besides promoting conservation, a less strict zone also permits more activities and easier access, enhancing place identity and cultural heritage felt by the users [80]. MPAs were able to provide three roles of CES: a place to encourage well-being, a place of spirituality, and a place of preservation and freedom [75,81], or as simple as being able to inspire a positive mindset on social media activity [17]. Based on these findings, managing and shifting the current MPA core zones seaward and increasing the other zones around the island could elevate the CES well-being dimensions of residents. Nevertheless, due to the intangible nature of CES, managers of MPAs may only be able to control the spiritual value of MPAs indirectly [75].
Distance to the coastline showed a negative effect on all well-being dimensions. In Gili Matra Islands, the coastlines were surrounded by a range of MPA zones. This means that those coastlines are regulated, restricting human activities, including core zones (strict, notake zone) and other zones (less strict zones), which consist of fisheries zones, rehabilitation zones, port zones, protection zones, and utilisation zones [23,25]. The coastline next to or close to the core zone has the most restrictions, as the core zone only allows very limited activities [64,65]. This restricted area might lead to a decrease in the well-being factors observed by people living next to this segment of the coastline. A high proportion (25.05%) of the coastline is adjacent to the core zone, which could also add to the effect on wellbeing. However, this was not the case for Spiritual well-being and the core zone. Higher proximity to the core zone significantly increased Spiritual well-being. By restricting human activities and preserving natural ecosystems in the core zone, residents may have felt more spiritual well-being through the cleaner water resource [82], quiet situations, open spaces, and natural sounds that the core zone provides [83]. Distance to beach spots and the coastline visibility also impacted spiritual well-being. The beach environments are believed to improve spirituality through outdoor experiences, scenic views, sounds, and the smell of the sea [75].
In Gili Matra Islands, the amenities included tourism commercial centres, schools, and others [30]. These artificial sites presumably detracted from natural ecosystems reducing the spirituality and place identity. Furthermore, artificial structures were considered to lessen well-being and were avoided by people, even though these structures could trigger future development [84]. The urban-like growth across the islands may also impact mental health [85]. Conversely, public amenities provide a space for residents to gather and socialise. This may have been reflected in the finding that visibility of any amenity to the resident could increase Social value. Overall, these results suggested that several features of marine tourism and protection significantly influenced some dimensions of CES wellbeing felt by the residents. These features, such as coral/dive spots, MPA's core zones and other features, are unique to Gili Matra Island, which can contribute to constructing local people's well-being and help establish societal goals and engagement [86]. Some other features and CES well-being dimensions did not show a significant relationship, providing an opportunity for future studies to improve the factor analysis model for assessing CES well-being in small islands.

Conclusions
This research explored the influence of marine tourism activities and marine protection on CES in small islands with increasing tourism growth by using the Gili Matra Islands as a case study. A hedonic model investigated the CES monetary state reflected in residential property prices. Across the three islands of Gili Matra, higher housing prices were found closer to the coastline, beach spots, and coastlines with a sunset view. When explored separately, each island demonstrated a particular characteristic. The multiple models showing multiple patterns across different islands suggested that although those islands are in the same location, each island's property prices were influenced by different marine tourism and protection features. These models revealed unique influence and thus specific management action to each island. To further inform the complex influence of any existing features on the island, future research can utilise other types of anthropogenic activities on the island, such as increasing fisheries activities.
The non-monetary state of CES explored in the eudaemonic subjective well-being suggested that residents of Gili Matra Islands experienced significant CES benefits from the island environment. These benefits contributed to several well-being dimensions, primarily to Spiritual value followed by Social bonds and Place identity. All well-being dimensions experienced by the residents showed significant relationships with several features of marine tourism and conservation. The zones in MPA that applied fewer restrictions (i.e., other zones), enabling water activities for residents, showed significant interrelations with all well-being dimensions explored.

Limitations and Future Research
The hedonic CES monetary model's temporal patterns and changes over time are yet to be explored. Another limitation was the small number of samples collected for tourism accommodation properties. Future research could follow up on the findings, such as the distribution of property prices and dive spots, incorporating them into vulnerability assessment. Future research could also utilise other anthropogenic activities that may influence the price of CES-related properties.
Future research should also explore additional factors relevant to small islands by using exploratory factor analysis in the CES eudaemonic well-being model. More samples and specified categories of respondents may improve outcomes, such as those who directly interact with the features. Lastly, adapting the same methods for future research may help assess temporal patterns and consistency in findings.

Implications for Management and Policy
This study provided several valuable and novel findings that could help inform management decisions. These impacts could be incorporated into the spatial planning and MPA zoning management of the Gili Matra Islands. For example, the correlations between high property price distribution with the features of interest could indicate a higher risk of coastal hazards. Therefore, some MPA management actions could be considered, like updating vulnerability assessment and construction of coastal protection structures. The significant influence of particular marine tourism features should also be viewed as a warning to prevent the depletion of those resources by mass utilisation. Management could revise zoning plans and update the recreation fees to prevent resource depletion. Situated mainly on the eastern coast of Gili Trawangan, high-priced properties could be included in the future cost-benefit analyses to provide quantitatively informed recommendations, such as beach nourishment to prevent erosion.
The positive influence of MPA's other zones (less strict zones) on CES well-being could be considered in MPA management by shifting the current core zones seaward and allocating more zones with less strict restrictions. Additionally, spatial planning for land use could be updated as the public amenities reportedly affected Gili Matra residents' spiritual value and place identity. However, CES eudaemonic well-being values are often subjective and intangible. MPA managers must understand that they could only influence such well-being outputs passively.
Supplementary Materials: The following supporting information can be downloaded at: https: //www.mdpi.com/article/10.3390/ijerph191912078/s1, Questionnaire S1: Hedonic and Eudaemonic Questionnaire; Equation (S1): Residential Property Prices Gili Matra (model HR1); Table S1: T-test result from standard multiple regression analysis for model HR1; Table S2: ANOVA test result from standard multiple regression analysis for model HR1; Equation (S2): Residential Property Prices Gili Trawangan (model HR2); Table S3: T-test result from standard multiple regression analysis for model HR2; Table S4: ANOVA test result from standard multiple regression analysis for model HR2; Equation (S3): Residential Property Prices Gili Air (model HR3); Table S5: T-test result from standard multiple regression analysis for model HR3; Table S6: ANOVA test result from standard multiple regression analysis for model HR3; Equation (S4): Residential Property Prices Gili Meno (model HR4); Table S7: T-test result from standard multiple regression analysis for model HR4; Table S8: ANOVA test result from standard multiple regression analysis for model HR4; Equation (S5): Tourism Accommodation Property Prices, after elimination (model HT1); Table S9: T-test result from standard multiple regression analysis for model HT1, after elimination; Table S10: ANOVA test result from standard multiple regression analysis for model HT1, after elimination; Equation (S6): "Skill" factor scores for Response #1; Table S11: T-test result from standard multiple regression analysis for model R1 "Skill"; Table S12: ANOVA test result from standard multiple regression analysis for model R1 "Skill"; Table S13: T-test result from standard multiple regression analysis for model R2 "Engagement"; Table S14: ANOVA test result from standard multiple regression analysis for model R2 "Engagement"; Table S15: T-test result from standard multiple regression analysis for model R3 "Memory"; Table S16: ANOVA test result from standard multiple regression analysis for model R3 "Memory"; Table S17: T-test result from standard multiple regression analysis for model R4 "Place Identity"; Table S18: ANOVA test result from standard multiple regression analysis for model R4 "Place Identity"; Table S19: T-test result from standard multiple regression analysis for model R5 "Social"; Table S20: ANOVA test result from standard multiple regression analysis for model R5 "Social"; Table S21: T-test result from standard multiple regression analysis for model R6 "Spiritual"; Table S22: ANOVA test result from standard multiple regression analysis for model R6 "Spiritual"; Table S23: T-test result from standard multiple regression analysis for model R7 "Therapeutic"; Table S24: ANOVA test result from standard multiple regression analysis for model R7 "Therapeutic"; Table S25: Communalities values for Model E1; Table S26: The goodness of fit indicators used for Model E1; Table S27: Factor score weight produced for Model E1 Gili Matra.