Growing-Season Temperature Change across Four Decades in an Arctic Tundra Pond

We examined temperature dynamics across a 42-year period in a low-centered tundra polygon pond on the Arctic Coastal Plain of northern Alaska to assess potential changes in thermal dynamics for ponds of this type. Using water temperature data from a pond near Barrow (now Utqiaġvik), Alaska, studied intensively during 1971 – 73 and again in 2007 – 12, we built an empirical model coupling historical air temperatures to measured pond temperatures for four summers. We then used the model to predict summer pond temperatures over a 42-year span, including 1974 – 2008, for which direct aquatic temperature records do not exist. Average pond temperatures during the growing season (1 May through 31 October) increased by 0.5 ̊C decade-1 or 2.2 ̊C over the 42-year period. Our simulations predicted the average date of spring thaw for the pond as 2 June (± 3 d), which did not change over the 42-year time period. However, average pond temperature during the first 30 days of the growing season increased from 1971 to 2012, suggesting that recently, ponds are warmer in early spring. The average date of pond sediment freeze over the 42 years shifted later by 15 days, from 28 September in 1971 to 13 October in 2012. These changes correspond to a growing season that has increased in length by 14 days, from 118 days in 1971 to 132 days in 2012. Contemporary temperature measurements in other shallow tundra ponds in northern Alaska show a high degree of temporal coherence (r = 0.93 – 0.99), which warrants the general conclusion that tundra ponds on Alaska’s Arctic Coastal Plain have undergone a significant change in thermal dynamics over the past four decades. Our results provide a means to incorporate these pond types into larger-scale simulations of Arctic climate change.


INTRODUCTION
Approximately 30% of Alaska's Arctic Coastal Plain consists of fresh surface water, with much of this total contributed by mosaics of small tundra ponds. A prominent pond type occupies low-centered polygons that form in drained lake basins because of thermokarst activity (Britton, 1966;Hinkel et al., 2005). The Arctic Coastal Plain is a landscape low in relief and underlain with continuous permafrost that restricts groundwater flow to levels near the surface. In some regions, tundra ponds are becoming desiccated as evaporation increases in response to a combination of increasing water temperatures and vapor pressure deficits, with consequent impacts on the number, size, distribution, and limnology of these ponds (Schindler and Smol, 2006;Lilijedahl et al., 2011). Tundra thaw ponds in ice-rich permafrost are a dynamic feature of the land surface on the Arctic Coastal Plain: ponds often increase in size through thermal erosion and coalescence before draining and reforming, in a cyclical pattern occurring over centuries (Smith et al., 2005;Jones et al., 2011).
In the Arctic, climate is changing at a faster rate than in temperate areas . Satellite imagery shows some of the highest temperature anomalies, relative to long-term averages, in the North American portion of the Arctic, with estimates of recent temperature increases averaging 0.60˚C per decade (Comiso and Hall, 2014). Earlier melt dates for sea ice, longer periods of abovefreezing temperatures, warmer average temperatures, greater temperature variability, and changes in precipitation patterns are all manifestations of climate change in the Arctic (Rouse et al., 1997). While long-term records exist for atmospheric (e.g., Brohan et al., 2006) and marine temperatures (e.g., Rayner et al., 2006) and for freshwater lakes (e.g., Schneider and Hook, 2010), fewer data exist for Arctic freshwater ponds. Because such small aquatic habitats represent more than 95% of total waterbody coverage on the tundra (Muster et al., 2013), this lack of data limits our understanding of how Arctic tundra ecosystems may be affected by climate change and may themselves influence climate dynamics.
There is a growing literature on tundra pond hydrology in the context of large-scale wetland complexes (e.g., Rovansek et al., 1996;Woo and Young, 2006;Assini and Young, 2012;Liljedahl et al., 2012), yet little is known about the thermal behavior of individual ponds. Historical temperature data for tundra ponds are sporadic-even for well-studied systems. Unlike lakes, which can be studied historically from satellite imagery, individual ponds on the landscape are too small for similar treatment given the resolution of most available imagery (Schneider and Hook, 2010;Muster et al., 2013). An indirect method to reconstruct thermal properties of ponds is necessary if data are missing. Air temperatures are readily available from national weather stations, and pond temperatures, like air temperatures, respond to solar radiation, cloud cover, and movement of air masses. Thus, we should expect a strong relationship between air temperatures and pond temperatures during the thaw period, especially in these shallow (< 1 m deep) freshwater ponds on the Arctic Coastal Plain that are well mixed by wind.
For this study, we built an empirical model of growingseason pond temperatures based on both historical (1971 -73) and recent (2007 -12) tundra pond temperature data to test the hypothesis that pond temperatures have warmed over the past four decades. Our model simulates pond temperature dynamics for the period 1974 -2006, for which no pond temperature data exist. Coupled with existing early and recent pond temperature data, model results help to characterize pond thermal behavior across a 42-year period. Understanding such long-term thermal dynamics is crucial to interpreting biological changes in these Arctic ponds (Lougheed et al., 2011;Andresen et al., 2016;Braegelman, 2016).

Study Site
We assessed temperature patterns in "Pond C," one of several tundra ponds near the village of Barrow (now Utqiaġvik), Alaska (Fig. 1), that were studied intensively during the International Biome Program (IBP) in 1971 -73 (Hobbie, 1980). Some of these IBP ponds, including Pond C, have also been studied more recently (Lougheed et al., 2011;Andresen and Lougheed, 2015;Andresen et al., 2016). Pond C (71˚17′40.43″ N, 156˚42′8.57″ W) is typical of ponds occupying low-centered polygons on the coastal plain tundra near Barrow. Water covers approximately 1500 m 2 of the Pond C polygon. Dense emergent Carex aquatilis occupies the peripheral 80% of this wetted area, with some Ranunculus sp. and Arctophila fulva at deeper depths. A central open-water region of approximately 300 m 2 , spanning depths of 25 -40 cm, is largely void of macrophytes. The exposed surficial sediments consist of coarse, unconsolidated particles with high organic content . Fish are absent, but during the growing season various waterbirds, including Red Phalaropes (Phalaropus fulicarius), Long-tailed Ducks (Clangula hyemalis), and Steller's Eiders (Polysticta stelleri), feed on the abundance of planktonic and benthic invertebrates.

Temperature Data Acquisition
All pond temperature data were collected using field loggers that recorded temperature at the pond-sediment surface every hour. We use "pond-sediment temperature" and "water temperature" as synonyms because previous research has shown that temperatures at the sediment surface and in the water column are closely related (Danks, 1971;Hobbie, 1980). Water temperatures were monitored in IBP Pond C during 1971 and 1973 and in adjacent Pond B in 1972 from shortly after pond thaw through the third week of August each year. We assume that Ponds B and C exhibit similar thermal behavior, given that the two ponds are within meters of each other and are of similar size and mean depth. Comparative monitoring of these two ponds during the early 1970s showed that temperatures in Pond B and Pond C differed by at most 0.5˚C for readings made at the same hour on the same date . These 1971 -73 temperatures were measured with a telethermometer placed on the pond sediment at water depth of 20 -30 cm, which recorded sedimentsurface temperatures hourly on a potentiometric recorder (Stanley, 1974). During 2007 -12, Hobo Pro ® waterproof temperature loggers with ± 0.2˚C accuracy collected hourly temperatures from the same location in Pond C year-round, beginning on 13 June 2007. Atmospheric temperatures taken every three hours were extracted from the online database of the National Oceanic and Atmospheric Administration for Barrow, Alaska (Station ID: PABR).
We consolidated pond temperatures to match time periods when air temperatures were obtained (i.e., every three hours) by deleting pond temperature observations without corresponding air temperatures (e.g., if pond temperatures were available at 0800 but air temperatures were not, we deleted the 0800 record from our pond temperature records). The result was that data were paired by daily 3 h intervals, which we then averaged into daily records. All time series were incomplete, but records for the later years (i.e., 2007 -12) were more nearly complete than those for the early years (1971 -73). Further, records were more nearly complete for the middle of the growing season (late June to late August), with most of the missing data occurring near thaw and freeze dates for the ponds. We analyzed only temperatures during the short growing season, spanning the time period from pond thaw to pond freeze. Thaw date for the pond was determined as the day during the year when pond-sediment temperature rose and remained above 0˚C, while freeze date was determined as the day when pond-sediment temperatures dropped below 0˚C and did not rise above freezing again that year.

Analytic Methods
Approaches to predicting water temperatures from air temperatures range from purely mechanistic (e.g., Mendez et al., 1998) to a variety of phenomenological models, including sinusoidal functions for long-term temperature patterns and air-water linkages through Markov chains (e.g., Cluis, 1972), Box-Jenkins transfer functions (e.g., Caissie et al., 1998) 2010), and Fourier functions (e.g., Kothandaraman, 1972). Mechanistic heat-budget type models can provide greater insight into temperature dynamics than do empirical models; however, these types of models can be complex and difficult to parameterize. Consequently, we did not take a mechanistic approach, but rather used an empirical approach related to that of Caissie et al. (1998) for post-thaw conditions in our focal pond.
For modeling purposes, we divided temperature data into three separate sets: (1) a calibration set that included years 1971 -72 and 2008 -09, which was used to fit the model; (2) a validation set that included years 1973 and 2007, which was used to verify the model; and (3) a prediction set that included years 2010 -12, which was used to determine how well our model could predict pond temperatures using only available air temperatures.
Using only the calibration data set (i.e., 1971 -72 and 2008 -09), pond-sediment temperatures (P T ) were separated into two components, a long-term annual component (P L ) and a short-term component (P R ) composed of the residuals, which together took the form: The long-term annual portion for the pond temperatures (P L ) was fit using a Fourier regression of the form: where β 0 , β 1 , A j , and B j are regression coefficients; φ is the number of days since 1 January 1971 = 1; j is the order of the Fourier series; f is the period, which was fixed at 365.25 days; and t is the Julian day within the year, where 1 January = 1 for each new year. The coefficient β 1 accounts for a linear increase in the pond temperatures between 1971 -72 and 2008 -09, while controlling for the regular seasonal rise and fall of temperatures. The number of Fourier components was determined to be the fewest parameters in the calibration model that led to a nonbiased fit to remaining data (i.e., the validation set : 1973, 2007; prediction set: 2010 -12). Bias was assessed by regressing observed temperature data against modeled temperature data generated by the calibration data function, then checking for an intercept not significantly different from zero and a slope not different from unity (Piñeiro et al., 2008).
The residuals (P R ) from Equation 2 were then regressed against residuals from a Fourier model fit to the air temperature, which were obtained in the same way as pond temperatures in Equation 1: where [Eq. 4] and α 0 , α 1 , C k , and D k are regression coefficients, k is the Fourier order, f is the period, and t is time in days. The extracted residuals from the fitted air model (i.e., A R (t)) were used to generate a Markov model (i.e., a series of lagged air temperature residuals) to fit the short-term residual pondsediment temperatures in the form [Eq. 5] where a 0 , b 0 , and b k are regression coefficients, and l is the order of the lag for air temperature residuals against pond-sediment residuals (Cluis, 1972). We then tested the hypothesis that β 1 = 0 (i.e., there is no linear increase in the temperature conditional on the season effects).
Once we had established an unbiased model, we generated predictions for summer pond temperatures during 1971 -2012 using Eq. 1. These predictions were compared to observed temperatures in 1971 -73 and 2007 -12 and were used to simulate pond temperatures during 1974 -2006 for which measured values were missing. We estimated thaw and freeze dates, duration of pond-sediment temperatures above 0˚C, and the annual cumulative temperature for the pond (available degree days). Degree-days (DD) are a proxy both for the potential for permafrost warming (Ling and Zhang, 2003) and for biological activity (Trudgill et al., 2005). DD is the product of the number of days a pond is thawed and the average temperature over that period and summarizes changes in heat balance over an entire growing season. We derived a measure for the rate of spring warming following pond thaw: the average daily temperature over the first 30 days post-thaw. Finally, we compared contemporary Pond C temperatures with temperatures in other local and regional tundra ponds that we have monitored comparably, using correlation tests to see how temporally coherent temperatures are among ponds on the Arctic Coastal Plain.
All statistical analyses were performed with R 2.12 program (R Core Team, 2010). Where appropriate, we report 95% confidence intervals instead of null hypothesis tests and associated probability values, since the former are more informative than the latter (Gardner and Altman, 1986).

RESULTS
Our calibration model (based on 1971 -72 and 2008 -09 data) was built using two Fourier coefficient pairs for pond temperature and a single Fourier coefficient pair for air temperature, and it produced a strong fit to observed data (R 2 = 0.81; RSME = 1.6; Fig. 2A). The best fit for pond temperature residuals against air temperature residuals was the air residual at time t plus one day lagged (i.e., t minus one day); the overall model coefficients are given in After fitting the model with the calibration data set, we applied the same model to the validation and prediction data and found good fit between predicted and observed data for all years. The model's fit to the validation set (R 2 = 0.78; RSME = 1.7; Fig. 2B) was quite similar to the calibration fit. There was no apparent bias for either the intercept (95% CI = −0.58, 0.61) or the slope (95% CI = 0.86, 1.01) when regressing the observed validation pond temperatures against values predicted by the calibration model (Fig. 2B). The model's fit to the prediction data set (2010 -12) was even stronger (R 2 = 0.89; RSME = 1.3; Fig. 2C). Overall, across all 12 years for which pond data exist, correlations between observed and predicted pond temperatures were high ( Fig. 3; r = 0.89, p < 0.001, n = 541).
We estimated dynamics of thaw/freeze events and DD over the period from 1971 to 2012. Despite an apparent trend toward earlier spring thaw, we found no significant   Fig. 5A). We also detected a positive relationship between this "spring" pond temperature and the date of pond thaw, with an increase of 0.1˚C in average temperature during the first 30 days post thaw for every day later that thaw occurred (p < 0. 01; Fig. 5B). There was a high temporal correlation in daily temperatures between Pond C at Barrow and other ponds on the Arctic Coastal Plain that were monitored over similar periods. These correlations ranged from 0.95 to 0.99 for three other ponds at Barrow, Alaska, and 0.93 for a single pond monitored near Atqasuk, approximately 100 km inland (Table 2).

DISCUSSION AND CONCLUSIONS
We successfully modeled missing pond water temperatures using existing pond and atmospheric temperature records. Model fits to data were strong: models explained more than 80% of the variability in pond temperatures. Given that we were not able to control for potential differences in precisely where the thermistors were placed within the ponds, or for use of different thermistors between the early and late periods, we consider the amount of variability explained by our pondtemperature model to be significant.
Pond sediments became warmer in the 42 years from 1971 to 2012, with modeled mean daily temperatures increasing by 2.2˚C (95% CI = 1.8, 2.6). The predicted thaw date of pond sediments did not change significantly over the 42-year period, while the date when pond sediments froze became later, so that the annual ice-free period has lengthened by more than 12 days. Taken together, these findings reflect a 33% increase in the number of DD and suggest an increase in the potential for greater pond-bed warming, enhanced permafrost thawing, and increased temperature-dependent biological activity each season. Barrow has continuous permafrost approximately 400 m deep, with an active layer of thaw ranging from 20 to 50 cm during the short summers. Deepest areas of thaw occur beneath ponds, with a shallower active layer in the surrounding grassy tundra (Brown and Sellmann, 1973;Hobbie, 1980). Recently, Streletskiy et al. (2017) have determined a deepening of thaw at Barrow over the period from 1962 to 2015, which would be consistent with overall pond warming.
The number of coefficients selected on the basis of parsimony for the Fourier portion of the model may provide insight into some of the mechanisms driving the seasonal temperature pattern. Because of solid-liquid phase shifts of water during melting at the beginning of the summer and freezing in the autumn, we should expect a truncated sinusoidal wave at the beginning of the year and a longer tail in temperatures at the end of the year, reflecting many false starts to the freezing of pond sediments. This is a pronounced feature of Arctic wetland soils and consistent with the so-called "zero curtain" effect described by Outcalt et al. (1990). Indeed, during years for which we had pond data during thaw and freeze events, we saw a rapid increase in pond temperature once the zero-degree threshold was crossed in May -June; at the end of the summer, the ponds stayed above freezing longer than might have been predicted using a simple harmonic function. In contrast, we would expect air temperatures to generate a more symmetric response, since no similar phase shift need be accounted for. Our most parsimonious air temperature model included a single Fourier coefficient, but two Fourier coefficients were selected in our best model of pondsediment temperature. This difference suggests that our models have captured the contrasting patterns between airtemperature and pond-temperature dynamics.
Over the period 1975 to 2016, the date for the disappearance of snow cover at Barrow advanced by approximately 0.27 days per year (p < 0.05). This rate, if applied to our 42-year period of study, equates to an average snowmelt that is 11.3 day earlier (Cox et al., 2017). Such advanced spring melt dates are consistent with patterns reported for other Arctic sites (e.g., Stone et al., 2002;Lynch and Brunner, 2007). Initial snow loss from the surface of pond ice and the surrounding tundra lowers the albedo of both aquatic and terrestrial habitats, allowing increased transmission and absorption of solar energy by dark sediments and early meltwater. Heat transfer via inflowing terrestrial meltwater further degrades remaining snow and ice in tundra pond basins. Exposed patches of organic sediments and stained water have a strong capacity to absorb additional heat, promoting rapid thaw of residual bedfast ice and subsequent warming of both water and sediments.
Spring thaw initiates a period of intense biological activity in these ponds. Date of thaw and subsequent water FIG. 4. Model-generated annual thermal dynamics in Pond C at Barrow, Alaska, over 42 years (1971 -2012). Shown are dates of (A) sediment thaw in spring and (B) freeze in autumn, (C) duration of unfrozen sediments, and (D) annual growing degree-days.
temperatures both influence the seasonality of ecological processes ranging from algal production rates and plant growth to zooplankton population dynamics and the timing of insect emergence (Stanley, 1974;Butler, 1980;Hobbie, 1980). Braegelman (2016) found an advance of adult emergence by chironomid midges in Barrow ponds during 2009 -13, relative to records from the mid-1970s. For one pond studied in both eras, the 14 most abundant midge species emerged an average of eight days earlier in 2009 -13 than in the 1970s. Braegelman attributed 70% of the variation in taxon-specific emergence timing during 2009 -13 to "pond" and "year" effects. These effects likely are driven by combinations of differing thaw dates and subsequent DD accumulation rates in different ponds and years.
To better understand the ecological implications of climate warming in these ponds, it would be useful to separate biological effects of variation in thaw timing from effects of subsequent water temperatures. Our pond temperature model failed to consistently predict actual thaw dates in Pond C when thaw dates were known from loggers, although the number of data points used to make the comparison is small (n = 5). This failure may reflect annual variation in local snow accumulation and melt, variables that are poorly correlated with air temperature at the scale of an individual pond.
We suspect that thaw date for pond sediments might be better predicted using some variable associated with snowmelt. Zhang and Jeffries (2000) found lake ice thickness to be much more strongly coupled to snow depth than to air temperatures. A strong connection between iceout dates and snow depth has been shown for lakes (Duguay et al., 2003;Surdu et al., 2014), and the same is probably true for ponds too, since early spring energy input must first be devoted to melting snow before pond sediments and water can be heated (Assini and Young, 2012). Also, we suspect that the color and composition of pond sediment have a strong influence on pond temperatures, but little research has included sediment variables as factors in pond temperature dynamics (MacKay et al., 2009).
Our model-generated trend suggesting earlier pond thaw at Barrow (Fig. 4A) is consistent both with reports of earlier ice-out dates in IBP ponds during 2010 -13 relative to the early 1970s (Andresen et al., 2016) and with the trend toward earlier snowmelt in Barrow cited by Stone et al. (2002) and Hinzman et al. (2005). Even without a significant change in melt date, our analysis reveals warmer water temperatures during the first month after pond thaw.
The ice-free season in these tundra ponds may also be extended by later freezing in the autumn. This possibility is consistent with warmer summer atmospheric conditions and decreases in precipitation, which should increase TABLE 2. Temporal correlations in average daily temperatures between our focal pond (Pond C) and other ponds on the landscape with similar morphometry (Area < 3000 m 2 ; Z max < 1 m) for periods during the year when pond temperatures were above freezing (i.e., > 0˚C). Distance represents the straight-line distance between the centroid of each pond and Pond C. The sample size (n) is the total number of hours used in each comparison, and the Pearson correlation coefficient (r) is provided as a measure of temporal coherence for pondsediment temperature relationships. evaporation-to-precipitation ratios, likely impacting solute concentrations in the ponds. Following the spring flush of meltwater, evaporation dominates water balance in ponds of this type, leading to increases in ion concentrations over the summer. In Arctic ponds such as Pond C, conductivity can more than triple from the beginning to the end of the growing season . When present, such solute effects on pond freezing do not carry over from the prior year to affect thaw date because of a resetting of hydrochemistry during the melt period (Thompson and Woo, 2009). Many Arctic ponds farther north on Ellesmere Island have shown a trend toward earlier seasonal drying and increased conductivity (Smol and Douglas, 2007). While altered pond hydrochemistry may play a role in the timing of freeze events, we suggest that pond thawing at Barrow is primarily influenced by solar energy. Solar inputs act both directly on pond ice and snow and indirectly through the release of warm meltwater from tundra uplands.
The effect of such meltwater inputs on ice melt varies among ponds, depending on landscape microtopography (Liljedah et al., 2012(Liljedah et al., , 2016, which contributes some local variability to pond thaw dates. Pond-sediment temperature dynamics over the period of this study are generally consistent with other climate data from the Arctic Coastal Plain, particularly at Barrow, Alaska. For example, atmospheric temperatures at Barrow have increased by 0.29˚C decade -1 from 1940 to 2000 . Our trend for pond-sediment temperatures showed nearly twice that rate, at 0.5˚C decade -1 from 1971 to 2012, which more closely matches the 0.60˚C decade -1 rate of temperature increase observed in Toolik Lake from 1975 to 2000 (i.e., 1.5˚C over 25 years at 1 m depth) .
Not only are the ponds likely influenced by climate change, but climate change is also likely influenced by tundra ponds. Albedo conditions provide feedbacks in the climate system, and sea ice dynamics are often cited for their role in the polar amplification of climate change, in large part through the influence of snow on albedo. However, given the abundance of freshwater lakes, ponds, and rivers on the Arctic Coastal Plain, the duration of the open-water season for melt ponds and other freshwater features on the landscape also has the potential to strongly influence albedo effects in feedback to the global climate system (Curry et al., 1995;Holland and Bitz, 2003;MacKay et al., 2009). Ponds of this type have long been considered potential sources of methane to the atmosphere as they warm; recent work on Pond C has shown that over the 42-year duration of this study, methane flux from individual ponds has increased by approximately 60% (Andresen et al., 2016). As the inclusion of pond effects may improve the accuracy of current global circulation models (Bowling and Lettenmaier, 2010), effects of climate change on these freshwater pond systems warrant further research.
In conclusion, our analysis of temperatures in a tundra pond on the Arctic Coastal Plain indicates a general warming trend throughout the growing season during the period from 1971 to 2012, as well as a lengthening of the ice-free period. The rate at which ponds, once thawed, warm in the early spring has increased over the study period. Our analysis of recent and historical data back to 1971 indicates that freeze dates have occurred later in the season in recent years. These findings agree in general with other research on surface air temperatures and other available findings for freshwater systems in the Arctic, a region shown to be more sensitive to climate forcers than lower latitudes.