Tracking the snow line: Responses to climate change by New Zealand alpine invertebrates

ABSTRACT

We review and test an ecological paradigm that asserts that alpine invertebrate communities may shift upslope with climate warming. Our model couples the end of summer snow line (EOSS) elevation with invertebrate populations in New Zealand’s Southern Alps, using a forty-year data set, from fifty index glaciers. We show the snow line has risen an average 3.7 m a⁻¹. This is equivalent to raising alpine isotherms by almost 150 m and presents alpine biotic populations with four possible scenarios: upslope tracking, stasis, horizontal dispersal, or local adaptation. We characterize the alpine invertebrate biota (AIB) and present two case studies that show that high-elevation taxa have tracked the snow line within a narrow range (<20 m), whereas lower elevation taxa have potentially shifted by tens of meters. Relationships between the EOSS and Southern Oscillation Index (SOI) are investigated because precipitation and temperature influence snow line elevation by 25 percent. We also highlight the utility of invertebrates for monitoring climate change impacts on alpine ecosystems with a proposal for alpine climate monitoring units (CMUs), complementing an existing network of ecological management units (EMUs). We include an annotated list of New Zealand alpine invertebrates as potential indicators of climate change.

KEYWORDS:

• Alpine ecosystems
• climate change
• snow line elevation
• invertebrate communities
• weta
• New Zealand glaciers

Introduction

A 220-year-old paradigm

In 1807 the Prussian naturalist Alexander von Humboldt produced an illustration that established alpine ecology as a credible discipline. Humboldt’s diagram (“Naturgemälde”) showed vertical zonation of plants on the slopes of 6,263-m Mt. Chimborazu in Ecuador. With Aimé Bonpland, Humboldt climbed the mountain, taking fastidious measurements of altitude, vegetation, and air temperature. On publishing that data, Humboldt forever changed our perception of mountains from the inhospitable to scientifically valuable “islands in the sky” (Humboldt 1808; Humboldt and Bonpland 1808; Wulf 2015).

Today, alpine ecosystems are recognized as sensitive indicators of climate change (R. J. Wilson et al. 2005; Gobbi et al. 2007; Finch and Löffler 2009; Forister et al. 2010). Global temperature increase may force alpine communities to shift upslope, with increased susceptibility to fragmentation, stochastic events, and competition (Figure 1). For example, Morueta-Holme et al. (2015) repeated Humboldt’s Chimborazu measurements, showing a 500 m ascent of alpine communities during the preceding 207 years (2.4 m a⁻¹).

Figure 1. The Diminishing Real Estate Model (DREM). The alpine biota (shaded band) undergoes two climate regimes. Upper left: At T1, cool adapted alpine biota occupy a broadly connected area of the mountain scape during cool climates. At T2, the climate has warmed, displacing and isolating the alpine habitat on separate peaks of the same mountain scape. The alpine habitats represent two frustums (conic sections) whose surface area can be calculated (center left column). Note that both frustum depths remain uniform while the radii reduce. Under the DREM, several ecological outcomes are possible across generation time (center right column). Within alpine habitats, community ecology indices are expected to have characteristic curves (colored lines) with increased elevation and, all things being equal, respond monotonically to climate warming (red) or cooling (blue). The model assumes typical species–area curves with the formula for a single mountain analogous to the island biogeography area curve: S= CA^z, where S is the number of species, A is the area of frustum, and c and z are fitted intercepts

Although solar radiation and lapse rates differ between equatorial and temperate latitudes (Chen and Chen 2013), New Zealand alpine ecosystems may respond similarly to climate warming. Lacking approximately 210 years of data, information on New Zealand alpine ecosystems remains patchy (Macfarlane et al. 2010; Hoare et al. 2017). Climate and spatial data on alpine communities need to be acquired for ecological reference and for managing the alpine estate in a warming climate.

Here, we review the current understanding of New Zealand’s alpine invertebrate biota (AIB) in the Southern Alps (Figure 2) and present an investigation testing upslope tracking of that fauna. We couple population densities with elevation and relate these to a more than forty-year data set that indicates a constant rise of the end of summer snow line (EOSS).

Figure 2. Map showing central Southern Alps of the South Island, Aotearoa/New Zealand, with additional mountain ranges labeled. Snow crystal icons show large regions of present-day glaciation. Dark blue dots are the end of summer snow line index glaciers (n = 50) including potential locations for establishing Climate Monitoring Units (CMUs), indicated in light blue. Brewster Glacier and French Ridge are marked in red

The Southern alps and its invertebrate biota

The South Island’s Southern Alps/Ka Tiritiri O Te Moana are New Zealand’s highest and longest mountain range, extending some 500 km (Figure 2). Between 11 to 13 percent of the island is above the timberline (O’Donnell, Weston, and Monks 2017; this study), the majority of that includes the central alps, with Aoraki/Mt. Cook the highest point at 3,724 m a.s.l. We define the alpine zone as a spatial band extending from the climatic timberline to the permanent snow line (Wardle 1964; Burrows 1967; H. D. Wilson 1976; Mark and Dickinson 1997; Mark, Dickinson, and Hofstede 2000; Marris, Hawke, and Glenny 2019). Above the alpine zone is the nival, a composite of permanent snow, ice, and bare rock.

Most alpine-adapted invertebrates complete their life cycle between the timberline and permanent snow line. Populations about this zone are subject to the same climatic variables as most glaciers, including temperature (−15°C to +25°C), relative humidity (40–100 percent), precipitation (0–700 mm/h), radiation (short and longwave; 0.30–50 μm), insolation (e.g., 140 W/m², direct and diffuse), wind (up to 180 kph), and air pressure (600–1,100 mb; N. J. Cullen and Conway 2015). In these conditions, invertebrates have shortened breeding, feeding, and developmental periods and usually overwinter under snow, avoiding predators and completing their life cycles (Sinclair, Lord, and Thompson 2001).

Many alpine invertebrates differ sufficiently from lower-elevation conspecifics to be considered separate species. These “true” alpine taxa share a suite of adaptive characters, notably, flightlessness, freeze physiology, melanism, omnivory, free-living life stages, hemimetaboly, dormancy, and sun basking or boundary layer retention behavior (Dillon, Frazier, and Dudley 2006). The orthoptera (weta and grasshoppers) represent alpine keystone species given their semimetabolous life cycle, physical mobility during all stages (except egg), omnivorous diet, freeze physiology, and diurnal activity (Field 2001) and as prey for alpine spiders (W. Chinn pers. obs.).

Climate change effects on the New Zealand alpine biota

Data from Antarctic ice cores, the Mauna Loa Observatory (Keeling curve), and New Zealand’s Baring Head Station show that CO₂ levels have risen over a fifty-year period (Brailsford et al. 2012) and are far greater than any values recorded for the late Pleistocene (Etheridge et al. 1996; Lüthi et al. 2008). Carbon dioxide is a significant greenhouse gas and a major factor driving current atmospheric warming (Intergovernmental Panel on Climate Change 2014). Between 1950 and 2018 the concentration of CO₂ rose from 300 to 412 ppm, a rate that is unlikely to cease within the next 100 years due to atmospheric inertia (C. D. Keeling and Whorf 2004; R. F. Keeling 2008; Sundquist and Keeling 2009). New Zealand’s average air temperature increased by 1.02°C between 1909 and 2017, a rate of 1°C per century and roughly one order of magnitude faster than the preceding 2000 years of global average (Mann and Jones 2003; Mullan et al. 2010; NIWA 2019a).

Between latitudes 40° and 50° S, the alpine biota experience persistent west-to-east circumpolar air currents, the so-called roaring forties (Pendlebury and Barnes-Koeghan 2007; Jiang, Griffiths, and Lorrey 2013; Strother et al. 2014). This climate regime maintains a four-way partitioning of humidity and temperature across the Southern Alps, specifically, warm + wet, warm + dry, cool + wet, and cool + dry. The effect has been correlated with speciation among alpine taxa (W. G. H. Chinn and Gemmell 2004; Heenan and McGlone 2012; Goldberg, Morgan-Richards, and Trewick 2015); however, the extent to which warming might alter this local climate scenario (and life history traits) remains unknown.

Alpine communities may respond to warming by (1) upslope habitat tracking, (2), in situ adaptation (via phenotypic plasticity or long-term selection), (3) horizontal migration, or (4) localized extinction of populations (Table 1). The upslope tracking model holds that populations will colonize higher substrates in lockstep with rising optimal (realized) niche parameters. We describe this scenario as the “diminishing real estate model” (DREM); see Figure 1 for a summary of effects.

Table 1. Potential responses of, and effects on, alpine invertebrate populations to a predicted rise in elevation of optimal habitat conditions, driven by a warming climate in the Southern Alps of New Zealand

Response	Predicted population effects
Able to track, unable to adapt	Population survives where habitat exists. Localized extinction on isolated mountain tops likely.
Able to adapt, unable to track	In situ adaptive radiation and speciation. Typical response for isolated mountaintops and ranges.
Able to track and adapt	High degree of intraspecific variation across leading and trailing edges of an actively tracking population. Potential for hybridization at contact zones.
Unable to track and adapt	Local extinction.

The snow line as an ecological datum

Climate forcing of New Zealand’s alpine biota has an established theoretical pedigree (Cockayne 1928; Zotov 1938; Willet 1950; Fleming 1963; Bigelow 1967; Dumbleton 1969). These authors assumed that the entire alpine ecosystem will shift upslope (with warming) by slow colonization, relative to an ecological datum. By introducing the moniker “bio-glacial control,” Willet (1950) established the fluctuating snow line paradigm. The snow line represents an isotherm with a hypothetical equator-to-pole thermal gradient (ETPTG), which clips the North Island volcanoes at 2,500 m, lowering to about 1,400 m in southern Fiordland, some 157 m/1° latitude (a 0.1° slope; Willet 1950).

The EOSS is the highest net elevation of the regional snow line following cold-season snow accumulation and warm-season (summer and early autumn) melt (T. J. H. Chinn 1995; T. J. H. Chinn, Heydenrych, and Salinger 2005; Willsman, Chinn, and Macara 2015). Snow lines are influenced by precipitation, evaporation, temperature, albedo, and energy balance throughout the year. The EOSS is a transient feature and can be viewed on glaciers as dust-impregnated snow (firn) where seasonal snow accumulation equals loss (Figure 3). Since 1978, an aerial photography program has documented the EOSS on fifty index glaciers (Willsman, Chinn, and Macara 2015). The snow line rarely, if ever, survives over bare rock (due to high-temperature radiative melting) and is necessarily measured on glaciers. In this capacity, the EOSS is important for approximating the equilibrium line altitude (ELA), a hypothetical elevation where glacier mass balance will be at equilibrium for a given location (Meier 1962; Nesje 1992). Currently, the Southern Alps EOSS is rising between 3.0 and 3.7 m a⁻¹ and is correlated with climate warming (Willsman, Chinn, and Macara 2015; this study). Notably, snow line recession is seasonal and dynamic, whereas the EOSS rise is decadal (protracted) and significantly coupled to long-term response to ambient temperatures.

Figure 3. A typical End of Summer Snow line (EOSS) on the Rolleston Glacier, a very small cirque-type glacier in Arthur’s Pass National Park, Southern Alps. The glacier albedo acts to thermally buffer winter snow, which ablates slowly from the glacier snout at the onset of summer. The highest elevation of the snow line represents the average of all ablation sources and provides a relative position for the upper boundary isotherm of the alpine zone. Two snow patches are visible either side of the massif indicating that even under high Infra–Red radiation (from scree and rock) the snow line remains detectable

Multiple climate factors influence the EOSS in New Zealand. In a study by Purdie et al. (2011), seven climate variables were correlated with snow ablation (water equivalent) on the Tasman Glacier and four were found to be statistically significant, including proximal sea surface temperature (0.579: p < .001), ambient air temperature (0.510: p < .001), the Southern Annular Mode (SAM; 0.344: p = .017), and the Southern Oscillation Index (SOI; 0.317: p = .026). The SAM represents a circumpolar ring of differential air pressure between the Antarctic and mid-latitude atmosphere and is closely linked to the SOI (Jiang, Griffiths, and Lorrey 2013). The EOSS may also fluctuate on decadal scales in response to latitudinal shifts of meridional airwaves, themselves a function of circumpolar troughs. Further climatological work by Mackintosh et al. (2017) showed a connection between low- to high-latitude climate modes, regional atmospheric circulation patterns, sea surface temperatures offshore of the Southern Alps, and air temperatures, which collectively impact glacier mass balance.

The SOI measures atmospheric pressure differentials between Tahiti and Darwin, Australia (NIWA 2019c), and accounts for at least 30 percent of New Zealand’s seasonal climate variability as part of the El Niño–Southern Oscillation (ENSO) phenomenon (Troup 1965; Clare et al. 2002; Bhaskaran and Mullan 2003). For our purposes, the SOI typically manifests as weather pattern reversals (El Niño/La Niña seasons; Table 2; Jiang, Griffiths, and Lorrey 2013). Because the SOI influences New Zealand weather by at least 25 percent (Purdie et al. 2011), we discuss the relationships between the SOI, the EOSS, and potential influences on alpine communities. It should be noted, however, that several modes of variability drive Southern Alps glacier response (Mackintosh et al. 2017), but for simplicity we evaluate New Zealand’s mid-latitude climate in response to tropical influences as represented by ENSO.

Table 2. Characteristics of the Southern Oscillation Index (SOI) and effects on New Zealand weather and snow line elevation in the Southern Alps

Scenario	New Zealand air pressure (sea level)	New Zealand sea surface temperature	SOI	New Zealand weather patterns	Correlated snow line
El Niño	Low	Cool	−ve	Latitude increase of west-to-east jet stream. More cyclonic westerlies (low-pressure systems). Dry in the east, higher precipitation on west coast. Winter southwesterlies, deep snow	Lowered
La Niña	High	Warm	+ve	Anticyclones (high-pressure systems). Warm summer northeasterlies. Little alpine precipitation, reduced snowfall	Elevated

Note. Data from Clare et al. (2002), Jiang, Griffiths, and Lorrey (2013), Ministry for the Environment (2018), and NIWA (2019a).

Testing the diminishing real estate model

Evidence supporting upslope tracking within the Southern Alps is limited (Halloy and Mark 2003; Harsch et al. 2009) and has met with some skepticism (McGlone et al. 2010; McGlone and Walker 2011). Studies relating timberline elevation with climate warming have also produced varying conclusions (Wardle 1985; Harsch et al. 2012). Beech forest seed production has been shown to increase with rising temperatures, although a considerable lag phase occurs (Wardle 1971, 2008; Whitehead, Leathwick, and Hobbs 1992). Efforts to detect the impact of atmospheric warming on alpine ecosystems include alpine grasshopper studies (White and Sedcole 1991; Koot 2018), timberline position (Wardle 1971; L. E. Cullen et al. 2001; Harsch et al. 2012; Case and Duncan 2014; Cieraad and McGlone 2014), plant community composition (Halloy and Mark 2003; Lord et al. 2017), and predator influxes (Christie et al. [2017] and references therein).

Ensemble climate modeling has shown promise for relating geographic range shifts among New Zealand insects (Bulgarella et al. 2014), including alpine weta and grasshoppers (Sivyer et al. 2018). For example, Sivyer et al. (2018) used MaxEnt (Elith et al. 2011) and other modeling approaches to suggest that mean annual temperature had the highest predictive power (>70 percent) for determining grasshopper range. Interestingly, they also found that isotherms had the lowest predictive power, yet our study relies on calculated isotherms in the alpine zone for population elevation. The EOSS, by contrast, is an actual elevation parameter, observable on the landscape. In all cases, correlating a warming atmosphere with a response by alpine communities requires at least two measures;

1. Multidecade climate trend data, specific to the alpine environment.

2. Biota physiologically coupled to the niche parameters within the alpine zone.

By coupling an appropriate environmental temperature lapse rate with a fixed elevation, we can expect the biota to experience warming at that location or to track upslope in order to maintain an optimal thermal envelope. For example, if the spring lapse rate at Mueller Hut in Aoraki/Mt. Cook National Park is about 6.6°C/1,000 m (Norton 1985), with a mean annual air temperature around 1.8°C at 1,760 m a.s.l. (NIWA 2019a, 2019b; Cliflo 2019), then we could expect an average temperature of 2.8°C in forty years, if scaled against the snow line rising at 3.7 m a⁻¹. The question is does the alpine biota adapt or shift in step with this rate function?

Conservation and alpine populations

The goal of alpine conservation is maintaining range, diversity, and pattern of high-elevation indigenous communities. The impact of climate warming on New Zealand’s alpine diversity remains unclear despite high endemism (up to 90 percent), with populations confined to single mountain ranges (W. G. H. Chinn and Gemmell 2004; Mark 2012). Greater predictability of climate effects would assist management while providing an ecological “early warning.” Ideally, alpine climate models should calculate extinction risk for select populations using three terms: which mountainous regions will be most affected (spatial), when climate warming effects might manifest (temporal), and at what rate alpine communities might respond (a constant).

Habitat area is an important consideration for extinction risk (Rosenzweig 1995). Rapid climate warming may fragment invertebrate populations, leading to critically small sizes. On the other hand, long-term (>1,000 years) segregation can drive speciation through genetic drift in isolation (Trewick 2001, 2005; Buckley and Simon 2007; Trewick and Morgan-Richards 2019; Trewick, Wallis, and Morgan-Richards 2011; Wallis et al. 2016), the so-called Pleistocene species pump (Schoville, Roderick, and Kavanaugh 2012). Assessing long-term persistence will involve a trade-off between the rate and extent of environmental change and the adaptive capacity of populations. Community-level indices of the alpine biota are critical because species richness is expected to decline as the climate warms (Peters 1992; Schneider 1992; Parmesan 2005). However, species abundance is less clear because the biomass of a keystone species may increase at specific elevations (Figure 1).

Current management of New Zealand alpine ecosystems includes bird and lizard monitoring, specifically Rock Wren (Xenicus gilviventris) and Kea (Nestor notabilis), and to some extent alpine gecko (O’Donnell, Weston, and Monks 2017). This work is framed within ecological management units (EMUs), areas in which community indices are repeatedly measured. The Southern Alps are conveniently scaled for EMU work. Cirque basins, long valleys, and spurs offer useful ecological gradients. Table mountains and sharp climatic gradients present natural limits to populations which can be employed to protect and monitor communities. The Zero Invasive Predators (ZIP 2018) initiative in the Perth Valley, west of the Main Divide is an example. We propose augmenting alpine EMUs with climate monitoring units (CMUs), potentially stimulating an alpine climate and conservation paradigm.

Aims

We investigate three themes: (1) Can upslope tracking be detected using the EOSS as a climate proxy? (2) What rates of upslope movement can we estimate using alpine invertebrates? (3) An assessment of the alpine invertebrate fauna for monitoring climate change effects on alpine ecosystems.

Methods

Study sites

We selected two alpine study sites Brewster Glacier (44°4′15″ S, 169°26′11″ E) and French Ridge (44°25′20″ S, 168°42′09″ E; Figure 2, Figure 4a-h). Brewster Glacier is a significant glaciology and climatology research area (George 2005; Willis et al. 2009; Anderson et al. 2010; Gillet and Cullen 2011; N. J. Cullen and Conway 2015; N. J. Cullen et al. 2016; Sirguey et al. 2016), and data from those studies may benefit future ecological work. French Ridge, located in Mt. Aspiring/Tititea National Park, is 70 km southwest of Brewster Glacier and represents an ecologically confined ridge system. French Ridge is situated between two EOSS index glaciers (Fog and Findlay; 14 and 23 km, respectively).

Figure 4. Composite images of the two sampling locations, Brewster Glacier (a-d) and French Ridge (e-h). Both sites support an ecologically isolated biota, confined to a ridge and extending from the timberline to the permanent snowline (or mountain summit). (a) Map of the Brewster Glacier and environs (Topomap 1:50 000 Series). Red line depicts invertebrate sampling transect, grid squares=1km. (b) aerial photograph of Brewster Glacier in 2016, view toward south west. (c) Typical snowline aerial photograph of Brewster Glacier showing firn line and End Of Summer Snow line (EOSS) with sampling transect (red line). (d) Ground view of Mt Brewster taken from transect line. Image insert is a female Pharmacus brewsterensis: Rhaphidophoridae, a resident of this habitat. (e) Map of French Ridge, Mt Aspiring/Tititea National Park (Topomap 1:50 000 Series). Red line depicts invertebrate sampling transect. Grid squares=1km. (f) Aerial photograph of French Ridge (showing transect) and Bonar Glacier. (g) View looking up to Mt French, from 1800 m a.s.l. Tapering red line represents transect path. Insert image shows the scree weta Deinacrida connectens feeding (at night) on a native Ranunculus herb. (h) Down-slope view of French Ridge showing hut location, sampling transect and example of alpine herb field/stone pavement

Invertebrate sampling

Invertebrates were sampled from transects between 1,400 and 2,200 m a.s.l. (800 m vertical, the maximum range for thermal effect). Sampling occurred in March 2016 (the last month of highest snow line in the Southern Alps) and coincided with the EOSS survey flights.

At 100-m contours, quadrat lines (N = 9) were established perpendicular to the main transect with five quadrats per line. Quadrats were 3 × 3 m with 15-m separation over 60 m. If terrain or obstacles prevented a complete line, the series was extended in the opposite direction. Invertebrates were counted per hour, per 100-m contour, and the highest density per elevation (the median) was used to infer the optimal thermal niche. Two people made counts by eye. Raw data were collected with counters and score cards, and elevation and distance were determined using a Garmin 62s Global Positioning System instrument and tape measures.

Invertebrate collecting was by hand, sweep net, and rock turning (day and dusk). Timed sweep netting of vegetation was employed for flying insects. Specimens were identified to species or lowest recognizable taxonomic unit. We omitted diptera from the study, a diverse group difficult to correlate with elevation. Representative specimens were preserved and lodged with the author’s collection at the Department of Conservation.

Population densities

Population density × Elevation was assessed using the violin plot function in R Studio and ggplot (R Core Team 2013). For these plots, invertebrate counts per contour were normalized against the highest value in order to visually accommodate the range (without resorting to a log scale or similar). Significance tests were performed in R (R Core Team 2013; Richards 1972).

The EOSS data set as a climate proxy

Long-term climate trend data are required to correlate alpine community responses to change.

In New Zealand, long-term in situ observations above 1,000 m are sparse (fewer than fifteen climate stations) and temporally incomplete, hindering ecological study. The EOSS survey provides surrogate long-term environmental data on the alpine zone with inherent linkages to multiple climate variables, particularly temperature.

Snow line elevations were derived from aerial photographs of fifty index glaciers (Table 3). The photographic series began in 1978 and continues today (Willsman, Chinn, and Macara 2015; SIRG 2019). Prior to the 1990s, the EOSS was estimated from photographs by eye, a method that incurred unknown error. Since 2001, the photographs were orthorectified in Photogeoref v 1.0 (Corripio 2004) using five correction variables: camera Global Positioning System position, focal length, image size, land digital elevation models, and ground control points (cross-referenced with topographical maps).

Table 3. Summary data for fifty snow line index glaciers in the Southern Alps, New Zealand. Glaciers are listed geographically from north to south

Glacier name	Glacier type^a	Latitude/longitude (decimal)	Névé elevation^b	Mean EOSS elevation (m a.s.l.) (SE)	Estimated rate of snow line rise (ma^−1)^c
Kaikoura Range	Rock	42.010/173.637	2,420	2,491 (7.5)	1.16
Mt. Ella	Glacierette	42.081/172.579	2,165	2,148 (13.9)	2.87
Mt. Faerie Queen	Glacierette	42.264/172.514	2,070	1,991 (17.3)	3.12
Mt. Wilson	Snow patch	42.932/171.678	1,885	1,908 (16.3)	3.16
Mt. Franklin	Cirque	42.870/171.661	1,900	1,825 (18.9)	3.15
Rolleston Gl.	Cirque	42.889/171.527	1,800	1,763 (10.3)	1.54
Mt. Carrington	Cirque	42.912/171.473	1,840	1,707 (23)	3.46
Mt. Avoca	Glacierette	43.034/171.404	1,960	1,964 (10.8)	2.82
Marmaduke Gl.	Mountain	42.986/171.382	1,860	1,837 (17.4)	2.51
Retreat Gl.	Mountain	42.961/171.304	1,800	1,737 (17.6)	3.17
Browning R.	Small cirque	42.920/171.269	1,640	1,547 (47)	1.87
South Cameron	Mountain	43.355/170.983	2,400	2,272 (17)	4.66
Mt. Butler	Shelf Gl.	43.249/170.931	1,840	1,824 (16.7)	2.00
Dainty Gl.	Mountain	43.232/170.889	1,980	1,942 (14.6)	2.82
Kea Gl.	Cirque	43.175/170.800	1,900	1,800 (24)	3.39
Jaspur Gl.	Glacierette	43.181/170.756	1,780	1,724 (21)	4.56
Seige Gl.^d	Valley	43.264/170.531	1,900	1,717 (37.7)	7.44
Vertebrate R. a	Mountain	43.318/170.614	1,900	1,865 (15.5)	2.97
Vertebrate R. b	Mountain	43.328/170.618	1,900	1,835 (11)	2.20
Ridge Gl.	Mountain	43.620/170.362	2,280	2,234 (15.4)	3.52
Langdale Gl.	Cirque	43.578/170.270	2,240	2,196 (28)	4.46
Tasman Gl.	Valley	43.509/170.337	2,400	1,806 (18)	3.30
Salisbury Gl.	Mountain	43.469/170.224	2,000	1,805 (17)	2.61
Jalf Gl.	Saddle Gl.	43.469/170.148	1,700	1,770 (26.7)	5.23
Chancellor dome	Cirque	43.510/170.135	1,900	1,737 (23)	1.37
Glenmary Gl.	Cirque	43.990/169.887	2,260	2,168 (11.3)	2.77
Blair Gl.	Mountain	43.986/169.764	2,140	1,940 (15)	3.86
Mt. MacKenzie	Mountain	44.800/167.710	1,660	1,891 (16)	2.27
Jackson Gl.	Mountain	43.889/169.790	2,060	2,060 (10.6)	2.13
Jack Gl.	Small cirque	43.809/169.639	1,824	1,901 (15.7)	2.56
Mt. St. Mary	Rock	44.260/169.654	2,080	1,917 (24.5)	6.00
Thurneyson Gl.	Mountain	44.163/169.600	2,200	1,967 (14.9)	4.20
Brewster Gl.^d	Mountain	44.070/169.435	1,980	1,924 (27)	6.25
Mt. Stuart	Mountain	44.108/169.261	1,840	1,671 (17.4)	3.54
Lindsay Gl.	Mountain/shelf	44.000/169.128	1,800	1,722 (18)	3.29
Fog Pk	Cirque	44.523/168.808	2,000	1,999 (16)	4.95
Snowy Gl.	Small mountain	44.648/168.746	1,700	2,089 (11.7)	2.84
Mt. Caria	Small cirque	44.380/168.523	1,700	1,460 (15.7)	2.61
Findlay Gl.	Cirque	44.321/168.446	1,800	1,690 (18.6)	4.45
Park Pass Gl.	Valley	44.584/168.232	2,000	1,828 (17.2)	2.88
Mt. Larkins	Glacierette	44.880/168.495	2,200	2,031 (28.4)	9.45
Bryant Gl.	Mountain	44.825/168.284	1,800	1,767 (21.2)	3.49
Ailsa Mts.	Cirque	44.910/168.241	1,700	1,631 (12)	1.97
Mt. Gunn	Cirque	44.758/168.087	1,680	1,591 (19.6)	3.66
Mt. Gendarme	Mountain	44.788/167.941	1,700	1,590 (17.6)	2.54
Llawrenny Pk	Cirque	44.658/167.808	1,800	1,464 (18)	3.48
Barrier Pk	Mountain	44.848/167.719	1,780	1,600 (27)	7.34
Mt. Irene	Glacierette	45.178/167.361	1,800	1,566 (22)	5.50
Merrie Ra.	Patchy	45.705/167.188	1,500	1,526 (28)	7.47
Caroline Pk	Glacierette	45.945/167.196	1,500	1,369 (24)	8.02
Mean (SD)					3.73 (1.8)

Notes. Glacier list and EOSS derived from Willsman, Chinn, and Macara (2015). Gl. = Glacier; R. = Range; Pk = Peak.

^aData from Willsman, Chinn, and Macara (2015).

^bTaken at the approximate geometric center of each neve.

^cRise derived from differential points of the regression slope [y = 2mx – b] – [y = 1mx – b].

Snow line rise = snow line elevation at T2 − T1/years elapsed.

^dPotential locations for CMUs.

Snow line elevations are derived values, estimated from plotting an expert guided snow line demarcation against a geodetic topographic map. Following digitization of the glacier snow lines, an ablation-to-accumulation ratio is calculated. Snow line elevations can then be accurately assessed from an area–altitude curve for each glacier or directly from a digital elevation model (Willsman, Chinn, and Macara 2015, 2017). These values were used for our analysis.

Snow line trend data

Snow line elevations were plotted (linear regression) against space and time to assess slope and trend throughout the Southern Alps. Linear distance was measured on longitude with the Kaikoura Mountains as a datum for all plots. Rates of annual snow line rise were derived by the equation (slope × T2 + intercept) − (slope × T1 + intercept). The relationship between the snow line and timber line was also investigated. Timberline elevations were derived from 1:50,000 scale maps (NZ Topo50 Series) within a 5-km orbit of index glaciers. Where visible, timberline elevations were cross-referenced with aerial photographs per glacier.

Southern oscillation index

Frequency and amplitude of the SOI were assessed as a composite of snow line elevation and, by association, invertebrate population elevation. The SOI was selected as an index to explain snow line fluctuation because it convolves several atmospheric and surface climate variables and has a strong relationship with Tasman sea surface temperatures, which in turn correlate strongly with glacier volumes and the ELA (Purdie et al. [2011]; Bowen et al. [2017] and references therein; Mackintosh et al. [2017]). The SOI was acquired from the U.S. National Oceanic and Atmospheric Administration and constrained for the period 1977 to 2016 prior to correlating and plotting annual snow line data (NOAA 2019).

The snow line tracking model

Our model (Figure 5) includes five assumptions: (1) Population densities are thermally coupled to elevation. (2) Population densities have remained thermally linked to elevation for the preceding forty years (the snow line data period). (3) Upslope movement is proportion to shifts of the snow line isotherm between periods T1 and T2. (4) For a given slope, the elevation differential between the snow line at T2 and population median provides a rough colonization rate (m a⁻¹). (5) If all populations are strictly coupled to the snow line elevation, we should expect faster rates of colonization across shallow slopes compared to steep slopes (varying rates for equivalent vertical gain and period). Conversely, if elevation varies across differing gradients following an equal time period, then abiotic isotherm tracking may be subordinate to biotic factors.

Figure 5. Model of upslope habitat tracking by an Alpine Invertebrate Biota (AIB). Colored temperature ellipses represent populations of different species and their optimal (realized) niche elevation. Densities are represented by three median width classes. The diagram shows an upslope shift of 120 m by the End of Summer Snow line (EOSS) between 1977 (T1) and 2017 (T2), a response correlated with atmospheric warming. T2 populations located at, near, or even extending into the previous snow line are assumed to have tracked up slope since T1. Populations can be directly sampled today for occupancy on snow- and ice-free slopes at T2. Nine populations protrude above the T1 snow line (blue dashes). For the remaining populations, upslope tracking has occurred if their T2 median densities are at elevations previously cooler than today’s apparent optimal temperatures, given a lapse rate of 6.6°C km⁻¹ relative to the snow line. The “red” population provides an example; at T2 the highest population density is located at 1,520 m a.s.l., with an isotherm of 3.6°C. Forty years ago, the same isotherm was 120 m lower

The model includes an estimate of the snow line isotherm. We derived an average snow line elevation for the Southern Alps at 1,600 m a.s.l. (excluding the ETPTG) across all years by converting the glacier ELA to meters above sea level using data from Willsman, Chinn, and Macara (2015, 2017). At 1,600 m a.s.l. the snow line isotherm was calculated using an eighty-year temperature record from near sea level (at Hokitika airport, Westland; NIWA 2019d). Coupling the mean sea level temperature at Hokitika (11.5°C) with a 6.6°C km⁻¹ lapse rate gives an estimated mean annual temperature of 0.94°C for the snow line, our isotherm of interest. Detailed studies of environmental lapse rates in New Zealand are limited, and we have not included the various uncertainties here. The 6.6°C km⁻¹ lapse rate conforms to Norton’s alpine spring maximum (Norton 1985; Tait and Macara 2014) and serves to test the fidelity of alpine invertebrate communities to narrow thermal zones since a high lapse rate may amplify any tracking signal (Equation (1)(1) $11.5^{\mathrm{o}}\mathrm{C}-\left(1600\times \left(6.6^{\mathrm{o}}\mathrm{C}\mathrm{k}{\mathrm{m}}^{-1}\right)\right).$(1) ):

(1) $11.5^{\mathrm{o}}\mathrm{C}-\left(1600\times \left(6.6^{\mathrm{o}}\mathrm{C}\mathrm{k}{\mathrm{m}}^{-1}\right)\right).$(1)

Hokitika is situated near the Southern Alps (42° S and 33 km from the main divide), approximately halfway along the Alps. A measurable shift of the snow line represents a repositioning of the 1°C isotherm and potential ecosystem driver. Populations of the AIB located adjacent to the EOSS should track at a rate equal (or similar) to the snow line shift (Equation (2)(2) $\mathrm{A}\mathrm{I}{\mathrm{B}}_{\mathrm{e}\mathrm{r}}=\frac{\mathrm{\Delta }{\mathrm{S}}_{\mathrm{l}}}{\mathrm{\Delta }\mathrm{y}\mathrm{r}}$(2) ):

(2) $\mathrm{A}\mathrm{I}{\mathrm{B}}_{\mathrm{e}\mathrm{r}}=\frac{\mathrm{\Delta }{\mathrm{S}}_{\mathrm{l}}}{\mathrm{\Delta }\mathrm{y}\mathrm{r}}$(2)

where AIB_er indicates alpine invertebrate biota elevation rise, ΔS_l is the change in elevation of the snow line, and Δyr is the number of years elapsed between snow line elevations.

For populations below the snow line, maximum densities are assumed to be static at specific elevations, whereas isotherms shift upslope with warming. Alternatively, those populations may keep pace with their optimal thermal envelope; that is, track. Because we cannot go back in time to observe shifts, we can measure the relative density of today’s population elevation and construct density elevation plots. Then we calculate the optimal isotherm (and elevation) for maximum population densities (Equation (3)(3) $\mathrm{C}\mathrm{a}\mathrm{l}\mathrm{c}\mathrm{u}\mathrm{l}\mathrm{a}\mathrm{t}\mathrm{e}\mathrm{d}\mathrm{e}\mathrm{l}\mathrm{e}\mathrm{v}\mathrm{a}\mathrm{t}\mathrm{i}\mathrm{o}\mathrm{n}={\mathrm{O}}_{\mathrm{p}\mathrm{o}\mathrm{p}}-\left[\frac{\mathrm{P}\mathrm{o}{\mathrm{p}}_{\mathrm{T}}-{\mathrm{\Delta }}_{\mathrm{T}\mathrm{S}\mathrm{L}}}{{ }^{\mathrm{o}}\mathrm{C}\mathrm{k}{\mathrm{m}}^{-1}\hspace{0.17em}}\right]$(3) ):

(3) $\mathrm{C}\mathrm{a}\mathrm{l}\mathrm{c}\mathrm{u}\mathrm{l}\mathrm{a}\mathrm{t}\mathrm{e}\mathrm{d}\mathrm{e}\mathrm{l}\mathrm{e}\mathrm{v}\mathrm{a}\mathrm{t}\mathrm{i}\mathrm{o}\mathrm{n}={\mathrm{O}}_{\mathrm{p}\mathrm{o}\mathrm{p}}-\left[\frac{\mathrm{P}\mathrm{o}{\mathrm{p}}_{\mathrm{T}}-{\mathrm{\Delta }}_{\mathrm{T}\mathrm{S}\mathrm{L}}}{{ }^{\mathrm{o}}\mathrm{C}\mathrm{k}{\mathrm{m}}^{-1}\hspace{0.17em}}\right]$(3)

where O_pop is the observed population elevation, Pop_T is the optimal temperature of population at T2, ∆_TSL is the temperature change of snow line rise, and ^oC km⁻¹ is the lapse rate (6.6°C/1,000 m).

Expanding the above equation terms gives

Pop_T = (Δ elevation_SL-CP) × (Lapse rate) + Temp_SL,

where Δ elevation_SL-CP is the elevation difference between T2 snow line and T2 population elevation, the lapse rate = 6.6°C/km⁻¹, Temp_SL = + 1°C, CP is the core population elevation, SL is the snow line elevation, and

∆_TSL = (T2 snow line elevation − T1 snow line elevation) × (6.6°C/km⁻¹).

Slope geometry

Snow lines and alpine communities represent positions on the hypotenuse of a right angle triangle (the mountainside). For an equivalent vertical movement, elevation shifts will appear faster on shallow slopes relative to steep slopes. We can exploit this property to roughly differentiate between biotic and abiotic drivers of elevation change. For example, if the AIB is coupled to an isotherm, a shift of that optimum should induce varying rates of colonization depending on the slope gradient. Conversely, if upslope colonization does not change across varying slopes, irrespective of an upward shift of the snow line (isotherm), biotic or other factors may be more significant than maintenance of a given thermal zone.

Under this scenario, population-level adaptation may be occurring, expressed as behavioral and physiological plasticity; for example, locating mates, predator avoidance, reduced competition (inter and intra), or exploiting an increase of organic matter as the snow line rises. In all cases, it will be necessary to calibrate tracking rates at each location as part of any model seeking to measure elevation change of the biota. The goal would be to plot varying rates of population shift throughout the Alps, highlighting regions with “slow” or “fast” tracking rates, if indeed they occur.

Climate monitoring units

We made a desktop assessment of locations within the Southern Alps suitable for establishing CMUs. Criteria included ecological isolation, contiguous habitat from subalpine to nival, conical mountains, cirques or broad ridges extending from valley floor to mountaintop, and close proximity to index glaciers with measurable snow lines (details in Table 4).

Table 4. Site selection criteria for establishing Climate Monitoring Units (CMUs) and measuring upslope tracking of the End of Summer Snow line (EOSS) by alpine invertebrates in the Southern Alps of New Zealand

Alpine site criteria	Optimal characteristics
Mountain characteristics	Shape: conical, ramp, ridge, or cirque basin with minimal “ecological escape” Relief: Simple with measurable frustum and defined surface area polygons. Slope angle: 15° to 20°
Index glacier criteria	Preferably an index glacier is present on study massif or within an 8-km radius
Snow line shift rate	Known and measurable gradient (>0.5 m a^−1) with low variance
Geology	Relatively stable geology. Avoid argillite, shales, and active screes. Geochemically benign (e.g., not ultramafic metal ore-bearing rock)
Latitude	Wide and evenly distributed range
Ecology	Indigenous alpine communities with measurable vertical zonation
Access	Accessible by air and foot (potentially all-weather foot access)

Results

The alpine invertebrate biota

Slope and total recognizable taxonomic unit counts for Brewster Glacier and French Ridge were respectively 11° and twenty-nine and 20° and twenty-five (see Table 5 for species checklist). Species abundance and richness curves were typical: few groups with numerous individuals and many taxa with few (Figure 6a and Figure 6b). Taxon richness at Brewster Glacier showed a negative relationship with elevation, where the highest proportion of taxa was at 1,600 m a.s.l. (R² = 0.965), between shrublands and alpine snowgrass. Values were derived from a best-fit, least-squares fourth-order polynomial (Figure 7a). Results from French Ridge also showed a negative richness relationship with elevation (linear best-fit model R² = 0.946; Figure 7a and Figure 7b).

Table 5. Checklist of alpine invertebrate taxa from Brewster Glacier and French Ridge, considered as members of the Alpine Invertebrate Biota (AIB). Maximum counts per elevation (at T2) are listed

				Max count and elevation (m a.s.l.)
Order	Family	Genus	Species	BG		FR
Punctoidea	Charopidae	Allodiscus	kakano?	18	1,600	–	–
	Punctidae	Indet sp.	sp.	–	–	14	1,400
Amphipoda	Talitridae	Talorchestia	sp.	10	1,400	48	1,400
Geophilomorpha	Chelenophilidae	Zelanion	sp.	10	1,400	–	–
Polydesmida	Dalodesmidae	Icosidesmus	sp.	–	–	8	1,400
Opiliones	Phalangidae	Pantopsalis	sp.	–	–	18	1,400
	Triaenonichidae	Prasma	sp.	10	1,400	15	1,400
Pseudoscorpiones	Cheliferidae	Nesidiochernes	sp.	8	2,000	–	–
Araneae	Araneidae	Cryptaranea	subalpina	15	1,400	–	–
	Agelenidae	Neoramia	mamoea	60	1,500	–	–
	Stiphidiidae	Cambridgea	obscura?	–	–	–	–
	Lycosidae	Anoteropsis	alpina	17	1,800	7	1,700
	Lycosidae	Anoteropsis	hilaris	32	1,500	–	–
	Lycosidae	Anoteropsis	Indet sp.	6	1,900	22	1,400
	Salticidae	Indet	sp.	10	1,900	12	1,800
	Zoropsidae	Uliodon	sp.	–	–	15	1,400
Blattodea	Blattidae	Celatoblatta	montana	18	1,400	14	1,600
	Blattidae	Celatoblatta	quinquaemaculata	50	1,800	–	–
Coleoptera	Byrrhidae	Liochora	sp.	38	1,900	32	1,800
	Curculionidae	Inophloeus	sp.	26	1,600	8	1,500
	Carabidae	Mecodema	sp.	6	1,500	4	1,500
Hemiptera	Cicadidae	Maoricicada	phaeoptera	32	2,100	25	1,800
	Cicadidae	Kikihia	subalpina	18	1,400	12	1,500
Lepidoptera	Crambidae	Orocrambus	sp.	36	1,400	22	1,400
	Geometridae	sp. 1		12	1,400	8	1,400
	Geometridae	sp. 2		4	1,500	–	–
	Lycaenidae	Boldenaria	boldenarum	15	1,400	2	1,500
	Nymphalidae	Argyrophenga	antipodum	15	1,400	–	–
	Nymphalidae	Percnodaimon	merula	16	2,000	10	1,800
Orthoptera	Acrididae	Sigaus	australis	138	1,400	43	1,500
	Anostostomatidae	Deinacrida	connectens	4	1,800	4	1,800
	Anostostomatidae	Hemiandrus	sp.	5	1,400	4	1,400
	Rhaphidophoridae	Pharmacus	brewsterensis	12	2,000	16	2,000
	Rhaphidophoridae	Pharmacus	chapmanae	28	2,100	18	1,900

Notes. – Indicates absent taxa (n = 34). BG = Brewster Glacier; FR = French Ridge.

Figure 6. Alpine invertebrate species abundance plots for (a) Brewster Glacier (n = 29) and (b) French Ridge (n = 25). Data are total individual counts per taxon; elevation range was from 1,400 to 2,200 m a.s.l. Both locations generated characteristic curves, with few taxa represented by many individuals and many with few counts. Grasshoppers represent the majority of alpine invertebrate biomass, probably reflecting the density of snow tussock (Chionochloa sp.)

Figure 7. Species richness curves (number of taxa) per 100-m contour for (a) Brewster Glacier and (b) French Ridge. Both locations show strong negative trends. For Brewster Glacier, a fourth-order polynomial curve provided the best fit (R² = 0.965). Taxon richness peaked at around 1,600 m a.s.l., declining by approximately five taxa per 100-m upslope shift. For French Ridge the best-fit trace was linear, with a strong negative trend (R² = 0.947). French Ridge taxon richness declined at 2.7 taxa per 100 m of elevation gain

Grasshoppers (Sigaus australis) were 30 percent more abundant at Brewster Glacier (n = 426; 1,400 m a.s.l.) than cockroaches, the next highest taxon group (n = 129). Grasshopper median density declined to less than five insects per hour by 1,900 m a.s.l. Conversely, cave weta (Pharmacus chapmanae), first appeared at 1,700 m a.s.l. (n = 6) and at 1,900 m a.s.l. the conspecific P. brewsterensis was detected (at two per hour). By 2,100 m a.s.l., P. chapmanae counts were highest for any taxon at that elevation (twenty-eight per hour). At Brewster Glacier species richness declined to six taxa by 2,100 m a.s.l., including cave weta (P. chapmanae and P. brewsterensis), a terrestrial micro snail (Charopidae), alpine cicada (Maoricicada sp.), soil centipedes (Geophilomorpha), and the mountain ringlet butterfly Percnodaimon merula.

At French Ridge Sigaus australis grasshoppers were less abundant (109 versus 426 at Brewster Glacier) and land hoppers (Talorchestia sp.) were common (forty-eight per hour at 1,400 m a.s.l., declining to one specimen at 1,600 m a.s.l.). Snow grass (Chionochloa sp.) is a major forb for alpine grasshoppers, and French Ridge supports considerably fewer than the Brewster Glacier environment. Punctid snails were most frequent at 1,400 m a.s.l. (n = 14/h) and at both locations snails were found within crevices, moss clumps, and short-stature herbs (Ranunculus sp.). Both sites showed a consistent decline in diversity and richness with increasing elevation. Taxa above 1,700 m a.s.l. became restricted to a high-altitude assemblage of moss beetles (Byrrhidae), cave weta (Rhaphidophoridae), mountain ringlets (Percnodaimon), spiders (Lycosidae and Salticidae), cicada, and snails (Punctoidea). Of cicada, four individuals were detected per hour at 2,200 m a.s.l. up to the snow line.

Population density and elevation

Elevation data are presented in Figure 8a and Figure 8b. Most invertebrates were between 1,400 and 2,000 m a.s.l., although three taxa at Brewster Glacier were present at 2,200 m a.s.l., higher than the local EOSS. Values for lower elevation taxa are likely to overestimate the median elevation because counts were not taken below 1,400 m a.s.l. (the minimum sampling elevation). Several medians clustered at the same elevation, a potential sampling bias.

Figure 8. Violin plots showing relative densities of invertebrate populations for T2 (March 2016) at (a) Brewster Glacier and (b) French Ridge. Black dots are population medians; horizontal blue lines are snow lines (T1 and T2). Any taxa extending above the snow line were present either on snow or basking on rock outcrops close to the 2016 snow line. Colors are representative of woody shrubland (green), snowgrass (beige), and rock/snow (blue). Above 1,800 m a.s.l. the alpine habitat comprised short-horned grasshoppers (Acrididae), cave weta (Rhaphidophoridae), and free-living alpine wolf spiders (Anoteropsis alpina)

At Brewster Glacier, 55 percent of taxa were classifiable as members of the AIB (i.e., their vertical distributions were confined by lower and upper boundaries). At French Ridge, 60 percent were classifiable alpine species (including populations with a distribution “taper”; Figure 8b). The remaining taxa show density tails extending below 1,400 m a.s.l., indicating a wide elevation range. For example, the ubiquitous New Zealand wolf spider Anoteropsis hilaris was sampled between 1,400 and 1,900 m a.s.l. despite being known from sea level to the alpine zone.

At least 72 percent of high-elevation taxa were shared between locations with cave weta (P. chapmanae and P. brewsterensis; Rhaphidorphoridae), the highest at both sites (from 1,600 to 2,200 m a.s.l. and higher). These weta species are among the highest free-living, flightless insects in New Zealand.

Snow line spatial trend

Summary data for the fifty index glaciers are provided in Table 3. For the years 1977 and 2016, the mean rate of snow line rise was 3.73 m a⁻¹ (SD = 1.8), with an acceptable difference of 0.7 m. The regression fit of all snow line data (n = 3,207) against latitude (linear distance) shows a 0.07° slope, approximately 130 m per degree of latitude (Figure 9). The timberline showed a consistent 800-m separation below the snow line, the intervening area representing the alpine zone. We interpret the total slope as the Southern Alps ETPTG.

Figure 9. Graph showing elevations for the end of summer snow line (upper blue line) and timberline (lower green) plotted against latitudinal distance from the Kaikoura Mountains to Fiordland. Data are from fifty index glaciers in the Southern Alps of New Zealand and monitored between 1977 and 2016. Timberlines derived from 1:50,000 topographic maps. The snow line and timberline descend at approximately 130 m per degree of latitude (data from Willsman, Chinn, and Macara 2015, 2017). Colored shading represents our definition of the alpine zone and the region inhabited by the Alpine Invertebrate Biota (AIB), indicated as invertebrate silhouettes. Skyline trace shows high point elevations of mountains at specific distances from origin

Snow line temporal trend

Figure 10 shows the total snow line (and SOI) data for all index glaciers from 1977 to 2016. The small positive correlation (0.144) was significant, t(3192) = 26.6057, p = .001, α = .05, x̄ = 1,841.4 (m a.s.l.), n = 1,597, σ = 235.28, σ_x̄ = 5.894. Slope calculations indicate a snow line rise of between 3 and 3.7 m a⁻¹ across all index glaciers. Data points are considerably scattered about the regression line although up to 85 percent of the error is a function of the index glacier latitude relative to the ETPTG slope. Snow line elevation data oscillated about the mean by 55 m a⁻¹ (σ = 27.8; n = 23). Eight seasons showed oscillations greater than 100 m at an average frequency of five years (σ = 3; blue line); a local climate response to regional conditions arising from the Southern Oscillation (red line; see below).

Figure 10. Combination graph showing a relationship between actual End of Summer Snow line (EOSS) data (blue line) and the Southern Oscillation Index (SOI), sea level air pressure anomaly as fluctuating red line. The ascending linear blue line traces the least-squares fit for the EOSS data. The red line trace shows an increase of the SOI over the same period. Here, the snow line rises at 3.10 m a⁻¹ with a 0.144 correlation coeffcient (significantly greater than zero). The EOSS was correlated with SOI fluctuations (Kolmogorov-Smirnov test: D = 0.5, p = 0.1678). Positive SOI values represent La Niña weather anomalies and negative values represent the El Niño system

Snow line trend for Brewster Glacier and French Ridge

Site snow line tends are shown in Figure 11a and Figure 11b. Brewster Glacier indicated a snow line rise of 5.73 m a⁻¹ between 1978 (1,800 m a.s.l.) and 2016 (2,010 m a.s.l.; R² = 0.1573). The Findlay Glacier snow line (French Ridge proxy) showed a 4.45 m a⁻¹ rise over thirty-four years (130 m; R² = 0.1531). This 29 percent difference between Brewster Glacier and Findlay Glacier illustrates the local variability of the snow line isotherm and the potential sensitivity of the alpine biota to localized mountain ranges.

Figure 11. Local End of Summer Snowline (EOSS) plots for each sampling site. Upper plot: Snowline rise for Brewster Glacier showing a rate of 5.73m a^-1 (y5.7438x -9549.8 = 2023m–1811m a^-38). The positive correlation (0.43) was significant: t(34) = 22.2715, p < 0.05, x̄ = 1923, σ = 157.88 (R² = 0.16). Lower plot: Snowline rise over 35 years for the Findlay Glacier (nearest analogue to French Ridge) showing a rate of 4.45m a^-1, between years 1981 and 2016. The positive correlation (0.40) was significant: t(32)=22.2715, p < 0.05, x̄ = 1690, σ = 101.82 (R² = 0.1531)

Southern Oscillation Index

Correlating the SOI and the EOSS data produced a slightly significant but ambiguous positive coefficient (0.520), t(2) = −116.8, p = .05 (Figure 10). Additional analysis using a nonparametric Kolmogorov-Smirnov test was significant, D(20) = 0.5, p = .1678 (i.e., the data were from the same distribution with no overall difference in shape).

Figure 10 illustrates the phase relationship between snow lines and the SOI. which peaks at a one- to two-year lead on the EOSS. Several features of the SOI plot are noteworthy: (a) over half of the SOI accounts for the EOSS variability per year; (b) amplitudes of the EOSS oscillations are tightly coupled to the four highest SOI positive values for the years 1989, 1999, 2007, and 2011 (ten, eight, and four years, respectively), but the reverse was not true; and (c) a positive SOI appears the greater driver of upslope snow line elevation (than negative values for lowering the snow line).

Positive SOI values are representative of the La Niña (warmer) air mass systems (anticyclones from the north) and lower snowfall. However, the year 2015 shows a decoupling between the negative SOI (−0.7500) and EOSS (1,900 m a.s.l.). In both cases, there was a positive trend (R² = 0.09 and 0.08).

Snow line tracking: Actual

In Figure 12a and Figure 12b we present a synthesis of our model. The plots combine the EOSS, the AIB, two time periods, and both sites. The graphs show upslope shifts of several invertebrate populations from the T1 snow line elevation. The weta Pharmacus brewsterensis and P. chapmanae tracked at a rate equivalent to, or greater than, the maximum snow line rise at both sites (5.73 m a⁻¹ at Brewster Glacier and 4.45 m a⁻¹ for Findlay Glacier/French Ridge). This indicates that the weta are nearly equivalent in colonization ability at each site.

Figure 12. A synthesis of the model. Each graph shows four variables: The End of Summer Snow line (EOSS), the Alpine Invertebrate Biota (AIB, on abscissa), two time periods, and the two sites (a) Brewster Glacier and (b) French Ridge. Blue and yellow lines represent observed and calculated population density–elevation values, respectively. Snow line elevations are shown by the two horizontal blue lines for T1 and T2. For Brewster Glacier the snow line has risen 220 m, with seven taxa occupying slopes previously under the end of summer snow line some thirty-eight years ago (actual tracking). Six of those taxa had medians located at 1,800 m a.s.l. and two were measured at about the T2 snowline (1,910 m a.s.l.). This suggests that those taxa have actively tracked the snow line but not necessarily at the same rate. By contrast, the modelled tracking rise shows a considerably greater increase (up to 600 m) in elevation of populations over the same time period. The French Ridge (Findlay Glacier) snow line elevation shift was 160 m between T1 and T2. Lower elevation taxa potentially shift further with a rising snow line, whereas high elevation taxa are confined to a narrow thermal zone. The area between the taxa elevations is suggestive of a hysteresis function. See Discussion for details

By contrast, several taxa presented slower tracking rates. At Brewster Glacier, the median for jumping spider and moss beetle populations was 1,900 m a.s.l., midway between the T1 and T2 snow line elevations, with a nominal upslope tracking rate of at least 2.7 m a⁻¹. At French Ridge, five taxa were detected midway between T1 and T2 snow lines, some 40 m above the (proxy) T1 snow line, indicating a tracking rate of 1 m a⁻¹. The French Ridge taxa comprised jumping spiders, moss beetles, weta (Deinacrida connectens), Maoricicada, and mountain ringlets (Percnodaimon sp.).

Evidence of long-term invertebrate presence at elevations previously above the snow line (for at least thirty years) included silken retreats, exoskeleton molts, frass accumulation, and abandoned occupancy cells showing persistent mechanical redistribution of particulates. Taxa of particular interest included moss beetles (Liochoria sp.: Byrrhidae), jumping spiders (Salticidae), and two weta groups—the scree weta Deinacrida connectens and the cave weta P. brewsterensis and P. chapmanae—at both sites. These high-elevation taxa have limited dispersal ability (most are flightless) and nearly all are detritivores or scavengers, except for the spiders, which presumably feed on weta, collembola (spring tails), mites, and possibly pscoids (Psocoptera).

Snow line tracking: Modeled

The model formula was applied to all T2 median values (actual elevations). Populations at lower elevation showed greater apparent vertical movement than those at higher elevation, particularly for populations at 1,400 and 1,500 m a.s.l., in contrast to those at higher elevations (1,600–2,000 m a.s.l.; Figure 12a and Figure 12b).

At Brewster Glacier there was an apparent difference of 543 m between actual and calculated population elevation for seven taxa (the landhoppers, day flying moths, spiders, cockroaches, weevils, grasshoppers, and grass cicada). The elevation differential reduced to 518 m (for Geometrid moths), 443 m for Agelenid spiders, 143 m for D. connectens weta, and as little as 43 m for P. chapmanae. French Ridge data showed a similar taper in elevation differential between T1 and T2, specifically 603 m (weevils), 403 m (Mecodema ground beetle), 303 m (cockroaches), 63 m (jumping spiders), and 0 m for P. brewsterensis. The maximum calculated upslope tracking rates for lower elevation populations was about 16 m a⁻¹, at least five times greater than that for the snow line taxa.

For taxa occupying elevations above 1,700 m a.s.l., their summer average was probably close to 0°C (and cooler during winter). Those taxa were apparently restricted by the 1°C isotherm (snow line), with smaller shifts within the temperature/elevation envelope, in contrast to lower populations. These “confined taxa” included the weta (D. connectens, B. brewsterensis, and P. chapmanae), spiders (A. alpina, Lycosidae, and jumping spiders) and Maoricicada sp., byrrhid beetles, and the mountain ringlet Percnodaimon sp.

Nevertheless, populations showed a time-dependent rising trend relative to the snow line rise. Although the snow line positions are stepped (220 m for Brewster Glacier and 160 m for Findlay Glacier), the invertebrate populations show a convergence between actual and calculated elevation, per unit of height gain.

Climate monitoring units

We selected fifteen candidate locations in the Southern Alps meeting our criteria for long-term CMUs (Table 6, Figure 2). Each site included (or was close to) an index glacier where the recorded snow line rise was at least 2.0 m a⁻¹. Most sites were sufficiently high for the snow line to encompass at least 70 percent of a mountain’s circumference, which is desirable for testing the DREM.

Table 6. List of Southern Alps locations with potential for establishing Climate Monitoring Units (CMUs). Each location has a corresponding index glacier (within 8 km of the monitoring site), which can provide End of Summer Snow line (EOSS) elevation data

Study site (CMU)	Location and coordinates (latitude/longitude)	Index glacier(s)	Glacier elevation (m a.s.l.)^a	Snow line rise rate (m a^−1)
Mt. Adams	43.260317S/170.503296E	Seige, Vertebrae Col	2,000	5.0
Mt. Annette	43.751282S/170.061851E	Tasman	2,000	3.3
Mt. Brewster	44.065618S/169.450130E	Brewster	2,000	3.7
Chancellor Dome	43.510155S/170.135808E	Chancellor Dome, Jalf	1,900	3.2
Mt. Dechen	43.797061S/169.755760E	Jack and Jackson	1,824	2.3
De La Beche Ridge	43.527594S/170.259468E	Tasman	2,400	3.3
Mt. Ella	42.078472S/172.578195E	Ella Glacierette	2,000	4.5
The Footstool	43.675516S/170.063506E	Stocking/Tewaewae	2,000	3.3
Intervention Ridge	44.432962S/168.351712E	Findlay, Park Pass	2,000	3.0
Little Unknown	43.327305S/170.623581E	Vertebrae Col, Seige	1,900	3.0
Mt. Murchison	43.003933S/171.376704E	Avoca, Retreat, and White	2,100	3.3
Park Dome	43.159277S/170.928687E	Ivory, Kea, Jaspur	2,200	3.9
Park Pass	44.585826S/168.232006E	Park Pass	2,000	3.3
Plunket Dome	44.454622S/168.627251E	Snowballs	9,100	4.0
Mt. Rolleston	42.911314S/171.512303E	Rolleston and Crow	2,100	3.9

Note. ^aThe arithmetic mean elevation of snow accumulation (at névé center).

Several sites fell within the Southern Alps “beech gap” (Orwin 2007) and others were within beech-forested mountains (Table 6). This geographic spread presents an opportunity to test beech and non-beech ecosystems under the same climate shifts. In our opinion, alpine invertebrates offer multiple advantages for investigating ecological responses to climate shifts. Consequently, we provide a list of candidate invertebrate taxa and their ecological information (Table 5 and Table 7). These taxa are common to the Southern Alps and are generally confined to subalpine and alpine zones. We assume that these taxa are responding to optimal niche parameters including life cycles that use all, or part, of the available substrate (Nagy and Grabherr 2009).

Table 7. Candidate invertebrate taxa for monitoring alpine ecosystem condition within Ecological Management Units (EMUs) and Climate Monitoring Units (CMUs). Taxa are common to the New Zealand Southern Alps

Candidate groups found within the alpine zone	Elevation range (m a.s.l.)	Habitat	Life history traits	Alpine adaptations	Adult behavior (summer)	Diet	Sampling methods
Insecta: Blattodea Blattidae (Cockroaches)	0–2,000	Subalpine shrublands (litter); stony pavements; talus; stable rocky assemblages	Paurometabolous Highly mobile nymphs and adults (r-selected breeding)	Flightless. Winter torpor (all life stages) Freeze tolerant^a	Nocturnal and colonial	Nymphs and adults Omnivorous and detritivorous	Quadrat counts Rock turning (day), photo records Pitfall trapping (where possible)
Coleoptera Byrrhidae (Pill/moss beetles)	0–1,800	Forest to alpine Rocks; crevices; mosses	Holometabolous Congregates	Adults (winter quiescence) Spherical, heat retention Dark coloration	Nocturnal	Bryophagous; fungivory (mosses/lichens)	Transects; rock turning; quadrat counts; photo records
Carabidae (Ground beetles)	0–1,800	Forest; shrubland; cushion plants; seepages; snowbank vegetation	Holometabolous Egg brood care (K selection)	Winter quiescence (or diapause?) as adults	Nocturnal: Carabidae Diurnal: Cicindelidae)	Predatory on small invertebrates	Pitfall traps (night activity); wooden discs; quadrats; rock turning
Curculionidae (Weevils)	0–2,000	Subalpine and alpine herb fields	Holometabolous Usually large bodied	Adults (winter quiescence) Spherical, heat retention Dark coloration Flightless	Crepuscular and diurnal	Phytophagy; seeds Aciphylla, Anisotome sp.	Pitfall traps; quadrats; transect counts
Scarabaeidae (Chafer beetles) Odontria Prodontria Scythrodes Sericospilus	600–2,000	Alpine herb fields; cushion plants; tussock lands; snowbank communities	Holometabolous	Volant Volant Flightless Flightless Flightless	Crepuscular (adults) Diurnal Diurnal Diurnal	Larvae feed underground (roots)	Night (light trapping) and pitfall traps
Silvanidae Bark beetles (e.g., Protodendrophagus)	1,450–2,000	High alpine Rock tors, outcrops, slabs	Holometabolous	Winter torpor, in all life stages	Heat retention behavior	Crustose lichens/fungivory	Crevice searching/rock turning
Tenebrionidae (Darkling beetles) Artystona spp. Mimopeus spp.	600–1,800	High alpine (Artystona sp.) Rock tors; outcrops; under slabs	Holometabolous	Flightless Winter hibernaculum and quiescence	Nocturnal	Lichens; fungivory; herbivory	Crevice searching; rock turning
Diptera Muscidae (e.g., Exul singularis) Tachinidae (bristle flies) Tipulidae (crane flies) Syrphidae (hover flies)	0–2,200	Subalpine Shrublands; tussocks; herb fields; scree; snowbanks and melt pools	Holometabolous Ovoviviparity (tachinid flies)	Over wintering as maggot and pupae Cold tolerant	Very active on hot summer days	Nectar (tachinids and hover flies); vegetable matter	Transects; timed sweep netting
Hemiptera Cicadidae Maoricicada spp.	900–3,000	Tussock grasslands; dry rocky slopes High-altitude scree	Hemimetabolous	Freeze tolerance; sun basking. melanic; pilose/pubescence	Active on hot sunny days Males produce high-pitched call	Alpine herbs; tussock grasses	Transect sweep netting; 5-min recordings; sonogram analysis
Lygaeidae (Seed bugs) Rhypodes spp.	0–2,000	Cushion plants; herb fields; riverbeds; recently colonized surfaces	Hemimetabolous	Wide thermal tolerance Occupy microclimates (wind avoidance behavior)	Active on hot sunny days Adult stage tightly coupled with plant phenology	Seeds (Celmisia and Raoulia) Some nectar	Quadrat sampling; photo points (of host plants); pooters/D-vac or equivalent
Pentatomidae (Alpine soldier bug) Cermatulus nasalis hudsoni	0–2,000?	Shrublands and tussock grasslands Endemic to South Island high country	Hemimetabolous; univoltine; volant; eggs November–February	Adult overwintering stage Melanism	Adults active November–May	First instar: Saprophagous adults: Prey on small invertebrates	Transects counts; beating shrubs; sweep netting tussock
Lepidoptera Crambidae: Orocrambus spp. Grass moths 25± alpine species	0–1,800+	Open tussock grasslands; fell field.	Holometabolous	Cold-tolerant Microclimate inhabiting (wind avoidance behavior)	In several species females are flightless Males active flying on warm sunny days (November–February)	Poa and Chionochloa grasses	Transect counts; timed sweep netting; light trapping (dusk)
Erebidae (Arctiidae) “Tiger moths” (e.g., Metacrias huttoni)	600–1,800	Cocoons under flat rocks Larvae on native plants Adults open tussock lands	Holometabolous Male volant (female wings reduced) Colony formation	Melanism; sun basking; larvae very pilose	Diurnal, males are day-flying	Larvae feed on flowers and seed of Raoulia sp. and Wahlenburgia sp.	Transect sweep netting; timed observations (day); larval counts
Geometridae Approximately 90 spp. Aponotoreas spp. Asaphodes spp. Notoreas spp. Paranotoreas spp.	200–2,200	Tussock grasslands; snowbanks; cushion; herb fields	Holometabolous Geometrids are a significant part of the alpine biota Wide diversity of life history features including flightlessness, diet, crypsis, and circadian activity	Flightless females; pupal overwintering; sun basking; microclimate inhabiting (wind avoidance behavior)	Diurnal flying	Larvae may feed on species of Aciphylla, Anisotome, Pimelia, and Kellaria	Transect sweep netting; timed observations (day); larval counts; light trapping
Nymphalidae Argyrophenga spp. Percnodaimon spp. Ringlet butterflies	Argyrophenga 600–1,800 Percnodaimon 1,200–2,800	Snow grass and tussock grasslands High alpine screes; rock out crops; tors; summit ridges	Holometabolous Strong fliers	Cold tolerance Melanism (Percnodaimon)	Argyrophenga spp. active day flyers in tall tussock habitat Percnodaimon spp. Nap-of-the-earth flight above warm scree and rock May cross snow patches	Argyrophenga caterpillars feed on Chionochloa Percnodaimon feed on Poa.	Transects with 5-min timed counts Percnodaimon strongly correlated with scree and rocky terrain
Orthoptera Acrididae (Short-horned grasshoppers)	0–2,100	Open tussock country, short and tall sward	Multivoltine Three-year free-living stage Highly mobile Large geographic range Four genera and up to 14 species in Southern Alps Sex dimorphism, r-selection strategy	Adult overwintering diapause and quiescence Freeze tolerant	Diurnal May establish prolific densities in summer months	Herbivory of all stages except egg Grasses and herbs	Day transect counts; timed sweep netting
Orthoptera Anostostomatidae Giant, ground, and tree weta: Deinacrida spp. Hemideina spp. Hemiandrus spp.	600–2,500	Slab rocks; stable talus; stone pavements; low alpine vegetation	Paurometabolous Dens; harems (giant and tree) K-selection strategy, five to seven molts; free-living stage 2–4 years Highly mobile.	Freeze tolerant Lower lethal temperature of all taxa at least −5°C to −10°C	Deinacrida and Hemideina taxa nocturnally active (feeding and mating) Hemiandrus spp. probably crepuscular	Wide omnivorous diet	Night searching; time-series tracking tunnels; nonlethal pitfall traps; day quadrat rock turning
Rhaphidophoridae (Cave weta; e.g., Pharmacus; Mark 2012)	1,200–3,200	High alpine Rock tors; outcrops; slabs; scree; felsenmeers	Paurometabolous Adults and nymphs highly mobile Five to seven molts; often clusters; 1- to 2-year free-living stage	Melanism (partial) Subnival activity Highly freeze tolerant at least to −10°C, possibly cooler	Nocturnal and diurnally very active (pers. obs.) Congregate among north-facing rock outcrops	Omnivory; lichens; fungivory; canabalism	Night transects (if feasible) and searches; daytime rock turning; crevices; tracking tunnels if feasible
Arachnida Araneae (spiders) Agelenidae/Clubionidae/Huttonidae (e.g., Neoramia spp.)	600–1,800	Subalpine fell field; tussock; herb fields; stone retreats	Egg–instars–adult All stages except egg highly mobile	Probably freeze avoidance proteins coupled with behavioral warmth seeking Adults overwinter in hibernaculum	Diurnal (wolf spiders) All others nocturnal	Predators Almost any small invertebrate prey, including cave weta nymphs	Night observation (eye shine); pitfall traps (where feasible); rock turning; quadrats; transects; timed counts
Lycosidae Wolf spiders (e.g., Anoteropsis alpina and A. montana)	0–2,200	Fell field; herb field; tussock grassland; open scree; boulder/felsenmeers	Young carried on female abdomen for first molt (analogous to a K-selection strategy)	Significant melanism Sun basking Probably freeze avoidance proteins and seek warm microclimates	Diurnal hunting Very active over scree and rocky surfaces in search of prey Fast runners	Almost any small invertebrate prey, including cave weta nymphs	Day observation; stone turning; pitfall traps (where feasible)
Salticidae Jumping spiders Several undescribed alpine taxa	1,400–3,000	Fell field; herb field; tussock grassland; open scree	Salticids balloon (wind dispersal), yet most alpine taxa local endemics Data poor	Extensively pilose (moisture retention) and cryptic	Active daytime hunters Sun basking Swift and mobile	Probably spring tails and other small invertebrates	Pitfall traps (where feasible); rock turning; quadrats; timed observations
Opiliones (Harvestman) Triaenonychidae (Short-legged harvestman) Neopilionidae (long-legged)	0–1,800	Subalpine shrubland to upper vegetation zone Among litter and stable rock fissures (freeze–thaw fractured) High relative humidity	Few studies Silken retreats	Unknown; probably freeze tolerance	Possibly crepuscular Active when located in retreats during the day	Predators of micro invertebrates (mites, collembolan, etc.)	Litter/regolith sampling (tulgren extraction); night observation; occasionally collected in pitfall traps This group warrants further study in the alpine zone
Mollusca Sigmurethra: Punctidae/Charopidae (Micro or “dot” snails; e.g., Thermia sp., Charopa sp.)	200–1,800	Very high-humidity microhabitats Wetter regions near main divide	Very small snails (<10 mm Ø) Poorly described, many knowledge gaps	Alpine snails have close affinities with lower elevation taxa, if not a continuation Little information on cold adaptation physiology	Nocturnal/twilight zones of cavities	Detritivores Animal and plant matter	Quadrat; transect and hand searching by day; litter sieving
Athoracophoridae (Leaf-veined slugs; e.g., Pseudoneitea sp.)	600–1,900	Damp, complex, often rocky fern/Anisotome “grottos”	Hermaphrodites Adult overwintering, large animals	Congregational behavior; possibly thermoregulatory	Nocturnal Occasionally diurnal in wet conditions	Mosses, algae, and possibly lichens	Quadrat searches These are slow moving, potentially useful indicators of predator control

Notes. Biogeographic and taxonomic information sourced from (but not limited to): Bigelow 1967; Buckley and Simon 2007; Gibbs 1999; Koot 2018; Mark 2012; Marris et al. 2019; Richards 1972 and Trewick and Morgan-Richards 2019, unpublished Department of Conservation surveys, and observations by W. Chinn. Classification is insects, spiders (and kin), others.

^aCan withstand ice formation in body tissues (Sinclair 1999).

Discussion

The snow line monitoring program shows a dramatic reduction of ice volume throughout New Zealand’s mountain glaciers (T. J. H. Chinn et al. 2012; Mackintosh et al. 2017; NIWA 2020). As time elapses, the data set becomes increasingly reliable, and today the study represents one of the longest proxy climate records from within the Southern Alps. In this article, we have explored the ecological utility of the snow line data as a complementary measure of climate forcing on alpine communities (Willet 1950; Fleming 1963; Wardle 1971; Burrows 1990; Mark 2012).

During the 1990s, the climate control paradigm saw renewed interest as DNA sequencing tools provided fresh estimates of population divergence times, particularly when coupled with molecular genetic clocks (Avise 2000). When applied to New Zealand’s alpine invertebrates and plants, genetic analysis has shown that long-term habitat tracking has been the rule, especially during the Pleistocene and early Holocene climate fluctuations (Gibbs 1999, 2016). However, few studies have provided rate functions for these elevation shifts, despite this being a necessary prerequisite for the tracking model to hold. More pressing is today’s current rate of warming (at least 1°C per century), driving an unprecedented rise of the snow line (by at least 3.7 m a⁻¹). At this rate, the potential to force alpine populations into local extinction under the DREM is considerable if habitat tracking holds true.

Did we answer our questions?

Our study showed that all taxa above 1,800 m a.s.l. had indeed tracked upslope in the preceding thirty-eight years of snow line rise (Figure 12a and Figure 12b). Below 1,800 m a.s.l. tracking results were calculated values only, and there was considerable rate variation, which is open to interpretation. Small groups of high-elevation taxa (cave weta in particular) tracked within a narrow isotherm band (or at least colonized refugia recently exposed). For these taxa, measured tracking rates were at least as fast as the upslope shift of the snow line (approximately 3 m a⁻¹). Meanwhile, at lower elevations the calculated shift rates appeared larger, showing higher vertical differentials with a maximum upslope rate around 16 m a⁻¹, shown in Figure 12a and Figure 12b as a time-dependant convergence between actual and calculated taxon elevations. We interpret the shape of these graphs as hysteresis curves, also known as “alternative state theory” (Lewontin 1969).

In our opinion, the convergence situation describes an ecological lag function (a hysteresis; Doering et al. 2007; Bader, Rietkerk, and Bregt 2018). Although speculative, we suggest that lower elevation populations have a wider thermal tolerance than those at higher elevations (near the snow line). In our case, the area between the actual and calculated population (median densities) indicates “work done” (Vertical movement × Number of taxa) by the entire community at specific elevations. The degree of population lag (the graphed area between past and present) apparently depends on the rate and extent of the vertical shift for a given isotherm. Comparing population trends between sites is perhaps best interpreted using integrals for the two time curves (T1 and T2), per site; for example,

${\int }_{T_1}^{T_2}\left[T2(x)-T1(x)\right]dx.$

We estimate that the alpine hysteresis curve will persist as the snow line continues to rise. With further warming, the area of the curve may become deep and wide as more taxa (including introduced) are able to make a greater elevation shift between any given T1 and T2 position. On that basis, alpine ecosystems may never return to the same population-elevation configuration as they were at T1 even in the unlikely event that temperatures were to again cool in the next few centuries.

The Southern Oscillation and CMUs

Our analysis showed a relationship between the SOI and EOSS elevation, consistent with the findings of Purdie et al. (2011) and Salinger, Fitzharris, and Chinn (2019). In these studies, ENSO-driven temperature was a stronger influence on snow line elevation than precipitation (Mackintosh et al. 2017). Although we are aware of the very close relationship between temperature and ice volume change in the Southern Alps, other modes of climate variability also exert significant influence on regional temperature patterns (Mackintosh et al. 2017). Therefore, we caution reliance on the SOI for anything other than investigating ecological oscillations in the alpine habitat and broader climate studies.

When combined with several regional climate indices, the SOI will be an important part of EMU and CMU data. For instance, during El Niño (negative SAM) conditions, alpine invertebrate communities may be depressed into lower elevations and potentially more prone to rodent predation. The reverse is true for La Niña years (positive SAM), because warmer temperatures and less precipitation favor a snow line rise (and a vertical “stretching” of the alpine zone). In all cases, it is advection of temperatures onto New Zealand from several thousand kilometers that appears to be the most important factor controlling Southern Alps air temperatures (Mackintosh et al. 2017). In this light, even during cool years (El Niño), Tasman Sea temperatures may be warmer than normal.

During the austral summer of 2017–2018, New Zealand experienced an unprecedented marine and atmospheric heat wave. The regional heating covered about 4 million km², with air temperatures across New Zealand higher than any recorded in over a century (BoM and NIWA 2018). A consequence of this warming saw the highest snow line elevation since records began (T. J. H. Chinn, pers. obs.). This illustrates the impact of extreme temperature swings in the alpine environment and highlights the urgency to establish baseline ecosystem information from a series of CMUs. We suggest that alpine CMUs should measure a suite of climate variables in conjunction with population elevations. Because we do not yet know the critical area and density values that may herald imminent population extinction we should, at least, establish a datum and monitoring system. The goal would be to elucidate the ratio between habitat area (frustum), invertebrate population densities, and their thermal parameters.

Similarly, CMUs could provide useful sites for investigating linkages between snow line isotherms and the “delta-T” model, which correlates summer threshold temperatures with mast seeding events (Kelly et al. 2012). Differential mast seeding effects between beech- and non-beech-forested mountains along the Southern Alps could also be investigated within CMUs. Finally, incorporating CMUs into alpine ecosystem management could enable cross-links between ecologists and alpine climatology groups, including the Snow and Ice Research Group, the Southern Alps End of Summer Snowline project (NIWA), The New Zealand Climate Change Research Institute (Victoria University, Wellington), and the Alpine Ecology Research Group (Otago University).

Caveats and considerations

We highlight several caveats with our proposal. We are not stating that snow line trends are a definitive predictor of alpine community response to climate change but rather that they present a measurable relationship between the two. We also caution the geographic scale of our method, given the minimal signal-to-noise ratio of the data from individual index glaciers: More is better, and any model needs to be a composite of several sites. Therefore, future sampling biases will require mitigation or correction using a random effects model, for instance (not included in our study).

Limitations of the EOSS

Our model was confined to alpine regions with actual EOSS data. For mountains located further north (warmer), model extrapolation is necessary. Snow line data for the Kaikoura Mountains could be coupled with the North Island volcanoes, providing a datum slope for modeling population elevations. High, steep-sided volcanoes are also prime candidates for future work on this topic (see Heine 1963). Mountainous areas below the snow line but above the timberline include the Paparoa Ranges, Tasman Mountains, and the Richmond Mountains.

At the single-mountain scale, snow lines are inherently convoluted, biased by aspect and often involute. The snow line is lower on shady slopes than on warmer slopes and often absent on sunny aspects. The resulting thermal niche-scape will describe an undulating surface embedded within the physical landscape that supports “pockets” of suitable habitat within otherwise unsuitable terrain of ice and snow.

Acquiring large-scale population-elevation data will therefore require two scales; individual mountains and regional averages combined, possibly using cluster analysis of snow line data using methods similar to Kidson (2000). Slope angle is important, too: The optimal is between 10° and 50°. Too shallow and the snow line elevation gain is negligible and insufficient for measuring biotic shifts; too steep and snow will avalanche off, sweeping any biota clear. At the regional scale (100 km²), snow line contours are less convoluted while displaying a distinct relationship with the west-to-east precipitation gradient (T. J. H. Chinn 1995).

Density-dependent effects

Our model assumed that population densities were stratified by temperature. As poikilotherms, most invertebrates have thermally cued life stages, and dispersal is limited by lower lethal temperatures. However, temperature is a density-independent variable, and other drivers of density may be operating in the alpine zone (including positive and negative effects). For example, in the resource-constrained alpine environment, the Allee effect may be significant.

The Allee effect is a positive correlation between population density and individual fitness; a good example is mate encounters—the more encounters, the higher the intrinsic growth rate and density—in contrast to the resource competition model (Allee and Bowen 1932). If the Allee effect is operating on the alpine biota (particularly free-living taxa like grasshoppers, weta, and cockroaches), the intrinsic rate of growth would be skewed against population density, and this should hold true for all elevations that alpine taxa occupy. Measuring population growth rates would be one method of determining which factor, biological or physical, has the greater influence on density.

The Allee effect may also combine with climate forcing to influence rates of population range expansion. For example, a climate-forced range expansion may incur increased genetic diversity as the population density rises within the colonizing leading edge (Tobin et al. 2007; Garnier, Roques, and Hamel 2012). Alternatively, models show that very rapid climate shifts may reduce genetic diversity at the leading edge of a colonizing population (Garnier and Lewis 2016), a probable situation in the Southern Alps as rapidly rising snow lines continue to force alpine invertebrate populations upslope.

A third consideration for alpine invertebrate populations is mammalian predator incursions (rodents, mustelids, and marsupials). Predators may target rare invertebrates in favor of common taxa, causing extinction of the rare compared with populations of common species, with potentially significant effect on the structure of alpine invertebrate communities. Target invertebrates may include energetically lucrative individuals such as large weta, a situation analogous to the “anthropogenic Allee effect” where rare taxa are disproportionately exploited, rendering them even rarer (Courchamp et al. 2006). These uncertainties (and contradictions) about population density highlight the infancy of alpine climate change ecology and the requirement for concerted research around this emergent discipline.

Finally, we note that high-elevation invertebrates apparently trade off beneficial traits (e.g., cold physiology) against various niche constraints. Thus, invertebrate taxa with hemimetabolous life histories, cold-tolerant physiology, omnivorous diets, and all-season mobility are well scaled to life at high altitude in New Zealand. The alpine habitat is also an environment absent of many predators, but the cost is a narrow envelope in which survival may be seriously compromised if climate warming continues at current rates.

Satellite populations

Our invertebrate sampling ignored populations on snow-free exposed rocky tors above, or adjacent to, the transect lines. In most cases peripheral populations are also subject to an upslope rise of isotherms and the snow line remains the only isotherm proxy in the alpine system. However, without establishing long-term alpine transects and data sets, information from any aspect is spurious. Our method remains useful because the upper boundary of the AIB is always coupled to a measurable thermal envelope with a known time series.

High mountains and snowbank ecology

In some cases, seasonal snow cover may actually afford survival benefits, especially for alpine plant communities (Lord et al. 2017). This is an important characteristic that should be included in the enumeration of herbivorous invertebrates associated with snowbank communities. Although a reduction of snowbank communities is another undesirable result of climate warming in alpine zones, the spatially fixed nature of these communities represents an in situ adaptation study and not a case of isotherm tracking.

Conclusion

We have shown that alpine invertebrate communities may be tracking upslope under a warming atmosphere in New Zealand. What is not certain, however, is the “rate and fate” of this movement for high alpine communities, and this alone is an important direction of future study. Predicting the effects of climate change on alpine ecosystems is complex and multidisciplinary. Data and concepts are required from meteorology, climatology, hydrology, glaciology, and geomorphology coupled with ecology, evolutionary biology, and physiology. Amalgamating these disciplines is worth the effort given the considerable conservation and scientific value.

The challenge is to combine reliable and quantitative indicators of climate change with alpine community indices, within the same location. Here, we have sought to establish a linkage between a layer of ambient seasonal temperature across the Southern Alps (the EOSS) and the elevation of alpine invertebrate populations. Though there may be other more complex and revealing dynamics between alpine ecology and climate, yet to be explored, we encourage their investigation under the CMU concept. Of the many pressing climate issues for conservation science to tackle, alpine ecology remains incomplete. Indeed, since Humboldt’s Chimborazu zonation study, alpine ecology seemingly languished, yet at no time since Humboldt has there been such an unprecedented rate of climatic forcing on mountain ecosystems.

Acknowledgments

We are grateful to the following people for their valuable comments and critique of the manuscript: Bob Frame, Andrew Lorrey, Jenny Christie, Kerry Weston, Colin O’Donnell, Petra Pearce, Barbara Chinn, and external reviewers. Many thanks to numerous Department of Conservation staff, including James Holborow, Mark Fitzpatrick, Dong Wang (statistical assistance), Andrew Grant (RStudio), Richard Earl (GIS expertise), Eric Edwards, and Joanne Monks. For field assistance and encouragement, we are grateful to Joanne Kane, Sylvia Thorpe, and Roger Parkyn. We acknowledge the New Zealand Department of Conservation for technical support.

Additional information

Funding

The authors acknowledge the National Institute of Water and Atmosphere (NIWA) for materials, data, and funding for the End of Summer Snowline program.