The role of teleconnection patterns in the variability and trends of growing season indices across Europe

Teleconnection patterns affect the weather and climate on both interannual and decadal timescales which in turn affects various socio‐economic sectors such as agriculture. We use three climate indices based on E‐OBS data from the INDECIS dataset (growing season onset [ogs10], growing season rainfall [gsr] and growing season temperature [ta_o]) to assess the interannual variability and trends over 1950–2017 associated with four teleconnection patterns (North Atlantic Oscillation [NAO], East Atlantic pattern [EA], Scandinavian pattern [SCA] and East Atlantic/West Russia pattern [EAWR]) using linear regression to extract the signal of each teleconnection pattern and their contribution to interannual variability. Trends towards an earlier growing season onset are found across most of Europe in low‐lying regions. The NAO dominates interannual variability in northwest Europe where an NAO index of 1 is associated with earlier ogs10 of about 10 days and the EA dominates the continent with a trend towards the positive EA phase driving an earlier growing season onset of 1.1–1.7 days·decade−1 in five regions (Great Britain and Ireland, France, Italy, Poland and North Germany, Hungary, Balkans). The EA and SCA gsr signals have north/south splits of orientation: positive EA is linked to increased gsr in northern regions and reduced gsr in southern Europe, and vice versa for SCA. Correlations between gsr interannual variability and the teleconnection contributions are strongest in the Mediterranean regions and south Scandinavia with maxima of 0.41 and 0.46, respectively. Decreasing ta_o trends in Romania are explained by poor data coverage causing problems with the EOBS gridding algorithm when new stations are incorporated from 1961. The net effect is that Romanian ta_o is about 1.5°C cooler than expected compared to trends from surrounding countries. Improved spatial and temporal data coverage will benefit the EOBS dataset and prevent such erroneous trends.


| INTRODUCTION
Agriculture is strongly affected by weather and climate with up to 80% of variability in crop production linked to weather conditions (Hoogenboom, 2000).For example, crop yields are sensitive to teleconnection patterns and provide advance notice to farmers regarding adaptation and mitigation for potential losses during the growing season (Gonsamu and Chen, 2015;Ceglar et al., 2017).The agriculture sector must therefore adapt to climate change since the type and yield of crops will be affected by the increased risk of heatwaves and droughts (Olesen et al., 2011;Zhu and Troy, 2018) and hydrological changes such as floods (Falloon and Betts, 2010).Increasing temperatures will also cause the northward movement of crop zones (Beck et al., 2018;Ceglar et al., 2019) with longer growing seasons in northern Europe due to warmer winters and the increased likelihood of droughts in southern Europe (Bindi and Olesen, 2011).Climate indices relevant to the agriculture sector are therefore likely to be very useful in understanding how the growing season has been affected by climate change over the past seven decades and how it is affected by teleconnection patterns.An improved understanding of the link between teleconnection patterns and agriculture indices will also be useful seasonal forecasting in the agriculture sector (Soares and Dessai, 2015;Nobre et al., 2019) which can make use of seasonal forecasting of teleconnection patterns (Lled o et al., 2020).Here we exploit a new dataset to evaluate changes and variability in growing season indices over Europe.
Climate indices are standardized metrics for assessing climate change based on various meteorological variables, often using daily data (Frich et al., 2002;Zhang et al., 2011).Some climate indices are used to detect extremes in temperature and precipitation, such as the 27 core indices developed by the Expert Team on Climate Change Detection and Indices (ETCCDI).However, many more indices are used for socio-economic sectors such as agriculture, tourism, health and energy.The INtegrated approach for the Development across Europe of user oriented Climate Indicators for GFCS high-priority Sectors (INDECIS) project compiled 125 sector-relevant indices for Europe (Dominguez-Castro et al., 2020;Peña-Angulo et al., 2020), including those from ETCCDI.
The trends in these indices can show where and how rapidly climate change is occurring.Peña-Angulo et al. (2020) have catalogued the climatologies and trends of the INDECIS climate indices in the European Climatology and Trend Atlas of Climate Indices (ECTACI).Linking the variability in climate indices to teleconnection patterns can also help us to understand how various aspects of the climate are changing in response to changes in atmospheric circulation.There have been several studies linking indices to teleconnection patterns in parts of Europe, such as temperature and precipitation indices in Serbia (Kneževi c et al., 2014;Arsenovi c et al., 2015), extreme precipitation in Finland (Irannezhad et al., 2017) and seasonal precipitation in northern Scandinavia (Marshall et al., 2020).These studies have tended to focus on particular countries or regions rather than all of Europe, but the release of the INDECIS dataset (Dominguez-Castro et al., 2020) has facilitated the study of climate indices on the pan-European scale and the role of teleconnection patterns in their trends and variability.
We have focused on three growing season metrics relevant to the agriculture sector in this study: growing season onset, growing season precipitation and mean growing season temperature.Previous studies have linked some aspects of the growing season to teleconnection patterns in various parts of Europe but have all focused on specific countries or regions.For example, Ylhäisi et al. (2010) and Irannezhad and Kløve (2015) looked at growing season parameters across Finland, Tomcyzk et al. (2019) and Tomczyk and Szyga-Pluta (2019) focused on Poland and Popotova et al. (2015) examined a specific region of the Czech Republic which produces a large amount of farmed vegetables.There have also been several studies applying growing season precipitation and mean growing season temperature to viticulture (wine production) in parts of Europe: the Iberian Peninsula (Blanco-Ward et al., 2007;Ramos et al., 2008;Santos et al., 2012;Moral et al., 2015;Blanco-Ward et al., 2017;Blanco-Ward et al., 2019), France (Neethling et al., 2012), England (Nesbitt et al., 2016) and Hungary (Kovacs et al., 2017).Only Santos et al. (2012), Irannezhad and Kløve (2015) and Tomcyzk et al. (2019) have linked the growing season indices to teleconnection patterns, although Cornes et al. (2019)  The aim of this study is to assess the role of four teleconnection patterns relevant to European climate in the interannual variability of three growing season indices and their contribution to trends in the indices since 1950 compared to the residual climate change trend.This is achieved using a linear regression method outlined in section 2.4 which can be used to extract the signal of each teleconnection from the growing season indices time series.The overall trends are split into the trends associated with each teleconnection (circulation changes) and the residual trends associated with climate change in various parts of Europe.Some of the trends are also compared to the underlying station data to assess their reliability in section 4.

| INDECIS dataset
The ERA4CS INDECIS project produced a dataset of 125 climate indices across eight broad categories (temperature, precipitation, bioclimate, wind, aridity/continentality, snow, cloud/radiation, drought) spanning Europe from 1950 to 2017 (Dominguez-Castro et al., 2020).The 27 core indices from the ETCCDI were included in this dataset (Frich et al., 2002).The data were constructed from the ECA&D E-OBS gridded dataset (Cornes et al., 2018) and the ERA5 reanalysis (Hersbach et al., 2020).The E-OBS data are taken from quality controlled meteorological records and ERA5 was used to supplement the observations.The indices were produced at a horizontal resolution of 0.25 to match the resolution of ERA5 and have been made available at monthly, seasonal and annual time scales (http://indecis.eu/indices.php).The domain of this dataset covers Europe including Iceland but not the European part of Turkey and western Russia due to a lack of data availability.Any blank regions in eastern and southern Europe shown in figures in this paper are a result of missing data.

| Growing season indices
INDECIS has various temperature-and precipitationbased indices related to the agriculture sector (Dominguez-Castro et al., 2020).For example, growing season onset (ogs6, ogs10), growing season length (gsl), growing season rainfall (gsr), non-growing season rainfall (ngsr), growing degree days (gdd) and mean growing season temperature (t_ao, t_ms) are all included in the dataset.In this study, we focus on three of these indices: growing season onset, growing season rainfall and mean growing season precipitation.It is important to note that these indices are intended to represent a range of crops and growing conditions since they are calculated from gridded datasets so cannot replicate exact conditions at field scale or incorporate local microscale features.
INDECIS provides two definitions of thermal growing season onset: one using the first 6 day span with each daily mean temperature above 5 C (ogs6) (Klein Tank et al., 2009;Zhang et al., 2011;Popotova et al., 2015) and another using the first 10 day span with each daily mean temperature above 5 C (ogs10).We have used ogs10 because synoptic scale temperature events which influence growing season onset typically last around 10 days (Cornes et al., 2019).This removes subsynoptic-scale variability which is problematic when using gridded datasets as 0 C may not cause frost damage across an entire grid cell.Using the longer span of days above 5 C also avoids any unusually early growing season onsets followed by late frosts to which the early growing season will become increasingly susceptible to as the onset gets earlier with climate change (Wypych et al., 2017;Liu et al., 2018;Ceglar et al., 2019).The units for ogs10 provided by INDECIS are simply the day of year (Dominguez-Castro et al., 2020) which is also referred to as Julian Days (Ati et al., 2002).
INDECIS provides the total growing season rainfall (gsr) as the cumulative total precipitation during the warm season between April and October (Ramos et al., 2008) using E-OBS data and beginning in 1950(Dominguez-Castro et al., 2020), which we use in this study.Some previous studies using European gsr that have focused on wine production in the Iberian Peninsula (e.g., Blanco-Ward et al., 2007;Ramos et al., 2008;Santos et al., 2012;Moral et al., 2015;Blanco-Ward et al., 2019) or Hungary (Kovacs et al., 2017) have used slightly different timespans over which to calculate gsr (April-October, April-September or May-September).
There are two versions of the mean growing season temperature in INDECIS: the April-to-October mean (ta_o) and the May-to-September mean (tm_s).Here, we focus on ta_o to maintain a consistent timespan with gsr.So far, there have been few studies using the mean growing season temperature for Europe and, like gsr, they are linked to viticulture in France (Neethling et al., 2012), England (Nesbitt et al., 2016) and Portugal (Blanco-Ward et al., 2017).

| Teleconnections
The interannual variability and trends in some of the indices are linked to four Northern Hemisphere atmospheric circulation patterns relevant to European weather and climate: the North Atlantic Oscillation (NAO), the Scandinavian pattern (SCA), the East Atlantic Pattern (EA) and the East Atlantic/Western Russia pattern (EAWR).The monthly time series of each teleconnection index were obtained from the National Oceanic and Atmospheric Administration's (NOAA) Climate Prediction Center (https://www.cpc.ncep.noaa.gov/data/teledoc/telecontents.shtml).These teleconnection patterns are identified using Rotated Principal Component Analysis (RCPA; Barnston and Livezey, 1987) applied to the monthly mean standardized 500 mb height anomalies between 20 N and 90 N. The teleconnection indices are calculated from the rotated leading modes of variability with a least squares solution.Composite anomalies of ERA5 500 hPa and 10 m winds are shown in Figures 1  and 2 over the 1950-2017 over which the indices discussed in section 2.2 are calculated.
The NAO is generally understood as the difference in sea level pressure (SLP) between the Azores and Iceland.In winter, the positive phase (NAO+) occurs when the pressure gradient is increased due to strengthening of the Icelandic Low and Azores High, and the negative phase (NAO−) occurs when the pressure gradient reduces due to a weakening of these centres of action (Hurrell et al., 2001).During NAO+ westerly winds are stronger (Figure 1a) with increased storm track activity and precipitation across northwest Europe, and during NAO− winters the storm track has a zonal orientation which causes increased precipitation over southern Europe.The NAO is weaker in summer (Figure 2a winter, but does have notable impacts on European precipitation, temperature and cloudiness (Folland et al., 2009).
The SCA consists of three centres of action with a strong positive height anomaly over Scandinavia and two weaker negative anomalies over eastern Russia and western Europe (Barnston and Livezey, 1987;Crasemann et al., 2018).The positive phase (SCA+) is associated with reduced precipitation over Scandinavia and increased precipitation over the Northeast Atlantic due to corresponding anticyclonic and cyclonic anomalies (Figures 1e,f and 2e,f) (Bueh and Nakamura, 2007;Marshall et al., 2020).In winter, negative temperature anomalies also occur over northern Europe during SCA+ (Sui et al., 2020) with warmer temperatures over western Europe during the negative phase (SCA−) (Vihma et al., 2020).
The EA is defined by a dipole across the North Atlantic with a centre of action at 55 N, 20-35 W (Figures 1c,d and 2c,d) and a band spanning the subtropics across the Mediterranean Sea (Barnston and Livezey, 1987;Bastos et al., 2016) although other studies have used SLP to identify an EA with a different pattern (Comas-Bru and McDermott, 2013;Comas-Bru and Hern andez, 2018;Mellado-Cano et al., 2019).We use the NOAA EA index based on the Barnston and Livezey (1987) method for consistency with the other teleconnection indices used since they are orthogonal.This EA index is associated with enhanced westerly winds from the Atlantic in the positive phase (EA+) (Figures 1c and 2c) and northerly wind anomalies over western Europe in its negative phase (EA−) (Figures 1d and 2d).Various studies have highlighted the relationship between the EA and NAO (Moore et al., 2011;Moore and Renfrew, 2012;Comas-Bru and McDermott, 2013;Comas-Bru et al., 2016) as they both interact to influence the latitude and strength of the eddy-driven jet (Woolings et al., 2010).Mellado-Cano et al. (2019) showed that EA+ shifts the jet stream north and EA− shifts it south.Hall and Hanna (2018) found positive correlations with United Kingdom summer and winter rainfall, and negative correlations with summer temperature.
The EAWR has four centres of action with two of the same sign over the North Sea (Figures 1g,h and 2g,  h) and Mongolia and two of the opposite sign over the central North Atlantic and western Russia north of the Caspian Sea (Barnston and Livezey, 1987;Lim, 2015).In the positive phase (EAWR+) the North Sea and Mongolian height anomalies are positive and the North Atlantic and western Russia height anomalies are negative.This corresponds to a meridional circulation over eastern Europe with southerly wind anomalies in winter (Figure 1g,h) and northerly wind anomalies between April and October (Figure 2g,h).Ionita (2014) showed that EAWR is strongest in winter with EAWR+ associated with dry conditions in southern Europe and stronger rainfall over Scandinavia, and increased temperatures across northern and central Europe but cooler temperatures over southern Europe.Lim (2015) also showed that the negative phase (EAWR−) is linked to cooler temperatures over western Europe and warmer temperatures over parts of eastern Europe.

| Regression method
To investigate the relationship between the growing season indices and four teleconnection patterns, the teleconnection signal (S) is removed from the growing season index time series (X) using the method described by Bhend and von Storch (2008) and Iles and Hegerl (2017).The detrended growing season index time series ( b X, where b Á denotes a detrending) is regressed on the detrended teleconnection index time series (b τ) at each grid point and the resulting slope (d b X=db τ) gives S. The teleconnection's contribution to the interannual variability in X is simply the product of the signal and the teleconnection index time series, defined: This produces a time series of C which can be compared to the time series of X.Note that we have not detrended τ when calculating C since our aim is to quantify how changes in teleconnection patterns influence both the variability and trend in each of the indices.Trends in teleconnection pattern could reflect multidecadal internal variability but potentially also a response of the atmospheric circulation to greenhouse gas and aerosol forcing.Therefore, it is difficult to disentangle climate change signals from trends in teleconnection patterns.The trend in the teleconnection contribution is represented by dC=dt, where t is time.The residual variability (R) represents the variability when the variability associated with each teleconnection pattern is removed and is simply calculated as the difference between the growing season index time series and the teleconnection's contribution to the variability, The residual trend, dR=dt, shows the remaining changes potentially associated with climate change after circulation changes associated with trends in the four teleconnection patterns have been removed.Statistical significance is calculated for each trend using a Wald test (Wald, 1943) with a t-distribution of the test statistic.Correlations between each C and X are also calculated with a Student's t test showing that all correlation coefficients greater than 0.239 are statistically significant for a sample size of 68.The statistically significant correlation coefficients presented are all positive even if their respective teleconnection signals are negative.This is because the teleconnection patterns are orthogonal (Barnston and Livezey, 1987) so there are no strong peaks in C of opposing sign to X.

| Growing season onset
Growing season onset (ogs10) is earliest in southern and western Europe, and later towards the north and east with mountainous regions also showing late onset (Figure 3a).Some of these regions are less relevant to agriculture as they are mainly occupied by forests (Rega et al., 2020) and do not have a suitable climate for growing crops (Metzger et al., 2005;Beck et al., 2018).The most variable regions are Ireland, south England, north France through Belgium and Holland into north Germany where standard deviation exceeds 30 days.Very low variability is seen in the rest of Europe and especially in Portugal and south Spain where the growing season is likely to last all year (Figure 3b).There is a general trend towards earlier onset over Europe which is strongest around the eastern end of the English Channel (Figure 3c) where ogs10 has moved earlier by at least 4 daysÁdecade −1 .Arnell and Freeman (2021) identified similar values for south England with the growing season onset advancing by about 30 days over the 21st century projections relative to the 1981-2010 average.Most of eastern Europe, which has late onset and low variability, has a weaker but statistically significant trend towards earlier onset of around 2 daysÁdecade −1 .
We focus on the JFM (January-February-March) NAO and EA teleconnections as the mean ogs10 for the main agricultural regions of Europe is generally within the first 90 days of the year (Figure 3a) and these two teleconnection patterns have the strongest signals with widespread statistical significance (Figure 4a,b).Positive phases in the JFM NAO and EA generally result in an earlier growing season (Figure 4) due to westerly wind anomalies advecting warmer air towards western Europe during NAO+ and southerly wind anomalies across the Mediterranean in EA+ (Figure 1a,b).The EA signal is strong across much of continental Europe with a statistically significant negative signal (earlier onset) over France, Germany, Italy and southeast Europe (Figure 4d).The NAO signal is strongest over southern England and northern France (Figure 4a).Over continental Europe, the NAO signal is much weaker than the EA signal although there is a region of statistically significant NAO signal over Hungary.In general, the influence of the NAO on ogs10 decreases further away from the Atlantic coast (Scheifinger et al., 2002;Wypych et al., 2017) as confirmed by the weaker signal (Figure 4a) and correlations (Figure 5b-g) across the continent.
The SCA has also been included as the signal is of opposite sign to the NAO and EA (Figure 4c) and has some statistically significant correlations with ogs10 interannual variability (Figure 5b,e).Easterly wind anomalies associated with the anticyclonic anomaly centred over northwest Russia during SCA+ (Figure 1e) advect cold air across Europe thus delaying the growing season onset.The SCA signal does not have widespread statistical significance but some isolated regions of significance are were selected (Figure 5a): Great Britain and Ireland (GBI); France; Italy; Poland and North Germany; Hungary; the Balkans.Hungary has been kept separate from the Balkans as it has strong and statistically significant signals from both the NAO and EA (Figure 4).The JFM teleconnection signals were averaged over these regions and compared to the area-averaged anomalies from the ogs10 climatology (Figure 5b-g).Years with particularly late or early growing season onsets can therefore be linked to one or more teleconnection component of interannual variability.Some of the interannual variability and anomalies in ogs10 (ogs10 0 ) can therefore be visibly linked to one of the teleconnection components (Figure 5b-g).
For example, the GBI NAO component closely follows the ogs10 variability in the late 1950s and from 2010 to 2013 (Figure 5b).The trend towards an earlier onset has been accompanied by a trend towards a positive NAO index with the extreme and persistent NAO− in 2010 (Cattiaux et al., 2010) causing a later onset (Figure 5h), although the NAO component exceeds ogs10 0 magnitude by about 5 days (Figure 5b).The SCA component has a statistically significant, albeit weak, correlation with ogs10 interannual variability over GBI.However, the contribution to ogs10 0 from the SCA component is generally weaker than the NAO component and any notable SCA contributions to ogs10 0 occur infrequently (e.g., 1989 and 1990).
Interannual variability in the EA components of the five continental regions corresponds well with the ogs10 variability, particularly for Italy and the Balkans (Figure 5d,g) which have the strongest correlations.The trend towards a positive EA index (Figure 5g) has also contributed to the trend towards an earlier onset.Many ogs10 0 are clearly linked to the EA component such as the negative ogs10 0 over Hungary in 1977 which is matched almost exactly by the EA component (Figure 5f) and the positive ogs10 0 in 2005 and 2006 over the Balkans (Figure 5g).However, the EA component often cannot explain the entire ogs10 0 .The Poland and north Germany region provides an extreme example of this in 2002 and 2007 with ogs10 0 of nearly −60 days compared to EA components of about −10 days in each year (Figure 5e).Both of these years had high EA indices (>1, Figure 5h) which are associated with ogs10 0 of at least 9 days in parts of this region (Figure 4d) but since ogs10 is normally around mid-March (Figure 1a) there may have been a role for weather patterns in December of the preceding years for these early onset events.
Five of the six regions in Figure 5a have statistically significant trends in ogs10, apart from the Balkans, with the NAO and EA dC=dt all negative or zero (Figure 3b and Table 1).The dC=dt maps are not shown for brevity since they essentially reflect the maps of signal (Figure 1a-c) multiplied by the trend in teleconnection but area-averages of dC=dt are presented in Table 1.The NAO and EA dC=dt therefore have the same spatial patterns as the signals (positive trends in NAO and EA) and SCA dC=dt has the opposite spatial pattern to its signal (negative trend in SCA).All six regions have statistically significant NAO and EA trends with most of the overall GBI trend coming from the NAO component and a negligible EA component.
Circulation changes therefore explain most of the trends towards an earlier ogs10 across Europe, but in the Balkans, circulation changes are opposed by dR=dt which implies a trend towards a later ogs10.This is particularly notable when removing the EA signal since dR=dt is positive (later onset) across most of the region and exceeds 3 daysÁdecade −1 in some places (not shown).The strongest overall trend is across the Poland and north Germany region (−3.0 daysÁdecade −1 ) which also has the strongest trend from the EA component (−1.7 daysÁdecade −1 ; Table 1).

| Growing season rainfall
The annual mean April-to-October gsr (Figure 6a) has a similar spatial pattern to the April-to-September gsr (Santos et al., 2012) with the highest rainfall totals in mountainous regions and the west coasts of the Iberian Peninsula, Scotland and Norway, and the driest regions (below 200 mm) in the south of the Iberian Peninsula, Italy and Greece.The trends since 1950 are mainly towards a wetter growing season of up to 30 mmÁdecade −1 in Scandinavia, greater than 40 mmÁdecade −1 in the Alps and less than 2 mmÁyear −1 across much of eastern Europe, and are statistically significant over much of Scandinavia, parts of Poland and the Balkans, the Eastern Alps, southern Portugal, Wales and West Scotland (Figure 6b).The coefficient of variation shows that although regions with the most variability in gsr have the wettest growing seasons, their standard deviation is only 15-20% of the mean gsr and the drier areas with low variability, such as southern Spain, have standard deviations around 50% of the mean gsr (Figure 6c).The coefficient of variation is shown in Figure 6c instead of standard deviation because precipitation data are not normally distributed so it is more informative to use the coefficient of variability to show the absolute variability relative to the mean.
The April-to-October average of the NAO, EA and SCA indices have statistically significant gsr signals in various parts of Europe (Figure 7).A more positive NAO is associated with less precipitation across continental Europe, GBI and Scandinavia but more precipitation along the northern fringes (e.g., Norwegian coast, northwest Scotland) and southern Europe (Folland et al., 2009;Bladé et al., 2012).The EA signal has a north/south split with statistically significant negative signal around the Mediterranean Sea and in eastern Europe around Moldova, with positive across parts of France and Germany, Denmark, GBI and Scandinavia.The SCA signal also has a north/south split of opposite sign to the EA signal: southern Europe has positive signal and some isolated regions of northern Europe have negative signal, which reflects the 500 hPa geopotential height anomalies with opposite cyclonic and anticyclonic anomalies between the extreme SCA phases (Figure 2e,f).The EAWR signal is not shown as it has only small, isolated regions of statistical significance in agricultural regions and its main region of statistical significance is in Finland and northern Sweden, where EAWR+ leads to a drier growing season (Irannezhad et al., 2016).
Six regions which have strong and statistically significant signals for at least one of the teleconnections shown in Figure 7 are GBI, south Scandinavia, France, the Iberian Peninsula, Poland and Italy and the Adriatic (ItAd) (Figure 8a).The ItAd region was chosen instead of separate regions for Italy and the Balkans (Figure 5a) because, for the Balkans, the SCA signal is only strong and statistically significant in the west along the coast of the Adriatic Sea (Figure 7g) and the gsr signals of all three teleconnections are generally of the same sign for the western Balkans and Italy.The interannual variability of the NAO, EA and SCA components are compared to the yearly gsr anomalies (gsr 0 ) in Figure 8b-g.
Few of the gsr 0 peaks can be explained by the teleconnection components, particularly when gsr 0 > 100 mm, and the statistically significant dC=dt correlation coefficients are all below 0.5 with south Scandinavian gsr correlation with EA dC=dt (0.46) the strongest (Figure 8).South Scandinavia is also the only region to have statistically significant correlations between gsr and three of the dC=dt.The strongest example of any of the teleconnection components matching gsr 0 is in 1960 when SCA > 1 (Figure 8h): the SCA component almost exactly matches the Iberian Peninsula gsr 0 (Figure 8e) and exceeds the ItAd gsr 0 (Figure 8g).Under such conditions gsr increases across these two regions, particularly in the northwest of each region (>100 mm; Figure 7e).This is caused by anomalous convergence of the surface winds across southern Europe (Figure 2e).Interannual variability of the SCA component corresponds fairly well with gsr 0 around 1990 and 2000, and it combines with the EA component in the 1970s for south Scandinavia and ItAd when the teleconnection indices are consistently of opposing sign (Figure 8h).The overall gsr trends are generally towards a wetter growing season across Europe (Figure 6b) although the Iberian Peninsula is an outlier with a weak area-averaged drying trend (−1.5 mmÁdecade −1 ; Table 2).South Scandinavia has the strongest increasing overall trend at 9.8 mmÁdecade −1 and is the only statistically significant overall gsr trend which is almost exactly matched by the trend in the EA component.However, the strong dR=dt (−7.0 mmÁdecade −1 ) offsets the smaller contributions from the other three components There is a clear north/ south split in each dC=dt with the EA (same spatial pattern as its signal) and SCA (opposite spatial pattern to its signal) trends contributing to increased gsr over some northern regions and reduced gsr further south.
Changes to gsr associated with changes to the atmospheric circulation have the opposite north/south split to the residual trends.The dR=dt are negative in the northernmost regions (−9.8 mmÁdecade −1 for GBI) but positive further south and particularly strong around the Mediterranean (Table 2).The ItAd region has the strongest dR=dt of 17.5 mmÁdecade −1 which mostly occurs in the northeastern Alpine areas and on the east Adriatic coast around Bosnia (Figure 6b).The opposing dC=dt and dR=dt lead to some cancellation and relatively weak overall gsr trends which are not statistically significant at the 95% level.

| Growing season temperature
Across much of southern Europe ta_o exceeds 20 C with milder temperatures further north and temperatures below 6 C at high latitudes and altitudes which are unsuitable for agriculture (Figure 9a).All of Europe shows a statistically significant increasing trend in ta_o apart from a striking negative trend in Romania which crosses the border into northern Bulgaria.(Figure 9b).There are also some small areas in Greece, north Macedonia, Norway, Sweden and northwest Spain which have weak decreasing trends (<0.1 CÁdecade −1 ).Romania also has a remarkably high standard deviation (>1 C for most of the country) compared to the rest of Europe which generally has standard deviations below 1 C (Figure 9c) and this is investigated in section 4.
The April-to-October averages of all four teleconnections have statistically significant ta_o signals in various parts of Europe (Figure 10a-d).The EA signal is positive across continental Europe and some of Scandinavia and generally has the strongest correlations with area-averaged ta_o in the six regions chosen in Figure 11.The Balkans region has the strongest statistical relationship with the EA of about 1 C increase in ta_o per unit of the EA index.Kneževi c et al. (2014) also found warmer April-to-September temperatures in Serbia associated with EA+.
Statistical significance in the EAWR signal is confined to eastern Europe (Figure 10d) close to the centre of action over Siberia and affected by northerly wind anomalies bringing cooler air (Figure 2g,h).The region of statistical significance spans from the high latitudes in Lapland to the Mediterranean on the south coasts of Italy and Greece, but north of the Mediterranean regions the strongest EAWR signal is found east of 20 E (Putnikovi c et al., 2018).In Romania, the EAWR signal has strong negative relationships which exceed −1 C per unit of the EAWR index.Arsenovi c et al. (2015) found that EAWR+ is linked to cooler summer and autumn temperatures at stations across Serbia, and Irannezhad et al. (2016) found negative correlations between the EAWR and growing season temperature at Finnish stations.
There is a clear north/south split in the SCA signal with positive and statistically significant relationships across Scandinavia and Britain, and negative relationships around the Mediterranean consistent with the SCA geopotential height anomalies (Figure 2e,f).The NAO has the weakest ta_o signal with statistical significance confined to a small region of weak, positive relationship in the northernmost region of Scandinavia and negative relationships in parts of southern Europe, mainly the Balkans.The weak role of the NAO in the variability of ta_o is likely a direct consequence of the summer NAO having a much weaker impact on European weather than its winter counterpart (Folland et al., 2009) and the April-to-October NAO index's tendency to remain in NAO+ or weak NAO− (Figure 8h) which give weak geopotential height and surface wind anomalies (Figure 2a,b).
The EA component most closely matches the interannual variability of ta_o for the Iberian Peninsula (r = 0.65) and ItAd (r = 0.71) regions, particularly in the 1970s (Figure 11d,e).The periods where the Iberian Peninsula and ItAd EA components most closely match the ta_o variability are also periods of strong EA− or EA+ (Figure 8h).The trend in the April-to-October EA index represents a transition from the northerly wind anomalies in the 1970s causing cooler temperature to southerly wind anomalies advecting warmer air towards southern Europe (Figure 2c,d).During the period of persistent EA− in the 1970s, there also appears to be some small contributions from the SCA component towards Iberian Peninsula ta_o 0 since the persistent SCA+ conditions cause cooler temperatures over this region (Figure 10c).
The EAWR component has very weak variability and does not contribute to any ta_o 0 peaks in western Europe but it has a very strong impact in eastern Europe (Putnikovi c et al., 2018).For example, the negative Romanian ta_o 0 peaks in 1985 and 1997 are substantially explained by the EAWR component during EAWR+ (Figure 11g,h).Ukraine and Belarus has the strongest correlation with the EAWR component with the negative ta_o 0 peaks in 1965 and 1997 receiving substantial contributions from the EAWR component.However, when ta_o 0 becomes positive the EAWR component struggles to match it despite the trend towards EAWR− since 1995 (Figure 11h).
In five of the six regions presented in Figure 11a the area-averaged ta_o are increasing and statistically significant (Table 3) with Romania's decreasing trend being the outlier (Figure 9b).There is a north/south split regarding the dominant contribution to the overall trend between dR=dt and EA dC=dt, which is the only dC=dt to have the same spatial pattern as its signal (not shown).In GBI and south Scandinavia, dR=dt dominates the actual ta_o trend with EA dC=dt balanced or exceeded by negative SCA dC=dt.The EA dC=dt dominates the actual ta_o trend in Ukraine and Belarus, the ItAd region and the Iberian Peninsula where each EA dC=dt exceeds 0.1 CÁdecade −1 .Therefore, temperature changes associated with changes to the atmospheric circulation dominate the ta_o trends in southern and eastern Europe, but temperature changes in northern regions are linked to other factors not directly included in the teleconnection indices.
Romania has the strongest EA trend (0.22 CÁdecade −1 ) and its dR=dt is larger and of opposite sign which would imply that, when circulation changes are removed, the trend in Romanian ta_o is negative.This is extremely unlikely since the rest of Europe has increasing ta_o in line with observed increases in global temperatures and is caused by the large positive ta_o 0 in the 1950s which mask the recent increase in ta_o (Figure 11g).This is discussed further in section 4. The EA-linked changes in southern and eastern Europe are therefore associated with a trend towards EA+ since the 1970s (Figure 8h).These parts of Europe have therefore experienced more positive geopotential height anomalies and southerly surface wind anomalies (Figure 2c).This in turn results in advection of warmer air towards Europe and anticyclonic conditions with reduced clouds cover which leads to increased temperatures in the April-October period.

| USING STATION DATA TO ASSESS ROMANIAN TEMPERATURE TRENDS
In Romania there is a striking feature in the ta_o trends where the trend is of opposite sign to the rest of Europe and statistically significant (Figures 9 and 10).To determine the cause of this, we have used the ECA&D station data for Romania and the surrounding countries (Figure 12).The daily station data are available from the ECA&D website (www.ecad.eu) in non-blended and blended formats.The blended time series are constructed by merging the data from two stations that are within 12.5 km of each other and have a height difference of no more than 25 m, and the non-blended data are simply the underlying series at each station with no merging applied.In Figure 12a, blended stations were chosen for Ukraine and Moldova because there are multiple time series for some stations.However, the Romanian, Hungarian and Serbian stations all have one time series for each station so the non-blended data are identical to the blended data.The ECA&D also does not make all the daily station data publicly available such as the Bulgarian stations highlighted in Figure 12a.
The most notable difference between Romanian ta_o 0 and other parts of Europe is the large positive anomalies in the 1950s which gradually decrease to negative anomalies similar to the other regions (Figure 11).Only three Romanian ECA&D stations have data in the 1950s (from north to south: Baia Mare, Targu Jiu and Turnu Magurele) but none of them show a steep decreasing trend (Figure 12b).Instead, the addition of several stations in 1961, including three cooler high altitude stations (from north to south: Ceahlau Toaca, Miercurea Ciuc and Varfu Omul; Figure 12b), coincides with a sudden negative shift when area-averaged ta_o for the surrounding countries (excluding Ukraine) all increase by about 1 C (Figure 12c).There appears to be no impact on areaaveraged Romanian ta_o when the three stations with data in the 1950s terminate in 1993 (Figure 12b).
T A B L E 3 Area-averaged trends for all six regions shown in Figure 11a Actual  The post-1961 time series for area-averaged Romanian ta_o is therefore about 1.5 C cooler than expected compared to Hungary, Bulgaria, Serbia and Moldova (Figure 12c).Although Romanian ta_o is cooler than expected, its year-to-year variability is consistent with Ukraine's after 1961.However, Ukrainian ta_o is likely to be heavily dependent on stations further east and therefore of little relevance to the temperature in Romania.The negative shift between 1960 and 1961 has two possible explanations: 1.The three cooler stations (Figure 12b) have too much influence on neighbouring stations in the EOBS gridding algorithm.2. Data sparsity before 1961 caused difficulties in training the 3D monthly thin-plate spline since altitude is a covariate in the EOBS gridding algorithm (Cornes et al., 2018).
Had it not occurred, area-averaged Romanian ta_o would be about 1.5 C warmer and would have a weak increasing trend across 1950-2017 like all the surrounding countries (Figure 9b).Romanian dR=dt should also therefore be much smaller and possibly increasing like some other parts of Europe (Table 3).Reduced Romanian EA and EAWR dC=dt would also be expected without the 1960-1961 temperature shift since the signals for both teleconnections (Figure 10b,d) are stronger in Romania than in surrounding countries and this is passed onto both dC=dt (not shown) which are notably larger than the areas-averaged dC=dt in the other regions considered (Table 3).
A lack of data in mainland Greece and in north Macedonia also contributes to the weak increasing and decreasing ta_o trends with no statistical significance across the southern region of the Balkans peninsula (Figure 9b).The ECA&D has only two stations which provide daily average temperature with data for only one (Prilep) made publicly available.Greece has nine stations within the region where the ta_o trend data is plotted in Figure 9b, but only Corfu (an island off the west coast) and Larissa (east coast at about 39 N) have data as far back as 1955.Improved spatial and temporal data coverage in ECA&D in this part of Europe is therefore required to facilitate more realistic trends in ta_o.
Overall, the region of southeast Europe focused on Romania appears to display artificial changes in surface temperature since the 1950s which can be linked with limitations in station coverage and regional representivity.Once these artefacts are accounted for, the remaining changes across Europe display consistency with the expectation of a warming of climate and changing atmospheric teleconnection patterns.

| SUMMARY
We have assessed three climate indices relevant to the agriculture sector, specifically growing season onset (ogs10), growing season rainfall (gsr) and mean growing season temperature (ta_o).Their variability and trends are related to four teleconnection patterns: North Atlantic Oscillation (NAO), East Atlantic pattern (EA), Scandinavian pattern (SCA) and East Atlantic/West Russia pattern (EAWR).This is the first study to link these indices to teleconnection patterns from a pan-European gridded dataset as previous studies have looked at specific regions in Europe (Irannezhad and Kløve, 2015;Tomcyzk et al., 2019;Tomczyk and Szyga-Pluta, 2019) or at the ECA&D station data (Cornes et al., 2019).
A linear regression method (Bhend and von Storch, 2008;Iles and Hegerl, 2017)  The ECA&D station data show a negative shift in area-averaged ta_o when many new stations are added in 1961 which is not seen in neighbouring countries.This coincides with the introduction of many new stations in Romania including three cooler high altitude stations that are not representative of the wider region.The negative shift in Romanian ta_o may either be a result of the high altitude stations exerting too much influence in the EOBS algorithm or data sparsity causing problems with the spline used in the algorithm.
Based on the findings in (7), data rescue efforts may be required to improve the spatial and temporal distribution of stations in regions with erroneous temperature trends in EOBS (Figure 9).Coll et al. (2019) recovered temperature data from various stations in the former Yugosloavian countries for EOBS, but there remain substantial gaps in coverage for some parts of Europe.The lack of data in the Carpathian region (44-50 N, 17-27 E) in the 1950s restricts the CARPATCLIM (Climate of the Carpathian Region) dataset to beginning in 1961 (Antolovi c et al., 2013), when many stations are added in Romania (Figure 12b).A more accurate understanding of the trends and variability of temperature-based climate indices could be gained by improving the spatial and temporal coverage of the underlying data from which they are calculated.The various socio-economic sectors that use these indices will therefore be better informed as to how to adapt to climate change.
Although different types of crops depend on multiple factors, such as sunlight and soil type and precursor conditions in the months prior to planting and growth, here we have focused on the basic direct relationship between growing season onset, rainfall and temperature without determining the direct impacts.For example, a warmer growing season also poses the increased risk of heat stress on plants, particularly in regions such as northern Britain which have started from a lower base temperature (Arnell et al., 2021), although these aspects are beyond the scope of the present study.Examination of a more comprehensive set of variables including meteorological conditions prior to the growing season are clearly of importance in future work and could be investigated in future work using the INDECIS datasets.The findings in this study may also provide a basis for seasonal forecasting in the agriculture sector regarding crop yields and types to suit the expected weather conditions.
,b) and influences a smaller region which is further north than in F I G U R E 1 Anomalies of 1950-2017 ERA5 JFM 500 hPa geopotential height and 10 m winds for extreme positive (left column) and negative (right column) phases of the JFM NAO (a, b), EA (c, d), SCA (e, f ) and EAWR (g, h).Geopotential height units are in metres and wind speeds are in mÁs −1 [Colour figure can be viewed at wileyonlinelibrary.com]

F
I G U R E 2 Anomalies of 1950-2017 ERA5 April-to-October 500 hPa geopotential height and 10 m winds for extreme positive (left column) and negative (right column) phases of the April-to-October or growing season (GS) NAO (a, b), EA (c, d), SCA (e, f) and EAWR (g, h).Geopotential height units are in metres and wind speeds are in mÁs −1 [Colour figure can be viewed at wileyonlinelibrary.com]

F
I G U R E 3 (a) The 1950-2017 mean ogs10, (b) the trend in ogs10 over 1950-2017 where stippling indicates statistical significance at the 95% level using a Wald test with a t-distribution of the test statistic, and (c) the standard deviation of ogs10 across 1950-2017 [Colour figure can be viewed at wileyonlinelibrary.com]found in Ireland, Great Britain, Denmark and Romania.The EAWR signal is not included as it is weaker across Europe and has limited statistical significance.To investigate the role of the teleconnections in ogs10 interannual variability, six regions with statistically significant signals for one or more of the teleconnections F I G U R E 4 Relationships between the 1950-2017 JFM (January-February-March) NAO, EA and SCA teleconnection indices and ogs10.The top row shows the teleconnection signals and the bottom row shows the residual trends when the teleconnection patterns are removed: (a, d) show the NAO signal and residual trend, (b, e) show the same for the EA, (c, f ) show the same for the SCA.The teleconnection signals have the units daysÁindex −1 and each dR=dT has the units daysÁdecade −1 .Green contours indicate a reduction in ogs10 (earlier onset) and brown contours indicate an increase in ogs10 (later onset).Stippling indicates statistical significance at the 95% level using a Wald test with a t-distribution of the test statistic [Colour figure can be viewed at wileyonlinelibrary.com]F I G U R E 5 (a) Regions over which the JFM NAO and EA ogs10 signals from Figure 4 are averaged; (b-g) yearly ogs10 anomalies (black lines) with the yearly NAO (red lines), EA (blue lines) and SCA (gold lines) components; and (h) the NAO, EA and SCA indices averaged across JFM.In (b-g) negative values represent an earlier growing season onset and positive values represent a later growing season onset.The coloured numbers in the bottom left corner are the statistically significant Pearson correlation coefficients between the teleconnection components and the ogs10 time series where statistical significance is calculated from a t-distribution [Colour figure can be viewed at wileyonlinelibrary.com] The actual trend is the ogs10 trend from Figure 3a, the NAO, EA, SCA and EAWR columns represent the trends in the teleconnection component (dC=dt) of ogs10, and dR=dt (Figure 4d-f) is the remaining trend when changes associated with the teleconnections are removed.Negative trends indicate an earlier growing season onset and positive trends indicate a later growing season onset.Units are daysÁdecade −1 and numbers in bold are statistically significant at the 95% level.F I G U R E 6 (a) The 1950-2017 annual mean gsr (cumulative total over April-October), (b) the trend in gsr over 1950-2016 where stippling indicates statistical significance at the 95% level using a Wald test with a t-distribution of the test statistic, and (c) the coefficient of variation of gsr across 1950-2017 [Colour figure can be viewed at wileyonlinelibrary.com]

F
I G U R E 7 Relationships between the 1950-2017 April-to-October NAO, EA and SCA teleconnection indices and gsr (cumulative total from April to October).The top row shows the teleconnection signals, and the bottom row shows the residual trends when the teleconnection patterns are removed.(a, d) show the NAO signal and residual trend, (b, e) show the same for the EA, (c, f) show the same for the SCA.The teleconnection signals have the units mmÁindex −1 and each dR=dT has the units mmÁdecade −1 .Blue contours indicate an increase in gsr and red contours indicate a decrease in gsr.Stippling indicates statistical significance at the 95% level using a Wald test with a t-distribution of the test-statistic [Colour figure can be viewed at wileyonlinelibrary.com]F I G U R E 8 (a) Regions over which the GS NAO, EA and SCA gsr signals from Figure 7 are averaged; (b-g) yearly gsr anomalies (cumulative totals from April to October, black lines) with the yearly NAO (red lines), EA (blue lines) and SCA (gold line) components; and (h) the NAO, EA and SCA indices averaged across the GS.In (b-g) negative values represent a decrease in gsr and positive values represent an increase in gsr.The coloured numbers in the bottom right corner are the statistically significant Pearson correlation coefficients between the teleconnection components and the gsr time series where statistical significance is calculated from a t-distribution [Colour figure can be viewed at wileyonlinelibrary.com] The actual trend is the gsr trend from Figure6a, the NAO, EA, SCA and EAWR columns represent the trends in the teleconnection component (dC=dtÞ of gsr, and dR=dt (Figure7d-f) is the remaining trend when changes associated with the teleconnections are removed.Negative trends indicate a decrease in gsr and positive trends indicate an increase in gsr.Units are mmÁdecade −1 and numbers in bold are statistically significant at the 95% level.

F
I G U R E 9 (a) The 1950-2017 annual mean ta_o, (b) the trend in ta_o over 1950-2017 where stippling indicates statistical significance at the 95% level using a Wald test with a t-distribution of the test statistic, and (c) the standard deviation of ta_o across 1950-2017 [Colour figure can be viewed at wileyonlinelibrary.com]

F
I G U R E 1 0 Relationships between the 1950-2017 April-to-October NAO, EA, SCA and EAWR teleconnection indices and ta_o.The top row shows the teleconnection signals and the bottom shows the residual trends when the teleconnection patterns are removed.(a, e) Show the NAO signal and residual trend, (b, f) show the same for the EA, (c, g) show the same for the SCA, and (d, h) for EAWR.The teleconnection signals have the units CÁindex −1 and each dR=dT has the units CÁdecade −1 .Red contours indicate an increase in ta_o and blue contours indicate a decrease in ta_o.Stippling indicates statistical significance at the 95% level using a Wald test with a t-distribution of the test statistic [Colour figure can be viewed at wileyonlinelibrary.com]F I G U R E 1 1 (a) Regions over which the GS NAO, EA, SCA and EAWR ta_o signals from Figure 10 are averaged; (b-g) yearly gsr anomalies (black lines) with the yearly NAO (red lines), EA (blue lines), SCA (gold line) and EAWR (pink line) components; and (h) the EAWR index averaged across the GS.The time series of the NAO, EA and SCA indices are shown in Figure 8h.In (b-g) positive values represent an increase in ta_o and positive values represent an increase in ta_o.The coloured numbers in the bottom right corner are the statistically significant Pearson correlation coefficients between the teleconnection components and the ta_o time series where statistical significance is calculated from a t-distribution [Colour figure can be viewed at wileyonlinelibrary.com]

Note:
The actual trend is the ta_o trend from Figure 9a, the NAO, EA, SCA and EAWR columns represent the trends in the teleconnection component (dC=dt) of ta_o, and dR dt (Figure 10e-h) is the remaining trend when changes associated with the teleconnections are removed.Negative trends indicate a decrease in ta_o and positive trends indicate an increase in ta_o.Units are CÁdecade −1 and numbers in bold are statistically significant at the 95% level.F I G U R E 1 2 (a) Map of the ECA&D stations in Romania and surrounding countries which provide non-blended (circles) or blended (triangles) data.The data from the Bulgarian stations (triangles) are not publicly available.The coloured circles represent the Romanian stations which have data from 1950 to 1993 and the red circles represent the three "cold" stations in Romania which have data from 1961 or later.Panel (b) highlights these stations and also shows the time series of yearly ta_o ( C) at the other Romanian stations in grey.Panel (c) shows the area-averaged yearly ta_o for Romania and the surrounding countries.The dashed vertical lines in (b, c) highlight the years where most of the data is added for Romania (1961) and in which year the time series of Romanian stations which have data from 1950 stop (1993) [Colour figure can be viewed at wileyonlinelibrary.com] linked the growing season onset and end to circulation patterns with the European Climate Assessment & Dataset (ECA&D) blended station data for all of Europe.
was used to remove the teleconnection signals from the time series of each growing season index.Our key findings are: 1. NAO+ and EA+ in JFM (January-February-March) are both associated with early ogs10.The EA influences areas on the main continent and the NAO affects Great Britain, Ireland and northern France most strongly by a magnitude of about 10 days.2. The EA component of interannual variability has a statistically significant correlation with ogs10 in five (France, Italy, Poland and north Germany, Hungary, the Balkans) of six regions European regions considered.The NAO component correlates significantly with ogs10 for all six regions although the EA dominates in continental regions.Hungary has the best combined statistical relationships between ogs10 and the JFM NAO and EA. 3. The April-to-October NAO signal for gsr is mostly negative (NAO+ causes a drier growing season) in a wide band between 45 N and 60 N with a wetter growing season along the Norwegian coast, western Scotland and Italy.The EA and SCA signals are confined to northern and southern Europe with EA+ associated with a drier growing season across much of southern Europe, but SCA+ is linked to a wetter growing season in the same regions.4. The strongest correlations between the teleconnection components and gsr are in south Scandinavia with the EA (0.46), the Iberian Peninsula with SCA (0.41), and Poland with the NAO (0.4).South Scandinavia is the only selected region which has statistically significant correlations with gsr for three teleconnection patterns.Comparing the gsr and C time series shows that the SCA component visibly explains some of the interannual variability in southern Europe-particularly in the 1970s.5.The April-to-October EA signal for ta_o is positive (warmer growing season for EA+) for almost all of Europe.A negative signal for EAWR (cooler growing season for EAWR+) is confined to eastern Europe and the NAO and EA teleconnections have far more localized impacts on the fringes of Europe.Romania appears as an outlier in both the overall and residual trends for ta_o as it consistently has a cooling trend in contrast to the warming trend everywhere else in Europe.6. EA variability correlates significantly with ta_o in six regions around Europe and is strongest in the Italy and Adriatic region and the Iberian Peninsula.Yearly ta_o anomalies are also clearly explained by interannual variability in the EA component with some notable influence from the EAWR component in eastern Europe.7. Romania has large positive ta_o anomalies in the 1950s which are not found elsewhere in Europe.