Tall willow thickets return to northern Yellowstone

. Northern Yellowstone National Park provides an example of passive restoration, as wetlands and riparian areas there lost most tall willows in the 20th century, due to intensive herbivory by elk ( Cervus canadensis ). Following large carnivore restoration in the late 1990s, elk numbers decreased, and some researchers reported willows growing taller with reductions in browsing, evidence of a shift toward willow recovery. Others questioned the extent and signi ﬁ cance of these changes. To investigate how willow heights have changed in northern Yellowstone since the 1990s, and to assess the importance of browsing as a driver of willow height and canopy cover, we compared data from three time periods: 1988 to 1993 when elk densities were high and most willows very short, 2001 to 2004 when willows may have begun to recover, and 2016 to 2018. We found a strong contrast between sites along streams, compared to wet meadows (meadow sites). Willows in meadow sites did not increase in height, but willows in stream sites increased signi ﬁ cantly, exceeding 200 cm mean height in summer by 2001 – 2004, and in spring by 2016, a height indicative of recovery. Where height did not increase, this was due to loss of annual growth to herbivory. Overall willow height distribution changed from mostly short to become clearly bimodal, with a new peak around 300 – 400 cm. Bison increased, and in some sites where bison congregate willows remained suppressed at heights below 120 cm, a condition strongly correlated with summer browsing. We also located and measured willow thickets (groups of 5 of more willows >200 cm spring height and <2 m apart) along streams in the study area. We found willow thickets in all stream reaches surveyed, a signi ﬁ cant change from past conditions. Thickets occupied >80% of willow patches in some sites, but as little as 22% in others. Tall willow thickets are an important habitat feature and an indicator of willow recovery. Our results demonstrate that the reduction of elk herbivory over the last two decades in northern Yellowstone has allowed willows to grow taller in many places, despite a warming and drying climate, while increased herbivory by bison continues to suppress willows in some locations.


INTRODUCTION
Stream banks and floodplains in the Rocky Mountain region of the western USA typically contain tall willows (Salix spp.), along with cottonwoods (Populus spp.), quaking aspen (Populus tremuloides), and other trees and shrubs. Willow bushes anchor against erosion, producing overhanging banks and shade that are important for fish. Branches drop leaf litter and wood into streams, providing food and habitat for invertebrates. Willows support beavers (Castor canadensis) that may build extensive dam complexes on small streams, or on side channels of large streams, expanding wetland habitat for willows as well as for amphibians, waterfowl, and many other species. Throughout the American West in the 20th century, many riparian areas lost willows and other woody plants due to livestock as well as wild ungulates browsing on streamside shrubs, and trampling fragile banks and wetlands (Ferguson and Ferguson 1983, Belsky et al. 1999, Baker 2003.
Many elk (Cervus canadensis), bison (Bison bison), and other native ungulates including deer (Odocoileus spp.) and moose (Alces alces) spend winters in the large valleys of the upper Yellowstone and Lamar Rivers of northwestern Wyoming, an area called the northern ungulate winter range, or northern range, of Yellowstone National Park (YNP). Willows in this area were strongly affected by herbivory for most of the 20th century (Chadde and Kay 1996, NRC 2002, Barmore 2003, Wagner 2006. Concern about damage to the range led to a culling program in the park to reduce the elk herd, but culling was not sufficient to reverse the suppression of willows. Culling ended in 1968, and a debate ensued concerning the role of herbivory in the loss of tall willows and aspen, and the degree to which this was mediated by elk density (Houston 1982, Wagner et al. 1995, Kay 1997, Singer et al. 1998, Keigley 2000, Wagner 2006). Researchers generally agreed that browsing was the proximate cause of willow suppression (Singer et al. 1994), but differed in their emphasis on the possible effects of drought. With high elk densities, browsing had strong effects on willow communities, limiting height and crown size and preventing development of thickets (Dobkin et al. 2002, Zeigenfuss et al. 2002, Singer et al. 2003, Anderson 2007. The northern Yellowstone elk winter range provided a rare opportunity to observe the effects of large carnivore restoration as a form of passive restoration of the range . Wolves, extirpated from the park by 1930, were reintroduced in 1995-1996, as predation by bears (Ursus spp.) and cougars also was increasing (Ruth 2004, Barber-Meyer et al. 2008. By 2003, browsing had decreased and willow height and cover increased in some sites, suggesting a trophic cascade, where predators indirectly benefit willows through effects on elk (Singer et al. 2003, Beschta and Ripple 2007, Beyer et al. 2007. Much of the research on northern range willows during this time was led by Francis Singer, who in 2003 observed that: ". . .willow stands have less browse pressure, are taller, and are being released from browsing suppression since wolf restoration in 1995" (Singer et al. 2003: 476). As elk densities declined further (Painter et al. 2015), willows grew taller in some locations (Tercek et al. 2010a, Baril et al. 2011, Painter and Ripple 2012, Beschta and Ripple 2016, exceeding an average height of 200 cm by 2010 in some study sites (Marshall et al. 2013), a threshold indicative of tall willow recovery (Dobkin et al. 2002, Singer andZeigenfuss 2003).
As willows have responded to reduced herbivory in the northern range, increases in willow height and canopy cover have been mediated by site conditions such as moisture availability and soil nutrients (Tercek et al. 2010a, Johnston et al. 2011. Stream incision or a low water table may slow or prevent willow growth or expansion in some sites (Wolf et al. 2007, Bilyeu et al. 2008, Marshall et al. 2013. Herbivory also continued at high rates in some places, particularly where bison were in high concentrations and browsed on willow in the summer, as in Lamar Valley Ripple 2012, Beschta et al. 2020).
The recent shift from herbivory as the primary factor limiting willow height and canopy cover, to greater importance of bottom-up factors, is evidence of a major ecological change in northern Yellowstone (Singer et al. 2003, Tercek et al. 2010a, Baril et al. 2011). This change is part of a pattern of increased size of willows, and also of young aspen, cottonwood, alder, and berry-producing shrubs, in association with a reduction in browsing (Ripple and Beschta 2012, Painter et al. 2014, Beschta and Ripple 2016, Klaptosky 2016. As willows grow tall, their canopy cover increases, and they may grow together into thickets, enhancing habitat structure and complexity. However, researchers using artificial beaver dams to study willow growth have suggested that much willow habitat has been affected by stream incision, lowering water tables and causing an alternative stable state in which willows cannot recover, even with reduced herbivory; that herbivory has not been significantly reduced in their study sites; and that the abundance of tall relative to short willows on the northern range has not changed since 1990 (Bilyeu et al. 2008, Marshall et al. 2013. Thus, the extent of tall willow recovery, and the mechanisms driving it, have remained unclear, and there was a need for a broader examination of willow growth in northern YNP. Our research compared data from the early 1990s, early 2000s, and 2016-2018, to investigate how heights of potentially tall willow species on the northern range have changed in the two decades since large carnivore restoration, and to assess the importance of browsing as a driver of willow height and canopy cover.

Study area
The study area included the central and eastern portions of the Yellowstone northern ungulate winter range (Houston 1982, Painter et al. 2015. Study sites ( Fig. 1) were chosen because of their use in previous willow studies from which we obtained data, or because of extensive areas with tall willow species (suppressed or not) or historical tall willow habitat. Study sites were dominated by tall willow species, primarily Geyer (Salix geyeriana), Bebb (S. bebbiana), Booth (S. boothii), and sometimes Drummond (S. drummondii) willows. On the larger streams-Lamar River, lower Slough Creek, and lower Soda Butte Creek-sandbar willow (also called coyote willow, S. exigua or melanopsis and hybrids) was also present. Sandbar willow is different from these other species in that it spreads vigorously through root sprouts (suckering). All of these tall willow species can commonly grow to heights above 400 cm. However, most willows on the northern range were less than 100 cm in height in the 1990s, and there was much concern about the suppression of willows due to intensive browsing by elk, primarily in fall and winter (Singer et al. 1998, NRC 2002. Wolf willow (S. wolfii) was present in some sites but this species typically does not grow taller than 150 cm and was not included in the analysis. Bison were present in high densities year-round in the eastern part of the northern range (Lamar Valley and Junction Butte vicinities) and browsed willows throughout the summer (Painter and Ripple 2012). Elk used the study sites primarily in winter, but near the town of Mammoth a resident herd was present year-round.

Willow data
An opportunity to measure changes in willow heights over time was provided by three previously collected datasets, obtained from the United States Geological Survey (USGS) and YNP: 1. The Singer Transects, measured in 1988-1993, which were established to monitor willows on the northern range (Singer et al. 1994 (Tercek et al. 2010b), hereafter referred to as the "willow surveys." Patches are groups of willows of any height, whereas thickets are groups of tall willows (here defined as 5 or more willows >200 cm spring height and <2 m apart). Thickets are located within willow patches.
Willow sampling took place in two phases: (1) measuring willows on the Singer Transects in July and August of 2016 and 2018, using methods consistent with previous height measurements, and (2) locating and measuring willow thickets in 2017 and 2018; these thicket surveys were compared to the 2003-2010 willow surveys by Tercek et al. (2010b) and used as an index of tall willow recovery.

Singer willow transects
The Singer Transects were located from approximate GPS coordinates, written descriptions, and prior visits to these study sites. These transects were originally designed to be representative of a range of browsing and height conditions (Singer et al. 1994). Our resampling ( Fig. 1) included West Blacktail Creek (one transect), lower Slough Creek (four transects), upper Slough Creek (one transect), Lamar Valley (two transects), Lamar-Soda Butte confluence (one transect), and three transects, called "willow belts," established by park researchers in 1958 in wet meadows near exclosures at Mammoth, Junction Butte, and in Lamar Valley (one transect each). Some of the original Singer Transects were located in headwaters outside the boundary of YNP, but we limited our study to the northern ungulate winter range within the park and did not include the additional transects in our analysis.
At each Singer Transect location, following methods consistent with Singer et al. (1994), we placed a 100-m transect through the centroid of the willow patch, along the long dimension of the patch, and following the curve of the patch or stream bank if necessary to stay in the patch. At intervals of 5 m, the closest willow within 5 m was chosen for sampling. At each sampling position, a canopy plot 2 m in radius was used to index tall willow cover for each site, by recording the number of plot quadrants that contained willow canopy taller than 200 cm. We also measured the height of the browsing line outside the Junction Butte and Mammoth exclosures, recording the lowest unbrowsed stem at intervals of 3 m (Mammoth) or 6 m (Junction Butte) along one side of the fence.
For each sampled willow, we recorded the species, previous spring height (height of the stem node from which new growth began in spring), and summer height (height of the top of new summer growth), for comparison to spring and summer heights in the Singer Transect data. We measured browsing intensity (browsed or not, in the current summer and previous year) on five of the most accessible stems (spaced around the perimeter of the bush) at two different heights: 100 cm or less (easily within the reach of bison) and 120-200 cm (above the reach of most bison but accessible to elk; Painter and Ripple 2012). This provided a measure of the browsing pressure on stems likely to be browsed but is not comparable to the percentage of annual production removed by browsing (browsing offtake) as previously measured on the Singer Transects (Singer et al. 1994, Bilyeu et al. 2007). The two measures are much different, as many stems on a willow bush may be protected inside the bush or above the reach of herbivores, and not easily accessible. Thus, browsing offtake as a measure of browsing pressure may be confounded by differences in plant morphology when comparing plants affected by different browsing regimes. Instead, we based our browsing estimates on accessible stems as an index of browsing pressure in the site.
The Singer Transects were designed to track individual plants marked with ID tags; thus, only the first year of sampling was random. Established plants were suppressed and not growing taller, and the transects were intended to study and track this process over several years. We analyzed the growth over a full year of individual plants in these data to test the hypothesis that plants were suppressed by browsing. In 2016-2018, the original tags could not be found and we initiated new random sampling of willows on the transects.

Willow thickets
We mapped willow thickets in 2017-2018 ( Fig. 1) along portions of West and East Blacktail Creeks, Oxbow, Geode, Elk, Crystal, Slough, and Soda Butte Creeks, and the Lamar River below the Soda Butte confluence to Lamar Canyon (Lamar Valley). The previous (2003-2010) willow surveys (Tercek et al. 2010b) were used to aid in identifying stream reaches with patches of tall willow species. Thickets were defined by Singer et al. (2005) as five or more willows >200 cm tall, and <2 m apart, when these newly emerging clumps of tall willows were observed in the early 2000s as willows began to grow taller (Baril 2009). We required the minimum five willows to be >200 cm in spring height, before the summer growing season, because willow stems can grow 50-100 cm in a season but this growth may be browsed and lost in winter. Our definition rules out identifying a thicket based on transient summer height. Each tall willow defining a thicket was separated from the others by at least 1 m, and two thickets separated by a gap of <4 m were considered part of the same thicket. We marked Fig. 2. Relative frequency height distributions of willows in the Singer Transects in three time periods, for all transects that were remeasured in 2016. In 1988-1993, most willows were short, with a few tall willows making a long tail. By 2001-2004 (n = 182), there was an overall shift toward taller willows, and by 2016-2018 (n = 200), height distribution had become bimodal, with a new peak in the 300-400 cm range. The many short willows in 2016-2018 reflect continued high rates of browsing that suppressed willow growth in some sites. the beginning and end of each thicket with GPS and recorded the maximum height in the thicket.
For both the Singer Transects and the thicket surveys, we recorded the presence of hydrophytic plants such as sedges (Carex spp.) and rushes (Juncus spp.) as indicators of moisture where willows were growing or obligate upland plants such as sagebrush (Artemesia spp.) as indicators of dry conditions. Shrubby cinquefoil (Potentilla fruticosa) was a common transitional species between wet and dry. We also noted the presence or absence in a site of surface water in late summer (the dry season), as well as signs of recent overbank flooding, beaver ponds, or stressed willows in dry sites.

Analysis
We first divided the Singer Transect data into spring and summer groups, as heights in the two seasons are not directly comparable. We then pooled the transect data for each season in the three time periods of data collection : 1988-1993, 2001-2004, and 2016-2018. Preliminary data inspection showed strongly divergent patterns of height and browsing between the willow belt sites in meadow seeps near the three exclosures (Mammoth, Junction Butte, and Lamar), compared to the rest of the study sites which were associated with streams. Therefore, we further separated the data into two groups as meadow and stream sites. Thus, for each season, there was a meadow and a stream dataset with data for the three time periods. Assumptions of normality and homoscedasticity were assessed using histograms (Fig. 2) and Bartlett's test for equal variance. As preliminary analysis showed the data were heteroscedastic, we used analysis of variance without the assumption of equal variance (Welch's ANOVA) to compare mean height between the three time periods within each dataset and used Dunnett's T3 test (R package PMCMRplus) to generate adjusted p values for multiple pairwise comparisons of mean heights (Dunnett 1980). We compared mean browsing offtake percentage between 1988-1993 and 2001-2004 using Welch's t-test for unequal variances, after applying an arcsine of square root transformation to the percentage data.
The Lamar willow belt transect was monitored in 1988-1993, but not 2001-2004. In 2016, it had no willows remaining and the site was severely trampled by bison. We replaced it with a new transect in a willow patch near the Lamar West exclosure to represent this area in the 2016 data but dropped this site from comparisons across the time periods. Therefore, the meadow dataset included only two meadow sites, the Mammoth and Junction Butte transects.
We tested the hypothesis that browsing was the main driver of differences in height increase, by comparing annual versus summer growth increments for all individual plants that were tracked for a complete year in the Singer Transects. For this comparison, we excluded plants >300 cm spring height, as they were too tall to be affected by elk browsing. If plants increased in height in the summer, but then lost the height increase by the next spring, this strongly indicates that browsing removed the annual height increase and suppressed the plant. If height increases were preserved to the next spring, then browsing was not suppressing height growth. We used linear regression to test the hypothesis that summer browsing (the percentage of stems browsed in summer) in 2016-2018 was driving differences in canopy cover (the percentage of canopy plot quadrants in a site with canopy 200 cm tall). To compare summer browsing across the three time periods, we calculated a summer browsing index as the percentage of plants on which any summer browsing was found on sampled stems, a measure that could be obtained from both past and current data.
For willow thickets found in our surveys, we compared maximum height in 2017-2018 to the previous maximum recorded for the corresponding patch in the willow surveys by Tercek et al. (2010b). We calculated (1) the percentage of thicket locations that fell outside of a previously mapped patch, (2) the percentage of patches that could have had a thicket in the previous willow surveys based on maximum height, and (3) the percentage of previously mapped patches in a site that contained a thicket in 2017-2018.

The Singer transects
Of the original 15 Singer Transects, we sampled 12 that were within the park boundary and the winter range (Appendix S1: Table S1, S2). Seven of the 12 transects were in Lamar Valley and lower Slough Creek, so this area of large streams and recent high bison density was strongly represented (Fig. 1). Lower Slough Creek was also the location of most of the very tall willows in 1988-1993 (Fig. 2), and these large, old clumps of willows persisted in places where they were protected from browsing by  (Dunnett 1980) for heteroscedastic data (*indicates comparisons that are significant with a minimum 95% confidence). Stream sites increased in height, while meadow sites did not change, or decreased. steep banks and logjams deposited by floods (L. Painter, personal observation). Mean heights varied significantly within each dataset, with all ANOVA P-values < 0.001 (Tables 1, 2). Stream sites showed significant height increase over time, while the three meadow sites (willow belts) did not, or even decreased (Tables 1, 2; Fig. 3-5). Data distributions of heights in 1988-1993 and 2001-2004 were non-normal but similar, with mostly short willows and a long tail of taller heights. However, the distribution of height overall (Fig. 2), and in stream sites in 2016-2018, was bimodal; thus, it is important to note that a comparison of means does not reflect the actual amount of change in height, as the distribution around the mean also changed (Fig. 4). Despite this increasingly bimodal height distribution, tall willows in stream sites exceeded 200 cm in mean height by 2001-2004 in summer, and by 2016 in spring (Fig. 3a), a remarkable change given the history of willows on the northern range.
Mean winter browsing offtake (percentage of annual production) decreased in stream sites between 1988-1993 and 2001-2004, but not in meadow sites (Fig. 3). Analysis of individual willows tracked for one year indicated a causal connection between browsing and plant height (Fig. 5). The percentage of stems browsed in summer explained 95% of the variation in the percentage of canopy cover taller than 200 cm in  Table 2  2016-2018 (Fig. 6a), and almost all of this relationship was explained by browsing at heights <100 cm, suggesting bison as the browsers. Summer browsing decreased after 1993, but then increased again after 2004 (Fig. 6b). Tall willows in the Lamar Valley area exhibited a browsing line at about 120 cm in height. The mean height of the lowest unbrowsed stem was significantly different (P < 0.001) at the Mammoth exclosure (mean 196 cm, n = 15) compared to Junction Butte (mean 145 cm, n = 20).

Thickets
We surveyed a total of 32 km of stream reaches for thickets, and found 225 thickets ranging in length from 10 to 110 m. Of these thickets, 25% were (in whole or in part) in willow patches not previously mapped. Of the previously mapped patches in the surveyed stream reaches, 37% (range 22-88%) had thickets in 2016. Forty-six thickets (20% of all) were dominated by sandbar willow, and 95% of thickets had max height >300 cm. Of those patches that had thickets in 2017-2018 and had been previously mapped in 2003 (on Blacktail, Geode, Elk, Crystal, and Slough Creeks), 57% had a max summer height >200 cm in 2003, and 39% were >300 cm in 2003. None had mean height >200 cm in 2003.

DISCUSSION
This study used three main approaches to assess willow recovery in northern Yellowstone: (1) comparing height and browsing offtake over time, using the Singer Transect data from 1988 to 1993, 2001 to 2004, and 2016 to 2018; (2) measurements of tall willow canopy cover and browsing intensity in 2016-2018 on these transects; and (3) surveys of tall willow thickets in 2017-2018 along streams in the study area. All three approaches indicated that willows have greatly increased in both height and tall canopy cover since the 1990s, but also that willows in some locations were still suppressed by herbivory. Canopy cover of tall willows in 2016-2018 was strongly associated with summer browsing intensity (Fig. 6a), pointing to herbivory as the primary driver of willow height and cover. Butte transect was representative of the three meadow sites, which decreased from an already short height. (c) Lower Slough Creek had some isolated tall willows in 1988-1993, resulting in a mean higher than other sites; the mean did not increase, but the distribution of height changed (Fig. 2), with many more willows >200 cm in 2016-2018, along with many short and suppressed.
Where tall willows had been only the long tail of a distribution with mostly short heights, a new peak emerged in the 300-400 cm range (Fig. 2), representing new growth of tall willows on the northern range since large carnivore restoration and subsequent reductions in elk browsing.

The Singer transects
This study marks the first time that a causal connection between browsing and willow height has been demonstrated statistically for the northern range, a result made possible by tracking the heights of individual plants on the Singer Transects. Willows in stream sites showed marked increases in height by the early 2000s, while there was little change in meadow sites (Figs. 3, 4). The significant height increases in the stream sites, to an average height >200 cm, were largely the result of reduced herbivory, and not changes in summer productivity, moisture, or other bottom-up factors. This conclusion was suggested by the reduction in browsing offtake relative to the 1990s (Fig. 3b), but, more to the point, a much greater proportion of summer growth remained through the winter (Fig. 5), demonstrating the mechanism by which height increased. In contrast, the lack of height gain in meadow sites was due to continued loss of annual growth to herbivory-the plants grew, but they were eaten. Meanwhile, the tall, dense growth of healthy willows inside the meadow exclosures (Fig. 8b) demonstrated that conditions were suitable for willow growth in the absence of browsing. The importance of browsing as a limitation on height was also emphasized by Barmore (2003) who studied these sites in 1962, as well as Singer et al. (1994) in 1989, and others (Chadde et al. 1988, Kay and Chadde 1992, Singer 1996. Field observations suggested that the lack of recent willow height increases in the meadow sites was primarily due to browsing in summer, but by different ungulate species. Browsing by elk likely suppressed willows at the Mammoth exclosure, while browsing and trampling by bison appear to have caused a similar result at Junction Butte and Lamar Valley (Painter and Ripple 2012). Although browsing by elk has decreased overall on the northern range, the Mammoth willow exclosure was frequented year-round by a local resident herd of elk, as indicated by fresh elk scat and frequent sightings of elk in the vicinity, but few bison. At Junction Butte and Lamar, the situation was different, with bison present in large numbers throughout the summer, and the sites severely trampled. Herds of bison were seen frequently, and bison scat and tracks were common on the transects, while no signs of elk were found. Furthermore, the height of the browsing line supported this inference, as it averaged nearly 200 cm along the fence at the Mammoth exclosure, typical of a browsing line caused by elk, but was significantly lower at Junction Butte (Fig. 8b). Bison have increased in recent years , White et al. 2015, Beschta et al. 2020, removing new growth and keeping short willows even shorter, browsing throughout the growing season as well as in winter.
In areas of the greatest bison concentrations (Lamar Valley, Lower Slough Creek, Soda Butte Creek, and Crystal Creek), we observed a horizontally suppressed condition, where willows had escaped browsing vertically, but were constricted horizontally by bison browsing and horning below about 120 cm. Some tall willows were broken down by horning, when bison twisted and broke larger stems with their horns. We documented some tall willows in Lamar Valley that grew tall in the early 2000s but were completely torn down by bison between 2010 and 2017 (Fig. 8c, d).

Thicket surveys
We found thickets along all streams surveyed in 2017-2018, an important component of habitat diversity and quality that was lacking on the northern range in the 20th century (Houston 1982, Chadde and Kay 1996, Wagner 2006, Baril et al. 2011. The percentage of previously mapped willow patches (Tercek et al. 2010b) that now contain thickets could function as an index of willow recovery (Fig. 7). Some stream reaches had thickets in more than 80% of mapped willow patches (Fig. 8a), suggesting a strong recovery of tall willows, while in other locations, such as Lamar Valley, thickets were rare. Overall, <40% of previously mapped willow patches contained thickets, suggesting that most stream reaches had not reached their potential for tall willow cover in existing willow habitat. Some thickets were very tall and expansive, while others were just emerging. Many (25%) fell outside of previously mapped willow patches, representing an expansion of cover of tall willows. Almost all thickets had willows >300 cm, but such tall willows were uncommon in previous willow surveys. Therefore, it is reasonable to conclude that most patches did not have thickets in 2003 when the surveys occurred, and reports by researchers at the time suggest that thickets were rare. The expansion of thickets is a strong indication that height and cover of willows have increased, as also was found by Baril et al. (2011).

Sandbar willow
An interesting result of the overall decrease in browsing along Slough and Soda Butte Creeks has been the recent expansion of sandbar willow on these large streams. This species has a singlestem growth form and spreads through rhizomes and root sprouts, forming dense patches on sandbars and banks. These willows grow quickly but previously were suppressed by herbivory. Sandbar willow was rarely mentioned in the data and publications related to willows on the northern range prior to 2003. Chadde et al. (1988) reported finding sandbar willows on lower Slough Creek, but noted that heavy browsing often limited height to <100 cm, and "may eventually eliminate this species." Similarly, some patches of short sandbar willow were found on lower Slough Creek and Soda Butte Creek in the willow mapping surveys in 2003 and 2007. Our thicket survey in 2017 found new tall thickets of sandbar willow in these locations and also in Lamar Valley. On these large streams, tall sandbar willows sometimes filled in the space between tall willows of other species. We assessed the growth history of some of the tallest of these sandbar willows and found they had been growing above an older base of browse-killed stems for about 4-6 yr (L. Painter, unpublished data). Thus, sandbar willow recently has become an important species on large streams in the northern range, where prior to the reduction in browsing it was suppressed and insignificant. As sandbar willows grow tall, they spread through root sprouts, creating dense thickets of tall willows of a type never before recorded on the northern range, but common along similar streams in the region.

Stream incision and willow habitat
Some researchers have emphasized the importance of stream incision as a factor limiting willow growth and recovery on Yellowstone's northern range (Wolf et al. 2007, Marshall et al. 2013; however, it is important to distinguish between recovery of tall willows in existing habitat, as we have shown, versus recovery of a greater extent of willow habitat and wetlands, as existed in the early days of the park. We agree that loss of riparian shrubs and beaver wetlands, and the stream incision that followed, has limited the extent of willow habitat to less than it was historically, and these changes may not be quickly reversible (Wolf et al. 2007). However, we found that stream incision was not preventing willow recovery in most existing willow habitat. For example, even though a portion of Elk Creek has become deeply incised (Wolf et al. 2007), lower Elk Creek in Yancey's Hole had shallow banks allowing seasonal flooding, with tall willow thickets emerging along the stream and in the adjacent meadow. Similarly, floodplains along West and East Blacktail Creeks were well-watered by groundwater seeps, with tall willows forming extensive thickets and overhanging the streams in 2016-2018 (Beschta and Ripple 2018). In these stream reaches, 88% of willow patches had thickets (Fig. 7), suggesting an extensive recovery of tall willows.
We found little evidence that existing willow habitat has been strongly affected by stream incision. For example, beavers built a series of dams on West Blacktail Creek in 2012 (L. Painter, Personal observation), flooding banks and further connecting the stream to the floodplain (Beschta and Ripple 2018). Many willow sites we surveyed were wet at the surface due to seeps of groundwater, which can support willows independently of stream incision (Johnston et al. 2011). Large streams-Lamar River, Soda Butte Creek, and Slough Creek-regularly overflow their banks in spring, and willows are supported by groundwater in old oxbow bends and toe-ofslope seeps. Most of our study sites, including the severely suppressed meadow sites, had moisture at or near the surface throughout the growing season, as indicated by the presence of water or the dominance of hydrophytic plants. Bilyeu et al. (2008) found that willows in experimental dry sites did not grow significantly taller after herbivory was removed; however, these dry sites were not representative of most willow habitat, as they had sagebrush growing alongside stunted willows (L. Painter, Personal observation). None of our study sites had willows growing with sagebrush, except perhaps at the transitional edges. In fact, willow recovery, like aspen recovery, began in the early 2000s during a time of drought (McMenamin et al. 2008, Painter et al. 2014) and has continued under a long-term warming and drying trend. Further expansion of willow habitat, on the other hand, will require beavers to flood wider areas (Marshall et al. 2013), restoring ecological services provided by beaver ponds and their associated wetlands in stream systems of the Rocky Mountains.

Reductions in browsing
Studies using experimental dams and exclosures on West Blacktail, East Blacktail, and Elk Creek have suggested that browsing offtake as a percentage of total annual production has not decreased significantly, based on measurements that began in 2003 (Bilyeu et al. 2008, Hobbs and). This appears to contradict our conclusion that browsing has been reduced, at least in some sites. One of these three experimental locations (Elk Creek) is on a trail used by ungulates (L. Painter, Personal observation), so it is not surprising that browsing has remained intensive there, but the others are in an area (Blacktail creeks) where tall willows have been recovering, as shown by our results. There are good reasons to think that browsing overall has decreased significantly in the northern range since the 1990s, including (1) our analysis showing a reduction in browsing offtake on the Singer Transects between 1988-1993and 2001-2004 (Fig. 3), using the same method employed at the experimental sites; (2) our analysis showing an increase in the amount of summer height gain retained over the next winter, demonstrating a biologically significant decrease in browsing (Fig. 5); (3) documented decreases in browsing pressure on aspen during the study period (Painter et al. 2015, Klaptosky 2016); (4) reductions in browsing and associated increases in willow growth documented by other studies (Singer et al. 2003, Beschta and Ripple 2007, Beyer et al. 2007, Beschta and Ripple 2018; and (5) the substantial, unprecedented, and sustained reduction in elk density in the northern range within the park since the 1990s (Painter et al. 2015(Painter et al. , 2018. How then to explain the lack of a significant decrease in browsing at the three experimental sites? One explanation is that average browsing offtake had already decreased by 2003, and willows were already beginning to grow taller by that time on some sites, as indicated by our results and these previous studies. Another consideration is that browsing offtake as a percentage of annual production may be a weak indicator of the change from an intensive to a moderate browsing regime, because plant morphology changes during the process of recovery. Suppressed plants have a hedged growth form, with much of the annual growth under the protection of dead stems. Even if all of the accessible growth is eaten annually by herbivores, much (perhaps 40%, for example) is protected by dead stems and will not be eaten. When browsing pressure is reduced, willows grow beyond the protection of dead stems and most annual production becomes accessible. If 60% of these taller, more accessible stems are eaten, that leaves 40% to grow taller. The browsing offtake percentage would be the same in both of these examples; however, one represents an intensive browsing regime with hedged, suppressed plants and the other a moderate browsing regime where plants grow taller. Thus, the ability of the browsing offtake method to detect changes in browsing intensity may be partially confounded by changes in plant morphology.
Willows began growing taller prior to 2005 when elk densities were reduced but still relatively high, suggesting that behavioral effects of predators may have played a role in reducing browsing (Singer et al. 2003, Ripple and Beschta 2006, Beschta and Ripple 2007, Beyer et al. 2007). Studies using tracking collars found that elk avoided riparian areas during times when wolves were active, suggesting a mechanism by which behavioral responses could reduce browsing of willows (Beyer 2006). Elk are adept at working around predators by avoiding risky places during risky times (Kohl et al. 2018), but responses to risk may change the amount of time spent foraging in those places, adding to the direct effects of predation. These behavioral studies, as well as studies of aspen, suggest that with high elk densities most plants will be browsed; however, with lower elk densities, behavioral responses to risk could reduce browsing in some places (White et al. 2003) and increase it in others (Creel and Christianson 2009), resulting in a patchy plant response.

CONCLUSIONS
Our results indicate that willows began to recover as elk densities decreased in the early 2000s (Figs. 2-4a), that height increases were due to a reduction in herbivory (Fig. 5), and also that some sites continued to be heavily browsed and have not recovered at all. These conclusions parallel those reported for aspen (Painter et al. 2015, Klaptosky 2016. Much recent herbivory of willows in some locations can be attributed to bison, based on the height and season of browsing (Fig. 6;Painter and Ripple 2012), and the presence of herds of bison. Since the late 1990s, bison have increased greatly in number on the northern range and are present year-round (Gates et al. 2005, Painter et al. 2015, while elk densities have decreased to historical lows. As with aspen, the exceptions to a trend of increased willow height and cover are often in places of highest bison density, particularly Lamar Valley. Even there, some willow thickets have emerged (Fig. 7), but many willow patches in the valley remain suppressed, and some tall willows have been torn down by horning (Fig. 8c,d). In some places such as Crystal Creek, the combination of beavers cutting tall willows and bison browsing them could cause retrogression of willows that recently grew tall when elk browsing was reduced (Baker 2003).
This analysis fills a gap in the story of the Yellowstone northern range ecosystem, putting exclosures and experimental sites into a landscape context. Most stream systems in the northern range have shown some willow recovery, with new thickets emerging, and beavers occupying some streams. Our results agree with those of researchers in the 1990s, who concluded that browsing was the cause of suppressed willow heights. Where browsing has been moderated, willows have begun to recover, creating a more normal and functional riparian ecosystem. The most important factor currently limiting recovery in existing willow habitat is the high densities of bison in portions of the northern range. While a drying climate and incised streams may slow or prevent willows and beavers from reoccupying some sites they formerly dominated (Johnston et al. 2011, Marshall et al. 2013, the general pattern in northern Yellowstone has been toward a moderate browsing regime that supports tall willows and development of willow thickets. This is a profound change in the ecology of the northern range, and a break from the past in which high elk densities prevented tall willow thickets. These new willow thickets will stabilize banks, shade streams, and provide habitat for birds, invertebrates, and beavers, increasing biodiversity and ecosystem resilience into the future.