On the spatial and temporal shift in the archetypal seasonal temperature cycle as driven by annual and semi-annual harmonics

03/15/2020 ∙ by Joshua S. North, et al. ∙ University of Missouri 0

Statistical methods are required to evaluate and quantify the uncertainty in environmental processes, such as land and sea surface temperature, in a changing climate. Typically, annual harmonics are used to characterize the variation in the seasonal temperature cycle. However, an often overlooked feature of the climate seasonal cycle is the semi-annual harmonic, which can account for a significant portion of the variance of the seasonal cycle and varies in amplitude and phase across space. Together, the spatial variation in the annual and semi-annual harmonics can play an important role in driving processes that are tied to seasonality (e.g., ecological and agricultural processes). We propose a multivariate spatio-temporal model to quantify the spatial and temporal change in minimum and maximum temperature seasonal cycles as a function of the annual and semi-annual harmonics. Our approach captures spatial dependence, temporal dynamics, and multivariate dependence of these harmonics through spatially and temporally-varying coefficients. We apply the model to minimum and maximum temperature over North American for the years 1979 to 2018. Formal model inference within the Bayesian paradigm enables the identification of regions experiencing significant changes in minimum and maximum temperature seasonal cycles due to the relative effects of changes in the two harmonics.

READ FULL TEXT VIEW PDF
POST COMMENT

Comments

There are no comments yet.

Authors

page 6

page 14

page 15

page 16

page 18

page 19

This week in AI

Get the week's most popular data science and artificial intelligence research sent straight to your inbox every Saturday.

1 Introduction

In ecology, “spatial synchrony” is a concept that describes coincident in time variations in an ecological process across geographically separated populations (Liebhold2004). In many cases, this synchrony leads to symbiotic relationships that are tied to environmental seasonal cycles. For example, an important climate-driven issue facing forest health concerns native bark beetle infestations, with current beetle outbreaks among the most severe in recorded history (Bentz2010). Historically, exposure to very cold temperatures is necessary to control the beetle population, but increasing seasonal minimum temperatures in northern latitudes has disrupted this historical synchrony, limiting beetle mortality, and increasing tree mortality. Other examples of spatial synchrony being disrupted by changes in environmental seasonal cycles include ocean primary productivity (Defriez2016), lake stratification (Kraemer2015), migration patterns (Usui2017), and flood hazards (Arnell2016), among others. The broad extent of these impacts highlight the importance of understanding how spatial patterns in environmental seasonal cycles vary through time.

The seasonal cycle in atmospheric variables is a direct response to the variation in solar insolation due to the Earth’s orbital path around the Sun. Specifically, the atmospheric response to the overhead sun crossing the equator twice a year suggests a more complicated seasonal variation that includes both an annual and a semi-annual harmonic. Harmonic analysis has been used by meteorologists and climatologists to characterize the connection between these harmonics and the observed seasonal cycle since the early 20th century (e.g., see Hsu1976; Hsu1976b for a review of this early work). Although the semi-annual harmonic typically contributes less variance to the seasonal cycle than the annual harmonic in the Northern mid-latitudes, its amplitude and phase vary considerably across space, and there are regions in which it can significantly affect the seasonal cycle (e.g., shifting the phase, strengthening the peak, flattening the minimum; see White1978). One of the first studies to discuss a specific dynamical mechanism behind the semi-annual cycle was VanLoon1967, in which he showed that the semi-annual cycle in the high Southern latitudes was the result of differential heating due to north–south land/sea contrasts between Antarctica and the Southern ocean. Unlike the north–south contrast exhibited in high latitudes, Wikle1996 showed evidence that the Northern hemisphere extratropical semi-annual cycle exhibits a strong east–west structure and is governed by the spatio-temporal asymmetries in the seasonal variation of the northern hemisphere stationary eddies (e.g., the wave structure in the atmospheric circulation). Because this variation in stationary eddies is due to east–west differential heating from land/sea contrasts, the fundamental mechanism for both the extra-tropical and high-latitude semi-annual cycles is due to land–sea contrasts and the impact this has on the atmospheric circulation.

It is well-known that atmospheric circulation patterns are varying due to differing responses of land and sea to climate forcing (e.g., Sutton2007). This suggests that the annual and semi-annual harmonics are also likely varying. Stine2009

investigated the change in the annual harmonic component of surface temperature between the years 1900-1953 and 1954-2007 by looking at the lag (difference between temperature and local solar insolation phases) and gain (ratio of temperature and insolation amplitudes). Based on simple t-test comparisons of the lag and gain for the different time periods, they showed that the annual temperature cycle has changed, but asymmetrically across space.

Dwyer2012 also showed that there has been heterogeneous variation in the annual amplitude and phase of the mean surface temperature cycle in response to greenhouse gases. These analyses did not explicitly consider the semi-annual component of the seasonal cycle, multivariate seasonal variation of atmospheric variables, nor did they consider a formal model-based uncertainty quantification framework that could accommodate spatial and temporal variation of the harmonics.

Dynamic spatio-temporal models (DSTMs) are well established in the literature for modeling complex spatial processes that evolve over time (see Cressie2011

for a collection of references and methods). Statistical DSTMs are able to capture spatial and temporal dependence in the process across different scales, while retaining the ability to capture uncertainty in parameter estimation and prediction. Surface temperatures over land generally can be decomposed into three components: a term to account for trend, a seasonal component, and a “weather” component. We would expect the first two of these components to vary somewhat slowly across time, with near-by locations experiencing similar temperatures and temperature variation through time, whereas the weather component corresponds to a dynamic process observed at finer spatial and temporal scales

(Wikle1998). The most natural way to accommodate slowly varying time variation in trend and seasonal parameters is via the dynamic linear model (DLM) paradigm (e.g., West2006). Such models are commonly extended to the dynamic evolution of parameters in spatio-temporal and multivariate settings by representing the parameters as spatial fields that vary in time according to a DLM (e.g., see the overviews in Gelfand2010; Cressie2011; Banerjee2015; Gelfand2017).

The main modeling contribution of this work is the development of a joint statistical framework for time-varying minimum and maximum temperature cycles, which are specified through the annual and semi-annual harmonics, while accounting for spatial and temporal dependence. This model is motivated by the fact that responses to changes in heating are asymmetric in space; thus, we expect that the annual and semi-annual harmonics in temperature are varying in time differently across space, leading to time-varying differences in seasonal cycles. By adopting a Bayesian framework for parameter estimation for the associated DSTM model, we are able to quantify the extent to which regions across North America are experiencing significant shifts in the minimum and maximum temperature cycles. Significant asymmetric shifts in minimum and maximum temperature seasonal cycles may seriously affect biological processes that are synchronously linked to such cycles.

The remainder of this manuscript is structured as follows. In Section 2 we describe the data used in the analysis and provide a brief exploratory analysis to motivate the work. Section 3 details the joint model specification using the annual and semi-annual harmonics and methods for model inference. Section 4 presents the findings of the analysis and Section 5 provides a discussion and directions for future work.

2 Data and Preliminary Analysis

For the analyses presented here, we consider air temperature (deg C) data at two meters above the surface obtained from the National Center for Environmental Prediction (NCEP) Reanalysis111https://www.esrl.noaa.gov/psd/. The data are available at three hour intervals for each day from January 1, 1979 to December 31, 2018, which we summarize as daily minimum and maximum temperature. The data are on a Northern Lambert Conformal Conic grid, with corners at approximately (1.000N, 145.500W), (0.898N, 68.320W), (46.354N, 2.570W), and (46.634N, 148.642E). All data exploration was conducted on a reduced spatial domain with corners at approximately (16.103N, 140.543W), (15.997N, 73.229W), (57.601N, 22.274W), and (57.856N, 168.499E).

Let denote temperature (say, minimum or maximum) on day where , and is the number of days in the year (365 or 366 for leap years). The Fourier series representation of the time series, expressed in amplitude-phase form, is given as

(2.1)

where and are the amplitude and phase, respectively, for the harmonic component. Reparameterizing Eqn. 2.1 in terms of its Fourier coefficients results in

(2.2)

where and are the Fourier coefficients for the harmonic, related to the amplitude by , and the phase by . We restrict our estimation to the first two harmonics () as discussed in the introduction, and refer to a “cycle” as the sum of the first and second harmonics hereafter.

Figure 1: Comparison of estimated seasonal cycle for minimum temperature using the first and second Fourier harmonics (solid line) compared to just the first Fourier harmonic (dashed line) in central Texas and Kings Canyon National Park, California. The top plot is the estimated cycle for 1979, the middle for 1999, and the bottom showing the difference (1999 year - 1979 year).

To investigate the possible relative importance of the semi-annual harmonic in North American minimum and maximum temperature cycles, Figure 1 shows the estimated seasonal cycles for minimum temperature when both the annual and semi-annual harmonic components are included (solid) compared to those in which only the annual component is considered (dashed). These estimated cycles are shown for two different years, 1979 and 1999, and two different locations, one in central Texas and the other in Kings Canyon National Park, California. The top panel shows the estimated cycle for the year 1979, the middle panel for the year 1999, and the bottom panel shows the difference between the two cycles, with the estimates from the 1979 cycle subtracted from the 1999 cycle. These plots illustrate the impact the semi-annual component can have on the temperature cycle, how the impact changes through time, and how these features vary across space. The estimated cycles for 1979 are very similar for both locations, suggesting the semi-annual harmonic had little influence on the minimum temperature cycle at these locations. Conversely, the temperature cycles for 1999 are much more dissonant for both locations, implying the semi-annual harmonic had a greater impact on the temperature cycle in 1999 than in 1979. The impact of the semi-annual harmonic of minimum temperature can be seen clearly in the bottom panel, where the difference between the estimated cycles between the two years using both annual and semi-annual components in central Texas shows a cyclical deviation from the difference in cycle estimates using only the annual component. At Kings Canyon, the cycle with the annual and semi-annual components oscillates about the difference using only the annual component. This suggests that the impact of the semi-annual component is spatially and temporally varying, and is important in capturing shifts in the temperature cycle through time.

Figure 2:

Computed t-statistics of the Fourier coefficients for the minimum temperature cycle. Locations with 95% point-wise confidence intervals not including 0 are shown, with the color corresponding to the t-statistic value.

Figure 3: Computed t-statistics of the Fourier coefficients for the maximum temperature cycle. Locations with 95% point-wise confidence intervals not including 0 are shown, with the color corresponding to the t-statistic value.

To identify differential change in space through time in the annual and semi-annual harmonics, we investigate daily temperature data for each location for two separated 15-year time periods; period one from 1979 - 1994 and period two from 2003 - 2018. From the annual and semi-annual harmonic estimates, we calculate the phase and amplitude for each year for each of the two time periods. As an exploratory comparison of the two periods at each location, we compute t-statistics of the difference (second period - first period) for both the annual and semi-annual phase and amplitude. Figures 2 and 3 show the t-statistics over the region for all four components (, , , ) for minimum and maximum temperature, respectively. T-statistics that are large (in magnitude) suggest possible shifts in temperature cycles. For example, in both the minimum and maximum temperature annual phase, a large area in the northeast of the continent has experienced a negative shift, implying the peak of the wave is occurring earlier in the recent years. Additionally, the shift in the annual amplitude for both the minimum and maximum cycles are alike, with areas over the Pacific, central United States, and northern Canada having similar spatial patterns and values. The shifts in the semi-annual amplitude and phase of both cycles have similar spatial patterns, with areas in the northwest United States and northern Canada most closely resembling each other. These results for the annual amplitude and phase closely resemble the results over North America reported in Stine2009, however our findings are purely exploratory as they do not account for any spatial, temporal, or process dependence, which could have important effects on reported regions of significant change.

These preliminary analyses identified possible shifts in the annual and semi-annual harmonic components of minimum and maximum temperature over North America, which suggest changes in the cycles themselves, and that these changes might vary considerably across space. Our aim is obtain full probabilistic inference with regard to these changes in minimum and maximum temperature cycles across the region. Specifically, we propose a multivariate statistical DSTM that captures the relationship between the minimum and maximum temperature cycles, as well as spatial and temporal dependence. This model will be used to quantify the uncertainty associated with potential seasonal cycle changes over space and time.

3 Bayesian Hierarchical Formulation

3.1 Data model

Let and denote the centered, by year, minimum and maximum temperature, respectively, for year , at location , where is the number of days in year . The linear model for temperature at location , year , and variable can be specified in terms of the Fourier coefficients (e.g., Eqn. 2.2) by

(3.1)

where , and , with the element of and equal to and , respectively, for , and . Here, is the variance for the variable at location , which is assumed constant over years for each cycle and location. Although the simplifying assumption of iid errors assumes no residual temporal autocorrelation as would be present in “weather” processes, preliminary analyses that accounted for this correlation through a daily random effect did not substantially impact parameter inference. As such, the model with daily random effects was not considered further due to the added computational complexity.

Spatial and temporal dependence is modeled using spatially-varying harmonic coefficients and with a random walk time structure. Specifically, the harmonic coefficients, , are spatial processes and the spatial field evolves according to a random walk. Letting denote the number of Fourier coefficients per cycle, the

-vector of spatially-varying coefficients,

, for year , location is modeled by

(3.2)

where with , and

where , is a unstructured covariance matrix, and denotes a Gaussian process over the spatial domain (Cressie1992). Each spatial process accounts for the residual spatial variation in the Fourier coefficients between year and . We assume an exponential covariance function, where , is the Euclidean distance between locations and , and consists of the spatial variance and decay parameters, respectively, for process . Note that the spatial covariance is process and parameter specific, but assumed constant across years.

Based on the similarities in the parameter estimates discussed in Section 2, the harmonic coefficients for both the minimum and maximum cycles are modeled jointly to borrow strength. Dependence between the minimum and maximum temperature cycles is captured through the covariance structure in the coefficients, . Lastly, we let , completing the model specification.

3.2 Predictive Process

Model inference presents computational challenges due to the number of spatial locations, years, harmonics, and processes. For example, a single draw from the conditional distribution of requires matrix operations on a matrix, which is computationally prohibitive for even modest sized data sets. Therefore, we propose using spatio-temporal predictive processes (Finley2012) to enhance computational efficiency.

Let be the locations where minimum and maximum temperature data are available. Next, define knot locations located inside the domain of interest where . For cycle , at location , year , for process , we define the predictive process (Finley2012) as

(3.3)

where , is a vector whose element is , and is the matrix with element given by . Let and , then the predictive process for the coefficients at year , , is predicated on all previous predictive processes,

where and is defined as in Eqn. 3.2. The resulting distribution of the coefficients is

and the data model (Eqn. 3.1) can be rewritten in the form of the predictive process.

where is defined the same as in Eqn. 3.1. Therefore, conditioned on the predictive process, the coefficient process is spatially independent, and draws from the full conditional of can be obtained univariately.

3.3 Parameter Models

To fully specify the Bayesian hierarchical model, we assign prior distributions to all remaining parameters. Conjugate, non-informative priors were chosen when available to ease computational burden. For , the variance for location that is shared across time is modeled . For the covariance matrix of the parameters we assign . Lastly, the spatial variance for the spatial process is modeled

. All hyperpriors were chosen such that the priors have finite first moments, specifically

, , and for all . Preliminary analyses with a uniform prior distribution for the spatial decay parameter, indicated that this parameter had little impact on the inference of the parameters of interest. Therefore, we set to a fixed value in our analysis.

The full hierarchical model can be written

(3.4)

where is a deterministic composition of , as shown in Eqn. 3.3.

3.4 Model Inference

Recall from Section 2 that the motivation for this modeling effort is purely inferential. Specifically, we are interested in inference with respect to the spatial processes of harmonic coefficients, for . We obtain samples from the joint posterior distribution using a Gibbs sampling algorithm, the details of which are given in Appendix A. Each of the parameters described above have a conjugate full-conditional distributions. To improve computational efficiency of our sampling algorithm, we also took advantage of parallel computation when possible. Specifically, conditioned on the predictive process, , the parameters are spatially independent and can be updated in parallel.

Posterior inference will focus on the comparison of the amplitude and phase of the annual and semi-annual harmonics across all years and spatial locations. Visualizing their spatial evolution over time provides insight into how the temperature cycle changes across years. Using samples from the posterior distribution of , we can obtain full posterior inference for the phase and amplitude of the annual and semi-annual harmonics using composition sampling. We can also compute important characteristics of these cycles, such as the day at which the cycle reached its peak or trough. These peak and trough days can be compared across time to quantify shifts in temperature cycles which may have important impacts on spatial synchrony. In addition, we can identify similarities and differences in shifts in minimum and maximum temperature cycles.

4 Results

We fitted the model to the NCEP Reanalysis temperature data introduced in Section 2. The spatial domain of interest spanned the continental United States and portions of Mexico and Canada, with corners at approximately (25.039N, 120.243W), (21.399N, 79.588W), (46.471N, 64.321W), and (52.7418N, 128.744W). We thinned the data spatially to reduce the overall dimension, keeping every other location in both the longitudinal and latitudinal directions. This resulted in daily minimum and maximum temperature values at 3621 spatial locations over 40 years, for a total of over data points. For the predictive process outlined in Section 3.2, we chose 144 knot locations evenly spaced across the domain of interest. All temperature time series were centered since the focus of inference is on the harmonics and change in harmonics as opposed to raw temperature.

Using MCMC and the Gibbs sampling algorithm (see Appendix A), we obtained 5000 samples from the joint posterior distribution. The first 1000 samples were discarded as burn-in and the remaining 4000 samples were retained for posterior inference. All computation and posterior inference was performed on a high performance computing infrastructure222Computation was performed on a Linux workstation using an Intel(R) Xeon(R) CPU E5-2680 v4 @ 2.40GHz processor utilizing 24 cores.

due to the dimensionality of the data and Bayesian inference output. Convergence of model parameters was assessed visually via trace plots, with no issues detected.

Figure 4: Posterior mean estimates of the minimum temperature annual and semi-annual harmonics for 2004. The left two panels show the annual and semi-annual amplitude (deg C), and the right tow panels show the annual and semi-annual phase (radians). Location A corresponds to the top plot in Figure 6; Location B corresponds to the middle plot in Figure 6; Location C corresponds to the bottom plot in Figure 6.

Posterior distributions of the annual and semi-annual phase and amplitude of minimum and maximum temperature cycles were obtained for each year using composition sampling. Posterior mean estimates of these four cycle quantities for minimum and maximum temperature for the year 2004 are shown in Figures 4 and 5, respectively. For both figures, estimates of the annual and semi-annual harmonics are shown on the top and bottom panels, respectively, while the amplitude and phase estimates are shown on the left and right, respectively. A prominent spatial feature of the temperature harmonics is the wave-like pattern that appears in the semi-annual amplitude and phase. This wave can been seen most prominently in the bottom-right panel of Figure 4 where a band of semi-annual phase values close to 0 (or ) spans from the north-west United States down through the center of the United States. The semi-annual phase estimates on either side of this band deviate from 0 (or ). The same spatial pattern appears in the semi-annual amplitude, where there is a band of smaller amplitudes following approximately the same path.

Figure 5: Posterior mean estimates of the maximum temperature annual and semi-annual harmonics for 2004. The left two panels show the annual and semi-annual amplitude (deg C), and the right tow panels show the annual and semi-annual phase (radians). Location A corresponds to the top plot in Figure 6; Location B corresponds to the middle plot in Figure 6; Location C corresponds to the bottom plot in Figure 6.

These same spatial patterns can are seen in Figure 5 for the semi-annual components of the maximum temperature cycle. In the bottom-right panel, the western United States have a contiguous area of lighter colored values close to 1.5 that are bordered by values close to 0 (or ). While less prominent than the minimum semi-annual amplitude, the maximum semi-annual amplitude has a band of smaller amplitudes that spans from the southern United States up through the center of the United States. These spatial patterns likely arise because of the east-west structure of the semi-annual harmonic due to the land/sea contrast as discussed in the introduction and by Wikle1996.

To investigate the temporal variation in minimum and maximum temperature cycles throughout the 40 year period, we produced an animation of the posterior mean estimates of amplitude and phase for the first and second harmonics333https://joshuanorth.shinyapps.io/harmonics_application/. The animation illustrates the slow evolution of the annual components and the temporal volatility of the semi-annual components. While the wave-like patterns seen in 2004 (Figures 4 and 5) are the most common, variations of these

Figure 6: The phase in radians (angular direction of the arrow/line), , and amplitude (height of the point) in degrees Celsius, , for the minimum (blue) and maximum (red) temperature cycles at three spatial locations across the United States. See Figure 4 corresponding to geographic locations, with the top plot corresponding to location A, middle plot to location B, and bottom plot to location C.

spatial patterns appear in both the maximum and minimum semi-annual amplitude and phase in other years. Similar wave-like spatial patterns have been detected for geopotential height (Wallace1993; Thompson1998; Thiebaux1986; Wikle1996) as discussed in the introduction.

To illustrate the component-wise difference in minimum and maximum temperature cycles, Figure 6 shows posterior estimates of the annual and semi-annual amplitudes (height of the point on the y-axis) and phase (angular direction of the arrow) simultaneously for minimum and maximum temperature at three different spatial locations. The locations, denoted “A”, “B”, and “C” in Figures 4 and 5, are in the west, north central, and northeast, respectively. The relationship between the minimum and maximum annual amplitudes differ across space, which could be attributed to climatological variations. Specifically, for locations “A” and “B”, the range between the minimum and maximum annual amplitudes is greater than for location “C”, and the interannual variability for the annual amplitudes is much greater for location “B” than locations “A” and “C”. The minimum and maximum semi-annual phase are the same for most years (i.e., phase locked), with relatively little year-to-year change in the semi-annual amplitude. Compared to the annual phase, the semi-annual phase is more volatile with each location experiencing differing degrees of variability. The semi-annual phase appears the most variable for location “C”.

The extent to which the temperature cycles have shifted (i.e., how the temperature cycle determined by the estimates of the first two harmonics has changed over time) over the 40 year period can be seen by comparing the day of the year at which the temperature cycles are at their peak and trough. We computed the average peak and trough day for the years 1979-1988 and 2009-2018 as well as the differences between these two time periods (computed as 2009-2018 minus 1979-1988). These posterior distributions can be used to identify spatial regions experiencing significant shifts in the minimum and maximum temperature cycles. We consider a shift to be significant if the 95% credible interval of the difference does not include zero.

Figure 7: Day at which the maximum temperature cycle reaches its peak. The left image is the average over the years 1979-1988, the center image is the average over the years 2009-2018, and the right image is the difference between the two images (2009-2018 minus 1979-1988), showing only locations where there is a significance difference.
Figure 8: Day at which the minimum temperature wave reaches its peak. The left image is the average over the years 1979-1988, the center image is the average over the years 2009-2018, and the right image is the difference between the two images (2009-2018 minus 1979-1988), showing only locations where there is a significance difference.

Figures 7 and 8 show the average day of the year in which the maximum and minimum temperature cycles, respectively, obtain their peak, as determined by the first two harmonics. The peak day for maximum temperature can be thought of as the hottest day of the year, and the peak day for minimum temperature is the day at which the warmest low temperature occurs. The differences in peak days between the two decades for both maximum and minimum temperature clearly identify regions experiencing seasonal shifts in temperature. The areas in red indicate locations for which the peak day is occurring later in the year for the 2009-2018 decade, whereas blue regions correspond to locations in which the peak day is occurring earlier in the year for the more recent decade. The spatial patterns in the shifts in maximum and minimum temperature appear similar across the domain. The northern regions (Montana, the Dakotas, Minnesota, and Canada) appear to be experiencing the greatest shift towards later seasonal peaks, with much of the western United States experiencing more moderate shifts. For both minimum and maximum temperature, the only two areas experiencing a shift to earlier seasonal peaks are located in the Midwest United States and along the Northwest coast.

Figure 9: Day at which the maximum temperature cycle reaches its lowest. The left image is the average over the years 1979-1988, the center image is the average over the years 2009-2018, and the right image is the difference between the two images (2009-2018 minus 1979-1988), showing only locations where there is a significance difference.

Figures 9 and 10 show the day for which each maximum and minimum temperature cycle reaches its trough as determined by the first two harmonics. For minimum temperature, the trough corresponds to the coldest day of the year, and for the maximum temperature, the trough captures the day of the coldest high temperature. Again, we see very similar patterns between minimum and maximum temperature cycles and shifts. In contrast to the spatial distribution of the shift for the peak day, the spatial distribution of the shift for the trough day has a strong north/south pattern. The northern half of the United States and Canada are experiencing a shift toward later seasonal troughs, whereas the southern half of the United States and Mexico are experiencing a shift toward earlier seasonal troughs.

Figure 10: Day at which the minimum temperature cycle reaches its lowest. The left image is the average over the years 1979-1988, the center image is the average over the years 2009-2018, and the right image is the difference between the two images (2009-2018 minus 1979-1988), showing only locations where there is a significance difference.

To highlight the contribution of the semi-annual component on the temperature shift, we obtained posterior distributions of the peak and trough days as well as the decadal shifts for both the minimum and maximum temperature cycles using only the annual components. We then computed the posterior difference in the decadal shifts between those obtained when both the annual and semi-annual component were included and those when only the annual component was included. The posterior mean of these differences are shown in Figure 11. We considered the contribution of the semi-annual component to be significant if the 95% point-wise credible interval of the posterior distribution of differences did not include zero. Based on these posterior credible intervals, locations in white correspond to locations where the semi-annual harmonic did not significantly contribute to the shift in the temperature cycle. All other locations indicate that the semi-annual component contributed significantly in capturing shifts in temperature cycles between the two decades. In the left panels of Figure 11, red (blue) indicates locations for which the peak day has shifted to later (earlier) in the year when the semi-annual harmonic is considered. Similarly, in the right panels of Figure 11, red (blue) indicate locations where the trough day occurs later (earlier) in the year when the semi-annual harmonic is considered. In each of these figures, the shading corresponds to the magnitude of these differences. The spatial distribution of significant semi-annual harmonic contributions are similar for both the minimum and maximum temperature cycles. The magnitude of the shift is higher for maximum temperature than for minimum, which could be attributed to the maximum temperature having more seasonal variation than the minimum. The semi-annual component significantly contributes to the later peak day (positive shift) in both the minimum and maximum temperature cycle in the North (North Dakota, South Dakota, and Minnesota), Southwest (New Mexico and Arizona), and the Gulf of Mexico. The semi-annual component significantly contributes to the earlier trough day (negative shift) in both minimum and maximum temperature cycles in the Gulf of Mexico and western United States.

Figure 11: Difference of the decadal shifts for the estimates considering the annual and semi-annual component with estimates considering only the annual component. Shading corresponds to the magnitude of the differences, with white areas corresponding to locations where the semi-annual component is not significant. Red indicates the peak/trough day has occurred later in the year when the semi-annual harmonic is considered, and blue indicates the peak/trough has occurred earlier in the year.

5 Discussion and Future Work

We proposed modeling minimum and maximum temperature cycles jointly through the components of the annual and semi-annual harmonics using a DSTM to detect temporal changes in the seasonal temperature cycle that may vary across space. Implementing our model in a Bayesian paradigm, we obtain estimates of the annual and semi-annual phase and amplitude through composition sampling. Spatial maps showing the difference in peak/trough days of the minimum and maximum temperature cycles for the years 2009-2018 relative to 1979-1988 identified regions experiencing seasonal shifts, as well as regions for which the semi-annual component contributed significantly to these shifts. These maps showed that the peak day for both minimum and maximum temperature cycles has shifted to later in the year in northern regions, and the trough day has shifted toward later in the year in the northern half and earlier in the year in the southern half of the United States.

The results of our analysis can be compared to those presented in previous research. For example, a similar east-west structure in the semi-annual cycle was reported by Wikle1996 and attributed to the land-sea contrast. Additionally, using only the annual harmonic, Stine2009 found an asymmetrical spatial pattern in the shift of the temperature cycle. However, since we considered both the annual and semiannual harmonic, our results differed from theirs in terms of the regions identified as experiencing asymmetric shifts in temperature cycles. Lastly, our model detected spatially-varying shifts in the peak/trough of the temperature cycle ranging between 15 days earlier to 15 days later in the year, whereas Dwyer2012 reported the annual phase in the temperature cycle is shifting to only later in the year.

While the results of our model have scientific merit of their own, they can also be used to detect changes in spatial synchrony between temperature cycles and other important environmental processes. For example, incorporating estimates of shifting temperature cycles in models for bird migration could identify regions for which the spatial synchrony between migration patterns of birds and temperature cycles have been disrupted. Similarly, we can investigate the effects temperature shifts on the occupancy or abundance of native bark beetles, which could lead to improved predictions of beetle spread or risk of invasion as well as aid in conservation efforts. While these are just two brief examples, a better understanding of the direction and magnitude in the shifts in temperature cycles over the last 40 years will motivate future scientific hypotheses with regard to the effects of these changes on important environmental processes.

1

References

Appendix A Gibbs sampling algorithm

All code can be found at https://github.com/jsnowynorth/Harmonics. This includes code for downloading and processing the data, the sampler (listed below), and processing the results.

Defining notation that will be used and restructuring equations to match those used in the sampler, let

Writing the ’s in block notation, which will be used in the update of the predictive process, let and , resulting in

where Block-Diag.

The Gibbs sampler algorithm for Eqn. 3.4 is initialized by setting all parameters to some starting values. Then, for each iteration of the Gibbs sampler, parameters are updated:

  1. For , update where

    where Block-Diag.

  2. For and , update , where

    where and are chosen hyperpriors, and denotes the predictive process for cycle for all knot locations.

  3. For , update where

  4. For and , update , where

  5. , where

  6. For and , update , where

Acknowledgements

JSN and EMS were funded in part by the National Science Foundation #EF-1638550. JSN was also partially funded by Missouri EPSCoR, supported by the National Science Foundation under Award #IIA‑1355406 and #IIA‑1430427. CKW was supported by NSF Award #DMS-1811745. The computation for this work was performed on the high performance computing infrastructure provided by Research Computing Support Services and in part by the National Science Foundation under grant number CNS-1429294 at the University of Missouri, Columbia Mo. NCEP Reanalysis data provided by the NOAA/OAR/ESRL PSD, Boulder, Colorado, USA, from their Web site at https://www.esrl.noaa.gov/psd/.