Resetting the baseline: using machine learning to find lost meadows

Mountain meadows occur in specific geomorphological conditions where low-gradient topography promotes fine sediment accumulation and high groundwater tables. Over 150 years of human-caused hydrological degradation of meadows along with fire suppression has resulted in decreased groundwater elevations and encroachment of upland vegetation, greatly diminishing the ecological value of meadows for water storage, baseflow, sediment capture, wildfire resistance, wildlife habitat, and carbon storage. We aimed to understand where and how frequently meadows historically occurred to reset the baseline condition and provide insight into their restoration potential. We trained machine learning algorithms to identify potential meadow areas with similar hydrogeomorphic conditions to extant meadows while ignoring their unique vegetative characteristics because we hypothesized that vegetation would change but geomorphology would remain. We used a publicly available dataset of over 11,000 hand-digitized meadow polygons occurring within a 25,300 km2, 60-watershed region in the Sierra Nevada, California USA to train random forest models to detect meadow-like hydrogeomorphic conditions. Predictor variables represented topographical position, flow accumulation, snowpack, and topographical relief at differing spatial scales. We assessed model performance and produced maps delineating high probability meadow polygons. Our findings showed that there is nearly three times more potential meadow habitat than currently documented. The predicted area includes a mixture of existing but undocumented meadows, non-meadowlands that may have converted from meadows due to lost function and forest encroachment, and areas with meadow-like geomorphology that may never have been meadow. The polygons encompassing predicted meadows often expanded existing meadows habitats into adjacent areas with continuous topography, but with upland vegetation and incised channels. Using readily available data and accessible statistical techniques, we demonstrated the accuracy of a tool to detect about three times more historical meadows than currently recognized within a complex, mountainous landscape. This “found” area greatly increased the potential area that could be subject to meadow restoration with benefits for biodiversity, wildfire management, carbon sequestration, and water storage.


Introduction
Meadows are receiving attention from land managers because they are ecologically important habitats and most of them have become degraded and need restoration. They retain water Hill and Mitchell-Bruker 2010;Ankenbauer and Loheide 2016), serve as carbon sinks (Norton et al. 2011;Reed et al. 2021), and provide habitat for a high diversity of plants and animals (Allen-Diaz 1991; Oles et al. 2017;Ziaja et al. 2018;Campos et al. 2020). These characteristics are increasingly important in landscapes experiencing the effects of climate change such as more extensive wildfires or protracted droughts. Yet, meadows are rarely considered during drought and wildfire resilience planning in forested landscapes, likely because meadows tend to be small, isolated, and degraded. Therefore, meadow restoration efforts are thought to be more beneficial at the site level rather than for their ability to affect broader landscapes. For example, debate exists whether meadow restoration can substantially increase water storage to affect downstream baseflow (Hunt et al. 2018Nash et al. 2018Nash et al. , 2020. The rationale is that meadows have too small of a capacity to store and release water to have a significant effect on hydrological processes (Nash et al. 2018). But if meadows were larger and more common, they could become important components of efforts to manage climate change and other landscape-scale threats to forests and water supply.
Many meadows historically supported or were created by beavers (Castor canadensis in North America and C. fiber in Europe) through their damming of unconfined, low-gradient channels (McComb et al. 1990;Pollock et al. 2003;Polvi and Wohl 2012). New evidence shows that beaversupported wetlands in forested landscapes can significantly increase groundwater and surface water storage and decrease water temperature (Dittbrenner et al. 2022). These changes can influence wildfire dynamics and provide important fire breaks in otherwise burned landscapes (Fairfax and Whittle 2020). Since beavers were much more common prior to the 19th century fur rush (Naiman et al. 1988), restoration of beavers and beaver habitats would likely increase these hydrologically resistant landscapes (Wohl 2021;Dittbrenner et al. 2022;Jordan and Fairfax 2022). Similarly, meadows were likely larger, wetter, and more common in the mountainous western United States and Canada prior to Euro-American settlement and subsequent meadow degradation (Loheide and Booth 2011;Celis et al. 2017;Lubetkin et al. 2017). Understanding how large and common wet meadows once were could provide insight into their capacity for affecting water storage, ecologically important summer baseflow, sediment capture, wildfire resistance, wildlife habitat, and carbon storage at both local and landscape scales.
Mountain meadows tend to occur in specific geomorphological conditions that promote fine sediment deposition and elevated groundwater tables (Allen-Diaz 1991; Miller et al. 2001;Shaw and Cooper 2008). In California's Sierra Nevada, they typically occupy low-gradient geologic benches at elevations from 1500 to 3000 m where they can be recharged annually by snowmelt (Wood 1975). Meadows support vegetation dominated by grasses, sedges, herbs, and shrubs such as willows and alder (Weixelman et al. 2011). These specialized plants have adapted to a high groundwater table and are sensitive to changes in groundwater elevations (Allen-Diaz 1991; Loheide and Gorelick 2007;Lowry et al. 2011).
This sensitivity of meadow plants leads to dramatic changes in the composition of vegetation when reductions in groundwater elevations occur because of loss of infiltration capacity, increased runoff rates, or channel incision (Loheide and Gorelick 2007;Hammersmark et al. 2009;Loheide et al. 2009). Humans simplified and confined flow paths to access timber and other natural resources and used meadows extensively for livestock grazing resulting in concentrated flows and incised channels that drained groundwater to the elevation of the new channel bed (Ratliff 1985;Loheide and Gorelick 2005). With the addition of fire suppression, surrounding upland vegetation encroached onto the newly dry meadow surfaces and, in many places, replaced meadow vegetation (Halpern et al. 2010;Celis et al. 2017;Lubetkin et al. 2017). Researchers recently quantified forest encroachment and closure in the Rocky Mountains of Canada using repeat photography and found that forest coverage increased by over 35% during the 20th century and that meadow losses were most acute at high elevations (Stockdale et al. 2019).
Given humanity's propensity for shifting the baseline of environmental change to more recent conditions (Essl et al. 2015;Fernández-Llamazares et al. 2015), we are likely underestimating the historical extent of meadows, and therefore, their importance for ecological restoration gains as we attempt to sustain biodiversity, promote wildfire resilience, and mitigate the effects of climate change on our forested landscapes. Here, we attempt to recognize the baseline historical extent in the Sierra Nevada by training machine learning algorithms to identify characteristic geomorphic, hydrologic, and climatic settings that supported stream-associated meadows. We hypothesized that, while the vegetation has likely changed due to degradation, the underlying hydrogeomorphologic signatures remain. We test if these signatures can be identified using machine learning techniques and available data (Jordan and Mitchell 2015) to reset the baseline for landscape and local restoration planning.
We anticipate this approach can be applied across mountain ranges, but focused our modeling efforts on the Sierra Nevada, California U.S.A. because the range supports over 20,000 extant meadows, most of which have been hand-digitized and are represented on a publicly accessible meadow polygon layer (UC Davis and USDA Forest Service 2017). In addition, the Sierra Nevada Range supports the freshwater needs of millions of people (Liu et al. 2021), yet increased warming and drought have led to reductions in snowpack and earlier snowmelt (Maina et al. 2022). Widespread declines in dry season baseflows and earlier peak flows have caused concern about the future of the water supply (Patterson et al. 2022). The Sierra Nevada also supports a large wildland urban interface that is increasingly vulnerable to wildland fires fueled by dense, drought-stressed, and diseased forests (Van Gunst et al. 2016;Restaino et al. 2019;Das et al. 2022). These important water and wildfire concerns are motivating land managers to identify opportunities to restore ecologically important meadow habitats to retain water and reduce wildfire stress. We anticipate the model results described here will support watershed restoration planning to improve the range's resilience to climate change and resistance to wildfire.

Study area
Our study area encompasses 60 Hydrologic Unit Code 10 (HUC10) watersheds covering 25,300 km 2 with elevations ranging from 246 m to 4413 m (Fig. 1). The area extends from the Upper North Yuba River watershed in the north to the Brush Creek-Kern River watershed in the south and includes watersheds both west and east of the Sierran Crest. We selected this study area because all but five of the watersheds each support at least 100 hand-digitized meadows included in the Sierra Nevada MultiSource Meadow Polygons Compilation Version 2 (SNMMP, http:/ meado ws. ucdav is. edu); the remaining five watersheds were included to create a contiguous geographic area (Online Resource 1: Model Planning).
The study area's latitudinal and elevational range supports a diversity of climatic conditions and associated forest types. The western slope of the Sierra Nevada receives substantial precipitation during fall and winter months, mostly as snow at elevations above 1800 m, while the eastern slope receives significantly less precipitation (Albano et al. 2019). The gentler western slope is dominated by mixed conifer forest while the steeper eastern slope is dominated by dry forest and shrubland. The study area has been subject to hydropower development, mining, agricultural land use, and urban development but consists primarily of public lands managed by the USDA Forest Service and National Park Service.

Model development
We applied random forest models (Breiman 2001) to 60 watersheds in the Sierra Nevada to uncover the potential historical extent of meadows (Fig. 1). Of machine learning algorithms, random forest models have proven to effectively classify land cover characteristics to improve our understanding of landform development across a range of scales (Liu et al. 2020;Tan et al. 2021;Donager et al. 2022). The predictive accuracy of the models relies on accurate remotely sensed topographic, ecological, and geophysical data used to train the machine learning algorithms (Houser et al. 2022, Fig. 2). Accurately classified training data allow the models to distinguish the landforms of interest (meadows) from complex, high relief background landscapes found within the Sierra Nevada.

Meadow training data
We used extant meadows included in the SNMMP dataset clipped to our study area to train our models (Fig. 2). We first subset the dataset to include stream-associated meadows (within 50-m of a DEM derived flowline) and exclude fens and discharge slope and depressional meadows unassociated with surface flow paths (Weixelman et al. 2011), leaving 11,127 meadow polygons covering 442.6 km 2 (1.75% of the study area). From these polygons, we sampled 1000 random points per watershed to create positive training data for the 60 "local" HUC10 watershed models. Similarly, we sampled 9000 random points outside of the meadow polygons per watershed to create negative training data. We chose a 1 to 9 ratio in the training datasets because meadows comprise a small proportion of the landscape, even historically, and non-meadow habitats can be extremely Fig. 1 The study area in the Sierra Nevada, California, USA. Red perimeter delineates the area of the "Sierra Nevada" model and black lines delineate the 60 Hydrologic Unit Code 10 watershed "local" models. Blue dots represent the extant meadow training data from the Sierra Nevada MultiSource Meadow Polygons Compilation, Version 2 (UC Davis and USDA Forest Service 2017) geomorphically diverse. We selected a maximum 10,000 training points as a tradeoff between accuracy and computing power needed to run each model. We used the same sampling methodology to create a 25%-sized testing dataset with 250 meadow points and 2250 non-meadow points.
We also sampled from the entire study area to create a "Sierra Nevada (SN) model" by selecting 1000 stratified random positive and 9000 stratified random negative training points across every watershed in proportion to the relative area of meadows for positive points and watershed areas for negative points. Individual watersheds contributed between 8 and 28 positive training points and 60-338 negative training points to the SN training dataset. We also created a 25%-sized testing dataset to assess performance of the SN model. We hypothesized that local models with abundant training data (lots of meadows within the watershed) would perform better than the SN model for those watersheds, but that the SN model would perform better in watersheds with sparse or incomplete training data and could be better applied to watersheds without training data outside the study area. To turn these sampled points into a response variable for the training data, we assigned a zero to points sampled outside existing meadows and a one to points sampled with existing meadows.

Predictor variables
We identified geomorphic, hydrologic and snowpack variables that describe conditions known or expected to support stream-associated meadows (Wood 1975;Weixelman et al. 2011;Albano et al. 2019). Streamassociated or "riparian" meadows usually occur in low slope, depositional zones with sufficient upslope area to deliver both surface and subsurface flow late into the dry season to maintain a high groundwater table (Weixelman et al. 2011). We therefore selected geomorphic variables that represent topographical position, relief, and flow accumulation.
Hydrogeomorphic variables were derived using elevation data from the 3D Elevation Program (3DEP) digital elevation model (DEM) at 1 arc second (~ 30 m) and 1/3 arc second resolutions (~ 10 m; U.S. Geological Survey 2019) (Table 1). Rasters were developed for relative local elevation, slope, topographic variability, topographic wetness, and distance to stream channels and included 5 × 5 relative elevation (elev_5 × 5_rel), 5 × 5 elevation standard deviation (elev_5 × 5_std_dev), 5 × 5 slope range (slope_5 × 5_range), 5 × 5 slope standard deviation (slope_5 × 5_std_dev), topographic wetness index (TWI) at 10 m scale (twi_10m), TWI at 100 m scale (twi_100m), horizontal flowline distance to stream channel (dd_h), vertical flowline distance to stream channel (dd_v) and surface flowline distance to stream channel (dd_s, Table 1). We used a 5 × 5 pixel 'moving window' to calculate the local variation surrounding the focal pixel including elevational variability, topographic roughness, and slope. Other than twi_100m (1 arc second DEM resampled to 100 m), we used 10 m pixels with 50 m × 50 m sampling windows. We used two scales for TWI because we hypothesized water accumulation patterns that occur in meadows may prove more relevant at the 100 m scale. The moving window approach has proven accurate for interpreting the characteristics of complex landscape structure (McGarigal and Cushman 2002;Zhu et al. 2020).
The flow distance and TWI variables were derived using Terrain Analysis Using Digital Elevation Models (TauDEM, Tarboton 2015), a suite of tools for analyzing hydrographic information from topography. Hydrologic sinks (low points in the elevation file that have no expected outflows) in the DEM were filled to define a drainage path from each grid cell. A flow accumulation raster was then derived from the filled DEM. A channel network raster was derived from the flow accumulation raster where channels had a minimum threshold of 350 upslope pixels (roughly matching the flowlines in the NHD dataset) to be considered a channel. TauDEM calculates the distance to a stream using a D-infinity model with each pixel having multiple flow directions where the outflow from each pixel is partitioned between up to two downslope pixels (Tarboton 1997). The TWI is calculated as a ratio of the natural log of the accumulation area to slope; ln( ∕tan ) where is the catchment area and is the slope. We offset by 0.0001 to retain any pixels where the slope is zero.
In addition to topographic variables, we incorporated the median snowpack on 1 April (2010-2020) to represent available snowpack for late season water delivery, which was previously found to be an  (Flint et al. 2021). We sampled the developed predictor raster datasets using the random positive and negative training dataset points to generate the training and testing datasets.

Random forest analysis
We used random forest models to assess the relationship between the predictor variables and existing meadows on the landscape. Random forests are an ensemble modeling technique that combines a large set of independent decision trees, each using a bootstrap sample of the data to reduce the overall variance of the individual predictions while also decreasing the correlation between each decision tree, and thus, increasing its predictive accuracy (Breiman 2001;Liaw and Wiener 2002). Combining multiple decision trees is achieved by randomly sampling from the original training dataset, with sample reposition, to generate different subsets of the data. The subset trees are split according to a series of if-then rules maximizing the variance among resulting "nodes" and "leaves". The algorithm outputs continuous values between zero and one with zero representing nonmeadow predictions, one representing meadow predictions and values between zero and one reflecting the model's relative uncertainty in that prediction.
We refined the models by varying the number of trees (ntrees, 100-2000) and the number of variables sampled at each split (mtry, 2-6) and the iteration with the lowest mean of squared residuals across models was used for the final models. Mean squared residuals generally remained static after 300 ntrees and most watersheds had the lowest mean square of residuals at mtry = 4 so these parameters were used in the final models.
To explore the nature of the relationship between the predictor variables and the response variable, we graphed density plots of raw data, ranked variable importance within the models, and created partial dependence plots of marginal effects (Friedman 1991). Density plots were generated using a subsample of the testing data that were filtered to remove rare outliers by excluding the bottom and top 1% of each variable. This step serves a practical purpose of making the graphs more interpretable at relevant scales.
We calculated variable importance ranks by calculating the mean decrease in accuracy at each decision tree when the focal variable was randomized. Partial dependence plots show the marginal effect of each predictor variable when all other variables are held to their mean value.
After models were trained, we created several geospatial datasets from the predictions of the models across the study area. First, 10-m raster grids of continuous predictions (0 to 1) were constructed for each watershed for both the local and SN models. We aggregated clusters of pixels with high predictions to generate polygons of potential meadows. One method for selecting a threshold to assign predictions as meadow or non-meadow is to identify the threshold that minimizes prediction errors in the testing dataset. We instead used a consistent, higher confidence threshold (0.5) for both local and SN models because accuracy in identifying areas that likely were historical meadows was more important to us than maximizing the identification of all potential meadow areas.
Once the two polygon datasets were created, we removed polygons that were less than 0.1 ha, removed areas that overlapped with ponds and lakes (California Department of Fish and Wildlife 2012), removed areas that overlapped with National Land Cover Dataset classes that would never support meadow (Perennial Ice/Snow and Barren Land; NLCD 2019), and filled holes in the polygons that were less than 400 m 2 . These changes made for more accurate estimates of total historical meadowlands and allowed for comparison of this dataset with other datasets. We then compared the polygon datasets to existing hand-digitized polygon datasets and to each other. We compared overlap and unique area predictions to determine conditions where models most succeeded and failed.

Random forest model performance
We took a three-stage approach to model validation and assessment: (1) predictive accuracy of models with held back training data, (2) comparing model predictions to an alternative dataset of extant meadows, and (3) examination of predicted locations with higher resolution LiDAR-derived digital elevation models (DEMs) to identify potential sources of degradation that may have caused conversion from meadow habitat to other habitat types.

Predictive accuracy
We used the area under the receiver operator curve (AUC) to evaluate local and SN models' fits at each watershed. The AUC represents the probability that true positives are ranked higher than true negatives in an aggregate measure of performance. Scores range from 0.0 (all predictions incorrect) to 1.0 (all predictions correct) with 0.5 representing a random model (Fielding and Bell 1997). We used the 25% testing datasets to calculate the AUC for each watershed.

Comparison with alternative dataset
The Sierra National Forest (Sierra NF) maintains an independent meadow dataset that includes some meadows not found in the SNMMP dataset and has different shapes of meadows shared in both datasets. The Sierra NF dataset was derived by local land managers using both remotely sensed and field discovery and validation and thereby provides an opportunity to validate model predictions at locations that were not used to train the model. We subset the Sierra NF dataset to match the spatial area used to derive the models and similarly limited it to meadows greater than 0.1 ha occurring within 50 m of flowlines. The resulting dataset included 6,566 meadow polygons compared to 3,713 meadow polygons used in the SNMMP training data. We compared the quantity and overlap of meadow polygons across the three datasets (predicted meadows, SNMMP meadows, and Sierra NF meadows) with the expectation that the predicted meadows would overlap with a greater proportion of unique Sierra NF meadows than by chance.

Comparison with imagery
We randomly selected 25 model-predicted meadow polygons from 4 watersheds where high resolution LiDAR data were available to visually determine the appropriateness of the designation, signs of degradation, and vegetative characteristics. We used aerial imagery and LiDAR-derived DEMs to document presence or absence of a flat floodplain, visibly incised channels, road crossings, meadow vegetation, and forest vegetation within the predicted meadow polygons.
Except for raster processing using TauDEM (described above), all spatial, statistical, and graphical analyses were conducted using R 4.0.5 (R Core Team 2021). Random forest models were trained using the package randomForest (Liaw and Wiener 2002). Spatial analyses were conducted using the terra and sf packages (Pebesma 2018: Hijmans 2022.

Overview of model results
The 60 local models and one SN model predicted areas with similar geomorphic and climatic characteristics to extant meadows in the Sierra Nevada. The SN model predicted 2.6 times greater potential meadow area. The local models predicted 1.2 to 7.7 times the area of mapped extant meadows with a mean of 2.7 times (Online Resource 1: Model Output).
The SN model predicted 42,812 total meadow polygons representing 82,096 ha or 3.2% of the study area. Most of the predicted meadow polygons were small (mean = 1.9 ha, median = 0.38 ha) compared to the training data (mean = 6.0 ha, median = 1.2 ha). However, the 6507 predicted polygons larger than 2 ha covered 64,689 ha compared to 40,124 ha for the 4219 SNMMP polygons larger than 2 ha. In many cases, the large, predicted polygons encompassed one or multiple SNMMP polygons (Fig. 2). The threshold that minimized prediction error in the testing dataset across the 60 local models averaged 0.13 (range: 0.06-0.25), much lower than the high-confidence threshold we selected.

Importance of predictor variables
The most important variables for predicting meadows included topographic wetness index at the 100 m scale (twi_100m), snowpack, standard deviation of local elevations in a 5-pixel x 5-pixel moving window (elev_5 × 5_std_dev), and the vertical distance to a stream channel along the flow line (dd_v) (Fig. 3). Meadows tended to occur in areas with relatively high TWI, low local relief, and minimal change in vertical flow distance to the stream channel compared to non-meadow areas (Fig. 3). Snowpack was important overall but varied across watersheds: it had low predictive power in low elevation watersheds but high predictive power in watersheds with high elevational relief where meadows occurred primarily in the upper watershed. Partial dependence plots show similar patterns of variable relationships (Supplemental Materials 2).

Predictive accuracy
For the watershed models, 59 of 60 local models had AUC values above 0.89 (one watershed scored 0.67) and the SN model scored above 0.80 when assessed using testing data for the individual watersheds (Online Resource 1: Model Validation). When the SN model was assessed using the composite testing data from all watersheds, it had an AUC value of 0.90. Local models outperformed the SN model at every watershed based on AUC rankings, although, 42 of 60 watersheds had SN AUC values within 0.05 of the local model's AUC. The local and SN model results varied depending on the quantity and accuracy of positive training data for a particular watershed. In general, the SN model predictions were similar to the local models except in watersheds where meadows occurred in relatively unique hydrogeomorphic conditions (Fig. 4).

Comparison with Sierra NF dataset
The SN model predicted 7,025 meadow polygons with a total area of 14,254 ha or 2.8% of the 509,168 ha Sierra NF study area. For comparison the SNMMP training dataset included 3,713 meadow polygons with a total area of 7,415 ha or 1.5% of the Sierra NF study area and the Sierra NF Fig. 3 Relationships between predictor variables and existing meadows and the importance of predictor variables in random forest modeling. A Relative density of predictor variable values for positive meadow training data (blue lines) and non-meadow training data (green lines). B Scaled model variable importance from the local (small black dots) and "Sierra Nevada" (blue triangle) random forest models predicting potential meadow. Higher values indicate increased variable importance. Yellow triangle indicates median variable importance value for 60 local models. Variable abbreviations and definitions are outlined in Table 1. dataset included 6,566 polygons with a total area of 8,007 ha (1.6%). Of the 1,878 ha (23%) of Sierra NF meadows that do not overlap the SNMMP dataset (i.e., meadows not used to train the model), the SN model successfully predicted 524 ha (28%). This represents a 10-fold increase in predictions covering existing testing meadows compared to portions of the landscape that do not have documented meadows (28% vs. 2.8%). Considering meadow counts instead of areas, of the 5627 (37%) of Sierra NF polygons that do not touch SNMMP polygons, the model identified 1846 (33%).

Comparison with LiDAR-derived DEMs and imagery
Review of a random selection of predicted meadows that overlapped with existing meadow areas using DEMs revealed that the models consistently extended meadow areas to encompass contiguous flat floodplains (Fig. 4). Review of aerial imagery showed both models predicted less area than included in the SNMMP training data (local: 0.9×, SN: 0.8×). Many of the training meadows in the Middle Fork San Joaquin River Watershed (HUC10: 1804000604) occurred along steep hillslopes with minimal flow paths so had hydrogeomorphic characteristics contrary to much of the rest of the training data. In this case, the local model did a better job of identifying these meadows than the SN model these areas often dominated by vegetation indiscernible from the surrounding upland vegetation (Fig. 4). A close review of aerial imagery at predicted polygons that do not overlap with existing meadow polygons showed evidence of meadow vegetation (35 of 67), although usually within a mosaic of upland forest. About 90% of both existing and predicted meadows showed evidence of incised channels. A review of the largest predicted meadow polygons revealed a few that were likely historically meadows and wetlands but had been dramatically altered for human uses such as drainage diversions for drinking water (Owens Valley) and airports. However, most of the large, predicted meadow polygons showed no visible human infrastructure to preclude restoration, but clearly had incised channels and upland vegetation (e.g., Gerle Creek shown in Fig. 2).

Discussion
We trained machine learning models to identify high probability potential meadow habitats for 60 watersheds in the Sierra Nevada, California. The local watershed-scale models identified one to eight times more potential meadow area than is currently mapped, and the "Sierra Nevada" (SN) model that encompasses all 60 HUC10 watersheds identified nearly 3 times more potential meadow area than is currently identified. Additionally, the SN model output located 42,812 new high probability potential meadow areas > 0.1 ha in size across the study area. Expanding the number and size of meadows substantively resets their ecological baseline importance by countering shifting baseline syndrome and resetting the baseline to historical conditions. These modelpredicted meadows represent areas with similar hydrogeomorphic characteristics of existing meadows and could therefore be considered during watershed restoration planning as areas for increasing groundwater storage and floodplain connectivity. The output maps are available for use in single meadow or watershed restoration planning. In a related paper, we conceptually showed how the model output can be incorporated into watershed restoration planning to help manage wildfire, increase base flow, and capture post-wildfire sediment pulses (Pope and Cummings 2023). By restoring hydrological processes in these low-gradient meadowlands, we can start to recover a catchment's natural storage potential. Many of the model-predicted meadows could be considered "lost meadows" because they historically may have been meadows that later transitioned to forest as natural processes and human disturbances caused forests to enclose meadows and other open spaces (Stockdale et al. 2019). Indeed, in our analysis of imagery in a subsample of the identified potential meadows, we found that most of the modelpredicted meadows support vegetation more similar to the surrounding upland habitat-types than meadows. The LiDAR-derived DEMs consistently showed incised channels that drain these habitats, which would allow for encroachment of upland vegetation once the water table had been lowered (Loheide and Gorelick 2007). It is important to note that during a stratigraphic assessment of southern Sierra meadows, Wood (1975) found no evidence of distinct channels within the meadows until the arrival of Euro American people. He also found that meadow habitats had been stable for about 2,000 years prior to rapid degradation presumably associated with extraction of natural resources and livestock grazing that started around the 1850s (Wood 1975).
We focused on stream-associated meadows that are supported by surface flow and groundwater. These "riparian meadows" tend to occur lower in the watershed than meadows only supported by snowmelt (Weixelman et al. 2011) and, therefore, may be more likely to benefit from meadow restoration activities. Restoring these groundwater-connected habitat-types could increase carbon storage (Norton et al. 2011;Reed et al. 2021Reed et al. , 2022, increase resistance to wildfire (Fairfax and Whittle 2020), and sustain ecologically important dry season baseflows that support downstream aquatic organisms and human uses Snyder et al. 2015;Johnson et al. 2017). With the current urgency to reduce greenhouse gas emissions (IPCC 2022) and mitigate the effects of climate change on freshwater resources (Parmesan et al. 2022), maximizing low risk and low-cost restorations to regain lost meadow and floodplain could provide a myriad of benefits when included in watershed restoration and forest resilience planning (Jordan and Fairfax 2022;Pope and Cummings 2023).
However, some of the model-predicted "lost meadows" will probably not return to actual meadows even if restoration is implemented given that meadows are among the most vulnerable ecosystems to climate change (Hauptfeld et al. 2014). Historical meadows that are highly sensitive to changes in snowpack are likely to continue converting to upland vegetation due to continued decreases in spring snowpack and snow melt (Drexler et al. 2013).
Our ability to visualize the historical baseline of the extent of meadows within forested landscapes and the potential for restoration gains has been increased dramatically by advances in remote sensing and machine learning. We selected random forest models for the analysis because they have been found to be powerful machine learning algorithms for geomorphic mapping and they have great capacity for data explanation (Houser et al. 2018;Siqueira et al. 2021). Random forest models can work with relatively few predictors, making the relationships and importance of the predictors more easily interpretable (Siqueira et al. 2021). This is important for understanding whether the associations are meaningful and could potentially be extended beyond the training area (Houser et al. 2022). Of the predictor variables used in the SN model, topographic wetness index at the 100 m scale (twi_100m), standard deviation of local elevations in a 5 pixel × 5 pixel moving window (elev_5 × 5_std_dev), the vertical distance to a stream channel along the flow line (dd_v), and snowpack had the highest relative importance. We interpret these predictors as highly meaningful and generalizable for finding locations of meadows. Basically, the four most important variables define areas with relatively smooth and low-slope surfaces with unimpeded flow paths where water accumulates from precipitation and snowmelt. We used two scales of TWI and found the 100 m scale more informative than the 10 m scale. We assume that the 100 m scale is more relevant for the hydrogeomorphic processes (e.g., sediment accumulation, high groundwater table) to result in meadow formation.
Overall, both local and SN models performed extremely well based on accepted machine learning validation metrics. We trained models at both scales because we wanted to maximize accuracy and applicability. We found local models most effective in watersheds with plentiful, accurate training data and where the geomorphic conditions of meadows are unique compared to surrounding watersheds (Fig. 4). The SN model performed well across the study area, including in watersheds with little training data, suggesting it could be generalizable to the entire Sierra Nevada and easily adapted to other ranges that support similar stream-associated meadows and have them mapped to use as training data.
In a few sub-watersheds, our selection of a high confidence threshold resulted in less meadow habitat predicted than was delineated in the SNMMP database. This counterintuitive result seems to have occurred in watersheds with poorly or uniquely delineated existing meadows. The negative training dataset was randomly selected from areas not defined as a meadow in the SNMMP dataset, so if a relatively large area of existing meadows was not included in the database, it resulted in a higher proportion of the negative training data resembling the positive training data, and ultimately, less predicted meadowlands. This occurred in the Middle Fork San Joaquin River Watershed (Fig. 4). In this case, the local model outperformed the SN model because the SNMMP dataset included unusual stream-associated meadows on steeply sloped hillsides. The local model received more positive training data matching these conditions and so did a better job of identifying the unique meadow areas. In addition, underpredicting existing meadows as seen in this watershed and in the validation using an alternative dataset may indicate the threshold of > 0.5 is conservative and the estimates of historical meadow presented here are likely conservative too. Regardless, when meadows were not clearly geomorphically distinct from the surrounding landscape matrix, it was more difficult for the models to detect potential meadow with high confidence, resulting in less area predicted and highlighting the importance of consistent and accurate training data.

Conclusions
The models developed in this study help identify possible meadows that once occurred across the Sierra Nevada. These predictions nearly tripled the amount of area that was previously mapped as meadow. In some instances, the predictions made by these models expand existing known meadows, while in others new potential meadows are identified. These results have strong implications for meadow restoration, particularly in demonstrating the large, previously unrecognized area for possible meadow restoration across the Sierra Nevada. More broadly, meadow restoration can help to mitigate the consequences of climate change by modulating stream flows and supplementing baseflow, storing carbon, increasing soil moisture resulting in reduced wildfire behavior, and serving as important habitats for a diverse range of plants and wildlife.
Author contributions AC and KP contributed to the study conception and design. Material preparation, model development, and analysis were performed by AC and GM. The first draft of the manuscript was written by AC and KP. All authors read and approved the final manuscript.
Funding Funding to support this work was provided by the USDA Forest Service Pacific Southwest Research Station.
Data availability Geospatial data developed in this research including prediction rasters and polygons are housed at the USFS data archive (https:// www. fs. usda. gov/ rds/ archi ve/). R script can be requested from the corresponding author.

Competing interests
The authors declare no competing interests.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http:// creat iveco mmons. org/ licen ses/ by/4. 0/.