A pitchfork-like relationship between reduced Barents-Kara sea ice and Ural atmospheric circulation

Atmospheric circulation monthly anomalies over the Ural region are key indicators of Eurasian climate anomalies. Here, whether there exists a one-to-two correspondence relationship that generally agrees with the supercritical pitchfork bifurcation model, referred to as a pitchfork-like relationship, between reduced sea ice concentration (SIC) in the Barents-Kara Seas in specific months and the lagging Ural circulation anomalies is explored. Based on the monthly observational SIC data and two reanalyses during 1979/1980 − 2020/2021, two typical examples are found by estimating the joint probability density function. Results show that when the gradually reduced SIC in September (January) passes a critical threshold, the preferred Ural circulation patterns in October (February) exhibit a regime transition from the flat zonal westerlies to wavy westerlies with a Ural trough and wavy westerlies with a Ural ridge. Because both the barotropic and baroclinic conversion of energy from the climatological-mean flow to Ural circulation anomalies exhibit a regime transition from one regime to two regimes. It might be associated with the increased both positive and negative shear vorticity of background westerly wind over the Ural region before the regime transition, contributed by the thermodynamic effect of the SIC reduction. After the regime transition, positive and negative anomaly events of Ural atmospheric circulation occur with equal probability under the same SIC. Our results suggest an increased incidence of both positive and negative anomalies of Ural atmospheric circulation and also the Siberian High, under the recent SIC reduction, which implies a low predictability of Eurasian climate anomalies in October and February.


Introduction
Atmospheric circulation monthly anomalies over the region around the Ural Mountains are key indicators of Eurasian climate anomalies (e.g., Wang and Lu 2017;Cheung et al. 2013). On the one hand, because the Ural region is one of the preferred geographical regions over the Northern Hemisphere for the occurrence of persistent 500-hPa geopotential height anomalies (Dole and Gordon 1983;Liu et al. 2018;Miller et al. 2020), monthly anomalies of the Ural atmospheric circulation are usually associated with the recurrent, persistent large-scale ridges (or blocking highs) and troughs (or cut-off lows), which have a lifetime of 5 ~ 25 days (e.g., Miller et al. 2020) and often cause severe regional warm and cold events over Eurasia. For example, in winter, Ural blocking highs or persistent Ural ridges induce cold advection downstream, enhance the Siberian High and East Asian winter monsoon, lead to frequent cold waves over East Asia (Takaya and Nakamura 2005;Zhou et al. 2009;Peng and Bueh 2012;Cheung et al. 2013; Wang and Lu 2017;Zhang Y et al. 2021;Sui et al. 2022;Yao et al. 2022). On the other hand, because one center of action for the Eurasian teleconnection pattern 1 3 is just around the Ural Mountains (Wallace and Gutzler 1981; Barnston and Livezey 1987;Liu et al. 2014), monthly anomalies of the Ural atmospheric circulation are also closely connected with the Eurasian teleconnection pattern, which is a dominant mode of intraseasonal atmospheric variability over Eurasia.
Over recent decades, a rapid reduction of sea ice cover and amplified surface warming have occurred in the Arctic, especially in the Barents-Kara Seas (BKS) (Comiso et al. 2008;Serreze et al. 2009;Cavalieri and Parkinson 2012;Dai et al. 2019;Isaksen et al. 2022). The surface net energy flux in the BKS is generally upward in autumn, winter and early spring (Serreze et al. 2007;Smedsrud et al. 2013); thus, reduced sea ice cover in the BKS in these seasons are often regarded as an external forcing to the surrounding atmospheric circulation (Deser et al. 2007;Petoukhov and Semenov 2010;Zhang P et al. 2018;Zhang and Screen 2021). Since the Ural region is just to the south of the BKS, to what extent the recent BKS sea ice loss influences the observed monthly circulation anomalies over the Ural region, hence influencing Eurasian climate anomalies, is of interest to researchers.
Numerous studies have already investigated the influences of the recent BKS sea ice loss and/or the Arctic warming on the Ural blocking, the Siberian High, and the recent surface cooling trend over central Eurasia (e.g., Cohen et al. 2014;Huang et al. 2017;Dai and Deng 2022;Luo et al. 2022), etc. For example, many studies suggested that the recent BKS sea ice loss and/or the Arctic warming favored more frequent and persistent Ural blocking highs (Yao et al. 2017;Huang et al. 2017;Luo et al. 2018Luo et al. , 2019Ma et al. 2018), as well as an anomalous local Arctic anticyclonic circulation (Xie et al. 2020(Xie et al. , 2022) and a strengthened Siberian High (Wu et al. 2011;Kug et al. 2015;Ma et al. 2018;Mori et al. 2019) in cold season, hence contributing to the recent increased cold events and surface cooling trend over central Eurasia. He et al. (2020) showed that Ural blockings and Eurasian below-average temperatures were likely more frequent in winters with deep Arctic warming than shallow, nearsurface warming over the BKS. Ma et al. (2022) found that SIC perturbations in the Greenland Sea, Barents Sea and Okhotsk Sea in winter were important for extended-range prediction of the strong and long-lasting Ural blocking formation, where negative (positive) SIC perturbations favored (inhibited) the Ural blocking formation four pentads later. Zhang R et al. (2018) suggested that heavy (light) SIC in the Barents Sea in spring-summer favored more frequent (rare) Ural blocking in summer.
Moreover, quite a few studies argued that the recent Arctic sea ice loss had little influence on the Siberian High and Eurasian cooling, they suggested that the recent increased cold events and Eurasian cooling probably were attributed to internal atmospheric variability (McCusker et al. 2016;Sun et al. 2016;Collow et al. 2018;Ogawa et al. 2018;Blackport et al. 2019). Dai and Song (2020) suggested that the greenhouse gases induced Arctic warming and sea ice loss were unlikely responsible for the recent Eurasian cooling.
Furthermore, several studies reported a nonlinear response of midlatitude atmospheric circulation to the BKS sea ice loss (i.e., a non-monotonic response of midlatitude atmospheric circulation to the monotonic BKS sea ice loss), they suggested that a specific amount of sea ice loss in winter or autumn considerably contributed to the severe winters over midlatitude Eurasia (Petoukhov and Semenov 2010;Semenov and Latif 2015;Zhang and Screen 2021). Chen et al. (2021) found that the frequency and local persistence of Ural blocking showed a non-monotonic change for a linear increased BKS warming in early winter.
Overall, the different and even inconsistent conclusions among previous studies mentioned above imply that the linkages between the BKS sea ice loss and the Ural circulation anomalies, and also Eurasian climate anomalies, are quite complex, which is essentially due to nonlinearities in the climate system (Overland et al. 2016;Cohen et al. 2020).
In fact, most previous studies drew conclusions by using common analysis approaches, such as linear correlation, linear regression, ensemble-mean model experiments, under an implicit assumption that the relationship between sea ice anomalies and atmospheric circulation anomalies was one-to-one correspondence. That is, positive and negative anomalies of atmospheric circulation were assumed to be related to opposite-sign anomalies of sea ice, respectively. Because of the different heat capacities between the ice-ocean system and the atmosphere, both positive and negative anomalies of the BKS SIC usually persist for several months and even several years (Fig. 1a), whereas both positive and negative anomalies of the Ural atmospheric circulation rarely persist for longer than a month (Fig. 1b;Dole and Gordon 1983;Miller et al. 2020). Consequently, during the period that the BKS SIC keeps a positive or negative anomaly, both positive and negative anomaly events of the Ural atmospheric circulation can reoccur several times. Hence, if one-to-one correspondence was assumed on the time scales such as subseasonal, seasonal and annual, then a considerable number of positive and negative anomaly events of the Ural atmospheric circulation were not thought to be related to the BKS SIC anomaly. For example, by calculating the correlation coefficients, one may conclude that the variation of BKS SIC in September, October and January are hardly related to the Ural circulation anomalies in the lower troposphere both in the same month and at one-month lag (Supplementary Figs. S1 ~ S4), although the BKS is very close to the Ural region. However, the one-to-one correspondence is not necessarily the unique assumption.
Recently Li et al. (2018Li et al. ( , 2020) (hereafter, L18 and L20) proposed the supercritical pitchfork bifurcation (SPB) model, which provided the example of one-totwo correspondence relationship. L18 found that in a low-order coupled land-atmosphere model, there may exist two flow equilibria under the same solar radiation and orographic forcing, when an external forcing was beyond a critical value (see details in Sect. 2.1). That is, there existed a one-to-two correspondence relationship between an external forcing and two flow equilibria in the model. L20 suggested that to some extent, positive and negative anomalies of the Ural atmospheric circulation in winter could be considered two equilibria under the same total solar irradiance (which was an external forcing), when the total solar irradiance was below the critical value 1360.9 W m −2 .
Apparently, compared with the common perspective of one-to-one correspondence, much more positive and negative anomaly events of the monthly Ural atmospheric circulation might be thought to be related to the BKS SIC anomaly from the perspective of oneto-two correspondence. However, L20 had not found that reduced SIC in the BKS had similar effect as the decreased total solar irradiance on the simultaneous Ural atmospheric circulation in winter. As the BKS SIC may be considered an external forcing in autumn, winter and early spring, we conjecture that in specific months of these seasons, there may exist a one-to-two correspondence relationship that generally agrees with the SPB model, hence referred to as a pitchfork-like relationship, between the BKS SIC and the lagging Ural circulation anomalies.
In this study, to validate the aforementioned conjecture, guided by the SPB model, we have examined the relationship between the observed BKS SIC and the lagging Ural atmospheric circulation on the monthly time scale during 1979/1980 − 2020/2021. Results verify that the conjecture is true. It is found that when the BKS SIC is below a critical value in specific months, most positive and negative anomaly events of the lagging Ural atmospheric circulation may be related to the same BKS SIC. This finding appears to be novel. It is apparently different from most previous studies based on the common perspective of one-to-one correspondence.
The remainder of this paper is organized as follows. In Sect. 2, we describe the SPB model, the data and methods used in this study. In Sect. 3, we provide examples that the relationship between the reduced BKS SIC and the Ural atmospheric circulation follows the pitchfork-like relationship. The characteristics of three preferred Ural circulation patterns are also presented. Section 4 is the discussion. In Sect. 4.1, we discuss the differences between the pitchfork-like relationship and the SPB model. In Sect. 4.2, we discuss the sensitivity of the pitchfork-like relationship to noise. In Sect. 4.3, possible mechanisms of the existence of the pitchfork-like relationship are explained. In Sect. 4.4, the possible influences of the regime transition of Ural atmospheric circulation on the predictability of Eurasian climate anomalies are discussed. The summary is presented in Sect. 5.

The SPB model
The SPB refers to a transition in the qualitative behavior of a dynamical system from a stable equilibrium to one unstable equilibrium and two stable equilibria, as some control condition is varied (Strogatz 2018). The fundamentals of SPB can be introduced by a simple example known as the buckling of the Euler rod (e.g., Kielhöfer 2011). Consider a vertical elastic rod subjected to a compressive axial force (Fig. 2a), the diagram of the relationship between the compressive axial force F and the amplitude of buckling x looks like a two-tine pitchfork ( Fig. 2b) (L20; Strogatz 2018). When F is gradually increased and less than the critical threshold F 0 , the rod remains straight state (x = 0), but the stability of the straight state is gradually decreased due to the stress accumulation in the rod. When the gradually increased F passes the critical threshold F 0 (termed a bifurcation point), the straight state of the rod becomes unstable, any small random perturbation leads to the left-side (x < 0) or right-side (x > 0) buckling. In other word, the neutral state converts to either the negative or positive anomaly state with equal probability. In the supercritical case (F > F 0 ), negative and positive anomaly states of the rod store the same elastic potential energy as a result of the work done by the compressive axial force and gravity. Therefore, the negative and positive anomaly states can be attributed to the same external force, indicating a one-to-two correspondence relationship. It is noteworthy that the increased stress in the rod in the subcritical case may be small and not statistically significant, but it should not be ignored. Because "not statistically significant" cannot explain why the rod finally buckles to either the left or the right. It is not a problem of the statistical significance but a problem of the dynamic stability here.
L18 considered a low-order coupled land-atmosphere model, which combined a two-layer quasi-geostrophic channel model and an energy balance model. Using highly truncated spectral expansions, multiple equilibrium solutions were obtained (Fig. 3a ~ c). The variables a , b and c are nondimensional spectral components of the lower layer streamfunction field ( 3 = a F 1 + b F 2 + c F 3 ). Here, F 1 = √ 2 cos y , F 2 = 2 cos(nx) sin y and F 3 = 2 sin(nx) sin y are nondimensional orthogonal functions, and n is the zonal wavenumber. The external energy input is the zonally symmetric solar radiation. The shortwave radiation flux anomaly absorbed by the land is R g = √ 2C g cos y (Fig. 4a), where the parameter C g is an indicator of meridional gradient of the shortwave radiation flux; thus, it indicates an external forcing to drive the zonal flow in the model. The topography in L18 is an ideal sinusoidal topography ( Fig. 3d ~ i).
In the SPB model of L18, the bifurcation point is at about C g = 50 W m −2 . For C g < 50 W m −2 , there is one stable equilibrium ( a = b = c = 0 ), called the neutraltype equilibrium or Hadley equilibrium. This equilibrium is characterized by zonally symmetric westerlies in the upper layer (Fig. 3d) and horizontally motionless in the lower layer (Fig. 3e). Probably due to the surface friction, zonally symmetric westerlies in the lower layer vanish; thus, the neutral-type equilibrium does not interact with the orography. The zonally symmetric westerlies in the upper layer correspond to positive shear vorticity ( − u∕ y > 0 ) over the northern channel and negative shear vorticity ( − u∕ y < 0 ) over the southern channel (Fig. 4b). When C g is increased from 40 to 48 W m −2 , the stability of the neutral-type equilibrium is supposed to be decreased, because both positive and negative shear vorticity of the westerlies are strengthened (Fig. 4c); meanwhile, vertical shear of the zonal wind is also strengthened (Fig. 4d). For 50 ≤ C g ≤ 54 W m −2 , the neutral-type equilibrium becomes unstable, and there exist two new stable equilibria, called the trough-type equilibrium ( a > 0, b < 0, c > 0 ) and the ridge-type equilibrium ( a < 0, b > 0, c < 0 ). For the trough-type equilibrium, there are wavy westerlies with troughs (cyclonic flows) located slightly west of the mountain crests in the lower layer (Fig. 3g), and wavy westerlies with troughs over the west side of the mountains in the upper layer (Fig. 3f). For the ridge-type equilibrium, there are wavy easterlies with inverted ridges (anticyclonic  (Kielhöfer 2011;Strogatz 2018). a A vertical elastic rod subjected to a compressive axial force F, with the amplitude of buckling x. b The diagram of the relationship between F and x. Solid (dotted) line indicates the stable (unstable) equilibrium flows) located slightly west of the mountain crests in the lower layer (Fig. 3i), and wavy westerlies with ridges over the west side of the mountains in the upper layer (Fig. 3h). The trough-type (ridge-type) equilibrium shows steady wavy westerlies (easterlies) in the lower layer probably is the result of a balance between the frictional torque and the mountain torque. As suggested by L18, the trough-type and ridge-type equilibria originate from the orographic instability of the neutral-type equilibrium. In addition, the ridge-type equilibrium becomes unstable when C g = 56 W m −2 and disappears when C g > 56 W m −2 . The troughtype equilibrium is always stable for C g ≥ 56 W m −2 . The ridge-type equilibrium is less stable than the trough-type equilibrium, probably because it has stronger baroclinicity than the latter (see Fig. 11e in L18).
The ideal sinusoidal topography in L18 resembles the real topography over the Ural region, where the latter is characterized by the elongated Ural Mountains lying between the East European Plain and the West Siberian Plain ( Fig. 1 in L20). L20 suggested that characteristics of the three flow equilibria (neutral-type, trough-type and ridge-type equilibria) were also very similar to those of the three Ural circulation patterns (neutral, negative geopotential height anomaly and positive geopotential height anomaly) on the monthly time scale. However, a diabatic forcing in observations, such as the BKS SIC, is much more complex than the model parameter C g . Firstly, an observed diabatic forcing is not necessarily independent of the Ural atmospheric circulation to be an external forcing. Secondly, an observed diabatic forcing varies not necessarily slowly enough to be considered a fixed forcing. Thirdly, an observed diabatic forcing in different months and seasons are not necessarily the dominant factor among multiple factors that affect the westerlies over the Ural region. These three reasons may be responsible for that L20 failed to detect the pitchfork-like relationship between the reduced BKS SIC Fig. 3 The supercritical pitchfork bifurcation model (replotted according to Li et al. 2018). a, b, c Multiple equilibrium solutions in the low-order coupled land-atmosphere system. The ordinate shows the nondimensional equilibrium solutions of a zonal component a , b wave component b and c wave component c of the streamfunction field in the lower layer atmosphere. Black, blue and red circles denote the stable neutral-type, trough-type and ridge-type equilibria, respectively. Crosses denote unstable equilibria. d, f, h The streamfunc-tion fields (contours, intervals are 2 × 10 7 m 2 s −1 ) in the upper layer atmosphere, e, g, i horizontal wind fields (vectors) in the lower layer atmosphere at C g = 50 W m −2 and n = 1.3 for the d, e neutral-type equilibrium, f, g trough-type equilibrium and h, i ridge-type equilibrium, respectively. Note that there is no wind vector in e. The colored background shows the topographic height, and the warm (cool) tone indicates positive (negative) regions. "H" and "L" represent high and low pressure, respectively and the simultaneous Ural atmospheric circulation on the seasonal time scale.

Data
We use the monthly SIC on a 1º × 1º latitude-longitude grid derived from the HadISST dataset (Rayner et al. 2003). The monthly mean horizontal wind, air temperature, sea level pressure and geopotential height on a 2.5º × 2.5º latitude-longitude grid are derived from the NECP-NCAR reanalysis (Kalnay et al. 1996), and they are the default data to analyze the atmospheric circulation. The monthly mean horizontal wind and air temperature on a 0.25º × 0.25º latitude-longitude grid derived from the ERA5 reanalysis (Hersbach et al. 2020) are used for comparison. The monthly mean surface net shortwave radiation flux, surface net longwave radiation flux, surface net sensible heat flux and surface net latent heat flux on a 94 × 192 latitude-longitude Gaussian grid derived from the NECP-NCAR reanalysis are also used. The sum of the four surface net fluxes is referred to as the surface net energy flux. In this study, we focus our analysis on the time period from 1979/1980 to 2020/2021.

Analysis methods
The relative vorticity (ζ) is a measure of the rotation of the horizontal wind. It is given by where u and v are the zonal and meridional wind, respectively; x and y are the zonal and meridional distance, respectively. The central difference method is used to calculate the relative vorticity at each latitude-longitude grid, except the North Pole and the South Pole. The grid intervals are Δx = a cos Δ and Δy = aΔ , where a = 6.371 × 10 6 m is the radius of Earth, and are latitude and longitude, respectively.
The change in the stability of the neutral-type equilibrium with respect to the parameter C g in the model of Li et al. (2018). a Shortwave radiation flux anomaly absorbed by the land ( R g = √ 2C g cos y , unit: W m −2 ) at C g = 40 W m −2 . b The shear vorticity ( U y = − u y , unit:10 -6 s −1 ) for the neutral-type equilibria at C g = 40 W m −2 . c The mean positive shear vorticity ( U + y , red dots) and mean negative shear vorticity ( U − y , blue dots) for the neu-tral-type equilibrium as a function of the parameter C g . The U + y and U − y are positive shear vorticity ( − u∕ y > 0 ) and negative shear vorticity ( − u∕ y < 0 ) averaged over the domain ( 0 ≤ x ≤ 2 , 0 ≤ y ≤ ) shown in (b), respectively. d The vertical shear of zonal wind ( − u∕ p , unit: 10 -5 m s −1 Pa −1 ) averaged over the domain ( 0 ≤ x ≤ 2 , 0 ≤ y ≤ ) for the neutral-type equilibrium as a function of the parameter C g According to previous studies (Simmons et al. 1983;Kosaka and Nakamura 2006), the local barotropic energy conversion (denoted as CK) and baroclinic energy conversion (denoted as CP) integrated from the p 0 = 1000 hPa level to the p 1 = 100 hPa level can be written as where g = 9.81 m s −2 is the acceleration of gravity, p is pressure, f is the Coriolis parameter, S p = RT (c p p) − T p is the static stability parameter, R = 287 J K −1 kg −1 is the gas constant for dry air, c p = 1004 J K −1 kg −1 is the specific heat at constant pressure, T is temperature. In this study, overbars and primes indicate the climatological-mean flow and the Ural circulation anomalies, respectively. A positive value of CK (CP) indicate that the conversion of kinetic energy (available potential energy) from the climatological-mean flow to the Ural circulation anomalies.
In general, the observed monthly mean atmospheric circulation is not in the "equilibrium state", and it is timevarying due to the existence of interference. Nevertheless, the atmospheric circulation may have a high frequency to remain quasi-stationary in the vicinity of a stable equilibrium, which is the so-called regime behavior (Hannachi et al. 2017). Examining the peaks in the probability density function (PDF) of a variate is a direct method to recognize the circulation regimes (Hannachi et al. 2017;L20). We calculate the joint PDF by using a MATLAB function named "ksdensity" (https:// www. mathw orks. com/ help/ stats/ ksden sity. html), which can return the smoothed probability density estimate for bivariate data, based on the specified bandwidth of the kernel smoothing window. The larger bandwidth, the smoother curve for the PDF estimated. The bandwidth (W) for different variates can be objectively estimated by the following formula (L20): where D is the range of the sample data, N is the sample size. The bandwidth can be rounded up to an appropriate integer or decimal.
To calculate the joint PDF, we need define two onedimensional variables to indicate the BKS sea ice and the Ural atmospheric circulation. The mean SIC is defined as the area-averaged monthly SIC over the domain (70° − 85°N, 20° − 80°E), which covers the area with large SIC variability in the BKS (Fig. 5a, b; Fig. 1a). Note that the climatological mean is retained in the mean SIC for each month. The ζ850 (4) W = D 10 log 10 N index is defined as the area-averaged monthly relative vorticity at 850 hPa, after subtracting the climatological mean for each month, over the domain (50° − 65°N, 40° − 80°E) which covers the Ural Mountains (Fig. 5c, d). It is seen that a 0 s −1 contour of the climatological-mean relative vorticity at 850 hPa crosses the Ural Mountains, indicating that the Ural region is a transition zone between the climatologicalmean positive relative vorticity over the high-latitudes and negative relative vorticity over the mid-latitudes (Fig. 5c, d).
We will estimate the joint PDF between the mean SIC and the ζ850 index to examine the pitchfork-like relationship, with the guidance of the SPB model. The BKS SIC can be significantly affected by the simultaneous Ural atmospheric circulation (e.g., Chen and Luo 2018); thus, the former is not independent of the latter. This problem can be addressed by estimating the joint PDF between the mean SIC and the lagging ζ850 index. Moreover, the variation of the mean SIC is relatively slow from January to May and from August to October, while it is fast from May to August and from October to December (Fig. 1a). In addition, the surface net energy flux in the BKS is generally upward in January to April and in September to December, while it is downward in May to August (Fig. 1a, Fig. 7 in Smedsrud et al. 2013). It implies that in January to April and in September to October that the mean SIC in the BKS may be considered a quasi-stationary diabatic heating factor. Furthermore, the variation of the ζ850 index from month to month is much faster than the mean SIC (Fig. 1b). Thus, the ζ850 index can be considered a time-dependent variable on the monthly time scale. We will estimate the joint PDF between the mean SIC and the lagging ζ850 index month by month. Therefore, the sample size N for the mean SIC and the ζ850 index are both always 42 (1979/1980 to 2020/2021). The sample data in different months have different ranges D, and thus the corresponding bandwidths W may be different.

Pitchfork-like relationship between the SIC in the BKS and the Ural atmospheric circulation
We have examined the relationship between the mean SIC and the ζ850 index at 0-month, 1-month, 2-months and 3-months lags (Supplementary Figs. S5 ~ S8, only the results for 0-month and 1-month lag are shown). We merely have found two typical examples of the pitchfork-like relationship between the mean SIC and the ζ850 index at 1-month lag. Figure 5e shows the scatter diagram for the mean SIC in September and the ζ850 index in October, and their joint PDF during 1979 − 2020. Although both the negative anomaly (ζ850 < 0) and positive anomaly (ζ850 > 0) events occur whatever the mean SIC is high or low, there exist different 1 3 circulation regimes with respect to the different mean SIC. As highlighted by the peak regions of the joint PDF with at least a total of 40% probability, the distribution of scattered points apparently shows a tilted pitchfork-like pattern, the critical value of mean SIC is about 16%. When the mean SIC in September is subcritical (SIC > 16%), there is only one peak region of the joint PDF, which is near the mean line of the ζ850 index in October, indicating that the "neutral state" is a preferred Ural circulation pattern or circulation regime. When the mean SIC in September is supercritical (SIC < 16%), the neutral state seems to become unstable, and the new peak region of the joint PDF looks like the less-than sign " < ", whose upper (lower) part is above (below) the mean line of the ζ850 index in October, indicating that both the positive and negative anomaly states are circulation regimes. In addition, when 5% < SIC < 16% in September, the frequencies of occurrence for the negative anomaly (ζ850 < 0) and positive anomaly (ζ850 > 0) events in October are both 14 (Table 1). It should be noted that when SIC < 5%, there are three points that are far away from the pitchfork-like pattern highlighted by the peak regions of the joint PDF ( Supplementary Fig. S6a). If these three points Examples of the pitchfork-like relationship. a The standard deviation of SIC (shadings, unit: %) and the 1% contour of climatological-mean SIC (cyan dotted lines) in September during 1979 − 2020. The black box indicates the region to define the mean SIC. b As in a, but for January during 1980 − 2021. c The climatological-mean relative vorticity at 850 hPa (ζ850, unit: 10 -6 s −1 ) in October during 1979 − 2020. Black lines denote 0 s −1 contours. The red box indicates the region to define the ζ850 index. d As in c, but for February during 1980 − 2021. Blue lines in a-d outline the Ural Mountains. e Scatter diagram for the mean SIC in September and the ζ850 index in October (dots) and their joint PDF (shadings, unit: 10 -4 ) during 1979 − 2020. Black contours highlight the peak probability density, and the areas enclosed by black contours are the probabilities of 20% and 40%. Horizontal blue line denotes the mean of the ζ850 indices. Vertical black dashed line indicates the critical value of the mean SIC. Note that the abscissa is reversed. Bandwidths of the kernel smoothing window for the mean SIC and ζ850 indices are 2% and 10 -6 s. −1 , respectively. f As in e, but for the mean SIC in January and the ζ850 index in February during 1980 − 2021 are taken into account, that is the SIC < 16% is considered, the frequencies of occurrence for the negative anomaly and positive anomaly events in October are 16 and 15, respectively (Table 1). These results indicate that the negative and positive anomaly events should occur with equal probability in the supercritical case.
Other example of the pitchfork-like pattern is seen in the relationship between the mean SIC in January and the ζ850 index in February during 1980-2021, the critical value of mean SIC is about 64% (Fig. 5f). When 50% < SIC < 64% in January, the frequencies of occurrence for the negative anomaly (ζ850 < 0) and positive anomaly (ζ850 > 0) events in February are 10 and 12, respectively ( Table 1). Note that when SIC < 50%, there are also several points that are far away from the pitchfork-like pattern highlighted by the peak regions of the joint PDF ( Supplementary Fig.  S6e). If these points are also taken into account, i.e., the SIC < 64% is considered, the frequencies of occurrence for the negative and positive anomaly events in February are 14 and 16, respectively (Table 1). These results indicate that the negative and positive anomaly events should be equally likely to occur in the supercritical case.
If the ζ850 index is calculated by using the ERA5 reanalysis, there are similar pitchfork-like patterns and the same critical values of the mean SIC ( Supplementary  Fig. S9a, b). In the supercritical case, the frequencies of occurrence for negative and positive anomaly events are still almost the same (Table 1).
Therefore, there is a robust relationship between the mean SIC in September (January) and the ζ850 index in October (February): (1) The peak regions of their joint PDF look like a two-tine pitchfork, resembling the bifurcation diagram of the SPB model (Fig. 3); (2) When the mean SIC is below a critical value, the frequencies of occurrence for the negative and positive anomaly events are nearly equal to each other during the same SIC range, agreeing with the SPB model who implies that the trough-type and ridge-type equilibria are equally likely to occur in the supercritical case. Hence, the term "pitchfork-like relationship" is apt to describe the relationship between the mean SIC in September (January) and the ζ850 index in October (February).
In addition, two examples of distorted or destroyed pitchfork-like pattern can be recognized in the relationship between the mean SIC in October and the ζ850 index in November ( Supplementary Fig. S6b), and between the mean SIC in February and the ζ850 index in March ( Supplementary Fig. S6f). It is noteworthy that the typical, distorted or destroyed pitchfork-like patterns are seen only for the mean SIC in September, October, January and February (Supplementary Figs. S5, S6). As we mentioned before, the mean SIC in these months can be considered a quasi-stationary diabatic heating factor.
Furthermore, if the simultaneous relationship between the mean SIC and the ζ850 index is considered, no pitchfork-like relationship is found in all months (Supplementary Figs. S7,S8). It is likely because the BKS SIC is not independent of the simultaneous Ural atmospheric circulation. Although a pitchfork-like pattern is seen in February (Fig. S8f), the negative and positive anomaly events are not equally likely to occur in the supercritical case: when 55% < SIC < 65% in February, the frequencies of occurrence for the negative and positive anomaly events are 6 and 11, respectively. Therefore, the simultaneous relationship between the mean SIC and the ζ850 index in February is not regarded as a pitchfork-like relationship.

Pitchfork-like relationship between the potential temperature over the BKS and the Ural atmospheric circulation
To further demonstrate the robustness of the pitchfork-like relationship shown in Fig. 5e and f, the relationship between a substitution variable of the mean SIC and the ζ850 index is examined in this subsection.
In September and January, there are large upward surface net energy fluxes in the Greenland Sea, Norwegian Sea and Barents Sea, while downward or small upward surface net energy fluxes in the Kara Sea and Eurasia (Fig. 6a, b). The  (16) 15 (16) Jan/Feb θ925 ≈ 265.5 K 265.5 ~ 270 K (> 265.5 K) 11 (15) 12 (13) regions with large upward surface net energy fluxes correspond to a warm tongue of potential temperature at 925 hPa there. Thermodynamic effect of the SIC reduction on the atmosphere is mainly through its influence on the surface net energy fluxes (Xie et al. 2022;Smedsrud et al. 2013), while the surface net energy fluxes are also influenced by the atmospheric processes (e.g., winds, clouds) and ocean processes (e.g., ocean currents) (Smedsrud et al. 2013). In other word, there is a coupled sea-ice-atmosphere system over the BKS. To indicate the thermodynamic property of this coupled sea-ice-atmosphere system, we define the areaaveraged monthly potential temperature at 925 hPa over the domain (70° − 80°N, 10° − 70°E) as the mean θ925. Figure 6c, d show the time series of the mean θ925 and the mean SIC in September during 1959 − 2020 and in January during 1960 − 2021, respectively. An out of phase relationship has been seen between the mean θ925 and the mean SIC since the late 1970s. The correlation coefficient between the two variables is − 0.62 (− 0.75) in September (January) during 1979 − 2020 (1980 − 2021), which is  Fig. 5e, f, respectively, but the abscissa is replaced by the mean θ925. Vertical black dashed line indicates the critical value of the mean θ925. Bandwidths of the kernel smoothing window for the mean θ925 and the ζ850 indices are always 0.6 K and 10 -6 s −1 , respectively statistically significant, passing the p < 0.01 t test. It indicates that the increased (decreased) potential temperature in the lower atmosphere over the BKS is highly associated with the reduced (increased) BKS SIC after the late 1970s. Hence, we can use the mean θ925 as a substitute for the mean SIC. However, the out of phase relationship is not seen before the late 1970s. The correlation coefficient between the two variables is − 0.28 (− 0.12) in September (January) during 1959 − 1978 (1960 − 1979), which does not pass the p < 0.1 t test. In addition, before the late 1970s, the mean SIC are high and almost subcritical in both September and January. It suggests that before the late 1970s, because the SIC in the BKS is considerably high, thermodynamic effect of the SIC change on the local lower atmosphere is quite weak. Moreover, uncertainty of the SIC data is smaller after 1979, because the satellite microwave data is available after October 1978 (Rayner et al. 2003). Thus, we have already focused our analysis on the time period after 1979.
Using the mean θ925 as a substitute for the mean SIC, the relationship between the mean θ925 in September and the ζ850 index in October shows a distorted pitchfork-like pattern, and the critical value of the mean θ925 is about 278.5 K (Fig. 6e). In the supercritical case (θ925 > 278.5 K), there are two separable peak regions of the joint PDF, the lower (upper) one is below (above) the mean line of the ζ850 index in October, indicating that both the negative and positive anomaly states are circulation regimes. The two separable peak regions do not look symmetrical. That is, although the negative and positive anomaly events can co-occur at the same θ925, the negative anomaly events can also occur at smaller θ925 than the positive anomaly events. Nevertheless, when 278.5 K < θ925 < 283 K in September, the frequencies of occurrence for the negative and positive anomaly events in October are both 15 (Table 1). If θ925 > 278.5 K is considered, the corresponding frequencies are 17 and 15, respectively. If the mean θ925 and ζ850 index are calculated by using the ERA5 reanalysis, the results are similar ( Supplementary Fig. S9c, Table 1). Thus, the relationship between the mean θ925 in September and the ζ850 index in October can be still regarded as a pitchforklike relationship.
In addition, a pitchfork-like pattern is obviously seen in the relationship between the mean θ925 in January and the ζ850 index in February, the critical value of the mean θ925 is about 265.5 K (Fig. 6f). When 265.5 K < θ925 < 270 K, the frequencies of occurrence for the negative and positive anomaly events are 9 and 12, respectively (Table 1). If θ925 > 265.5 K is considered, the corresponding frequencies are 15 and 13, respectively. If the ERA5 reanalysis is used, we obtain similar results ( Supplementary Fig. S9d, Table 1). Thus, the relationship between the mean θ925 in January and the ζ850 index in February is also considered a pitchfork-like relationship.
Because the mean θ925 is a substitution variable of the mean SIC, the above results indicate that the pitchfork-like relationships shown in Fig. 5e, f are robust.

The characteristics of three preferred Ural circulation patterns
All of above results indicate that there exist three preferred circulation patterns or circulation regimes over the Ural regions in October and February. The characteristics of three preferred Ural circulation patterns are shown by the composites of circulation fields (Figs. 7, 8). The neutral state is the composite of 42 samples derived from 1979 to 2020 in October and 1980 − 2021 in February; hence, it also represents the climatological-mean state. The positive (negative) anomaly state is the composite of samples when the ζ850 index is greater (less) than or equal to its + 1 (− 1) standard deviation, with a sample size of 8 (6) in October and 6 (6) in February.
In October, when the Ural atmospheric circulation is in the neutral state, there are flat zonal westerlies at 850 hPa and 500 hPa over Eurasia (Fig. 7a, d); the negative and positive relative vorticity are dominant over the southern and northern Ural region in the troposphere, respectively (Fig. 7g, Fig. 5c); the Siberian High, as indicated by the 1020 hPa contour of sea level pressure, covers the Central Asia, Mongolia and northern China (Fig. 7a). The positive anomaly state exhibits positive relative vorticity anomalies and cyclonic circulation anomalies at 850 hPa over the Ural Mountains (Fig. 7b), and wavy westerlies with a Ural trough at 500 hPa (Fig. 7e); positive relative vorticity anomalies over the Ural Mountains extend from surface to the lower stratosphere, and the maximum anomalies occur at 250 ~ 300 hPa (Fig. 7h); meanwhile, the equivalent barotropic structure of positive relative vorticity anomalies is accompanied by a cold core in the troposphere and a warm core in the lower stratosphere (Fig. 7h); the Siberian High is weak-than-normal, with a southward shrinkage (Fig. 7b); there is a surface cold anomaly in the Ural region and a surface warm anomaly in the region south of Lake Baikal (Fig. 7e). By contrast, the negative anomaly state is presented by negative relative vorticity anomalies and anticyclonic circulation anomalies at 850 hPa over the Ural Mountains (Fig. 7c), and wavy westerlies with a Ural ridge at 500 hPa (Fig. 7f); negative relative vorticity anomalies over the Ural Mountains also extend from surface to the lower stratosphere, and the maximum anomalies occur near the tropopause (Fig. 7i); the equivalent barotropic structure of negative relative vorticity anomalies corresponds to a warm core in the troposphere and a cold core in the lower stratosphere (Fig. 7i); the Siberian High is stronger-than-normal, with a northward expansion (Fig. 7c); the surface temperature is warmer-than-normal in the northern Ural region and colder-than-normal in the region south of Lake Baikal and in Western Europe (Fig. 7f).
Similar characteristics are seen in February (Fig. 8), but all of the extents of Siberian High in the three states (Fig. 8a, b, c) are larger than those in October. Besides, the intensities of both Ural trough (Fig. 8e) and Ural ridge (Fig. 8f) in the two anomaly states are stronger than those in October. For the neutral state, the negative relative vorticity is dominant over the center of Ural Mountains, while the positive relative vorticity is dominant over the south and north of Ural Mountains (Fig. 8g, Fig. 5d). For the negative anomaly state, the characteristics of temperature field are different from those in October. In February, the warm core in the troposphere is located over the northern Ural region and the BKS (Fig. 8i); the surface warm anomaly occurs in the BKS and is not statistically significant, while the significant surface cold anomaly occurs in a belt area of 40° − 60°N (Fig. 8f).
It is notable that characteristics of the horizontal wind field for the three flow equilibria in the SPB model (Fig. 3 d ~ i) are similar to those of the three preferred Ural circulation patterns, namely, the neutral state, positive and negative anomaly states (Figs. 7 a ~ f, 8 a ~ f; also see Fig. 2 in L20). Of course, the two-layer model in L18 is not sufficient Fig. 7 Composites of circulation fields for three preferred Ural circulation patterns in October during 1979 − 2020: a, d, g neutral state; b, e, h positive anomaly state; c, f, i negative anomaly state. a 1020 hPa contour of the sea surface pressure (SLP) and wind (vectors) at 850 hPa. b, c 1020 hPa contour of the SLP, wind anomaly (vectors) and relative vorticity anomaly (shadings, unit: 10 -6 s −1 ) at 850 hPa. Wind vectors in b-c are plotted if they are significant in at least one direction (zonal or meridional). d Geopotential height at 500 hPa (Z500, black contours, unit: gpm). e, f Z500 and temperature anomaly at 1000 hPa (shadings, unit: ℃). Pressure-latitude cross section along the 60°E for (g) the relative vorticity (shadings, unit: 10 -6 s −1 ; black lines denote 0 s −1 contours), h, i relative vorticity anomaly (shadings, unit: 10 -6 s −1 ) and temperature anomaly (contours, unit: °C; negative contours are dashed). Black stippled areas in b, c, e, f, h, i indicate the anomalies are statistically significant exceeding the 95% confidence level (p < 0.05) based on two-sided Student's t-test. Blue lines in a-f outline the Ural Mountains. Two vertical blue dashed lines in g-i indicate the latitudinal extent of the Ural Mountains  Positive and negative anomaly events occur nearly at the same frequencies to depict more realistic vertical structure of the three flow equilibria.
Overall, our results show that there exists a pitchforklike relationship between the BKS SIC in September (January) and the Ural atmospheric circulation in October (February) during 1979 − 2020 (1980 − 2021). It indicates that when the gradually reduced BKS SIC in September (January) passes a critical threshold, a regime transition of Ural atmospheric circulation from the neutral regime to two anomaly regimes occurs in October (February). The regime transition corresponds to a transition of preferred Ural circulation patterns from the flat zonal westerlies to wavy westerlies with a Ural trough and wavy westerlies with a Ural ridge. After the regime transition, there exists a one-to-two correspondence relationship between the BKS SIC and the Ural circulation anomalies. In addition, as summarized in Table 2, there are many similarities between the pitchfork-like relationship and the SPB model.

Differences between the pitchfork-like relationship and the SPB model
We use the term "pitchfork-like relationship" instead of the term "supercritical pitchfork bifurcation (SPB)" to describe the relationship between the BKS SIC in September (January) and the Ural atmospheric circulation in October (February), because of some differences between the two concepts. Firstly, data of the mean SIC and the ζ850 index, as well as the mean θ925 and the ζ850 index, do not satisfy the Eqs. (19) − (27) in L18 that describes the SPB model. Clearly, the mean SIC and mean θ925 are not the model parameter C g , and the ζ850 index is not the geostrophic streamfunction ( 1 , 3 ) in the model. Hence, the pitchfork-like pattern does not necessarily look perfectly symmetrical (e.g., Fig. 6e). In addition, we use the statement "critical value of the mean SIC or θ925" instead of the term "bifurcation point". Because in the SPB model, the bifurcation point refers to the critical threshold of the parameter C g , at which new, non-trivial solutions (i.e., trough-type and ridge-type equilibria) of the Eqs. (19) − (27) are generated. In the joint PDF diagrams, a critical value of the mean SIC or θ925 is the value at which the peak regions that represent the positive anomaly and/or negative anomaly are generated. Hence, a critical value of the mean SIC or θ925 is determined by eye and is not necessarily accurate enough.
Secondly, the SPB model indicates the absence of neutral-type equilibrium (trough-type and ridge-type equilibria) in the supercritical (subcritical) case. However, the monthly mean Ural atmospheric circulation usually is time-varying and not in the "equilibrium state". The monthly mean Ural atmospheric circulation can be influenced by multiple factors, such as North Atlantic sea surface temperature (e.g., Li 2004;Sato et al. 2014), Eurasian snow cover (e.g. Han and Sun 2018), Arctic Oscillation (Thompson and Wallace 2000) besides the BKS SIC. Probably due to the existence of interference, the scatter diagrams in Fig. 5e, f show that the neutral, positive anomaly and negative anomaly events can occur in both the subcritical and supercritical cases. Nevertheless, in the subcritical (supercritical) case, the monthly mean Ural atmospheric circulation exhibits a high frequency to remain in the vicinity of one (two) "equilibrium state(s)", which is the so-called regime behavior and highlighted by the peak regions of the joint PDF.
In addition, on the daily time scale, the Ural atmospheric circulation usually is far away from the equilibrium state, because of the major interference from the high-frequency traveling Rossby waves. The highfrequency traveling Rossby waves can be partially filtered out by averaging the daily ζ850 indices, such as 10-day, 20-day and 25-day average. In fact, the joint PDF between the monthly-mean SIC in September (January) and the daily-mean ζ850 index in October (February) does not show a pitchfork-like pattern ( Supplementary Fig. S10), and so does not the joint PDF between the monthly-mean SIC and the 10-day or 20-day average ζ850 index ( Fig. S11 a ~ d). By contrast, the joint PDF between the monthlymean SIC in September (January) and the 25-day average ζ850 index in October (February) shows a distorted pitchfork-like pattern ( Supplementary Fig. S11e, f). It implies that the pitchfork-like relationship (Fig. 5e, f) may be a phenomenon on the time scale of one month.
Thirdly, the relationship between the external forcing parameter C g and three flow equilibria in the SPB model is a driver-and-response or cause-and-effect relationship. However, just by analyzing the scatter and joint PDF diagrams, whether the pitchfork-like relationship between the BKS SIC and three preferred Ural circulation patterns is a cause-and-effect relationship is unclear.

The sensitivity of the pitchfork-like relationship to noise
We have used the mean θ925 as a substitute for the mean SIC to demonstrate the robustness of the pitchfork-like relationship. One may guess that if the abscissa (ordinate) in Fig. 5e and f were replaced by any other variables that were highly correlated to the mean SIC (the ζ850 index), there were still recognizable pitchfork-like patterns; hence, it would suggest that the pitchfork-like relationship may be just a correlation relationship. Some unexpected, failed examples are given below to illustrate that the pitchfork-like relationship is very sensitive to noise, and thus the guess of correlation relationship has not been confirmed. We have added some white noise to the time series of the mean SIC (Fig. 9a, b). Correlation coefficients between the time series interfered with noise and the original time series are both 0.97 in September and January (passing the p < 0.01 t test), which are apparently higher than the correlation coefficients between the mean θ925 and the mean SIC (Fig. 6c, d). However, peak regions of the joint PDF between the mean SIC interfered with noise and the ζ850 index do not look like a two-tine pitchfork (Fig. 9c, d), suggesting that white noise in the time series of the mean SIC causes signal distortion. By contrast, a pitchfork-like pattern is apparently recognizable in the relationship between the mean θ925 and the ζ850 index (Fig. 6c, d), suggesting that the time series of mean θ925 has a higher signal-noise ratio. It indicates that the pitchfork-like relationship between the BKS SIC and the Ural atmospheric circulation is very sensitive to noise.
We also have defined and calculated other variables, based on the surface net shortwave radiation flux, surface net longwave radiation flux, surface net sensible and latent heat flux, and surface net energy flux. Some variables are significantly correlated to the mean SIC or the mean θ925 (Supplementary Table S1). However, the corresponding joint PDF diagrams do not show any pitchfork-like pattern ( Supplementary Fig. S12), except that the joint PDF between the mean shortwave flux in September and the ζ850 index in October shows a plausible distorted pitchfork-like pattern (Fig. S12a), but the mean shortwave flux is not significantly correlated to both the mean SIC and the mean θ925 (Table S1). These results also suggest that the pitchfork-like relationship is very sensitive to noise.
In addition, similar to the definition of ζ850 index, we have defined four new monthly relative vorticity indices over the domain (50° − 65°N, 40° − 80°E) but at 925 hPa, 700 hPa, 500 hPa and 300 hPa, respectively. The correlation coefficients between the ζ850 index and four new variables, namely, ζ925 index, ζ700 index, ζ500 index and ζ300 index are respectively 0.99, 0.98, 0.93 and 0.87 in October, and are respectively 0.99, 0.99, 0.93 and 0.87 in February (passing the p < 0.01 t test). Comparing Fig. 10 with Fig. 5e, f, we find that if the relative vorticity index is calculated at 925 hPa and 850 hPa, a typical pitchfork-like pattern can be seen in the relationship between the mean Fig. 9 Time series of the original mean SIC (blue lines) and the mean SIC interfered with white noise (red lines) in a September during 1979 − 2020 and b January during 1980 − 2021. c, d As in Fig. 5e, f, respectively, but the abscissa is replaced by the mean SIC interfered with white noise. Bandwidths of the kernel smoothing window for the mean SIC and ζ850 indices are always 2% and 10 -6 s −1 , respectively 1 3 SIC and the relative vorticity index (Fig. 10a, b; Fig. 5e, f). By contrast, if the relative vorticity index is calculated at lower pressure levels (higher altitude layers) such as 700 hpa, 500 hPa and 300 hPa, the pitchfork-like patterns are distorted and even destroyed (Fig. 10 c ~ h). It may because the signal source of the pitchfork-like relationship is located at the surface and lower troposphere, and the signal becomes weaker at lower pressure levels than 850 hPa. These results also indicate that the pitchfork-like relationship is highly sensitive to noise.

Possible mechanisms of the existence of the pitchfork-like relationship
In this subsection, we will investigate the possible mechanisms of the existence of the pitchfork-like relationship between the BKS SIC in September (January) and the Ural atmospheric circulation in October (February).
The key point is to explain that with respect to the BKS SIC reduction in September (January), why there exists the regime transition of monthly Ural atmospheric circulation from one regime to two regimes in October (February).
The SPB model can provide guidance, because of the similarities between the pitchfork-like relationship and the SPB model (Table 2). In the SPB model, in the subcritical case, an increase in the parameter C g causes an increase in both the positive and negative shear vorticity of the westerlies, as well as the vertical shear of the zonal wind (Fig. 4), hence decreasing the stability of the neutral-type equilibrium until losing stability. Losing stability of the neutraltype equilibrium is the cause of the bifurcation, i.e., equilibrium transition. In the supercritical case, the trough-type and ridge-type equilibria can coexist (Fig. 3), suggesting that both the troughs and ridges near the mountains can extract Fig. 11 Linear change in the temperature at 1000 hPa (unit:°C; left panels), zonal wind at 850 hPa (unit: m s −1 ; middle panels; black lines denote 0 m s −1 contours) and shear vorticity of zonal wind at 850 hPa (unit: 10 -6 s −1 ; right panels; black lines denote 0 s −1 contours) when the mean SIC decreases by 5%. The regressed atmospheric variables against the simultaneously decreased mean SIC in (a, b, c) September during 1979-2020, (g, h, i) January during 1980-2021. (d, e, f) As in (a, b, c), respectively, but for atmospheric variables in October against the SIC in September. (j, k, l) As in (g, h, i), respectively, but for atmospheric variables in February against the SIC in January. Cyan dotted lines denote the 1% contour of climatological-mean SIC in a September, d October, g January, j February. Cyan lines (crosses) in middle panels denote the climatological-mean 0 m s −1 contours of zonal wind (easterly wind regions). Black stippled areas indicate the regression is statistically significant exceeding the 95% confidence level (p < 0.05) based on two-sided Student's t-test. Blue lines outline the Ural Mountains energy from the basic westerly flow to maintain themselves. This is responsible for the existence of one-to-two correspondence relationship between the parameter C g and the two anomaly equilibria in the supercritical case. We find that the existence of the pitchfork-like relationship between the BKS SIC and the Ural atmospheric circulation has similar mechanisms. The results are presented as follows.

Linear change in the background atmospheric conditions with respect to the BKS SIC reduction
In this subsection, the role of the BKS SIC reduction in change the background atmospheric conditions are explained. In Fig. 11, when the mean SIC decreases by 5%, the change in the temperature at 1000 hPa, the zonal wind and the shear vorticity of zonal wind at 850 hPa are estimated by using the linear regression. With the decrease in the mean SIC in September and January, in the areas north of the 1% contour of climatological-mean SIC, there is a significant surface warming in the same month (Fig. 11a, g) and one month later (Fig. 11d, j). It suggests that because the BKS SIC anomaly in September (January) can persist into October (February) (Fig. 1a), the thermodynamic effect of the SIC reduction contributes to a similar local surface warming in September and October (January and February). It is consistent with Dai et al. (2019) who pointed out that the amplified surface warming only occurred over the areas with significant sea ice loss. The surface warming in the Arctic causes a reduced meridional temperature gradient, hence giving rise to a weakened westerly wind at 850 hPa over north Eurasia (Fig. 11b, e, h, k). Focus on the Ural region, it is seen that the linear change in the westerly wind at 850 hPa corresponds to an increase in the shear vorticity over the Ural Mountains, and a decrease in the shear vorticity over the regions south and north of Ural Mountains (Fig. 11c, f, i, l). Thus, the SIC reduction in the BKS contributes to similar background shear vorticity of the westerly wind over the Ural region in September, October, January and February.
Next, we focus on the linear change in the shear vorticity of zonal wind in October (February) when the mean SIC in September (January) decreases by 5% in the subcritical case. With the decrease in the mean SIC in previous month in the subcritical case, in October there is still an increase in the shear vorticity over the Ural Mountains, and a decrease in the shear vorticity over the regions south and north of Ural Mountains (Fig. 12a). These wind shear components correspond to an equivalent barotropic structure in the vertical direction (Fig. 12c), except that over the southern Ural Mountains, the decreased shear vorticity is confined to the lower troposphere from surface to 700 hPa, while the increased shear vorticity extends from 700 hPa to lower stratosphere. By contrast, in February there is an increase in the shear vorticity over the northern Ural Mountains and the region to the south of Ural Mountains, while a decrease in the shear vorticity over the southern Ural Mountains and the region to the north of Ural Mountains (Fig. 12b). These wind shear components extend from surface to lower stratosphere (Fig. 12d). Note that when the Ural atmospheric circulation is in the neutral state (i.e., climatological-mean state) in both October and February, there are respectively negative and positive shear vorticity over the southern and northern Ural Mountains in the troposphere, which generally resembles Fig. 4b. Thus, with the decrease in BKS SIC in previous month in the subcritical case, in both October and February the spatial distribution of the linear change in shear vorticity in the lower troposphere over the Ural region are roughly in phase with those in the neutral state, indicating the increased both negative and positive shear vorticity of the westerly wind over the southern and northern Ural Mountains, respectively.
To further investigate the relationship between the BKS SIC reduction and the change of shear vorticity over the Ural region, we define the area-averaged monthly positive shear vorticity ( − u∕ y > 0 ) and negative shear vorticity ( − u∕ y < 0 ) of zonal wind at 850 hPa over the domain (50° − 65°N, 40° − 80°E) as the mean U + y and mean U − y , respectively. As shown in Fig. 12e and f, a reduction in the mean SIC in September (January) corresponds to an increase in both the mean U + y and mean U − y in October (February). This result is similar to Fig. 4c.
In addition, with the decrease in the mean SIC in previous month, the linear change in the vertical shear of zonal wind over the Ural region does not show consistent results in October and February ( Supplementary Fig. S13). Hence, we do not discuss them.
Although the increased both positive and negative shear vorticity related to the reduced BKS SIC shown in Fig. 12 are small and not statistically significant, they should not be ignored. According to the SPB model, the increased both positive and negative shear vorticity in the subcritical case might decrease the stability of the flat zonal westerlies until losing stability; if dynamic instability of the flat zonal westerlies was triggered in the presence of any random initial perturbation, the Ural trough and ridge might be produced with equal probability, which might lead to a regime transition from the neutral regime to two anomaly regimes.
However, the observed atmosphere is more complex than the SPB model. Because the regressed shear vorticity is not statistically significant, it is difficult to distinguish the small change in the background shear vorticity linked to the reduced SIC from that linked to noise. In addition, because both the Ural trough and ridge can occur whatever the BKS SIC is high or low (Fig. 5e, f), the possible explanation on the production of the Ural trough and ridge above is not sufficient to explain the occurrence of the regime transition. Fig. 12 a-d Linear change in the shear vorticity of zonal wind (U y , unit: 10 -6 s −1 ) when the mean SIC decreases by 5% in the subcritical case. a The regressed U y at 850 hPa in October against the decreased mean SIC in September when SIC > 16% during 1979-2020. Black lines denote 0 s −1 contours. Cyan lines (crosses) denote the climatological-mean 0 s −1 contours of U y (negative U y regions). The red box indicates the region to define the mean positive shear vorticity ( U + y ) and mean negative shear vorticity ( U − y ). c As in a, but for the pressure-latitude cross section of the regressed U y along the 60°E. (b, d) As in (a, c), respectively, but for the regressed U y in February against the decreased mean SIC in January when SIC > 64% during 1980-2021. Black stippled areas indicate the regression is statistically significant exceeding the 95% confidence level (p < 0.05) based on two-sided Student's t-test. Blue lines in a-b outline the Ural Mountains. Two vertical blue dashed lines in c-d indicate the latitudinal extent of the Ural Mountains. e Scatter diagram for the mean SIC in September and the mean shear vorticity (red dots for U + y , blue dots for U − y ) in October in the subcritical case during 1979-2020. The red and blue lines are the corresponding least squares fitting straight lines, and s denotes its slope. The r denotes the correlation coefficient between the corresponding two variables. Note that the abscissa is reversed. f As in e, but for the mean SIC in January and the mean shear vorticity in February in the subcritical case during 1980-2021. All fitting lines (correlation coefficients) do not statistically significant at the 90% confidence level for the F-test (t test)

Changes in the transfer of energy to the Ural circulation anomalies with respect to the BKS SIC reduction
In this subsection, with respect to the BKS SIC reduction, changes in the transfer of energy from the climatologicalmean flow to the Ural circulation anomalies are presented. Figure 13 shows the vertically integrated local barotropic and baroclinic energy conversion (denoted as CK and CP, respectively) for the composited positive and negative anomaly events of Ural atmospheric circulation in October and February. The positive and negative anomaly events selected are the same as those in Figs. 7 and 8.
In October, for the positive anomaly state characterized by wavy westerlies with a Ural trough at 500 hPa, it mainly extracts kinetic energy from the climatological-mean flow over the regions west of the Ural Mountains and west of Lake Baikal, and it mainly releases kinetic energy into the Fig. 13 (a, b) The vertically integrated local barotropic energy conversion (CK, shadings, unit: W m −2 ) and the geopotential height at 500 hPa (contours, unit: gpm), (c, d) the vertically integrated local baroclinic energy conversion (CP, shadings, unit: W m −2 ) and the geopotential height at 500 hPa (contours) for the composited (a, c) positive anomaly and (b, d) negative anomaly events in October dur-ing 1979 − 2020. e-h As in a-d, respectively, but for February during 1980 − 2021. All of the contour intervals are 100 gpm. The contour ranges are (a-d) 5300 ~ 5700 gpm and (e-h) 5100 ~ 5500 gpm. Cyan boxes indicate the region to define the mean CK + and mean CP + . Blue lines outline the Ural Mountains climatological-mean flow over the northern Europe, and the regions northeast and southeast of the Ural Mountains (Fig. 13a); meanwhile, it mainly gains (releases) available potential energy from (into) the climatological-mean flow over the region west of the Ural Mountains (northern Europe and the region east of the Ural Mountains) (Fig. 13c). Besides, the baroclinic conversion of energy has contributed more than the barotropic conversion of energy to maintain the positive anomaly state. For the negative anomaly state characterized by wavy westerlies with a Ural ridge at 500 hPa, the results are generally similar, except that the kinetic energy and available potential energy are mainly extracted from the climatologicalmean flow over the region northwest of the Ural Mountains (Fig. 13b, d). In fact, for the cyclonic circulation anomalies accompanied by a cold core (Fig. 7b, e, h) and the anticyclonic circulation anomalies accompanied by a warm core (Fig. 7c, f, i), they correspond to the same sign of the terms (v �2 − u �2 ) 2 and u ′ v ′ ( v ′ T ′ and u ′ T ′ ) over the Ural region, hence leading to similar spatial distributions of the CK (CP) for both the positive and negative anomaly states. In addition, the spatial distributions of the CK and the CP for the negative anomaly state (Fig. 13b, d) are roughly consistent with Fig. 7 in Shi et al. (2022), who analyzed the maintenance mechanism of the winter blocking highs around the Ural Mountains.
In February, the spatial distributions of the CK for both the positive (Fig. 13e) and negative (Fig. 13f) anomaly states, and the CP for the positive anomaly state (Fig. 13g) are roughly similar to those in October, except that more kinetic energy and available potential energy are extracted from the climatological-mean flow to maintain the Ural circulation anomalies with a larger wave amplitude in February. However, for the negative anomaly state in February, because the anticyclonic circulation anomalies are not accompanied by a warm core in the troposphere over the central Ural region (Fig. 8f, i), the available potential energy that extracted from the climatological-mean flow over the region west of the Ural Mountains is less than that over the region east of the Ural Mountains (Fig. 13h).
The above results suggest that to maintain both positive and negative anomalies of Ural atmospheric circulation, the region west of the Ural Mountains is the common region that kinetic energy and available potential energy are  Fig. 5e, f, respectively, but the ordinate is replaced by (a, b) the mean CK + , (c, d) the mean CP + . Red and cyan dots indicate the positive anomaly events (ζ850 > 0) and negative anomaly events (ζ850 < 0), respectively. Horizontal blue lines indicate 0 W m −2 . Note that the ordinates are non-negative, which indicate the conversion of (kinetic or available potential) energy from the climatological-mean flow to the positive or negative anomaly of Ural atmospheric circulation. Vertical black dotted lines in (a, c) indicate the SIC = 16%, in (b, d) indicate the SIC = 64%. Bandwidths of the kernel smoothing window for the mean CK + , CP + and SIC are always 0.1 W m −2 , 0.6 W m −2 and 2%, respectively extracted from the climatological-mean flow. We define the area-averaged monthly positive CK (CP) over the domain (45° − 75°N, 30° − 60°E) as the mean CK + (CP + ). Figure 14a (Fig. 14c) shows the scatter diagram for the mean SIC in September and the mean CK + (CP + ) in October, and their joint PDF. When the mean SIC in September is gradually reduced, although the positive and negative anomaly events in October occur whatever the mean SIC is high or low, both the mean CK + and mean CP + for the positive and negative anomaly events exhibit a regime transition from one regime to two regimes. In the subcritical case (SIC > 16%), the unipolar peak region of the joint PDF shown in Fig. 14a, c indicate that a weak positive anomaly is the preferred Ural circulation pattern that extracts kinetic energy and available potential energy from the climatological-mean flow to maintain itself, which corresponds to the "neutral" regime shown in Fig. 5e. By contrast, in the supercritical case (SIC < 16%), the dipole peak regions of the joint PDF indicate that both positive and negative anomalies are preferred Ural circulation patterns that extract kinetic energy and available potential energy from the climatological-mean flow to maintain themselves, which correspond to the positive and negative anomaly regimes shown in Fig. 5e.
Similarly, with respect to the SIC reduction in January, both the mean CK + and mean CP + for the positive and negative anomaly events in February also exhibit a regime transition from one regime to two regimes (Fig. 14b, d), which correspond to the regime transition of preferred Ural circulation patterns shown in Fig. 5f. However, for the mean CP + , the lower part of the "dipole" peak regions of the joint PDF in the supercritical case (SIC < 64%) is partially destroyed (Fig. 14d). As mentioned above, it may be ascribed to the lack of a warm core in the troposphere over the central Ural region for the composited anticyclonic circulation anomalies in February (Fig. 8f, i).
Moreover, in other months of autumn and winter, both the mean CK + and mean CP + for the positive and negative anomaly events do not show a regime transition from one regime to two regimes (Supplementary Figs. S14, S15).
Therefore, it is suggested that when the gradually reduced BKS SIC in September (January) passes a critical threshold, both the local barotropic and baroclinic energy conversion for Ural circulation anomalies in October (February) exhibit a regime transition from one regime to two regimes, leading to the corresponding transition of preferred Ural circulation patterns. After the regime transition, under the same BKS SIC, both positive and negative anomalies are preferred Ural circulation patterns that extract kinetic energy and available potential energy from the climatological-mean flow to maintain themselves; hence, there exists a one-totwo correspondence relationship between the BKS SIC and the Ural circulation anomalies (Fig. 5e, f). According to the SPB model, the regime transition of the local barotropic and baroclinic energy conversion might be associated with the increased both positive and negative shear vorticity of background westerly wind over the Ural region before the regime transition.
To further investigate the mechanisms of the regime transition of monthly Ural atmospheric circulation, we need to investigate that what have changed for the Ural atmospheric circulation on the time scale shorter than one month and why the changes happen. They are beyond the scope of this paper. Luo et al. (2018Luo et al. ( , 2019 suggested that the BKS warming associated with the sea ice loss led to more persistent Ural blocking in winter, through reduced background potential vorticity gradients over mid-to-highlatitude Eurasia. Hence, the connections between the regime transition of monthly Ural atmospheric circulation and the changes in the duration of blocking highs and cut-off lows deserve an investigation in the future.

Possible influences of the regime transition on the predictability of Eurasian climate anomalies
The regime transition of Ural atmospheric circulation found in this study appears to be new. Next, we briefly discuss the possible influences of the regime transition on the predictability of Eurasian climate anomalies. The Ural atmospheric circulation anomalies on monthly time scale are key indicators of regional warm and cold events over Eurasia (e.g., Zhou et al. 2009;Wang and Lu 2017). We have confirmed that when the mean SIC is supercritical, the positive and negative anomaly events of Ural atmospheric circulation should occur with equal probability (Table 1). It implies that the occurrence of Ural atmospheric circulation anomalies in the supercritical case may be related to not only the sea ice anomalies, but also the random initial perturbations in the atmosphere. The latter is unpredictable. Noteworthily, the mean SIC and mean θ925 in September and January since 2005 are nearly both supercritical (Fig. 6c, d). It implies that the predictabilities of Ural atmospheric circulation anomalies in October and February, under the BKS SIC reduction since 2005, may be low.
In addition, similar to the definition of ζ850 index, we define the area-averaged geopotential height anomaly at 500 hPa over the domain (50° − 70°N, 40° − 80°E) as the Z500 index. Although the pitchfork-like pattern between the mean SIC in September and the Z500 index in October (Fig. 15a), and between the mean SIC in January and the Z500 index in February (Fig. 15b), are both quite distorted (sensitivity to noise), the composite of 500 hPa geopotential height field in the subcritical case (Fig. 15c, d) and supercritical case (Fig. 15e, f) both show flat zonal westerlies over the Ural region, with weak and not statistically significant geopotential height anomalies there. It implies that the predictabilities of 500 hPa geopotential height anomalies over the Ural region, related to the recent BKS SIC reduction, may be low in October and February if the ensemble-mean prediction is used.
The Siberian High (SH) anomalies on monthly time scale are also key indicators of Eurasian climate anomalies (e.g., Gong and Ho 2002;Wang and Lu 2017). We define the areaaveraged sea level pressure over the domain (50° − 65°N, 45° − 85°E) as the SH intensity, which represents the intensity variation of the SH associated with the Ural circulation variability. We also define the total number of the 2.5º × 2.5º latitude-longitude grid points enclosed by the 1020 hPa contour of sea level pressure over the domain (35° − 75°N, 0° − 150°E) as the SH extent, which represents the extent variation of the SH over Eurasia in the mid-high latitudes.
Interestingly, a distorted pitchfork-like pattern is seen in the relationship between the mean SIC in September and the SH intensity in October (Fig. 16a), but no pitchfork-like pattern can be recognized in the relationship between the mean SIC in January and the SH intensity in February (Fig. 16b). By contrast, no pitchfork-like pattern can be recognized in the relationship between the mean SIC in September and the  Fig. 5e, f, respectively, but the ordinate is replaced by the Z500 index. Vertical black dotted lines in a indicate the SIC = 16%, in b indicate the SIC = 64%. Bandwidths of the kernel smoothing window for the mean SIC and the Z500 indices are always 2% and 15 gpm, respectively. Composites of the geopotential height at 500 hPa (black contours, unit: gpm) and its anomaly (shadings) in October in c the subcritical case (SIC > 16% in Septem-ber) and e supercritical case (5% < SIC < 16% in September) during 1979 − 2020, in February in d the subcritical case (SIC > 64% in January) and f supercritical case (50% < SIC < 64% in January) during 1980 − 2021. Black stippled areas indicate the anomalies are statistically significant exceeding the 95% confidence level (p < 0.05) based on two-sided Student's t-test. Blue lines outline the Ural Mountains SH extent in October (Fig. 16c), but an apparent pitchforklike pattern is seen in the relationship between the mean SIC in January and the SH extent in February (Fig. 16d). It suggests that the regime transition of preferred Ural circulation patterns mainly influences the intensity of the SH in October and the extent of the SH in February. Thus, when the mean SIC is supercritical, the intensity (Fig. 16a) or the extent (Fig. 16d) of the SH is not one-to-one correspondence with the BKS SIC. It implies that the predictabilities of the SH anomalies in October and February, under the BKS SIC reduction since 2005, may be limited if the ensemble-mean prediction is used.

Summary
In this study, guided by the SPB model, whether there exist a one-to-two correspondence relationship between the reduced BKS SIC in specific months and the lagging Ural atmospheric circulation is explored. Based on the monthly observational SIC data and two reanalyses during 1979/1980 − 2020/2021, two typical examples that generally agree with the SPB model are found: (1) The joint PDF of the mean SIC in the BKS in September (January) and the ζ850 index over the Ural region in October (February) shows multiple peaks and looks like a two-tine pitchfork (Fig. 5e,f;Supplementary Fig. S9a,b), resembling the bifurcation diagram of the SPB model (Fig. 3). It suggests that over the Ural region in October and February, there exist three preferred circulation patterns (or circulation regimes), namely, the neutral, positive anomaly and negative anomaly states, which are characterized by the flat zonal westerlies, wavy westerlies with a Ural trough, and wavy westerlies with a Ural ridge in the troposphere, respectively (Figs. 7, 8). The characteristics of three preferred Ural circulation patterns are similar to three flow equilibria (namely, the neutral-type, trough-type and ridge-type equilibria) in the SPB model (Fig. 3).
(2) When the mean SIC is below a critical value in September (January), referred to as the supercritical case, the positive and negative anomaly events of Ural atmospheric circulation occur nearly at the same frequencies in October (February) ( Table 1), indicating that they should occur with equal probability. It agrees with the Fig. 16 a, b As in Fig. 5e, f, respectively, but the ordinate is replaced by the SH intensity. c, d As in Fig. 5e, f, respectively, but the ordinate is replaced by the SH extent. Horizontal blue lines denote the mean of a-b the SH intensity and c-d the SH extent. Bandwidths of the kernel smoothing window for the SH intensity are 1.0 hPa in October and 1.5 hPa in February, for the SH extent and the mean SIC are always 20 and 2%, respectively SPB model who implies that the trough-type and ridgetype equilibria are equally likely to occur.
Hence, we refer to the relationship between the BKS SIC in September (January) and the Ural atmospheric circulation in October (February) as a pitchfork-like relationship. The mean θ925 over the BKS, as a substitution variable of the mean SIC, has demonstrated the robustness of the pitchfork-like relationship (Fig. 6, Supplementary Fig.  S9c, d; Table 1).The similarities between the pitchfork-like relationship and the SPB model have been summarized in Table 2, and their differences have been discussed in Sect. 4.1. The pitchfork-like relationship suggests that: (1) in the subcritical case (high SIC), the SIC change in the BKS may have little influence on the Ural atmospheric circulation; (2) in the supercritical case (low SIC), both the positive and negative anomaly regimes, as well as most positive and negative anomaly events, of Ural atmospheric circulation may be related to the BKS SIC, which is a one-to-two correspondence (note that it does not mean that except for the BKS SIC, other factors are not related to the Ural circulation anomalies). We have illustrated that the pitchfork-like relationship is very sensitive to noise. However, whether the pitchfork-like relationship is a cause-and-effect relationship is still unclear and need be investigated in the future.
Inspired by the SPB model, we have provided a possible explanation for the mechanisms of the existence of the pitchfork-like relationship. After the late 1970s, the BKS SIC in September, October, January and February can be considered a quasi-stationary diabatic heating factor. The thermodynamic effect of the SIC reduction in September (January) may contribute to similar changes in the background atmospheric conditions in September and October (January and February), especially the increased both positive and negative shear vorticity of background westerly wind over the Ural region, which might indicate the decreased stability of the flat zonal westerlies in the subcritical case. When the gradually reduced BKS SIC in September (January) passes a critical threshold, both the barotropic and baroclinic conversion of energy from the climatological-mean flow to the Ural circulation anomalies in October (February) exhibit a regime transition from one regime to two regimes, leading to the transition of preferred Ural circulation patterns from the flat zonal westerlies to wavy westerlies with a Ural trough and wavy westerlies with a Ural ridge. After the regime transition (i.e., in the supercritical case), under the same BKS SIC, both positive and negative anomalies are preferred Ural circulation patterns that extract kinetic energy and available potential energy from the climatological-mean flow to maintain themselves; hence, there exists a one-to-two correspondence relationship between the BKS SIC and the Ural circulation anomalies. Furthermore, the regime transition of Ural atmospheric circulation mainly influences the intensity of Siberian High in October and the extent of Siberian High in February. Our results suggest that the BKS sea ice loss since 2005 leads to an increased incidence of both positive and negative anomalies of the Ural atmospheric circulation, as well as the Siberian High, in October and February. It implies that the predictabilities of Eurasian climate anomalies in October and February, which are induced by the Ural atmospheric circulation anomalies and the Siberian High anomalies, may be low if the ensemble-mean prediction is used.
Overall, guided by the SPB model and from a one-totwo correspondence perspective, our work provides a new understanding on the possible influences of recent BKS sea ice loss on monthly anomalies of the Ural atmospheric circulation, and predictabilities of Eurasian climate anomalies. The two typical examples of the pitchfork-like relationship in September/October and January/February also provide evidences for the existence of multiple flow equilibria in the atmosphere (Charney et al. 1979(Charney et al. , 1980 over the Ural region.
It should be pointed out that the observed atmosphere is more complex than the SPB model. The SPB model just provides a heuristic guidance for detecting the pitchforklike relationship between a diabatic heating factor and the Ural atmospheric circulation in observations, and exploring the possible reasons that why this relationship exists. It does not directly tell us the properties of the pitchfork-like relationship detected, such as whether it is a cause-andeffect relationship and through what physical processes it is established. These problems need the concrete analysis by using multiple approaches.
In addition, although we focus our analysis on the monthly time scale and we suggest that the pitchfork-like relationship may be a phenomenon on the time scale of one month, it does not mean that the nonlinear atmospheric processes associated with transient waves (e.g., nonlinear wave-mean flow interaction and eddy-mean energy transfer) contribute little to monthly anomalies of the Ural atmospheric circulation. Associated with the regime transition of monthly Ural atmospheric circulation, what have changed for the nonlinear atmospheric processes on the time scale shorter than one month, such as blocking highs, cut-off lows, the Eurasian teleconnection pattern and eddy-mean energy transfer, and why the changes happen, deserve careful investigation in the future. Moreover, as the monthly Ural atmospheric circulation can be influenced by multiple factors (e.g., North Atlantic sea surface temperature, Eurasian snow cover) besides the BKS SIC, what roles have other factors played in the existence of the pitchfork-like relationship also need a study.

Conflict of interest
The authors have no relevant financial or non-financial interests to disclose.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http:// creat iveco mmons. org/ licen ses/ by/4. 0/.