Study area

This research was conducted in the Pixquiac and the Gavilanes catchments in the upper Antigua river watershed in central Veracruz, Mexico. The drainages of the Pixquiac (~100 km²; 1,024–3,759 meters above sea level [m.a.s.l.]) and the Gavilanes (~41 km²; 1,090–2,960 m.a.s.l.) rivers are located on the eastern (windward) slopes of the Sierra Madre Oriental mountain range (19°28’, –97°01’; Figure 1). These two adjacent catchments are important water supplies for the city of Xalapa (providing 38% of the water supply for 500,000 inhabitants) and Coatepec (providing 90% of the water supply for 92,000 inhabitants). Dominant land covers in the catchments consist of forests (76%), pasture and agricultural lands (22%), and urban areas (2%) (Von Thaden et al. 2019). The forests are primarily composed of TMCF (>50%), followed by mixed temperate forest (~25%) (ESA 2015). The climate is temperate humid with abundant summer rains (García 1988). About 80% of the annual rainfall falls as convective storms during the wet season (May–October), followed by a prolonged dry season (November–April). Maximum groundwater recharge and runoff also occurs during the wet season (Muñoz-Villers and McDonnell 2013).

The local climate varies markedly with elevation. Data from climate stations in or near the catchments indicate that mean annual rainfall ranges from 1,120 mm to 3,185 mm to 855 mm as elevation increases from 1,200 m.a.s.l. to 2,100 m.a.s.l. to 3,000 m.a.s.l., respectively. Mean daily temperatures (5°C to 19°C) and mean annual evapotranspiration (855–1,215 mm, estimated as reference evapotranspiration) also vary considerably with elevation (Holwerda et al. 2013; Muñoz-Villers et al. 2012). Although these spatial trends in climate have been recognized at the regional scale, climatic data in the two catchments is sparse, which has limited the characterization of critical hydrologic processes controlling runoff. As a result, the impacts of land use on hydrologic cycling and the effectiveness of conservation programs focusing on the linkages between forest conservation and hydrologic services are poorly understood.

The Pixquiac and Gavilanes catchments meet several of the enrollment criteria for Mexico’s national PHS program, including the presence of priority ecosystems, proximity to downstream cities, and water scarcity issues, as well as active local PHS programs operating in the region. Figure 1 shows the areas receiving payments from PHS programs in the study area between 2003 and 2014, covering 26% of both catchments. Climatic data is especially limited in the mid-range elevation portions of these catchments, i.e., 1800–2500 m.a.s.l., which, if regional spatial trends were to hold true, is where most of the rainfall occurs. The goal of the Quiahua network was to concentrate collection of rainfall measurements at mid-range elevations (1,500–2,500 m.a.s.l.) across our study catchments and thus to contribute the information toward efficient targeting of PHS payments.

Data collection

Modeling watershed hydrologic services

The point-based daily rainfall data summed to 5,083 observations over 621 days and was interpolated over the study area with ArcGIS to map dry season (November–April), wet season (May–October), and annual precipitation. Three different datasets of rainfall information were compared: (1) the national climate stations alone, (2) the citizen scientist volunteer data plus the national climate stations, and (3) the automatic gauge data plus the national climate stations. Thiessen polygons were used to create the rainfall distribution maps for the national climate stations alone, since this procedure is appropriate for cases where only a few observation points are available. The inverse-distance-weighted technique was used to generate the precipitation maps for the rainfall monitoring data plus the national climate stations as well as the automatic gauge data plus the national climate stations.

A hydrologic model, the Soil Water Assessment Tool (SWAT) (Watson and Philip 1985), was parameterized for the study area and calibrated against available streamflow data (López-Ramirez et al. 2020). SWAT is a semi-distributed rainfall-runoff model that has been used to simulate watersheds in tropical montane regions similar to our study area (Arnold et al. 2012; Kim et al. 2017; López-Ramírez et al. 2020; Plesca et al. 2012; Schmaltz and Fohrer 2009). Soil-type and hydrophysical properties maps (Strauch et al. 2017), land-use/land-cover maps (INEGI, 2013a), watershed boundaries, and digital elevation models (CONANP et al. 2015) were obtained from government sources and used as input to the SWAT model with the ArcSWAT tool (INEGI, 2013b). Calibration was based on minimizing the difference between simulated and observed daily stream flows at the two watershed outlets, as is typical of hydrologic modeling efforts (Arnold and Allen 1999; Gupta, Sorooshian, and Yapo, 1999; Winchell et al. 2013). SWAT streamflow output includes spatial distributions of surface runoff and baseflow contributions. In this case, the spatial distributions of baseflow contributions were used to indicate spatial locations of hydrologic significance, since baseflow is the critical element for maintaining streamflow during the dry season and buffering against droughts. We used the annual average baseflow contributed over the simulation period as a simple metric representing baseflow magnitude.

Volunteer motivations and barriers

We had a 100% response rate prior to the workshop and two weeks after the workshop, and 75% (n = 9) six months after the workshop. Volunteers were in their late thirties (mean 39.2 ± 17.8 yrs), had lived most of their lives in their community (21.7 yrs), were highly educated (67% [n = 8] had completed more than high school education), and 58.3% (n = 7) reported earning above-average incomes (>8,170 pesos/month) (Table 1). Volunteers were equally split between men and women who had an average of 3.83 family members, with few having children less than 12 years old (Table 1).

The factor analysis with non-orthogonal rotations listed in Table 2 shows motive-related questions on responsibility (4), on learning (4), on being part of something (3), on socializing (2) and on escape (3) were loaded onto five factors. Factor 1, conservation values, is a recombination of motivation themes that explains 4.4% of the variance, and results from the recombination of motivations to be part of a conservation team and to be in nature (Ryan, Kaplan, and Grese 2001). Factor 2, motivations for learning, reflects the same interests in learning new things as conceptualized in (Ryan, Kaplan, and Grese 2001) and explains 4.2% of the variation. Factor 3, feelings of responsibility for taking care of nature, explains 3.6% of the variation and largely reflects the same themes described in (Ryan, Kaplan, and Grese 2001) while integrating one escape factor in the form of, “you like to explore the environment.” Factor 4, public service, is a theme that reflects a recombination of social and learning motivations, explaining 2.6% of the variation. Finally, Factor 5, community, is also a theme that reflects a recombination of social motivations and motivations to be part of something, explaining 2.1% of the variation.

Before the training, volunteer citizen scientists were highly motivated by learning (4.81; mean value) and conservation values (4.65; mean value). Motivations of responsibility for the environment, public service, and community were also high (Table 3). After two weeks, the highest motivations continued to be learning (4.81; mean value), while public service became more important (4.63; mean value). Six months after training, the highest stated motivations for those that responded were still conservation values (4.60; mean value) and learning (4.50; mean value), although they were lower than originally stated before the training.

A repeated measures ANOVA for the five factors yielded no significant variation for Factor 1 (p-value = 0.61), Factor 3 (0.459), or Factor 4 (p-value = 0.807) (Supplemental File 2). However, the repeated measures ANOVA demonstrated significant changes over time for Factor 5 (0.026) and Factor 2 (p-value = 0.049), where a post hoc Tukey-Kramer test for Factor 5 showed significant differences before and two-weeks after training (p = 0.031), but not before and six months after training (p = 0.937). The post hoc Tukey-Kramer test failed to show significant differences at p < 0.05 between any time points for Factor 2.

Qualitative data from the open-ended survey question on the greatest motivation for participation in rainfall monitoring were consistent with survey data, showing that volunteers were overwhelmingly interested in learning more generally, and learning about water specifically. Of the nine who responded six months after the survey, six responded on their motivation to learn, for example, “I’d like to know more.” Another important theme mentioned by four volunteers was their interest in helping research, for example, “I’d like to help research on precipitation.” Also mentioned by two people was their interest in contributing to the community, and one person mentioned their interest in learning about the region they live in, which we coded as “sense of place.”

Data from the open-ended question on the greatest barriers to participation in rainfall monitoring show that the frequently cited barriers for volunteers were insufficient time (n = 4), difficulty organizing others to collect rainfall data (n = 3), or traveling to the site (n = 2), after six months, where two people mentioned more than one barrier. However, the majority (n = 5) indicated that there were no barriers to participation.

Evaluating accuracy of manual versus automatic rainfall gauges

On average, volunteers measured precipitation 565 times (91% of days) between June 2017 and February 2019 (621 days), for a total of 5,083 observations. Table 4 shows the results of the comparison between the automatic gauges and the volunteer rainfall data annually and for wet and dry seasons. The results indicated that, overall, data collected by volunteers underestimated daily rainfall by almost 1 mm (12% error), and t-tests indicated significant differences between the automatic gauge and volunteer’s rainfall data (p < 0.0001). Most of the underestimation occurred during the wet season, where volunteer data underestimated the rainfall by 16%, and the t-test results indicated a high probability that the means were different (p < 0.0001).

The mean missing data error and misread data error were –4.5 mm/day (506 missing observations) and –1.0 mm/day, respectively. The missing data and misread errors were responsible for 29% and 71% of the mean error of –0.86 mm/day, respectively. The highest fraction of missing days for a volunteer was 29%; with this particular volunteer removed, the mean fraction of missing days was less than 7%. Of the days when readings were taken by volunteers, approximately one-third (32%) occurred on dry days and where both the automatic and manual gauges registered zero rainfall and a misread error of zero. Figure 2 shows the cumulative frequency of the misread errors without the days both the automatic and manual gauges registered zero rainfall, showing that the distribution of errors is roughly symmetric around zero, with a slight bias toward negative numbers. One-third (33%) and more than three-quarters (82%) of misread errors in daily rainfall fall between –1 mm and +1 mm and –10 and +10 mm, respectively.

Although volunteers’ measurements tended to underestimate total rainfall, the errors did not translate to substantial differences in spatial patterns when compared with patterns obtained from the automatic gauges alone, as shown in Figure 3. The results in Figure 3 show dramatic differences between the rainfall distributions using the national climate stations alone and distributions using volunteer rainfall data and automatic gauge data together with the national gauge data. Especially important is that the annual and wet-season rainfall in the central portions of the study area were substantially higher (more than 300 mm and 400 mm for the wet season and annual rainfall) when looking at volunteer or automatic gauge data combined with the national gauge data, versus with the national gauge data alone. Volunteer rainfall data and national gauge data underestimate annual and wet-season rainfall by 100 mm to 200 mm in the central portion of the watershed, in agreement with the results in Table 4. In particular, the gauge located in the center of the watersheds (ZAP_P2) was responsible for much of the underestimation by volunteers at that location in the central part of the watershed.

Incorporating rainfall datasets into hydrologic models

The three sets of rainfall data were used to generate annual baseflow contribution maps with the SWAT model in Figure 4. As expected, since the rainfall in the central-to-lower portion of the Pixquiac and lower portion of the Gavilanes watersheds was higher when volunteer rainfall data or automatic gauge data were combined with the national gauge data, the baseflow contribution rates were substantially higher in these maps (Figure 4b and 4c) than for baseflow contribution rates simulated with the national gauge data alone (Figure 4a). The increase in simulated baseflow contribution rates in the middle portion of the study area for volunteer rainfall and automatic gauge data (Figure 4b), and with the national gauge data, was not as large as one might expect, given the large differences in rainfall with the national gauges alone (Figure 3). This result could be explained by the relatively steep slopes in most of the middle portion of the study area (see Figure 1), which tend to increase runoff and reduce soil infiltration, ultimately decreasing the water available for baseflow. Similar to rainfall distributions, the simulated baseflow contribution was substantially lower when using combined volunteer rainfall data with the national gauge data as compared with combined automatic gauge and national gauge data.

The average baseflow contributions for the entire study area and for the areas receiving PHS payments alone are shown in Figure 5. The results in Figure 5 (white and gray bars) indicate that the average baseflow contributions predicted by the SWAT model increase substantially when comparing the national gauges alone with volunteer rainfall data plus the national gauges, and with the automatic plus the national gauges, matching the spatially distributed results in Figure 4. The results in Figure 5 also indicate that the current targeting results in 13%–18% greater contributions to baseflow in the areas receiving payments, compared with the mean baseflow contributions across the watershed.

To further assess the efficiency of the payments with respect to capturing baseflow, the areas with the highest baseflow contributions corresponding to the same total area (3,990 ha) currently receiving payments were identified. These results (black bars in Figure 5) indicate that, as expected, baseflow capture increases when the PHS payments target areas that maximize baseflow contribution, compared with the mean baseflow contribution across the watershed. The value of volunteer rainfall data toward increasing targeting efficiency is demonstrated by looking at the improvement in baseflow contributions for the national gauges alone (30%) versus the national gauges with volunteer rainfall data (54%) when maximizing baseflow contributions.

Designing citizen science to fulfill motivations and reduce missing and misread data

Understanding who signs up to participate in citizen science, why they stay, and what barriers they face are essential questions to address when engaging volunteers in citizen science projects. This study supports previous research that shows that citizen science typically attracts educated, financially stable individuals who are motivated by learning about the projects’ topics (Deutsch et al. 2009; Deutsch and Ruiz-Córdova 2015; Moriasi et al. 2015) and are interested in learning generally (Land-Zandstra et al. 2016; Newman et al. 2012; Brossard et al. 2005). In this study, we found that citizen scientists were also motivated by their conservation values. However, motivations did not change significantly over time; rainfall volunteers remained engaged, motivated by learning and conservation values, and very few listed barriers to volunteering beyond constraints to time. The results suggest that recruitment should target individuals who are interested in engaging intellectually with the topic, and provide support for ongoing learning opportunities.

Like many citizen science projects, understanding how to maintain a consistent volunteer effort in data collection is only part of the puzzle. Ensuring that the data itself is accurate is another critical issue. Our results show that volunteers consistently underestimated rainfall compared with automatic tipping bucket data, particularly in the wet season. There is a possibility that the automatic gauges were biased; however, most authors have noted that automated tipping bucket gauges are subject to underestimation, especially during high-intensity rainfall events (Crall et al. 2013; Molini, Lanza, and Barbera 2005; Upton and Rahimi 2003; Lanza and Vuerich 2009), rather than overestimation. In contrast to our findings, Reges et al. (2016) reported that manual gauge measurements by CoCoRaHS volunteers tended to be higher than automated, tipping bucket gauges. Parsing the time series error revealed that missing data and misread errors contributed to roughly one-third (29%) and two-thirds (71%) of the mean error of –0.86 mm/day, respectively. While strategies such as collecting data forms from volunteers or providing other opportunities for interactions more frequently may decrease the incidence of missing data, reducing the misread error may be difficult. However, analysis of the distribution of misread errors, illustrated in Figure 2, indicates that eliminating larger misread errors in daily rainfall, using < –10 mm and > +10 mm as thresholds, would decrease the misread error from –1.0 mm/day to –0.2 mm/day overall. Thus, applying quality control protocols (similar to those used to analyze automated gauge data) to eliminate outlier errors in volunteer data could improve volunteer estimates substantially.

Incorporating citizen science into hydrological models to better inform payment for hydrologic services programs

Evaluations of PHS effectiveness typically use changes in forest cover as a surrogate to measure ecosystem service provision (Shedekar et al. 2016). In Mexico, CONAFOR targets payments to areas of high deforestation risk, with primary forest cover receiving priority, particularly cloud forests, and where policy makers assume that greater forest cover increases hydrological services (Wunder 2007). However, other metrics of success for PHS programs include additionality, which measures program-caused forest conservation. (Muñoz-Piña et al. 2008). Recent studies have shown that the majority of payments in Mexico (61%) go to areas with low to medium deforestation risk (Von Thaden et al. 2019). bringing into question the additionality of such programs. Adding to efficacy concerns, studies have shown that PHS programs are often lacking in areas identified as priority areas for hydrological service provisioning (Mokondoko et al. 2018). Recent work (Costedoat et al. 2015; Berry et al. In review) suggests that hydrologic metrics, such as measures of dry-season streamflows, can be valuable indicators of ecosystems services provided by forests and converted land uses. The results in Figure 5 show that including a baseflow metric for prioritizing PHS can substantially increase baseflow in the study area. The reliability of these hydrologic metrics depends on the availability of accurate, frequently collected hydroclimatic data, such as the data collected in the Quiahua citizen science rainfall effort. The results in Figure 5 also show that, despite the fact that the citizen scientists consistently underestimated rainfall, particularly during the wet season, data contributed by the citizen scientists could lead to substantially greater targeting of baseflow contributions from the PHS policy design at a much lower cost than automated tipping bucket gauges.