Long-term patterns of mass stranding of the colonial cnidarian Velella velella: influence of environmental forcing

Velella velella is a pleustonic cnidarian noted worldwide for mass stranding of the colonial phase. Utilizing a 20 yr dataset (2000−2019; 23 265 surveys) collected by the COASST citizen science program, we examined the spatio-temporal occurrence of mass strandings of V. velella along the Pacific Northwest coast from Washington to northern California, USA. V. velella mass strandings were documented in 14 years, with expansive events in 2003−2006 and 2014− 2019. Events predominantly occurred in spring and were synchronous (April) among years, concurrent with shifts to prevailing onshore winds. Autumn mass stranding events occurred infrequently, with no consistent phenology (2005: November; 2014: August). In stranding years, reports of V. velella were mostly synchronous throughout the surveyed area, and events consistently spanned >400 km of coastline, with highest reporting rates in the vicinity of the Columbia River plume, collectively suggesting extensive V. velella blooms throughout the northern California Current system in some years. Annual metrics of spring V. velella reporting rate (proportion of beaches; January−June) were modeled as a function of indices representing sea surface temperature anomaly (SSTa), easterly (onshore) wind speed, and regional upwelling. The best models (based on Akaike’s information criterion corrected for small sample size) indicated that SSTa averaged over the preceding winter (December−February) was positively correlated with spring reporting rate, suggesting that mass strandings of V. velella may be more prevalent in warmer years. As planetary warming continues, and V. velella strandings are easily recorded by citizen science programs globally, we suggest that stranding prevalence may be one relatively easy measure providing evidence for epipelagic ecosystem response.


INTRODUCTION
Increasingly dramatic responses of natural systems to human-assisted top-down and bottom-up forcing are evidenced by massive shifts in the distribution (Perry et al. 2005, Cheung et al. 2012, abundance (Last et al. 2011), population health (Burge et al. 2014), biodiversity (Jones & Cheung 2015), reproductive output (Foo & Byrne 2017), and competitive viability (Sumaila et al. 2011) of marine species. Gelatinous zooplankton, including the larger scyphozoan jellyfish, have received attention worldwide due to an apparent increase in documented population explosions or blooms. Explanations vary, and include decadal oscillations which may be forced by climate (Condon et al. 2013) and/or specifically a warming ocean (Lucas et al. 2014); release from top-down and competitive forcing by global overfishing (e.g. Daskalov et al. 2007 but see Opdal et al. 2019); and/or anthropogenic degradation of nearshore systems (Purcell 2012). In addition to recorded blooms across a wide range of large marine ecosystems, species present in coastal systems, or transported to coastal systems, have also been recorded in mass beaching or landfall events (Purcell et al. 2007, Betti et al. 2019. With respect to gelatinous taxa, blooms of Velella velella have been recorded worldwide across temperate and tropical systems, frequently evidenced by mass stranding events. V. velella, hereafter Velella, is a cosmopolitan, tropical to temperate pleustonic species (Pires et al. 2018) often found in association with coastal fronts (Betti et al. 2019), Langmuir cells, and other convergences . A cnidarian in the family Porpitidae, the Velella life cycle phase most often encountered are colonies of asexual polyps, small pleustonic organisms with a ring of purple tentacles and a stiff sail-like dorsal extension (maximum disk size = 80 mm; Bieri 1977). Colonies produce pelagic medusae (the sexual phase) which sink 600−1000 m before releasing pelagic larvae which develop as they ascend to the surface, culminating in the pleustonic colonial phase. Colonies are thought to be produced over deep water (Purcell et al. 2015). Bieri (1977) suggested 2 full life cycles annually, with colonies extant in spring and autumn (see also McGrath et al. 1994, Pires et al. 2018, al though Evans (1986) noted that in the North Pacific, shipbased records of Velella were restricted to the spring and early summer (late April through very early July).
Unusually large beaching events or strandings of millions of Velella have been reported worldwide (Portugal: Pires et al. 2018;Mediterranean: Betti et al. 2019;Oregon: Kemp 1986;Ireland: McGrath et al. 1994;South America: Carrera et al. 2019;Pacific: Purcell et al. 2015). Mass strandings have been associated with warm waters (Pires et al. 2018, Zeman et al. 2018 including during El Niño events (Carrera et al. 2019), as well as with onshore winds (Pires et al. 2018, Betti et al. 2019). In part because Velella mass beaching events are highly visible, public engagement in reporting them is high. Several programs worldwide, including Jellyfish Watch (Purcell et al. 2015) and GelAvista (Pires et al. 2018) have specifically recruited observers to report Velella beachings, while other studies have incorporated public reports into a wider dataset (e.g. McGrath et al. 1994).
Observations collected by citizen science programs that employ a standard data collection protocol can provide a means of collating information over broad spatial and temporal scales in a cost-effective manner. Although marine programs are scarce relative to terrestrial ones (Theobald et al. 2015), several biodiversity monitoring programs collect data at the scale of large marine ecosystems on marine mammals (Moore et al. 2009), marine birds (Jones et al. 2018), finfish (Thorson et al. 2014), benthic marine invertebrates (Schultz et al. 2016), and gelatinous zooplankton (Purcell et al. 2015). Collectively, these programs may detect and monitor, for example, species range extensions (Champion et al. 2018), invasions (Scyphers et al. 2015), harmful algal blooms and their effects (Trainer & Hardy 2015), fishery impacts (Hamel et al. 2009), pollution (Keil et al. 2011), and population responses to a shifting climate (Parrish et al. 2007). Identification and documentation of ways in which data could be collected by hundreds to hundreds of thousands of people could greatly aid in the creation of rigorous longitudinal datasets needed to document cyclic, catastrophic, and chronic ecosystem change.
In this paper, we used data from a citizen science program, the Coastal Observation and Seabird Survey Team (COASST; www.coasst.org) to document mass strandings of Velella along the northern half of the California Current large marine ecosystem (CCLME). Although COASST is primarily focused on marine birds (Parrish et al. 2017) and marine debris, participants are invited to record unusual sightings of any nature during their monthly surveys in a 'Comments' box located on their datasheet. Mass strandings of by-the-wind sailors, as Velella are commonly known, are one example of an unusual sighting of which participants are shown a photograph during their training. Because COASST covers ~200 beach sites annually in the northern CCLME, spread from Cape Flattery at the tip of the Olympic Peninsula in Washington, to Elk, California (39.12°− 48.34°N: Fig. 1), we used these data to examine the phenology and spatial distribution of mass stranding events across the ecosystem, as well as to model annual patterns in reporting rate as a function of environmental forcing. Given previous literature, we hypothesized that mass strandings of Velella would be associated with short-term onshore wind events, strength of upwelling, and a warming ocean.

Data collection
COASST is a citizen science program in which expert-trained observers monitor a particular length of coastline on a monthly basis, searching for beachcast bird carcasses or marine debris. Although survey sites are different lengths, all are fixed in space (i.e. permanently marked edges along the long axis of the beach) such that the exact same section of coastline can be repeatedly monitored for years (see Parrish et al. 2017 for an extensive description of the program, and Parrish et al. 2019 for a description of the participants). For this study, we used all surveys from 2000 through 2019 conducted along the outer coast of Washington, Oregon, and northern California (N = 23 265; Fig. 1). This spatio-temporal range encompasses a total of 293 survey sites, with annual coverage varying from 10 sites in 2000 to 158 sites in 2019, and with >100 sites monitored annually since 2006 when the COASST program expanded into northern California. The majority of survey sites are >1 km in length (75%; median length = 1.2 km, range = 0.2−4 km). In addition to their primary data, participants are encouraged to report abnormal occurrences, including but not limited to: unusual singular events such as stranded marine mammals, Humboldt squid Dosidicus gigas, sharks, or finfish; and non-bird mass stranding events. All abnormal occurrence data are recorded in a 'Comments' (text string) field, and may be accompanied by photographs. Because abnormal occurrence data are not required, their presence in the COASST dataset represents certain occurrence whereas their lack does not necessarily represent absence. We note that as reporting of Velella is not mandatory, or subject to a rigorous search protocol, low-level deposition may be missed or not documented in survey comments. Consequently, our presence data likely represent colony presence in quantities above some unknown threshold. Therefore, our analyses should be interpreted within the context of massive stranding events, as our data may not be representative of shortlived/less abundant deposition.
To create a quantitative dataset, all survey comments were searched for character strings associated with Velella sightings (using the fuzzy matching search 'grep' in R). We searched for strings representative of alternate spellings of Velella (e.g. 'velel', 'vellel', 'velell', 'valel', 'vallel', 'vella'), as well as 'sailor', based on the common name 'by-the-wind sailor'. Matches were screened to verify reports were of genuine Velella sightings.
As our base data derive from opportunistic reports, we explored the possibility that variable rates of Velella reporting among years may be influenced by changes in the participant pool (i.e. due to recruitment/retirement of program participants). Specifically, we compared reporting rates in April, the month with the highest overall Velella reporting rate, among paired years (constrained to 5 years before and after the focal year) and holding participant pool constant. These calculations resulted in similar reporting rates to the overall (total participant population within the focal year) reporting rate throughout. Most notably, the transition from years with no reports of Velella to those with numerous reports, and vice versa, were maintained, indicating that overall reporting rate is minimally affected by changes in participant population (see Table S1 in the Supplement at www. int-res. com/ articles/ suppl/ m662 p069 _ supp .pdf).
In sum, our basic data consisted of survey date, beach location (lat, long), and Velella report presence/presumed absence. We split our analyses into: (1) determination of spatio-temporal patterns in Velella reporting rates; and (2) examination of relationships between environmental forcing factors and Velella reporting rates expressed on an annual scale.

Spatial pattern
To represent spatial patterns in Velella reporting, we calculated beach-specific reporting rates (number of Velella reports out of the total number of sur- study area, with major landmarks mentioned in the text and reference locations for wind (diamonds) and upwelling (circle) data highlighted. Light blue diamonds in panel B indicate alternate locations for which wind data were obtained from the North American Regional Reanalysis dataset to test for spatial variability in wind speed and direction veys for each beach) from our base survey data (i.e. 23 265 surveys from 293 beach sites spanning 2000 through 2019). To visualize larger-scale spatial patterns across our study system, we applied a kernel smoothing function (Gaussian, σ = 0.45° latitude, or 50 km) to beach reporting rate as a function of latitude. Due to apparent differences in reporting before (overall fewer reports) and after (extensive reports on an almost annual basis) the start of the Northeast Pacific marine heatwave (2014−2016), the smoothing function was applied to beach-specific reporting rates calculated across all surveys pre-2014 and post-2014 separately, to investigate whether the spatial distribution in reporting rate differed between these 2 periods.

Spatio-temporal pattern
To examine the extent and duration of mass stranding events through time, we processed our base survey data to provide presence/assumed absence of Velella reports each month from 2000 through 2019 within 50 km latitudinal segments of coastline (segment length = 0.45° latitude; 40.35°−48.5°N). Month by 50 km segments with fewer than 4 surveys (i.e. <1 survey per week on average within that segment), as well as those with no surveys, are differentiated from the remainder (i.e. month by 50 km segments with ≥ 4 surveys) to indicate time-locations with limited to no survey information, respectively.

Seasonal pattern
Although each beach is only surveyed on average once per month, there is sufficient coverage to examine patterns at a finer temporal resolution (i.e. 14 d) due to the asynchrony of surveys (i.e. surveys are not all performed on the same day each month, but are spread out dependent on participant preference) and the number and extent of surveyed beaches. Therefore, we split each year into 14 d windows (n = 26) based on day of year, and calculated window-specific reporting rate (number of Velella reports as a percentage of the total number of surveys performed within that time period) by year for the entire coastline. We performed these calculations separately for Years with low survey coverage, defined by 14 d windows with fewer than 5 surveys, were omitted from this examination of seasonality due to potential unreliability of inferred reporting rate at lower sample sizes.

Short-term effects of wind
Given that Velella beachings have been documented following periods of onshore winds (Pires et al. 2018, Betti et al. 2019, we examined the timing of Velella reports north of 44°N (i.e. our longest dataset) in relation to prevailing wind forcing. Wind data were obtained from the North American Regional Reanalysis (NARR) dataset, from 2002−2019, consisting of 3-hourly measures of easterly (zonal, u) and northerly (meridional, v) winds on a grid (0.3° ≈ 32 km) arrayed across North America (Mesinger et al. 2006). We extracted data from a grid point (46.20°N, 124.64°W; Fig. 1) situated ~45 km offshore and near the center of the latitudinal range encompassed by our survey data (44°−48.4°N). For each year, we applied a kernel smoothing function (Gaussian, σ = 7 d) to the easterly wind speed in order to provide a longer-term representation of onshore wind propensity throughout the year. We also calculated a 14 d moving proportion of easterly wind observations which represented onshore winds (u > 0), and marked intervals where this proportion exceeded 0.6 (i.e. the start to end of all contiguous 14 d windows) as a measure of predominantly onshore wind forcing. Based on literature describing spring and autumn blooms (e.g. Bieri 1977, Pires et al. 2018, we split the year into 2 halves (spring: January to June, autumn: July to December; offset by 1 mo relative to the bloom cycles presented by Bieri 1977 andPires et al. 2018 to account for the COASST monthly survey cycle) and overlaid the earliest reports of Velella strandings, as well as the 14 d period in which reporting rate was highest, for each half-year period. Although wind-speed information was taken from a single location (46.20°N, 124.64°W; Fig. 1), timeseries of wind-speed components were highly correlated (Pearson's correlation coefficient: 0.66−0.99; Table S2) between our chosen location and 3 other grid points at alternate latitudes (45−48°N; distance offshore maintained at 45−55 km; Fig. 1), suggesting that patterns of prevailing wind-speed and direction are broadly similar throughout the study region (Table S2).

Interannual variation in reporting rates
To identify whether Velella mass strandings were related to interannual variability in environmental forcing factors, we processed our base data for the northern region (north of 44°N) from 2002 onwards to provide annual measures of Velella reporting rate for spring and autumn periods defined as above, as a proxy for prevalence and extent (more beaches reporting Velella indicative of more extensive/continuous mass stranding events) of mass stranding events. For each half-year, we found the number of beaches where Velella was reported, as well as the total number of beaches surveyed within that period, which enabled us to calculate reporting rate as the proportion of surveyed beaches that reported Velella. We use beaches rather than surveys to calculate reporting rate to avoid double reporting (i.e. participants may report 're-found' Velella in a subsequent survey), and to more closely reflect spatial prevalence and extent of stranding events (survey-and beachbased reporting rates were tightly correlated among years; Table S3). Although we present annual reporting rates for spring and autumn periods, we restrict our subsequent analyses described in the following section to the spring period, as autumn strandings were only reported in 5 of the 18 years examined (Table S3).

Modeling environmental forcing
To examine the influence of environmental forcing factors on Velella beaching during the spring time period, we constructed a series of models linking indices of specific environmental conditions thought to influence the extent and/or intensity of Velella mass strandings (i.e. those extant prior to beaching) to spring annual reporting rates calculated across beaches north of 44°N (see Section 2.3.2). Given that spring blooms are thought to be associated with a growth cycle lasting from December (prior year) to May (e.g. Bieri 1977, Pires et al. 2018, and that earliest reports (mid-March to mid-April) and peak reporting rates during our spring window occurred prior to the end of April (Table S4), we constrained our consideration of environmental conditions to the December to April window to represent environmental conditions during the growth phase, and prior to landfall, respectively. Indices were constructed as follows: Temperature. Sea surface temperature anomaly (SSTa) data were obtained from the NOAA OI SST V2 High Resolution Dataset (global, daily, resolution = 0.25°; Reynolds et al. 2002) for 2002− 2019. Temperature anomalies were extracted for December to March of each year for all SSTa grid cells within 200 km of the coastline from 44° to 48.4°N, matching the extent of coastline represented in beach surveys. As initial stranding reports predominantly occurred in late March (Table S4), we excluded April SST values, as they would represent oceanic conditions following landfall in most years. Annual temperature indices were created by averaging SSTa data for 2 mo (n = 3; December−January, January−February, February− March), and 3 mo (n = 2; December−February, January−March) windows spanning December to March for each year, representing alternate timing, and time-scales, of importance.
Upwelling. Daily upwelling data were downloaded from the NOAA Southwest Fisheries Science Center: Environmental Research Division data portal (https://oceanview.pfeg.noaa.gov/products/upwellin g/dnld) between 2002 and 2019, for a location offshore of central Oregon (45°N, 125°W; Fig. 1B). The Bakun index (Bakun 1973) was used, representing the daily average of wind-driven cross-shore transport computed from 6-hourly surface pressure fields. While this measure does not represent realized upwelling, it does provide a measure of upwelling propensity off the coast of Washington and Oregon (García-Reyes et al. 2014). We created 3 indices of upwelling: (1) average upwelling intensity (UI), (2) prevalence of positive upwelling calculated as the proportion of time with positive upwelling (UI > 0), and (3) cumulative positive upwelling (Σu>0UI). Each of these indices was created for the same time periods as described for SSTa.
Wind. Wind data from the NARR dataset (see Section 2.3.1) were processed as above to obtain annualized measures (2002−2019) representative of prevailing easterly wind conditions. In order to represent alternate features of wind forcing, we created 3 indices: (1) average easterly wind speed, (2) prevalence of onshore winds calculated as the proportion of time where winds were onshore (u > 0), and (3) cumulative onshore wind speeds (Σu>0u). Each of these indices was calculated for the March, April, and March−April windows, respectively, as we were only concerned with onshore winds during the period immediately prior to the predominant landfall window (Table S4). We calculated the same indices for 3 points in space (Fig. 1), and the resulting indices were highly correlated among locations (Table S5), suggesting that results would be unaffected by choice of location.
Annual reporting rates of Velella were analyzed using generalized additive mixed models (GAMMs, fitted in R via the 'gamm4' package; Wood & Scheipl 2020) to account for potentially non-linear relationships between reporting rates and environmental factors. Our response of annual re porting rate is represented as a binomial variable, where the number of trials was given as the total number of beaches surveyed and the number of successes was the number of beaches reporting Velella. Each model contained a random effect of survey year to account for pseudo-replication/non-independence within our response data (i.e. among beaches within year; Kéry & Royle 2015).
As we had alternate metrics for each of the 3 environmental forcing factors (SSTa: n = 5; Upwelling: n = 15; Wind: n = 9) we firstly constructed all possible single ( where each forcing factor was represented in a specific model by only 1 (or none) of the candidate indices. For each model, we calculated Akaike's information criterion corrected for small sample size (AICc) as a measure of model fit (Burnham & Anderson 2002). We also calculated the maximum variance inflation factor (VIF) among model predictors (VIF = 1 / 1 -R i 2 , where R i 2 is the coefficient of determination of predictor i, as a function of all other predictors contained within the model) for each model as a measure of multicollinearity, and excluded models with max VIF > 2.5 to avoid the selection of models with collinear predictors (Akinwande et al. 2015). Following this step, we calculated Akaike weights for each model (WAICc = e -ΔAICc/2 ; where ΔAICc was measured relative to the best overall model), which represents the likelihood of a particular model being the best model given the data and the candidate model set (Burnham & Anderson 2002, Wagenmakers & Farrell 2004. As a measure of predictor importance, we then calculated the summed Akaike weight for each predictor across models in which it was included. The proximal index for each environmental forcing factor was then selected as the index that had the highest summed Akaike weight across the set of alternate indices for each environmental forcing factor. Models containing only those proximal indices for each environmental forcing factor were then compared based on AICc to identify the best overall model of reporting rate. In addition, we identified the best overall single-predictor models for each environmental forcing factor (i.e. among alternate indices) to examine the degree to which any individual factor (i.e. wind, temperature, upwelling) could explain interannual differences in reporting rate. For a summary of the spatial and temporal extent of each analysis, refer to Table S6.

RESULTS
Our searches returned 475 matches for Velella search terms, or 2% of all surveys. Of these, 465 were genuine reports of Velella occurrence. Velella was reported on 53% of all beaches within the study region (n = 293).

Spatio-temporal patterns of Velella reporting
Stranding events were reported in 14 of 20 years. Interannual patterns in Velella reporting rate appeared to show stanzas of mass strandings (e.g. 2003−2007, 2015−2019) interspersed by years where Velella mass strandings were apparently absent (2002, 2008−2014), a pattern particularly evident during spring (Fig. 2). For the region north of 44°N where we have the most complete data, spring strandings were typically first reported from mid-March to mid-April (Fig. S1, Table S4). Autumn stranding was less common and more temporally variable than spring, with notable events in 2005 (December) and 2014 (August) (Fig. S1).
The majority of spring events (2003−2006, 2015− 2018) were extensive, spanning > 400 km of coastline, and occurred synchronously up and down the coast (Fig. 3). In 2014 and 2015, Velella reports extended to the southern edge of COASST coverage at Cape Mendocino in northern California, essentially covering the entire 900 km latitudinal extent of the survey region (Fig. 3). Insufficient spatial coverage in earlier years precludes understanding whether mass stranding events in 2003−2006 also extended south of central Oregon (Fig. 3). Within the stranding years, reporting rates were generally higher north of Cape Blanco (Fig. 4). The highest reporting rates occurred within the region of the Columbia River mouth, broken only by beaches within the immediate vicinity of the river mouth (Fig. 4A), although this pattern appears to be driven by data from 2014 onwards (Fig. 4B).

Short-term effects of wind
During winter (November to February), prevailing winds within our study region tended to be directed north, taking on more of an onshore component (i.e. directed east) in March to April (Fig. S2). From April onwards, the prevailing winds were directed south to southeast, a pattern most clearly demonstrated throughout the summer months (Fig. S2). Based on previous literature (e.g. Pires et al. 2018, Betti et al. 2019, we expected Velella transport to be directly wind driven, thus onshore during/following the spring transition even though Ekman transport would be offshore during the spring and summer. Spring Velella strandings were first reported following the transition from prevailing south-easterlies in winter (wind-driven travel = offshore) to north-westerlies in spring (wind-driven travel = onshore; Fig. 5), supporting the hypothesis that strandings are directly tied to onshore transport conditions.

Environmental forcing associated with Velella mass strandings
When models were constricted to only a single predictor, the best metric for each of the environmental forcing factors included average winter SSTa (December− February), spring onshore wind prevalence (March−April), and positive upwelling prevalence during winter (December−January; Table S7). Overall, winter SST was the best predictor of Velella reporting rate during the ensuing spring, and the best overall model contained average winter SSTa (December−February) and late spring onshore wind prevalence (April) as predictors (Table 1). Winter SSTa was included in all of the highest ranked models, and was the only predictor included in the sec-ond ranked model, which was almost equivalent to the best model when judged on AICc (Table 1). By contrast, upwelling metrics were only included in lower ranked models (Table 1).
Fitted relationships for the best overall multiplepredictor model ( Table 1) were suggestive of a positive relationship between Velella reporting rate and winter SSTa (Fig. 6). For winters where average SSTa was < 0, annual springtime reporting rates were at or near 0, whereas warmer winters (SSTa > 0) displayed a step-like transition in Velella occurrence (Fig. 6A), a pattern accentuated in the single-predictor model (Fig. 6C). The best model also suggested that lower spring onshore wind prevalence was associated with higher springtime Velella reporting rates when assessed among years (Fig. 6B), a pattern which became modal in the single-predictor model (Fig. 6D). The fitted relationship for upwelling suggested a negative relationship between Velella re porting rate and upwelling prevalence in the preceding winter months (Fig. 6E). However, models con taining either upwelling or wind as the only predictor of reporting rate were considerably worse than models containing winter SSTa (wind: ΔAICc = 5.9, upwelling: ΔAICc = 8.0; Table 1), suggesting that prevailing patterns are best described by winter SSTa (Table 1). Because Velella reporting rate and average December−February SSTa were both ex treme in the spring of 2015, we re-ran the model selection and fitting process excluding this year. Resultant models (Tables S7 & S8) and fitted relationships ( Fig. 6) were not manifestly different.
Not surprisingly, winter SST, upwelling, and spring onshore winds were collinear, as indicated by VIF values (Table 1). In particular, winter SSTa (December− February) was negatively correlated (assessed via Pearson's correlation coefficient) with positive upwelling prevalence in December−January (ρ = −0.69) and with onshore wind prevalence in spring (March− April: ρ = −0.60; April: ρ = −0.27; Fig. S3). As such, disentangling the degree to which reporting rates were affected by a combination of all 3 environmental forcing factors, beyond the effect of winter SSTa, is problematic.

DISCUSSION
Our data demonstrate that Velella beaching can be extensive in the northern CCLME, arrayed along the coastline for >1000 km in some events. Beachings are associated with predictable annual shifts in onshore wind direction (Fig. 5), pushing massive colony aggregations to shore synchronously throughout our study system (Fig. 3). The occurrence of Velella mass  This model had a VIF max value exceeding the 2.5 cut-off and is excluded from the calculation of Akaike weight due to multicollinearity among included predictors Table 1. Model selection table for generalized additive models of spring Velella reporting rate. Multiple-predictor models were compared among all permutations of models constructed including average December to February sea surface temperature anomaly (SSTa), onshore wind speed prevalence in April, and positive upwelling prevalence in February to March, which were identified as the best representations of each environmental forcing factor based on summed Akaike weight. Best possible models consisting of only a single predictor of SSTa, wind, or upwelling are given in the latter half of the table. For each part of the table (multiple, single) ΔAICc (where AICc is Akaike's information criterion corrected for small sample size) is given relative to the best possible model in that set, and WAICc is the Akaike weight (WAICc = e -ΔAICc/2 ) as a measure of the evidence in support of that model being the best model given the data and the candidate model set. VIF max is the maximum variance inflation factor statistic calculated among predictors included in that model as a measure of model multicollinearity stranding events in spring appears to be related to warmer winter water temperature (Table 1), with stranding rate appearing to increase as a step function of SSTa (i.e. > 0°C; Fig. 6A,C). Mass strandings of Velella may also be influenced by local oceanographic features (Fig. 4), herein most prominently demonstrated by the elevation of reporting in the region of the Columbia River plume, and potentially by the diminution of reporting south of Cape Blanco, a known upwelling domain boundary (Huyer et al. 2005). Finally, although some studies (e.g. Pires et al. 2018) have hypothesized that years in which young colonies receive a nutrient boost may result in greater blooms, we found no evidence that upwelling strength was associated with an increase in Velella reporting rate after controlling for the effect of temperature (Table 1).

Stranding events
The COASST data provide a window into the geographic scale of Velella strandings, with re -peated annual strandings of over 500 km in spatial extent (Fig. 3). In fact, whether Velella strandings extended north and south of our study system is an open question. McGrath et al. (1994) also noted that Velella strandings can occur almost simultaneously over large stretches of coastline, in their case a single stranding event over at least 400 km of the Irish coastline in 1992. Other studies have noted smaller spatial extents (Pires et al. 2018, Betti et al. 2019, mostly due to sampling constraints. It is notable that while our data clearly suggest a prominent spring bloom, autumn beachings were much rarer and occurred at different times (2005: December; 2014: August) in the years in which they were reported (Fig. S1). Although both spring and autumn strandings of mature colonies have been reported (e.g. Bieri 1977, Purcell et al. 2015, Purcell et al. (2012) suggested that 2 generations per year, as proposed by Bieri (1977), is unlikely in cooler water systems where declining autumn−winter temperatures, lower production, and increasingly rough seas likely curtail second-generation growth and survival. 78 Fig. 6. Spring Velella reporting rate as a function of physical forcing factors, with fitted GAM predictions (mean ± 95% confidence interval) of reporting rate. Rows denote differences in model construction called out in headings; columns denote different predictors, called out in the upper right of each panel. Best overall model: Reporting rate ~ f(SSTa + Wind, where SSTa: sea surface temperature anomaly), with fitted regression lines between predicted reporting rate and (A) winter SSTa (proportion onshore winds held constant at 0.66) and (B) spring onshore winds (SSTa held constant at 0) for models with and without the inclusion of 2015. Single-predictor models: panels show the best fitting single-predictor models for each of (C) SSTa, (D) onshore winds, and (E) upwelling for models with and without the inclusion of 2015 The paucity of autumn data on Velella strandings in our study system prevented the use of modeling, precluding definitive conclusions about the influence of temperature or wind. However, both years of significant autumn Velella reporting (2005,2014) were also climatologically anomalous, with delayed coastal upwelling in the northern CCLME in 2005 (Schwing et al. 2006, Barth et al. 2007, Parrish et al. 2007, and the formation of the northeast Pacific Marine Heatwave in winter 2013/14 (Bond et al. 2015, Jones et al. 2018. Zeman et al. (2018) noted that the availability of fish eggs as a high-quality prey resource during winter/spring may explain the apparent disparity in Velella abundance between spring and autumn peaks in general, and that anomalous climate conditions, such as in 2014/15, may lead to prolonged availability of fish eggs due to alterations in spawning phenology and duration. The extension, or alteration, of fish spawning behavior during anomalous years provides one possible explanation for these sporadic autumn mass strandings, but other explanations related to these climatological anomalies (i.e. concurrent alterations to the physical environment) cannot be ruled out.

Role of wind
Our results clearly indicate that wind direction plays a crucial role in Velella beaching (Fig. 5). Given the pleustonic nature of Velella colonies, this is not surprising. Previous studies have also identified local wind speed and direction as important determinants of Velella strandings. Betti et al. (2019) found that spring strandings along the Ligurian coastline were associated with onshore winds. Pires et al. (2018) noted the importance of wind-driven surface waters in determining the shoreward distribution of Velella, and pointed to upwelling relaxation events as the primary factor influencing the immediate occurrence of Velella on Portuguese beaches, as colonies aggregated along the coastal front were transported inshore.
Our data indicate that spring reports of Velella stranding occurred following the transition to prevailing onshore winds, which typically occurred in March (Fig. 5). However, this transition was observed in all years, including those without reports of Velella stranding. Taken together, these results suggest that onshore winds are a necessary but not sufficient precondition for Velella stranding. Our modeling exercise revealed a negative (Fig. 6B) or modal (Fig. 6D) influence of onshore wind prevalence, with years of greater prevalence in March−April associated with few to no reports of Velella strandings. Although this may seem contradictory, onshore wind prevalence values included in our models were greater than 0.5 across all years (Fig. 6B,D), indicating that onshore transport conditions are present each year, as concluded above. Whether strandings occur may thus be governed by bloom formation, size, and/or persistence prior to the transport window. Our limited autumn data suggest that Velella strandings also occur during/following periods of onshore winds, but similar to spring, such conditions were evident across all years (Fig. 5), suggesting that autumn stranding events may also be related to factors other than solely favorable wind direction.

Role of temperature
Our data suggest multi-year stanzas of Velella mass beaching events, most prominently in 2014− 2018, interspersed by years where mass beaching events were seemingly absent, inferred from the lack of reports in those years (Figs. 2 & 3). Purcell et al. (2015) also reported on the anomalous occurrence of Velella on beaches in the North American Pacific during 2014. Their study corroborates the larger pattern we found, as only in 2014 did their data suggest that the relative occurrence (Velella sightings over all gelatinous zooplankton reports) was significant (52.8% or 141 out of 267 reports). Prior to 2014, Velella comprised no more than 3% of all reported beachings. Our model suggests that one obvious factor underlying these stanzas is winter SST (Table 1, Fig. 6), which was significantly warmer than climatological normal during/following the onset of the northeast Pacific marine heatwave  (Fig. 6) and lacked reports of Velella, suggesting that Velella blooms are either restricted (i.e. lower abundance), or prevented from beaching following colder/ stormier winters. Zeman et al. (2018) found that during the heatwave, Velella colonies along the northern CCLME (42.5°−46.5° N) had high ingestion rates of northern anchovy Engraulis mordax eggs, particularly in the region of the Columbia River plume where this fish species is known to spawn. They hypothesized that during the heatwave years, the relatively warmer winter−spring months allowed temporal expansion of spawning northern anchovy (Auth et al. 2018), which may have facilitated spring blooms of Velella via the additional food resource of anchovy eggs. Our spatial data support this hypothesis, as the highest reporting rates of Velella in the COASST dataset were observed in the region of the Columbia River plume during the spring months of the heatwave years (Fig. 4).
Our results, and those of Purcell et al. (2015) and Zeman et al. (2018), are demonstrative of Velella mass beaching associated with warmer temperatures within the northern CCLME. Other studies have also pointed to the association between warmer than normal SST and the occurrence of Velella elsewhere, including mass aggregations along the coast of Portugal (Pires et al. 2018), and anomalous sightings along the Pacific coast of South America (Carrera et al. 2019), the latter tied to El Niño− Southern Oscillation. More generally, recent state space modeling by Bellido et al. (2020) suggests an influence of positive winter SSTa on subsequent jellyfish swarms and mass beaching (Pelagia noctiluca along the coast of Malaga, Spain), and where data were collected in part by the digital citizen science application 'Info medusa App.' Collectively, this work suggests that further exploration of the role of warming ocean temperatures on the occurrence and relative abundance of jellyfish, including Vellela, is warranted.

Citizen science
Several studies of the distribution and interannual abundance of Velella have relied on citizen science data (e.g. METEOMEDUSE, EcoJel, JellyWatch Project: Purcell et al. 2015;GelAvista: Pires et al. 2018), although the majority of these programs are occurrence-only and often unassociated with systematic sampling. Purcell et al. (2015) suggested that Velella occurrence, as proxied by beaching reports, is particularly suited to citizen science, as 'such programs can provide data over large regions at relatively low cost' (p. 1064). The COASST dataset clearly demonstrates that even when Velella reporting is not re quired, a systematic monthly beach sampling program can return Velella data at volume and of a quality suitable for quantitative modeling. In part this is because Velella mass strandings are difficult to miss. At the same time, program participants often demonstrate the ability to make astute ob servations in addition to simple presence. As one COASST observer noted: 'The majority of the beach was covered with Velella Velella and they appeared to be smaller than what we have seen in the past as well as in larger numbers so I have attached a couple of photos just to document. Wrack is virtually nonexistent now as there is a major shortage of kelp. CM; 28 April 2016' Many authors have noted the necessity to move from Velella as a curiosity to Velella as a bona fide object of study , Pires et al. 2018, Zeman et al. 2018. Beach sampling programs such as COASST, designed to sample organisms or objects with a monthly/seasonal signal, may be in valuable to extending the study of Velella, as they can direct participants to sample both presence and absence, and could ask for basic information on disk size, colony freshness, and even colony density and percent cover along prescribed transects or quadrats (e.g. the information volunteered qualitatively in the quote above). That these latter measures could be collected digitally (i.e. with properly scaled images) makes the realistic possibility of large-scale, long-term data collection even easier.

Conclusion
Current remote sensing technology is able to detect physical signals in the marine environment (e.g. sea ice cover, SST, surface roughness) and lower trophic level response (e.g. chlorophyll), as well as indicators of human activity (e.g. ship tracking via vessel monitoring system), making calculations of global forcing relatively straightforward (Pettorelli et al. 2018, Werdell et al. 2019). However, tracking and understanding ecosystem response to these changes is lacking , which inhibits our ability to monitor for sudden shifts in ecosystems pushed beyond resistance and resilience thresholds (Tam et al. 2017. Given the persistence of spring onshore transport conditions in this region across all years such that Velella blooms, if they existed offshore, would be transported to the beach, our study suggests that winters with warmer than climatological average SST appeared most favorable to the appearance and persistence of large-scale blooms (Fig. 6A). Given the broad spatial extent of the blooms intuited here (e.g. Fig. 3), the predicted increase in heatwave events (Meehl & Tebaldi 2004, the presumed importance of Velella as an epipelagic predator  with potential impacts on forage fish populations (Zeman et al. 2018), and the massive transport of biomass from the coastal fronts to the nearshore represented by Velella strandings (e.g. Kemp 1986, Betti et al. 2019, we suggest that the widespread occurrence of Velella strandings may signal shifts in pelagic ecosystems, at least in eastern boundary current systems such as the CCLME. When combined with other ecosystem indicators of the dramatic, nonlinear effects of marine heatwaves in these systems, including harmful algal blooms (McCabe et al. 2016), multi-trophic level shifts in community composition and productivity (Cavole et al. 2016, Peterson et al. 2017, and seabird mass mortality events (Jones et al. 2018), Velella strandings fit into a larger picture of sudden and persistent shifts in trophic energy pathways, community biodiversity, and carrying capacity.