the Creative Commons Attribution 4.0 License.
the Creative Commons Attribution 4.0 License.

A methodology for tracking cold spells in space and time: development, evaluation and applications
Charles Chemel
Amanda Maycock
Paul Field
Cold spells, identified as periods of prolonged extreme low temperatures, are often analysed in an Eulerian framework or through the use of case studies. However, this restricts information about their spatio-temporal evolution and limits the ability to compare analogous events that share similar developments. This study identifies cold spells as a series of geographical objects that are connected across subsequent time-steps. These objects are characterised by persistent low-temperature anomalies and can be grouped into the same event by using a connected component method, previously applied to heatwaves. This work extends this method further by taking into account advection by the tropospheric mean (weighted vertical average) horizontal wind. We also extend the methodology to filter quasi-stationary events that may have different drivers to transient events. Once catalogued, the cold spells are easily accessed based on their properties, location, or time period for further study. This study applies the cold spell identification methodology to the European Centre for ECMWF Reanalysis v5 (ERA5) dataset to develop a climatology of cold spells in the Northern hemisphere, establishing their seasonal variations and associated patterns of atmospheric circulation. We compare the results to an existing Eulerian based methodology and explain some reasons for differences in the cold spells identified. This study also demonstrates the typical pathways by which cold spells reach regions of East Asia, Europe and North America.
- Article
(30731 KB) - Full-text XML
- BibTeX
- EndNote
Cold spells (CS) are prolonged periods of extreme low temperatures. They have been linked with exacerbated respiratory, cardiac and cerebrovascular related deaths (Keatinge et al., 1984; Donaldson et al., 1999; Chen et al., 2013). In a study analysing causes of premature deaths over the world, Gasparrini et al. (2015) found that around 5 million deaths from the 74 million analysed were linked with low temperature conditions in the 3 weeks preceding death. Physiological changes are not only related to the severity of the low temperature, but also the duration of the CS event (Bull and Morton, 1978), the relative change in temperature (Ryti et al., 2016), and the climatology – with climatologically warmer regions being at higher risk than climatologically colder regions (The Eurowinter Group, 1997). Extreme low temperatures can also cause severe ecological damage, by affecting coral reefs, mangrove forests and animal migration rates (Boucek et al., 2016), and economical damage (Sun et al., 2022).
Methods for identifying CSs in the literature vary, but many methods use a local temperature threshold (e.g., absolute or percentile based) that must be exceeded for a specified number of consecutive days (see for instance Vavrus et al., 2006; Kolstad et al., 2010; Peings and Magnusdottir, 2014; Liu et al., 2023). However, this classification does not explicitly retain information about spatio-temporal coherence of CS events, which is essential for examining their dynamics, severity, energy budget, growth rate, and pathways. To gain such insights, case studies must be conducted. Other methods of identifying CSs include extracting spatial information about temperature anomalies using self organising maps (Nygård et al., 2023), or preselecting specific regions to study anomalously low temperatures and their precursors (Huang et al., 2021; Rantanen et al., 2023). These studies are then constrained by the geographical region of choice, which can provide an incomplete representation of the extent of CSs. The need for a Lagrangian approach to examine CS dynamics has been expressed in the literature (e.g., Tuel and Martius, 2024), and whilst a spatio-temporal approach is available to study cold air mass (Iwasaki et al., 2014) and its thermodynamics, it does not investigate directly anomalous near-surface temperatures.
A method for event-based spatio-temporal tracking of temperature anomalies was developed for heatwaves, where connected geographical objects have to meet specific criteria (Lo et al., 2021). These criteria can be percentile based thresholds such as the 95th percentile of daily mean wet-bulb temperature with an absolute minimum threshold (Jackson et al., 2025). Similar connected object algorithms have been previously used to study precipitation (Sellars et al., 2015), ocean dynamics (Xue et al., 2023), storms (Dixon and Wiener, 1993), dust events (Yu et al., 2018), floods (Debusscher and Van Coillie, 2019) and more recently a cold spell event (Stone et al., 2025).
The present study adapts the connected object tracking algorithm and applies it to CSs. Additionally, it presents an enhanced version of the algorithm that incorporates advection, identifying regions of anomalously low temperatures that are related to the same large-scale dynamical event, but that are disjointed in space and time. The disconnection can occur because CSs exhibit different dynamics when over land, ocean, and sea regions, with varying surface and atmospheric conditions influencing their behaviour. Hence, this modified Lagrangian approach offers a more comprehensive identification of CSs compared to some other methods. Furthermore, a filter is developed to identify quasi-stationary events that might be related to slowly evolving processes (e.g., in response to persistent sea surface temperature anomalies). These new options allow the identification of CSs, based on a common set of thresholds, to be more refined. The tracking framework used to identify CS events is described in Sect. 2 and compared to previous methods for characterising regional CSs in Sect. 3. Section 4 presents applications of this framework to CSs, while Sect. 5 discusses the findings. Conclusions are given in Sect. 6.
2.1 Input data
The algorithm requires the following gridded input variables at a temporal resolution of 24 hours or less (see Sect. 2.7): 2 m temperature, zonal and meridional wind components on pressure levels for the advection option, and optionally the land-sea mask for additional filtering. For illustration purposes, the present study uses 6-hourly data, averaged to daily means, from the Medium-Range Weather Forecasts Reanalysis v5 (ERA5) dataset (Hersbach et al., 2023a, b), gridded at 0.25° × 0.25° horizontal resolution, over the period 1940–2022.
2.2 Pre-processing
The 2 m temperature is deseasonalised and detrended to remove the long-term increase in temperature associated with global warming. The deseasonalisation involves Fourier transforming the climatological annual cycle of temperature at each grid cell, and retaining the six lowest frequencies. The inverse Fourier transform of these frequencies forms a smooth seasonal cycle, which is removed from each year to leave the residual temperature anomalies from the seasonal cycle. The data is then linearly detrended at each grid cell. Since this study focuses on the Northern hemisphere, the November-March extended winter temperature anomalies are extracted, to capture the coldest days of the year. A threshold is applied on temperature anomalies to extract only those that fall below the 5th percentile for 3 or more consecutive days. These are then grouped using the contouring algorithm, as described in Sect. 2.3.
There is an extra step in pre-processing if using the advection-based method, where a pressure-weighted average of wind fields between 800 and 250 hPa is calculated. This approach disregards surface-level winds, which are influenced by local topography, and assumes that CSs move as deep tropospheric air masses. The advection method neglects the vertical component of the wind (as it is much smaller than the horizontal component) and does not consider the effects of other processes such as diabatic heating and interaction with the surface that may affect the near-surface temperature.
2.3 Contouring
The contouring stage evaluates the 2 m temperature anomaly field at each time-step, masking all grid points that do not meet the threshold criterion (see Sect. 2.2). It then groups all contiguous grid cells that meet the criterion, including diagonal neighbours into objects. The contour of an object is defined as the outer boundary of grid cells with at least one adjacent masked neighbour and where the grid cell itself is not fully enclosed within the contour. A diagram illustrating examples of contoured objects is provided in Fig. A1. For each object, the contour and time of occurrence is identified and recorded, along with its location, temperature anomaly, wind field, and any other relevant fields inside the contour. In this study, a track is defined as a group of such contoured extreme temperature anomalies (CETAs) that can evolve coherently over time, forming a CS event. There are two ways that the algorithm can form an event, either through the overlap-based method or by the advection-based method (see the following Sect. 2.4). When performing advection, the wind field is linearly interpolated to the midpoint of the grid cell edges that lie on the contour of each CETA. The wind field along the contour is then stored along with the other variables mentioned above for each CETA.
2.4 Tracking
The tracking algorithm using the overlap method follows the connectivity algorithm presented by Lo et al. (2021). This algorithm is based on whether a CETA at the next time-step overlaps with a CETA from the previous time-step. The algorithm then performs the following logical steps (see the schematics in Fig. 1): if a CETA overlaps another CETA (by at least one grid cell) after one time-step, then this CETA is a continuation of the original CETA, and hence grouped with all its information under the same event. However, if no CETA is found to be overlapping a CETA on the previous time-step, then the CETA is regarded as decayed. The algorithm repeats this process for all CETAs. Any new CETAs that do not have a parent, from a previous time-step, will form a new event.

Figure 1Algorithm for tracking cold spell (CS) events. The tracks consist of contoured extreme temperature anomalies (CETAs) that meet the minimum area threshold Amin. When considering the advection-based method, the CETAs are transported by the mean tropospheric wind, capped at a value set by Umax. The tracks that do not contain at least one CETA of area greater than Alft are discarded. QS stands for quasi-stationary filter.
The tracking algorithm using advection assumes that the CS is advected by the mean tropospheric wind. CETAs that intersect with the advected contour edge at the next time-step are then considered to be part of the same event as the original CETA. To perform this advection, the algorithm first solves a direct geodesic equation at each midpoint edge of grid cells that lie along the contour of the CETA. This uses the wind interpolated at the contour edge and the time-step duration to calculate the distance travelled, as well as the azimuthal angle from the meridional and zonal wind components to determine direction of propagation. A convex hull problem is solved by taking the starting point and all the points along the geodesic trajectory until the endpoint, from all advected contour edge points. The solution forms the smallest possible polygon that encompasses all points, which describes the advected contour. The same rules as for the overlap algorithm then apply for cataloguing the CS events, as shown in Fig. 1, except that the overlap between two contours at consecutive time-step is now determined by the intersection between the convex hull polygon and the CETAs at the next time-step.
In general, the overlap-based method can neglect discontinuities and fragmentation of CETAs from one time iteration to the next, while the advection-based method can merge distant CETAs. Figure 2 illustrates two examples in which the constituent CETAs in the CS events differ between the methods, using the well-documented winter of 2009/10 in Europe (e.g., Christiansen et al., 2018; Prior and Kendon, 2011) and East Asia, which caused heavy snow fall resulting in the closure of schools, highways, and airports (NASA Earth Observatory, 2010). Note that in these examples although the average temperature anomaly over Eurasia was low throughout the CS period, the CETAs do not cover the entire region due to the constraints in relation to the 5th percentile and consecutiveness.

Figure 2Comparison between the overlap-based and advection-based algorithms. The left column shows two different CS events identified through the advection-based algorithm and the right column shows the corresponding CSs identified through the overlap-based algorithm. Filled contours show the composite area coverage of all CETAs during the event. Hatched contours in the left column highlight CETAs that were not found through the overlap-based method. The single CS event in panel (c) is represented by two CS events in panel (d), shown using two different filled colors. The events (a–d) occur between 7–25 December 2009, 11–25 December 2009, 22 December 2009–12 January 2010, and 24 December 2009–11 January 2010 and 1–12 January 2010, respectively. Grey shading shows the average 2 m temperature anomalies during the CS days. The thin black lines show the orography, with contours every 500 m.
In the first example (see Fig. 2a and b), the overlap-based method, fails to identify several CETAs belonging to the CS event, when compared to the advection-based method. The second example (see Fig. 2c and d), shows that the CETAs located in Europe and East Asia are grouped together when using the advection-based method, but when using the overlap-based method they are considered as two separate CS events. The advection-based algorithm is more likely to categorise two CETAs that are in proximity to each other as the same CS event. Hence in this case, the advection-based method finds that the CS occurring in Europe in early January of 2010 is related to the CS in East Asia. This approach allows one to research interconnected CS events. Furthermore, the advection-based method identifies more nearby CETAs, including those over the ocean which tend to be more discontinuous. However, both approaches show similar results for these case studies. In the first example, the CS first emerges in the Ural region and propagates westward, severely affecting the European continent. In the second example, the CS emerged also in the Ural region, but affected both north-west Europe and central Asia. This temporal evolution was identified through a day-by-day examination of the CS track (not shown here).
To summarise, while the advection-based and overlap methods identify the same CETAs, they differ in how they group these into CS tracks. The advection-based method tracks the evolution of CETAs in space and time more coherently, capturing the effects of winds that cause advection of cold air masses. In contrast, the overlap-based method tends to form many separate tracks of CETAs. This can be further seen by the hatched contours in Fig. 2c, which will be classed as separate CS tracks in the overlap-based method despite being part of a larger-scale event. As a result, some tracks identified by the overlap method may lack physical meaning. The differences between the two methods will be further discussed in Sect. 3.
2.5 Merging and splitting of tracks
All CS events in the catalogue are made of branches, which can merge with or split from the event. The treatment of merging and splitting of CS events in the algorithm is illustrated in the schematics in Fig. 3. Merging occurs when a CETA from an existing CS event forms a connection to a different CS event. In the example shown in Fig. 3a, two events overlap the same CETA at a given time, with the later event merged to the earlier event by incorporating it as a new branch; see Fig. 3b. Each branch within a CS event retains information about its merge location. The CETA causing the merge and its subsequent behaviour is then described within the earlier (unmerged) branch.

Figure 3Illustration of how merges and splittings of CETAs are handled in the algorithm for a hypothetical CS track. (a) The tracks at time-step t+9. At this time-step, CS Event 1 (blue) and CS Event 2 (pink) are seen to be intersecting the same CETA. (b) After this time-step, Event 2 merges into Event 1, becoming Branch 3 (blue). The new Event 1 decays entirely at time-step t+11. In this case, t is an arbitary time-step in the extended winter period. Black filled circles represent the formation of a branch. Merges are shown using a bold cross, with an arrow showing the merging. Splittings are shown through a diamond shape, and a decayed CETA is represented by a cross.
As a result, CS events can contain multiple branches depending on the number of merges. Merges do not always represent multiple cold air sources; more commonly, the algorithm identifies multiple branches due to fragmentation of a single CS event. For example, at the start of a CS track, CETAs can be small in scale and fragmented and therefore may be, at first, identified as separate events. In subsequent loops of the algorithm (see Fig. 1), these fragmented events are merged into a larger, more coherent cold spell, and the initial CETAs are seen as branches. Whether the CSs branches are due to multiple cold air sources or fragmentations can be easily confirmed by plotting the track of the cold spell.
Splitting occurs when a CS event overlaps with two (or more) CETAs on the next time iteration, as shown for Event 2 from time-step t+6 to time-step t+7 (see Fig. 3a). In such cases, the splittings are a continuation of the original event. There are also cases of self-merges when splitting occurs and the associated CETAs that are tracked following the splitting merge with CETAs that are part of the same event, as seen in Event 1 at t+7 (see Fig. 3a). However, a branch is only created when two separate CS events, that are currently not identified as the same event, are joined by the algorithm. Therefore, this self-merge will still be part of the same branch as it is already part of the CS event it is joining.
2.6 Filtering
Further criteria can be applied to filter the CS events. These criteria can be applied to the CETA themselves, filtering them before they are tracked, or to entire event as a post process. Some example criteria and the stage at which they are applied are highlighted by black rectangles in Fig. 1. For example, a minimum area threshold (Amin) can be set on the CETAs to reduce the tracking of small fragments that can persist and merge with other CETAs. Such CETAs are typically remnants of larger-scale events and may be considered less significant from an impact perspective. It is also recommended to apply a speed threshold Umax when using the advection-based method. This threshold will prevent very high wind speeds at polar jet locations from merging CETAs together over unrealistically large distances. Furthermore, advection using very high wind speeds may not be accurate as air parcel properties can change dramatically over long distances due to other physical processes (e.g., radiation driven cooling). In this study, Umax was set to 20 m s−1, a value at the high-end of observed cold front speeds (Reeder and Smith, 1986; Smith and Reeder, 1988; Gohm et al., 2010). Other thresholds can also be used to filter CETAs (e.g., one can consider only the CETAs with average temperatures below freezing).
At the post-processing stage, a lifetime area threshold (Alft) can be applied to the CS events. This ensures that at some point during the event's lifespan, a CETA must reach a specified minimum size. This enables the study of large-scale CS events, that are likely to originate from large-scale dynamical drivers rather than from local conditions. Some large-scale events may be quasi-stationary, primarily those related to persistent temperature anomalies resulting from ocean–atmosphere interactions (e.g., La Niña events). These quasi-stationary CSs can merge with CETAs that are driven by synoptic-scale dynamical processes. A quasi-stationary filter was designed to locate cases for which in a single branch a CETA persists over any grid cell for longer than a specified period. If that branch is predominantly over the ocean (e.g., 80 % of the area of the region it encompasses is over the ocean), then the branch can either be removed or kept and tagged as a different type of CS event in the catalogue. The filter then re-evaluates all branches in the affected event, amends their connections or deletes them if they no longer meet any previous filtering criteria (e.g., Alft threshold). An example of where the quasi-stationary filter is effective is demonstrated using a CS from the winter 1998/99; see Fig. 4a. This CS lasted 139 d and is predominately driven by recurring CETAs over the Pacific Ocean. These CETAs emerged over the ocean due to a particularly strong La Niña event that persisted for the entire winter (see Shabbar and Yu, 2009; Dong et al., 2000). However, additional CETAs originating from the Arctic merged with this event toward the end of their lifetime due to their proximity to the Pacific Ocean. One could interpret this as distinct events driven by different mechanisms and choose to separate them, using the quasi-stationary filter, as shown in Fig. 4b to d. This results in the CETAs emerging from the Arctic being considered a separate CS event.

Figure 4(a) Composite of CETAs from a CS event identified from 1 November 1998 to 19 March 1999 (lasting 139 d) using the advection-based algorithm, with Amin = 70 000 km2 and Umax=20 m s−1 and no quasi-stationary filter applied. b-d) shows the composite CETAs belonging to the three independent CS events resulting from the event shown in panel (a), after a quasi-stationary filter of 30 d was applied, with the quasi-stationary CS event shown in panel (b) removed from the catalogue of events. The CSs in panels (b)–(d) occur between 1 November 1998–19 March 1999, 22 November–2 December 1998 and 23 January–23 February 1999, respectively. The thin black lines show the orography, with contours every 500 m.
2.7 Sensitivity tests
The present study uses a daily time-step, and while a shorter time-step could yield more accurate results, using a significantly longer time-step is not recommended, as it would result in inaccurate representation of events. Indeed, when using the advection option, the convex hull polygon would grow progressively larger for longer time-steps, increasing the likelihood of non-physical and distant connections between CETAs. As for the overlap-based method, the percentage overlap between non-stationary CETAs would decrease with longer time-steps, potentially disjoining otherwise connected CSs.

Figure 5(a) Number of CS tracks per year as a function of Amin. Inset in panel (a) shows a close examination of the behaviour below 40 tracks per year. (b) Number of CS tracks per year as a function of Alft; Alft = 0.8 × 106 km2 is highlighted. (c) Total area coverage across all CETAs at various Amin thresholds, averaged over all years. (d) Number of CS tracks that are altered as a function of quasi-stationary persistence threshold (QS), averaged over all years; QS = 30 d is highlighted. Green lines indicate the advection-based method and red lines indicate the overlap-based method.
This work focuses on large-scale extreme CS events, and thresholds are applied to extract only these events. Sensitivity tests are performed to evaluate the impact of these thresholds and ensure robustness in identifying the most significant events. To extract these CSs, Alft was set to be of the typical size of a mid-latitude weather system of 1000 km across, that is 0.8 × 106 km2. This threshold is appropriate to study CSs which likely originate from large-scale atmospheric structures, such as blocking events, rather than from local conditions (e.g., frost holes). Applying this threshold reduces the annual total cumulative area occupied by CETAs by 13 % (25 %) for the advection-based (overlap-based) method (see Fig. 5c). More CETAs are retained when using advection rather than overlap as the CETA overlap polygon is contained within the convex hull polygon. While both methods initialy include the same CETAs, Alft filters out smaller events, leading to a sharp decrease in the number of events, by 98 % (99 %) for the advection-based (overlap-based) method (see Fig. 5b). As a result, the remaining events include different sets of CETAs and hence occupy different total areas.
To assess the sensitivity of Amin, we begin by examining variations observed in the absence of other applied filters. Setting the minimum area threshold Amin to 10 000 km2 reduces the annual average total cumulative area occupied by CETAs by 3.3 %, when compared with no area threshold, for both the advection-based and overlap-based methods (see Fig. 5c). Increasing Amin further from 10 000 to 30 000 km2 and from 30 000 to 70 000 km2 leads to further reductions: 3.5 % and 4.5 %, respectively, for the overlap-based and advection-based methods. While the area occupied by CETAs decreases steadily as Amin increases, the number of events decreases sharply when CETAs smaller than 10 000 km2 are filtered (see Fig. 5a). The reduction is 76 % for the advection-based method and 79 % for the overlap-based method. Increasing Amin further from 10 000 to 30 000 km2 and from 30 000 to 70 000 km2 yields further reductions: 34 % and 47 %, and 28 % and 39 %, for the advection-based and overlap-based methods, respectively. Once Alft is applied, the average number of events increases by 30 % (7 %) from Amin of 0 to 70 000 km2, for the advection-based (overlap-based) method (see Fig. 5a), as less branches are now merging from these CETAs. The advection-based method is more sensitive to this filtering, as the resulting CS events are more reliant on CETAs to merge events together. For this study, Amin was set to 70 000 km2. While this value was chosen arbitrarily, it removes small-scale CETAs responsible for merges between large-scale events, with the number of events becoming more stable to slight changes in Amin, suggesting that the major CS events are identified.
The impact of removing quasi-stationary CETAs on the number of CS events is shown in Fig. 5d. The number of events altered (i.e., amended or removed) decreases rapidly as the threshold for persistence, the number of days a CETA must persist over a single grid point for its branch to be identified as quasi-stationary, increases. For thresholds of 10 and 30 d, 9.3 % and 1.0 % of events are altered for the advection-based option, and 8.7 % and 0.8 % for the overlap-based method, respectively.
A threshold of 30 d was chosen for this study, as a similar number of branches are filtered for both methods, suggesting that similar events are being identified. Additionally, using this threshold, 96 % of the total area of the filtered branches lies within the tropics (below or at 30° latitude) for both methods. The quasi-stationary behaviour identified is likely associated with slowly evolving tropical processes rather than faster-evolving extratropical synoptic-scale weather events. Applying this threshold reduces the annual average total cumulative area occupied by CETAs by 9 % for both methods, when also applying the other thresholds mentioned above. An intuitive understanding of the changes in area coverage between these filters is illustrated in Fig. 6. Alft and Amin reduce the total CETAs identified per year, whereas the quasi-stationary filter reduces the CETAs identified over the tropical oceans.

Figure 6Average number of days per year when a CS event is detected at each location using the advection-based method, when (a) there is no thresholds applied, (b) Alft = 0.8 × 106 km2, (c) Alft = 0.8 × 106 km2 and Amin = 70 000 km2, (d) Alft = 0.8 × 106 km2, Amin = 70 000 km2 and a quasi-stationary (QS) filter of 30 d. The regions defined by black polygons in panel (a) are used for analysis in Sects. 3, 4.2 and 4.3. The thin black lines show the orography, with contours every 500 m.
Unless stated otherwise, the results presented in the remainder of this study are obtained by using the following thresholds: Amin = 70 000 km2, Alft = 0.8 × 106 km2, and by filtering and discarding quasi-stationary CSs longer than 30 d. This catalogue contains 2371 and 3171 CS events for the advection-based and overlap-based methods, respectively. It is worth noting that the methodology presented in this section is not restricted to CS events; any feature with a detectable contour and an associated wind field could be analysed using this approach. It is recommended that if the methodology is applied to a different type of event, the thresholds are either motivated by physical characteristics and behaviour of the event or selected through a similar sensitivity analysis. This ensures that the thresholds are tailored to the specific application.
In order to assess the effectiveness of the tracking methodology, the CSs identified were compared with those following the Eulerian based method (EBM) used by Huang et al. (2021) for regions extending in longitude λ from λW to λE and in latitude φ from φS to φN. In Huang et al. (2021), CSs are identified according to the criterion:
represents the area of CSs relative to that of the region considered. Hence, a CS is identified if it occupies 30 % of the region. In this equation, H is the Heaviside function, being equal to one when , or zero otherwise. Here, is the daily 2 m temperature anomaly TS at time t, is the local standard deviation of TS, and α is a parameter that sets the severity of the CSs. Huang et al. (2021) selected α=0.8 for moderate CSs and α=1.2 for severe CSs. Assuming a normal distribution the 5th percentile corresponds to α=1.645. The same regions of study were chosen: the southern (35–55° N) and northern (55–70° N) regions of Europe (0–60° E), North America (120–60° W) and East Asia (90–150° E); these regions are illustrated in Fig. 6a. For a fair comparison with Huang et al. (2021), a constraint was placed in the tracking algorithm such that 30 % of the geographical region must be covered by CETAs on a given day and an additional criterion was placed on Eq. (1) such that the CS must occur for three or more consecutive days (see Sect. 2.2).
Table 1Comparison of CSs found using the methodology presented in this study and the Eulerian based method (EBM) of Huang et al. (2021). The percentages represent the number of CS days identified relative to all days sampled as a fraction of 100. The intersection (∩) between EBM (with persistence) CS days and the tracking methods CS days is calculated as a percentage based on the total number of CS days identified by the tracking methods. The tracking methods have an additional threshold that the area of all CSs occupying the region, on a specific day, must cover at least 30 % of the total area of the region. Persistence requires Eq. (1) to be valid for 3 consecutive days. Results involving the EBM are displayed in two rows, with those of the upper row based on Eq. (1) and those of the lower row based on this equation altered such that is equal to the 5th percentile.

The percentage of winter days where the specified criteria are met are similar for the overlap-based and advection-based methods (to within 0.1 %; see Table 1), with the percentage ranging from 0.5 % in South America to 1.9 % in Northern Europe. The EBM shows a significantly wider range of CS occurrence across regions, varying from 0.9 % in North America to 5.0 % in Northern Europe. The intersection between the tracking algorithm days and the EBM days is greater than 99 % except for the Northern America region where it is 81 % and 80 % when using the advection-based and overlap-based methods, respectively (see Table 1). This difference was found to result from the 2 m temperature anomalies themselves, as they did not exhibit a normal distribution, and so the value of α did not correspond to the 5th percentile. To correct for this, α was set to vary such that is equal to the 5th percentile. The introduction of this correction led to a comparable intersection of days between the methods, greater than 95 % (see Table 1). However, the EBM identifies more CS days than the tracking algorithm. There are two reasons for this: (i) the persistence criterion on regional scale used in the EBM is a weaker constraint than that used in the tracking algorithm at a grid cell scale, and (ii) the EBM does not consider the area thresholds Amin and Alft or the filtering of quasi-stationary CSs. Moreover, whilst the tracking algorithm does not capture all EBM CS days in a direct comparison, this does not imply that these CSs days are entirely unaccounted for. Most of these days are still represented in the catalogue, despite not meeting the 30 % area threshold as seen in Fig. A2.
Overall, the intersection of days found between the methods is greater or equal to 96 %, with the remaining days identified through the EBM also found in the tracking method at a smaller coverage ratio. This provides strong evidence that similar CSs are being identified using these methods. Moreover, the value of using the tracking algorithm is that a spatio-temporal reconstruction can now be be made for any CS event in the catalogue, as shown in Sect. 4.3.
4.1 Case studies of individual cold spell events
A CS case study is presented to highlight the information produced from the tracking algorithm. The case is that of a Eurasian cold air outbreak – causing the 5th coldest event in Japan on record (Tokyo Climate Center, 2018). The CS originated from eastern Siberia on 12 January 2018 (see Fig. 7), with the average 2 m temperature anomalies over the CETAs starting at −10 K and decreasing to as low as −19 K in the next 8 d (see Fig. 8). As the CS moved equatorward, temperature anomalies weakened. Once the cold air reached the Himalayas, the topography caused the cold air flow to diverge (see Fig. 7). Some of the cold air remained trapped in the Tarim Basin, due to its basin shaped topography, persisting from 28 January 2018 until it decayed (i.e., no longer met the tracking criteria) on 8 February 2018. This region typically experiences limited air movement due to the surrounding mountains, creating what is known as a “stagnation zone” (Wang et al., 2023). The higher concentrations of air pollution observed in this region serve as an indirect indicator of this stagnation, as the reduced airflow limits the dispersion of pollutants (Wang et al., 2019).

Figure 7Example of a CS event occurring between 12 January–11 February 2018, in 3 d intervals, shown in panels (a)–(k) respectively, using the advection-based method. Grey shading shows the 2 m temperature anomaly for each day. The blue filled contours represent the CETAs identified for each day. The dates denote the first day of each 3 d period. The thin black lines show the orography, with contours every 500 m.

Figure 82 m temperature anomaly averaged over each CETA (points) for the CS event shown in Fig. 7. The line represents the average of these points for each day.
The cold air westward of the Himalayas propagated to Central Asia, but had a very short lifetime and decayed on 31 January 2018. In contrast the cold air that was advected eastward, propagated equatorward for longer once it passed the Himalayan region, and decayed over the warm Pacific Ocean and the Maritime Continent. A case study of this specific CS was reported by Tokyo Climate Center (2018), showing that prior to the cold air intrusion into Japan on 23 January 2018, the polar jet's position was shifted equatorward in east Siberia, thereby drawing cold air equatorward and towards Japan, as shown in Fig. 7.
4.2 Climatology of cold spell events
The catalogue can be examined to characterise the spatio-temporal distribution of CSs over the entire Northern hemisphere. CS events are detected predominantly over land, especially over the Rocky mountains and Eurasia where they occur more than 5 d per year on average, with the least CSs located over the Atlantic and Pacific Oceans, as shown in Fig. 6d. CSs are also found to have preferential distributions depending on the specific month within the winter period (see Fig. 9). CSs are more prominent in northern East Asia (see Fig. 6a) in November, December and January (when 21 %, 23 % and 26 % of the CS days occur, respectively), with less CSs during February and March (18 % and 12 % respectively). In northern North America CSs occur predominately in December and January (when 24 % and 25 % of the CS days occur, respectively), with a large proportion of CSs affecting the Rocky mountains. Northern Europe experiences most CSs in January (when 29 % of the CS days occur), with fewer CSs during November and March (when 12 % and 13 % of the CS days occur, respectively); see Table A1 for more details. Midlatitudes experience least CSs during the month of March, when fewer CSs occur. The spatial distribution of the local standard deviation of the number of CS days per month for each extended winter month, corresponding to Fig. 9, shows that year-to-year variability is largest where there are more CS events (see Fig. A3).

Figure 9Average number of days when a CS event is identified at each location for each extended winter month with panels (a)–(e) representing November to March, respectively. Results are shown for the advection-based method, however, the overlap-based method yields similar results. The thin black lines show the orography, with contours every 500 m.
Figure 10 shows the direction and components of the tropospheric average wind during CS events (calculated using the wind field over the area occupied by each CETA) and of the average prevailing tropospheric wind across all winter days (climatological wind), respectively. This figure indicates that CS events in Europe are generally driven by tropospheric winds from the Norwegian Sea, the Barents Sea and Siberia (see Fig. 10a, c and e), with westward winds in Europe, in contrast to the eastward climatological winds (see Fig. 10c and d). In North America, the mean wind direction associated with CS events is more equatorward in the Rocky mountains and Bering Strait than that seen for the climatological winds. CSs tend to pass into North America through Canada, following the typical circulation pattern of the wind climatology (see Fig. 10a and b). However, during CS events, equatorward winds do tend to be stronger and more far reaching toward the Gulf of Mexico (see Fig. 10e and f). Eastward winds also appear stronger in the Pacific and Atlantic Oceans. In Greenland and Siberia, the regions of poleward wind become more localised (compared to the climatology) during CS events, occurring over a smaller geographical area (see Fig. 10a and b). Atmospheric conditions associated with positive Arctic Oscillation and Scandinavian blocking are known to cause a more east-equatorward flow from Greenland leading to more cold days in that region (Rantanen et al., 2023), whilst an Aleutian low can cause an equatorward flow in East Asia (Abdillah et al., 2017).

Figure 10(a, c, e) Average tropospheric wind from areas occupied by CS events. (b, d, f) Average climatological tropospheric wind during the extended winter period. (a, b) Wind direction on a cyclical color-map. (c, d) Average zonal wind speed. (e, f) Average meridional wind speed. The thin black lines show the orography, with contours every 500 m.
Overall, during a CS event equatorward winds intensify in southern Greenland, the Bering Sea, the Rocky Mountains, and Eurasia. In contrast, regions such as North America, southern East Asia and parts of Siberia experience equatorward tropospheric winds during CS events, that resemble those of the climatology. A monthly analysis provides a picture similar to that of the extended winter (not shown), namely there are preferential pathways for Arctic air to intrude into mid-latitudes.
4.3 Dynamics of regional cold spell events
Section 4.2 highlights that regions particularly affected by CSs are East Asia, Europe and North America. We now analyse the dynamics of CS events for each of these regions in turn. It is important to note that such a broad analysis was not possible before and has previously been limited to case studies. The same criteria as in Sect. 3 is used to illustrate this point, except that now a single event must occupy more than 30 % of the defined region at some stage in its lifetime (as opposed to previously where all the CETAs in the catalogue over a single day must cover more than 30 % of the defined region). This simple modification increases the sample size while allowing to focus on CSs that affect a specific region on a sufficiently large scale.
In southern East Asia (see Fig. 11), CSs first emerge in the Ural region at days −8 to −4. From day −4, the composite of CETAs shifts south-eastward. Day 0 marks the first day when CS event cover 30 % of southern East Asia. The CS composite appear partially blocked by the Himalayan topography while the composite of CETAs continues to shift mostly eastward. In the next 4 d, CS events are transported into the Pacific Ocean and southern China, where they decay from day 2 onward. A similar picture was provided by Yao et al. (2022) for the case study of the winter 2020/21, when a pool of cold air originating from Siberia progressed equatorward; this behaviour was associated with the development of a Ural blocking event. Likewise Shoji et al. (2014) and Abdillah et al. (2017) showed that a Siberian high and Aleutian low create a pressure gradient that guides the cold air mass from the Siberian region towards East Asia. The present study similarly shows that very large-scale CS events that occur in southern East Asia are associated with extreme low temperatures in the Ural region 8 d beforehand and on average decay 2 to 8 d after reaching southern East Asia.

Figure 11Subset composite of CSs events affecting East Asia and which cover at least 30 % of the region shown in the black box at some stage in their lifetime, identifed using the advection-based algorithm. (a)–(i) Composite lead-lag evolution of CETAs at 2 d intervals from day −8 to 8; the days shown represent the first day of each interval. Day 0 is taken as the first day the 30 % threshold is met. Vectors show the average wind direction at CETA locations on each day. In total, 24 CS events were identified. The thin black lines show the orography, with contours every 500 m.
In southern Europe (see Fig. 12), CSs appear over the Eastern European Plain (the Western Russian Plain), from day −6 to day −2, with some CSs also originating from the Barents Sea. The CS composite moves towards Western Europe while it is entrained in the zonal flow as it approaches the equator. The CSs typically disperse over the European continent from day 0 onwards, decaying as the CS shifts into the tropical regions and the Atlantic Ocean between day 0 to day 10. The CSs are very stationary in central Europe, but expand in area between day 0 to day 4, until the composite CS shifts toward the Atlantic Ocean after day 4. While some CSs are transported toward Eurasia throughout this period, they are small-scale, infrequent, and decay quickly compared to those in Europe. The picture presented above is consistent with known CS pathways during the negative phase of the North Atlantic Oscillation (NAO-) or Scandanavian blocking in Europe (Rantanen et al., 2023).

Figure 12Same as Fig. 11, but for southern Europe with panels (a)–(j) showing day −6 to 12. In total, 32 CS events were identified.
In southern North America (see Fig. 13), CSs first appear over north-west Canada and Hudson Bay about 4 to 6 d prior to reaching the region. Between day −6 and day −2, the CS composite increases in size in central North America. From day −2 to day 2, it becomes large-scale over North America, with the Rocky mountains channeling CETAs either towards the Atlantic Ocean or equatorwards. The CSs decay over the next several days in the Atlantic Ocean, off the Rocky mountains in the Pacific ocean, and in the Gulf of Mexico. It is known that cold air surges tend to be advected along and east of the Rocky mountains toward Mexico (Colle and Mass, 1995), as shown here. A similar behaviour was found by Walsh et al. (2001), when tracking air parcels during six major CS events in North America and Europe, with cold air in North America being advected equatorward, but in Europe being advected from Siberia westward. Further research is needed to better understand the role of synoptic-scale events, such as Ural blocking or the Aleutian Low, in the development of these CS events. The algorithm presented here offers a useful tool for isolating pathways of interest, and it has been shown that these pathways are consistent with those reported in the literature.

Figure 13Same as Fig. 11, but for southern North America region with panels (a)–(h) showing day −6 to 8. In total, 22 CS events were identified.
4.4 Cold spell events over the Arctic Ocean
In this section a criterion is applied to retain only CS events where at least one CETA extends over the Arctic Ocean, north of 60° N. These CS events are related to marine CSs that are typically identified using sea surface temperature anomalies (Schlegel et al., 2017). A heatmap of CETAs from all the CS events that meet this criterion (see Fig. 14a) shows regions of frequent CS occurrence, namely the Bering Sea, Hudson Bay, Labrador Sea, Norwegian Sea, and the Barents and Kara Seas. All of these regions are typically associated with seasonal variations in sea-ice cover (Danielson et al., 2011; Kowal et al., 2017; Wang et al., 1994; Germe et al., 2011; Kumar et al., 2021).

Figure 14(a) Average number of days per year when a CS is present at each location and at some stage in its lifetime is located within the Arctic Ocean (north of 60° N). (b–f) As in panel (a) but for subsets of CS events that intersect specific bounded regions: Hudson Bay, Labrador Sea, Norwegian Sea, Barents and Kara Sea, and Bering Sea, respectively (shown as white boxes in the panels). Note that no Alft threshold was applied, to understand the distribution of CS all events, not just the large-scale ones. The thin black lines show the orography, with contours every 500 m.
Marine CSs are associated with cold air masses that move from ice covered regions to the open ocean, thereby causing an increase in sensible heat flux, potentially causing baroclinic disturbances (Mansfield, 1974), leading to polar lows. These polar lows are then associated with extreme weather conditions (Landgren et al., 2019). Regions like the Labrador, Greenland and Barents and Kara Seas have been strongly associated with marine cold air outbreaks and have been extensively studied in the literature (see for instance Kolstad et al., 2009; Narizhnaya et al., 2020; Polkova et al., 2021; Dahlke et al., 2022). While not much research is available on CSs over the Hudson Bay and Bering Sea, (Fletcher et al., 2016) indicated that the strongest marine cold air outbreaks occur over the Bering Sea. Furthermore, temperatures in the Hudson Bay have been shown to be related to the temperatures experienced in North America and Canada (Rouse, 1991; Lochte et al., 2019), with CS dynamics there being different to that over the open ocean, due to being partially enclosed by land.
To further understand the distribution of these Arctic CSs, Fig. 14b to f show the composite of CETAs for all the CS events that pass through the Arctic Ocean and the regions: Hudson Bay (60–70° N, 70–100° W), Labrador Sea (60–80° N, 40–70° W), Norwegian Sea (60–80° N, 40–0° W), Barents and Kara Sea (60–80° N, 10–100° E) and Bering Sea (60–70° N, 160–210° W), respectively. CSs over the Norwegian Sea, and Barents and Kara Seas, have a larger impact over Europe. Those from the Barents and Kara Seas, in particular, are typically large-scale events, predominantly influencing Central and Eastern Europe, while extending their impact across Eurasia and parts of Western Europe. CSs from the Norwegian Sea are associated with temperature extremes occurring in the Labrador Sea and Barents and Kara Seas. These events form large-scale extreme events in north-western Europe. CSs from the Labrador and Bering Seas are typically linked to lower temperatures over the north pole. Additionally, CS events related to the Hudson Bay are connected to those from the Labrador Sea and contribute to extreme low temperatures in North America. In general, the application presented in this section can inform the study of future changes in the distribution of CSs over the Arctic Ocean, and their impact on continental regions.
Section 4 demonstrates how the CS identification methodology can be used to analyse both individual and composite CS events. To date, CSs have been mostly studied using Eulerian methods, limiting the analysis of their temporal evolution, unless case studies were considered. The framework developed here allows research on CSs to be more generalised than with case studies, with the benefit of enabling the analysis of thousands of CS events on a hemispheric scale. The CSs identified through the tracking algorithm are consistent with current knowledge of CS distribution. Yet, these CS events hold the additional information about their pathways (e.g., Fig. 7, showing the CS transported from East Siberia to central Eurasia and then toward the Maritime continent). An alternative approach to this methodology would be a Lagrangian analysis of back-trajectories of air parcels relating to a CS. This would trace the pathways of air masses responsible for the CS, thereby providing a deeper understanding of its origin, dynamics and spatio-temporal distribution. However, this approach is computationally demanding, and hence making it difficult to apply to large datasets. The method used in the present work, which tracks extreme temperature anomalies from their onset point, provides a simpler means to understand the pathways of CS events and captures key dynamics without the computational complexity of a full Lagrangian trajectory analysis.
The CS tracks will be sensitive to the definition of CS. In this work, a 5th percentile threshold was used to characterise the 2 m temperature anomalies. However, a fixed percentile threshold varies significantly with geographical location. The polar climate exhibits smaller temperature variations compared to mid-latitudes. On one hand, this is advantageous, as cold air moving equatorward warms up, but can still be detected as a 5th percentile event in that region. On the other hand, it complicates the identification of the origin of the air mass. Some CSs may originate closer to the extratropics simply because they did not meet the extreme cold thresholds near the pole, even though the air mass originated from that region.
The cooling that occurs as the air mass moves toward lower latitudes is related to the rates of advective and radiative cooling, which vary depending on the location of the temperature anomaly (Nygård et al., 2023). In the Southern hemisphere, the decay of cold air mass was shown to be dominated by diabatic processes made possible through the development of mesocyclones from baroclinic instabilities. The latent heat released from the ocean was found to be a key driver of the decay, along with the sensible heat flux (Papritz and Pfahl, 2016). In the Northern hemisphere, adiabatic warming can also be an important process, depending on the trajectory of the air parcel (Nygård et al., 2023). In this study, when CSs reach warm tropical oceans (Atlantic and Pacific Oceans and Maritime Continent), they decay rapidly (see Figs. 11, 12 and 13). However, to fully understand the dominant processes occurring across the CSs identified, in different regions of the Northern hemisphere, further exploration is needed regarding the respective role of adiabatic and diabatic processes.
This work addresses the need to characterise cold spell (CS) events, defined here as periods when the local near-surface temperature anomaly with respect to its climatological average is below the 5th percentile for a minimum of 3 d. For this purpose, a methodology is developed to identify the spatio-temporal behaviour of CSs. A novel element of the methodology is that regions of extreme temperature anomalies can be connected through space and time by advection using mean tropospheric winds. The output is a catalogue of CS events, in the form of tracks, which can be filtered for further analysis. Advection offers the benefit of identifying temperature anomalies under the same CS event that may not appear to be spatially connected at a given time. This is particularly useful over regions influenced by the same dynamical drivers but inhomogeneous surface forcing where extreme temperature anomalies can appear discontinuous. Constraints are applied to filter the CS events, in order to reduce fragmentations or false connections, resulting in fewer CS days being identified compared to Eulerian-based methods. However, the CS events identified from the tracking algorithm allow an enhanced insights into their spatio-temporal behaviour. The methodology is applied to the ECMWF Reanalysis v5 (ERA5) dataset, for the Northern hemisphere, over the period 1940–2022 to evaluate it and demonstrate its use for applications. The main conclusions are summarised below.
The CS events identified using this methodology are found to match those identified through the Eulerian-based approach used by Huang et al. (2021). While the Eulerian-based approach identifies a larger number of CS days based on certain thresholds, the tracking algorithm still captures the associated CS events, albeit at smaller scale. The differences in the CS days arise due to variations in the definitions of CS and the constraints applied to the CS events. Besides, the tracking algorithm provides a unique approach to retaining the spatio-temporal evolution of CSs, rather than simply identifying the days when CSs occur in a region.
Pathways for the transport of high latitude air to midlatitudes producing large-scale CSs, are concentrated in specific regions, as pointed out in the literature. CSs in East Asia are caused primarily by transport of colder air from the Ural region, with the topography playing a key role in channelling the air flow. In this case, low temperatures are detected in the Ural region 8 d before the cold air reaches southern Asia, and decays over a period of a couple of weeks. In North America, CSs are directed primarily to the east of the Rocky mountains, affecting multiple surrounding regions as they are transported toward the Gulf of Mexico and the Atlantic Ocean. However, some CSs are transported to the west of the Rocky mountains and decay in the Pacific Ocean. In Europe, CSs tend to be longer lasting when compared with the other regions and whilst the large-scale CSs originate from Siberia, other well-defined routes for the transport of colder Arctic air include along the Norwegian Sea and along the southern coast of Greenland.
The study also indicates that these preferred pathways result from anomalies in the large-scale circulation, facilitating the transport of colder Arctic air to midlatitudes. Specifically, the largest wind anomalies occur in the Bering Sea, the Rocky mountains, Greenland and the higher latitude coastal regions of Europe and Asia. The occurrence of CSs in midlatitudes is found to vary with time within the extended winter period. CSs in the ocean at high latitudes are found to be related to anomalous temperatures over the central Arctic. Further research is required to unravel the dynamics of CSs forming near these polar coastal regions, and over the nearby oceans, to understand their impact on terrestrial regions.
The applications presented in this work can be extended to examine spatial and temporal correlations between CS events or families of CS events. The catalogue of CS events can also be stratified according to metrics describing the large-scale circulation and the state of the climate. This will be reported in a future study.

Figure A1Illustration of two examples of contoured extreme temperature anomalies (CETAs). Grid cells that meet the threshold criteria and are part of CETA 1 (CETA 2) are shaded in pink (green) colour. Contours are highlighted using a thick black (red) outer line for CETA 1 (CETA 2). The locations where the wind field is interpolated along the contour of CETA 1 (CETA 2) are shown using black (red) dots.
Table A1Percentage of CS days for each extended winter month by region, based on the total CS days. Here, a CS day is defined as any day when at least one CETA intersects that region.


Figure A2Fraction of area covered by CETAs in northern (a, c, e) and southern (b, d, f) Europe (a, b), North America (c, d) and East Asia (e, f), for the set of days when the Eulerian-based method meets the 30 % threshold but the advection-based tracking method does not.

Figure A3Spatial distribution of the local standard deviation of the number of CS days per month for each extended winter month, corresponding to Fig. 9.
The tracking code is available to download through GitHub (see https://doi.org/10.5281/zenodo.14960646, Osmolska, 2025a), and the code used to make the analysis is available upon request.
The data is available to download through Zenodo (see https://doi.org/10.5281/zenodo.17378095, Osmolska, 2025b).
WO, CC and ACM designed the study. WO performed the analysis and produced the figures under the guidance of CC and ACM. WO, CC and ACM wrote the manuscript. PF provided suggestions to improve the algorithm and comments on the manuscript.
The contact author has declared that none of the authors has any competing interests.
Publisher's note: Copernicus Publications remains neutral with regard to jurisdictional claims made in the text, published maps, institutional affiliations, or any other geographical representation in this paper. While Copernicus Publications makes every effort to include appropriate place names, the final responsibility lies with the authors. Views expressed in the text are those of the authors and do not necessarily reflect the views of the publisher.
The authors would like to thank Alan Blyth for his help with this study, as well as the two anonymous reviewers for their feedback and constructive suggestions, which helped improve this manuscript. Copernicus Climate Change Service is acknowledged for providing the ERA5 dataset. This work used JASMIN, the UK's collaborative data analysis environment (https://www.jasmin.ac.uk, last access: 17 October 2025) for data processing (Lawrence et al., 2013). The colour palettes used to create some of the figures in this study were taken from Crameri et al. (2020).
This work was possible through the support from the Leeds-York-Hull Natural Environment Research Council (NERC) Doctoral Training Partnership (DTP) Panorama under grant NE/S007458/1. This work benefitted from and contributed to the Climate change in the Arctic–North Atlantic Region and Impacts on the UK (CANARI) programme funded by NERC, under grant NE/W004984/1. WO acknowledges additional funding provided by the UK Met Office through a CASE studentship. ACM was funded by the NERC StratClust project (NE/X011933/1) and the Leverhulme Trust.
This research has been supported by the Natural Environment Research Council (grant nos. NE/S007458/1, NE/X011933/1, and NE/W004984/1).
This paper was edited by Yang Zhang and reviewed by two anonymous referees.
Abdillah, M. R., Kanno, Y., and Iwasaki, T.: Tropical–Extratropical Interactions Associated with East Asian Cold Air Outbreaks. Part I: Interannual Variability, J. Climate, 30, 2989–3007, https://doi.org/10.1175/JCLI-D-16-0152.1, 2017. a, b
Boucek, R. E., Gaiser, E. E., Liu, H., and Rehage, J. S.: A review of subtropical community resistance and resilience to extreme cold spells, Ecosphere, 7, e01455, https://doi.org/10.1002/ecs2.1455, 2016. a
Bull, G. and Morton, J.: Environment, temperature and death rates, Age Ageing, 7, 210–224, 1978. a
Chen, R., Wang, C., Meng, X., Chen, H., Thach, T. Q., Wong, C.-M., and Kan, H.: Both low and high temperature may increase the risk of stroke mortality, Neurology, 81, 1064–1070, 2013. a
Christiansen, B., Alvarez-Castro, C., Christidis, N., Ciavarella, A., Colfescu, I., Cowan, T., Eden, J., Hauser, M., Hempelmann, N., Klehmet, K., Lott, F., Nangini, C., van Oldenborgh, G. J., Orth, R., Stott, P., Tett, S., Vautard, R., Wilcox, L., and Yiou, P.: Was the Cold European Winter of 2009/10 Modified by Anthropogenic Climate Change? An Attribution Study, J. Climate, 31, 3387–3410, https://doi.org/10.1175/JCLI-D-17-0589.1, 2018. a
Colle, B. A. and Mass, C. F.: The Structure and Evolution of Cold Surges East of the Rocky Mountains, Mon. Weather Rev., 123, 2577–2610, https://doi.org/10.1175/1520-0493(1995)123<2577:TSAEOC>2.0.CO;2, 1995. a
Crameri, F., Shephard, G. E., and Heron, P. J.: The misuse of colour in science communication, Nat. Commun., 11, 5444, https://doi.org/10.1038/s41467-020-19160-7, 2020. a
Dahlke, S., Solbès, A., and Maturilli, M.: Cold air outbreaks in Fram Strait: Climatology, trends, and observations during an extreme season in 2020, J. Geophys. Res.-Atmos, 127, e2021JD035741, https://doi.org/10.1029/2021JD035741, 2022. a
Danielson, S., Curchitser, E., Hedstrom, K., Weingartner, T., and Stabeno, P.: On ocean and sea ice modes of variability in the Bering Sea, J. Geophys. Res.-Oceans, 116, https://doi.org/10.1029/2011JC007389, 2011. a
Debusscher, B. and Van Coillie, F.: Object-based flood analysis using a graph-based representation, Remote Sens.-Basel, 11, 1883, https://doi.org/10.3390/rs11161883, 2019. a
Dixon, M. and Wiener, G.: TITAN: Thunderstorm identification, tracking, analysis, and nowcasting-A radar-based methodology, J. Atmos. Ocean Tech., 10, 785–785, https://doi.org/10.1175/1520-0426(1993)010<0785:TTITAA>2.0.CO;2, 1993. a
Donaldson, G., Seemungal, T., Jeffries, D., and Wedzicha, J.: Effect of temperature on lung function and symptoms in chronic obstructive pulmonary disease, Eur. Respir. J., 13, 844–849, https://doi.org/10.1034/j.1399-3003.1999.13d25.x, 1999. a
Dong, B. W., Sutton, R. T., Jewson, S. P., O'Neill, A., and Slingo, J. M.: Predictable winter climate in the North Atlantic sector during the 1997–1999 ENSO cycle, Geophysical Research Letters, 27, 985–988, https://doi.org/10.1029/1999GL010994, 2000. a
Fletcher, J., Mason, S., and Jakob, C.: The Climatology, Meteorology, and Boundary Layer Structure of Marine Cold Air Outbreaks in Both Hemispheres, J. Climate, 29, 1999–2014, https://doi.org/10.1175/JCLI-D-15-0268.1, 2016. a
Gasparrini, A., Guo, Y., Hashizume, M., Lavigne, E., Zanobetti, A., Schwartz, J., Tobias, A., Tong, S., Rocklöv, J., Forsberg, B., Leone, M., De Sario, M., Bell, M. L., Guo, Y.-L. L., Wu, C.-f., Kan, H., Yi, S.-M., de Sousa Zanotti Stagliorio Coelho, M., Saldiva, P. H. N., Honda, Y., Kim, H., and Armstrong, B.: Mortality risk attributable to high and low ambient temperature: a multicountry observational study, The Lancet, 386, 369–375, https://doi.org/10.1016/S0140-6736(14)62114-0, 2015. a
Germe, A., Houssais, M.-N., Herbaut, C., and Cassou, C.: Greenland Sea sea ice variability over 1979–2007 and its link to the surface atmosphere, J. Geophys. Res.-Oceans, 116, https://doi.org/10.1029/2011JC006960, 2011. a
Gohm, A., Mayr, G. J., Darby, L. S., and Banta, R. M.: Evolution and structure of a cold front in an Alpine valley as revealed by a Doppler lidar, Q. J. Roy. Meteor. Soc., 136, 962–977, https://doi.org/10.1002/qj.609, 2010. a
Hersbach, H., Bell, B., Berrisford, P., Biavati, G., Horányi, A., Muñoz Sabater, J., Nicolas, J., Peubey, C., Radu, R., Rozum, I., Schepers, D., Simmons, A., Soci, C., Dee, D., and Thépaut, J.-N.: ERA5 hourly data on pressure levels from 1940 to present, Copernicus Climate Change Service (C3S) Climate Data Store (CDS), https://doi.org/10.24381/cds.bd0915c6, 2023a. a
Hersbach, H., Bell, B., Berrisford, P., Biavati, G., Horányi, A., Muñoz Sabater, J., Nicolas, J., Peubey, C., Radu, R., Rozum, I., Schepers, D., Simmons, A., Soci, C., Dee, D., and Thépaut, J.-N.: ERA5 hourly data on single levels from 1940 to present, Copernicus Climate Change Service (C3S) Climate Data Store (CDS), https://doi.org/10.24381/cds.adbb2d47, 2023b. a
Huang, J., Hitchcock, P., Maycock, A. C., McKenna, C. M., and Tian, W.: Northern hemisphere cold air outbreaks are more likely to be severe during weak polar vortex conditions, Communications Earth & Environment, 2, 147, https://doi.org/10.1038/s43247-021-00215-6, 2021. a, b, c, d, e, f, g
Iwasaki, T., Shoji, T., Kanno, Y., Sawada, M., Ujiie, M., and Takaya, K.: Isentropic Analysis of Polar Cold Airmass Streams in the Northern Hemispheric Winter, J. Atmos. Sci., 71, 2230–2243, https://doi.org/10.1175/JAS-D-13-058.1, 2014. a
Jackson, L. S., Birch, C. E., Chagnaud, G., Marsham, J. H., and Taylor, C. M.: Daily rainfall variability controls humid heatwaves in the global tropics and subtropics, Nature Communications, 16, 3461, https://doi.org/10.1038/s41467-025-58694-6, 2025. a
Keatinge, W., Coleshaw, S., Cotter, F., Mattock, M., Murphy, M., and Chelliah, R.: Increases in platelet and red cell counts, blood viscosity, and arterial pressure during mild surface cooling: factors in mortality from coronary and cerebral thrombosis in winter., Brit. Med. J. (Clinical Research Ed.), 289, 1405–1408, https://doi.org/10.1136/bmj.289.6456.1405, 1984. a
Kolstad, E. W., Bracegirdle, T. J., and Seierstad, I. A.: Marine cold-air outbreaks in the North Atlantic: Temporal distribution and associations with large-scale atmospheric circulation, Clim. Dynam., 33, 187–197, https://doi.org/10.1007/s00382-008-0431-5, 2009. a
Kolstad, E. W., Breiteig, T., and Scaife, A. A.: The association between stratospheric weak polar vortex events and cold air outbreaks in the Northern Hemisphere, Q. J. Roy. Meteor. Soc., 136, 886–893, https://doi.org/10.1002/qj.620, 2010. a
Kowal, S., Gough, W. A., and Butler, K.: Temporal evolution of Hudson Bay sea ice (1971–2011), Theor. Appl. Climatol., 127, 753–760, https://doi.org/10.1007/s00704-015-1666-9, 2017. a
Kumar, A., Yadav, J., and Mohan, R.: Spatio-temporal change and variability of Barents-Kara sea ice, in the Arctic: Ocean and atmospheric implications, Sci. Total Environ., 753, 142046, https://doi.org/10.1016/j.scitotenv.2020.142046, 2021. a
Landgren, O. A., Seierstad, I. A., and Iversen, T.: Projected future changes in marine cold-air outbreaks associated with polar lows in the northern North-Atlantic Ocean, Clim. Dynam., 53, 2573–2585, https://doi.org/10.1007/s00382-019-04642-2, 2019. a
Lawrence, B. N., Bennett, V. L., Churchill, J., Juckes, M., Kershaw, P., Pascoe, S., Pepler, S., Pritchard, M., and Stephens, A.: Storing and manipulating environmental big data with JASMIN, in: 2013 IEEE international conference on big data, IEEE, 68–75, https://doi.org/10.1109/BigData.2013.6691556, 2013. a
Liu, M., Hu, D., and Guan, Z.: Modulation of a long-lasting extreme cold event in Siberia by a minor sudden stratospheric warming and the dynamical mechanism involved, Clim. Dynam., 60, 797–811, https://doi.org/10.1007/s00382-022-06353-7, 2023. a
Lo, S.-H., Chen, C.-T., Russo, S., Huang, W.-R., and Shih, M.-F.: Tracking heatwave extremes from an event perspective, Weather and Climate Extremes, 34, 100371, https://doi.org/10.1016/j.wace.2021.100371, 2021. a, b
Lochte, A. A., Repschläger, J., Kienast, M., Garbe-Schönberg, D., Andersen, N., Hamann, C., and Schneider, R.: Labrador Sea freshening at 8.5 ka BP caused by Hudson Bay Ice Saddle collapse, Nat. Commun., 10, 586, https://doi.org/10.1038/s41467-019-08408-6, 2019. a
Mansfield, D.: Polar lows: The development of baroclinic disturbances in cold air outbreaks, Q. J. Roy. Meteor. Soc., 100, 541–554, https://doi.org/10.1002/qj.49710042604, 1974. a
Narizhnaya, A., Chernokulsky, A., Akperov, M., Chechin, D., Esau, I., and Timazhev, A.: Marine cold air outbreaks in the Russian Arctic: climatology, interannual variability, dependence on sea-ice concentration, in: IOP C. Ser. Earth Env., IOP Publishing, vol. 606, 012039, https://doi.org/10.1088/1755-1315/606/1/012039, 2020. a
NASA Earth Observatory: Heavy Snow in Eastern China, https://earthobservatory.nasa.gov/images/42172/heavy-snow-in-eastern-china (last access: 25 November 2024), 2010. a
Nygård, T., Papritz, L., Naakka, T., and Vihma, T.: Cold wintertime air masses over Europe: where do they come from and how do they form?, Weather Clim. Dynam., 4, 943–961, https://doi.org/10.5194/wcd-4-943-2023, 2023. a, b, c
Osmolska, W.: coldsnap (v1.0.0), Zenodo [code], https://doi.org/10.5281/zenodo.14960646, 2025a. a
Osmolska, W.: Cold spell tracks from ERA5 dataset 1940–2022 (v1.0.0), Zenodo [data set], https://doi.org/10.5281/zenodo.17378095, 2025b. a
Papritz, L. and Pfahl, S.: Importance of Latent Heating in Mesocyclones for the Decay of Cold Air Outbreaks: A Numerical Process Study from the Pacific Sector of the Southern Ocean, Mon. Weather Rev., 144, 315–336, https://doi.org/10.1175/MWR-D-15-0268.1, 2016. a
Peings, Y. and Magnusdottir, G.: Response of the Wintertime Northern Hemisphere Atmospheric Circulation to Current and Projected Arctic Sea Ice Decline: A Numerical Study with CAM5, J. Climate, 27, 244–264, https://doi.org/10.1175/JCLI-D-13-00272.1, 2014. a
Polkova, I., Afargan-Gerstman, H., Domeisen, D. I., King, M. P., Ruggieri, P., Athanasiadis, P., Dobrynin, M., Aarnes, Ø., Kretschmer, M., and Baehr, J.: Predictors and prediction skill for marine cold-air outbreaks over the Barents Sea, Q. J. Roy. Meteor. Soc., 147, 2638–2656, https://doi.org/10.1002/qj.4038, 2021. a
Prior, J. and Kendon, M.: The UK winter of 2009/2010 compared with severe winters of the last 100 years, Weather, 66, 4–10, https://doi.org/10.1002/wea.735, 2011. a
Rantanen, M., Lee, S. H., and Aalto, J.: Asymmetric warming rates between warm and cold weather regimes in Europe, Atmos. Sci. Lett., 24, e1178, https://doi.org/10.1002/asl.1178, 2023. a, b, c
Reeder, M. J. and Smith, R. K.: A comparison between frontogenesis in the two-dimensional Eady model of baroclinic instability and summertime cold fronts in the Australian region, Q. J. Roy. Meteor. Soc., 112, 293–313, https://doi.org/10.1002/qj.49711247202, 1986. a
Rouse, W. R.: Impacts of Hudson Bay on the terrestrial climate of the Hudson Bay Lowlands, Arctic Alpine Res., 23, 24–30, 1991. a
Ryti, N. R., Guo, Y., and Jaakkola, J. J.: Global association of cold spells and adverse health effects: a systematic review and meta-analysis, Environ. Health Persp., 124, 12–22, https://doi.org/10.1289/ehp.1408104, 2016. a
Schlegel, R. W., Oliver, E. C., Wernberg, T., and Smit, A. J.: Nearshore and offshore co-occurrence of marine heatwaves and cold-spells, Prog. Oceanogr., 151, 189–205, https://doi.org/10.1016/j.pocean.2017.01.004, 2017. a
Sellars, S. L., Gao, X., and Sorooshian, S.: An Object-Oriented Approach to Investigate Impacts of Climate Oscillations on Precipitation: A Western United States Case Study, J. Hydrometeorol., 16, 830–842, https://doi.org/10.1175/JHM-D-14-0101.1, 2015. a
Shabbar, A. and Yu, B.: The 1998–2000 La Niña in the context of historically strong La Niña events, Journal of Geophysical Research: Atmospheres, 114, https://doi.org/10.1029/2008JD011185, 2009. a
Shoji, T., Kanno, Y., Iwasaki, T., and Takaya, K.: An Isentropic Analysis of the Temporal Evolution of East Asian Cold Air Outbreaks, J. Climate, 27, 9337–9348, https://doi.org/10.1175/JCLI-D-14-00307.1, 2014. a
Smith, R. K. and Reeder, M. J.: On the movement and low-level structure of cold fronts, Mon. Weather Rev., 116, 1927–1944, https://doi.org/10.1175/1520-0493(1988)116<1927:OTMALL>2.0.CO;2, 1988. a
Stone, J., Gervais, M., Bowley, K., and Zarczyki, C.: Identifying, Tracking, and Evaluating Mechanisms of North American Cold Air Outbreaks (CAOs) Using a Feature Tracking Approach, Mon. Weather Rev., https://doi.org/10.1175/MWR-D-23-0265.1, 2025. a
Sun, Q., Sun, Z., Chen, C., Yan, M., Zhong, Y., Huang, Z., He, L., and Li, T.: Health risks and economic losses from cold spells in China, Sci. Total Environ., 821, 153478, https://doi.org/10.1016/j.scitotenv.2022.153478, 2022. a
The Eurowinter Group: Cold exposure and winter mortality from ischaemic heart disease, cerebrovascular disease, respiratory disease, and all causes in warm and cold regions of Europe, The Lancet, 349, 1341–1346, https://doi.org/10.1016/S0140-6736(96)12338-2, 1997. a
Tokyo Climate Center: 51st issue of the TCC News – Japan Meteorological Agency, https://wmo.int/media/news-from-members/51st-issue-of-tcc-news-japan-meteorological-agency (last access: 30 August 2024), 2018. a, b
Tuel, A. and Martius, O.: Persistent warm and cold spells in the Northern Hemisphere extratropics: regionalisation, synoptic-scale dynamics and temperature budget, Weather Clim. Dynam., 5, 263–292, https://doi.org/10.5194/wcd-5-263-2024, 2024. a
Vavrus, S., Walsh, J., Chapman, W., and Portis, D.: The behavior of extreme cold air outbreaks under greenhouse warming, Int. J. Climatol.: A Journal of the Royal Meteorological Society, 26, 1133–1147, https://doi.org/10.1002/joc.1301, 2006. a
Walsh, J. E., Phillips, A. S., Portis, D. H., and Chapman, W. L.: Extreme Cold Outbreaks in the United States and Europe, 1948–99, J. Climate, 14, 2642–2658, https://doi.org/10.1175/1520-0442(2001)014<2642:ECOITU>2.0.CO;2, 2001. a
Wang, J., Mysak, L. A., and Ingram, R. G.: Interannual variability of sea-ice cover in Hudson Bay, Baffin Bay and the Labrador Sea, Atmos. Ocean, 32, 421–447, https://doi.org/10.1080/07055900.1994.9649505, 1994. a
Wang, J., Ju, T., Lei, S., Li, B., and Niu, X.: Study on Characteristics, Influencing Factors and Health Benefits of Atmospheric Multi-Pollutants in Southern Xinjiang, Atmosphere-Basel, 14, https://doi.org/10.3390/atmos14111681, 2023. a
Wang, Y., Zhang, J., Bai, Z., Yang, W., Zhang, H., Mao, J., Sun, Y., Ma, Z., Xiao, J., Gao, S., and Chen, L.: Background concentrations of PMs in Xinjiang, West China: An estimation based on meteorological filter method and Eckhardt algorithm, Atmos. Res., 215, 141–148, https://doi.org/10.1016/j.atmosres.2018.09.008, 2019. a
Xue, C., Niu, C., Xu, Y., and Su, F.: A Process-Oriented Exploration of the Evolutionary Structures of Ocean Dynamics with Time Series of a Remote Sensing Dataset, Remote Sens-Basel, 15, 348, https://doi.org/10.3390/rs15020348, 2023. a
Yao, Y., Zhang, W., Luo, D., Zhong, L., and Pei, L.: Seasonal Cumulative Effect of Ural Blocking Episodes on the Frequent Cold events in China during the Early Winter of 2020/21, Adv. Atmos. Sci., 39, 609–624, https://doi.org/10.1007/s00376-021-1100-4, 2022. a
Yu, M., Yang, C., and Jin, B.: A framework for natural phenomena movement tracking – Using 4D dust simulation as an example, Computers & Geosciences, 121, 53–66, https://doi.org/10.1016/j.cageo.2018.10.003, 2018. a
- Abstract
- Introduction
- Methodology
- Comparison with an established method to identify cold spells
- Examples of utilisation of the catalogue for studying cold spells
- Discussion
- Conclusions
- Appendix A
- Code availability
- Data availability
- Author contributions
- Competing interests
- Disclaimer
- Acknowledgements
- Financial support
- Review statement
- References
- Abstract
- Introduction
- Methodology
- Comparison with an established method to identify cold spells
- Examples of utilisation of the catalogue for studying cold spells
- Discussion
- Conclusions
- Appendix A
- Code availability
- Data availability
- Author contributions
- Competing interests
- Disclaimer
- Acknowledgements
- Financial support
- Review statement
- References