OR/14/004 Worst case scenario research at BGS

Thomson, A W P (Editor), Beggan, C, Kelly, G, Baillie, O, Viljanen, A, and Ngwira, C. 2014. Project EURISGIC: worst case scenarios (Technical note D5.1). British Geological Survey Open Report, OR/14/004.

Recent research at BGS has focussed on applying the technique of extreme value statistics (EVS) to geomagnetic and, more recently, to geoelectric field data. EVS analysis involves the fitting of a generic function — in this case the generalised Pareto distribution (GPD) — to the tail of a distribution of data. From the best functional fit, such characteristics as the return level and return period for extreme events can be deduced, along with uncertainties.

Estimates of extreme values and their return intervals can then be used to investigate the extremes in geomagnetically induced currents (GIC) in, for example, the UK or European high voltage transmission grid, by means of a model of each grid and a model of the surface electric field generated by the extreme magnetic variations.

Some BGS results were previously reported in EURISGIC technical note D1.3 (see Discussion and conclusions), given in BGS report number CR/12/024 (Beggan and Thomson, 2012). These results are not repeated here.

There are 4 specific results we do report here:

• EVS as applied to EVS and European geomagnetic observatory data
• EVS as applied to EVS and Global geomagnetic observatory data
• EVS as applied to EVS and Nagycenk geoelectric observatory data
• Extreme GIC in the UK grid resulting from extreme European magnetic variations

EVS and European geomagnetic observatory data

In the paper by Thomson et al. (2011), EVS was applied to a number of decades of one-minute mean magnetic data from 28 magnetic observatories in Europe, to provide a preliminary exploration of the extremes in magnetic field variations and their one-minute rates of change. These extremes were expressed in terms of the variations that might be observed every 100 and 200 years in the horizontal strength (H) and in the declination (D) of the field.

Before fitting a GPD to each data set, Thomson et al. first determined a threshold to mark the onset of extreme behaviour. The ideal threshold should be low enough to allow for a meaningful number of samples, but high enough that the modified ‘scale’ parameter of the GPD is approximately constant and the GPD ‘shape’ parameter approximately linear (within error-margins), above the chosen threshold. It was found that setting the threshold at the 99.97th percentile proved reasonable for each variable at most observatories.

Clusters of extreme values occur during geomagnetic storms. These can be where the storm maximum is accompanied closely in time with other near-maxima. Including these near-maxima in the data skews the extreme value statistics in the following sense: we wish to identify the return period for each major event, i.e. single magnetic storm, in the sense that one would describe 13th March 1989, or 30th October 2003 as a single event, even though there was much complex structure in the activity that occurred on each of these days. By experimenting with a de-clustering algorithm with a variable de-cluster length (from three hours to one week), Thomson et al. showed that the return level is generally only weakly dependent on the de-cluster length for most observatories, but that de-clustering increases the extrapolated return level, compared to not applying de-clustering.

For each observatory Thomson et al. computed the peak residual and peak rate-of-change predicted by the observatory GPD to be exceeded once within periods of 100 and 200 years, from examination of the return-level statistics. The results for the four time-series are summarised in Figures 2.1 (H residual), 2.2 (dH/dt), 2.3 (D residual) and 2.4 (dD/dt). Note that the H residual is computed from the north (X) and east (Y) component as the square root of the sum of squares in X and Y. dH/dt is then computed through the change in H, minute to minute, and not from dX/dt and dY/dt. Thus dH/dt may underestimate the strength of changes in the vector horizontal field. (Future work will investigate the vector dH/dt.)

In Figures 2.1–2.4 Thomson et al. do not consider any local time, or longitude, dependence as significant as the latitude dependence. During any individual storm, when the peak electrojet activity occurs will define where the extreme is localised in longitude (if it is localised). However at the time scales of years between single storms and with de-clustering of 12 hours or more, such small time dependency of a few hours is probably not significant.

Thomson et al. found that both measured and extrapolated extreme values generally increase with geomagnetic latitude though there is a marked maximum in estimated extreme levels between about 53 and 62 degrees north. At typical mid-latitude European observatories (55–60 degrees geomagnetic latitude) compass variations may reach approximately 3–8 degrees/minute, and horizontal field changes may reach 1000-4000 nT/minute, in one magnetic storm once every 100 years. For storm return periods of 200 years the equivalent figures are 4–11 degrees/minute and 1000–6000 nT/minute. More details are given in Table 5.1.

From the GPD fits to each observatory data set Thomson et al. found that, for nearly all observatories and data types, the GPD shape parameter is positive, but generally less than one, and exceeds its standard error. This suggests that extreme geomagnetic activity is not bounded by some maximum, but follows either a Gumbel or Frechet distribution, i.e. the probability of increasingly higher values diminishes exponentially or polynomially with value. The consequence of this is that over even longer return periods, compared for example with 200 years, higher extreme values can be expected to occur.

Table 1    Estimated 100 and 200 year maxima in each of H, dH/dt, D and dD/dt between 55 and 60 geomagnetic degrees north, summarised from Figures 2–1 to 2–4. Figures in parentheses apply where Valentia (Ireland) observatory data are excluded, on the basis that these data seem anomalously low, compared to others of similar latitude.

	H (nT)	dH/dt (nT/min)	D (degrees)	dD/dt (deg/min)
100 Year Return	2000–5000 (4000–5000)	1000–4000 (1500–4000)	4–15 (7–15)	3–8 (5–8)
200 Year Return	3000–6500 (4500–6500)	1000–6000 (2000–6000)	5–20 (10–20)	4–11 (7–11)

Figure 2.1    The measured maximum for each observatory (top), and estimated 100-year (middle) and 200-year (bottom) return-levels, for H, in nT, as a function of geomagnetic latitude. 95% confidence limits are shown.

Figure 2.2    The measured maximum for each observatory (top), and estimated 100-year (middle) and 200-year (bottom) return-levels, for dH/dt, in nT/min, as a function of geomagnetic latitude. 95% confidence limits are shown.

Figure 2.3    The measured maximum for each observatory (top), and estimated 100-year (middle) and 200-year (bottom) return-levels, for D, in degrees, as a function of geomagnetic latitude. 95% confidence limits are shown.

Figure 2.4    The measured maximum for each observatory (top), and estimated 100-year (middle) and 200-year (bottom) return-levels, for dD/dt, in degrees/min, as a function of geomagnetic latitude. 95% confidence limits are shown.

EVS and Global geomagnetic observatory data

Although the results from the European wide EVS analysis (Thomson et al., 2011, and as summarised in the previous section) show some consistency between extremes for observatories at similar latitudes, there are instances of large differences in return levels, which are unexplained (see Figure 2.2 for an example). For this reason it is useful to put the European results into a global context.

In Figure 2.5 we therefore show the distribution of global observatories, colour coded according to the length of digital data sets each observatory holds (data from INTERMAGNET and the Edinburgh WDC for Geomagnetism). There is clearly an uneven distribution, favouring the northern hemisphere, although the three longest data sets are actually in the southern hemisphere.

Figure 2.5    Digital data colour coded by duration, for global magnetic observatories.

Figure 2.6    Digital data for 30 global observatories, selected for EVS analysis.

Figure 2.6 then shows the locations of 30 observatories provisionally chosen for an EVS analysis. This has a reasonably even distribution in both latitude and longitude, and includes some European observatories for a consistency check with the earlier work. However we note that to obtain this even distribution we find that some observatories have relatively short digital histories. This may impact on the statistical robustness, as was the case in the European study.

EVS has been applied to some of these observatory data sets (this is a work in progress as it requires some laborious quality assurance to remove spikes and outliers). We have analysed the residual minute means — and their rates of change — in the D, H, X (north) and Y (east) components, as well as daily block maxima for each component, for context. De-clustering lengths and threshold selections were set as in the previous paper.

In Table 2 we show EVS results for Eskdalemuir from a repeat analysis, carried out by a completely independent researcher, to assess the confidence in the results of Thomson et al (2011).

Table 2    Estimated 100 and 200 year maxima in each of H, dH/dt, D and dD/dt, X, dX/dt, Y, dY/dt for Eskdalemuir. Data used to compute these extremes are the daily maxima in the one minute observatory data and the maxima per block of 15 minutes. The results from the Thomson et al. (2011) paper are shown for comparison with the 15 minute block maxima.

Esk		D	dD/dt	H	dH/dt	X	dX/dt	Y	dY/dt
100 yr	Daily max	7.11	6.79	3351	1491	3340	1415	2208	1397
	15 Min max	6.62	4.65	4050	4044	3953	3577	3268	1643
	2011 paper	~7.5	~5.1	~5000	~3700
200 yr	Daily max	9.49	10.93	3938	1892	3949	1781	2908	2014
	15 Min max	8.95	6.94	5195	6668	5051	5798	4916	2500
	2011 paper	~10.0	~8.0	~6000	~6100

The comparison with the earlier paper is made between the rows labelled ‘minute’ (the new results) and ‘2011 paper’ (published results). Given the wide 95% confidence limits in the earlier work these new results are broadly consistent, which provides confidence in the statistical robustness of our approach. Interestingly comparison with the daily block maximum results (‘daily max’ row) shows some large differences particularly in dH/dt and dX/dt. This will be investigated further.

Eight observatories have been assessed so far and provisional results are shown in Figure 2.7 for the dH/dt component. Even allowing for the relative lack of data (yet) from the southern hemisphere we see a degree of hemispheric symmetry, e.g. as reported by Pulkkinen et al. (2012) and Ngwira et al (2013a), as well as the ‘auroral bulge’ of Thomson et al. (2011).

Figure 2.7    dH/dt 100 and 200 year return levels using the minute max values, a threshold of 99.97% and a de-clustering run length of 12 hours (plotted by geomagnetic latitude).

One issue that has slowed progress on processing the global data set is the difficulty in selecting a consistent extreme level threshold across all observatories. For some observatories, e.g. Fort Churchill (FCC), we find a high number of occurrences per year (~46) for the rate of change in each component, given a 99.97% threshold. This would suggest that the threshold is maybe not high enough. However for other observatories, e.g. Eskdalemuir (ESK), this is not the case. For this reason we have investigated the trade-off between clustering length, threshold and return level. This is shown in Figure 2.8 for FCC and Belsk (BEL).

From Figure 2.8, for both FCC and BEL (apart from at the 99.97% threshold at BEL) the de-clustering doesn’t make much difference to the return levels (i.e. estimates are well within the 95% confidence intervals) and, as the threshold increases, the occurrences per year obviously come down. For FCC the return levels and the size of the confidence interval (range in the plots) increase as the threshold increases. For BEL the opposite is true: both the return levels and the size of the interval decrease with threshold.

What this seems to suggest is that we may need to use a different threshold for B and dB/dt, in terms of each component and different thresholds between observatories to get a consistent picture. In other words to make sure the GPD fit is appropriate for each observatory separately. This is quite different from the earlier Europe-only study, for which de-clustering length and threshold was fixed for all observatories and all components. For example, we note that in general there are many more occurrences per year for the rate of change of each component compared to the component itself.

Figure 2.8    The 100 year return levels (left hand scale and blue/red marks and error bars) for a range of thresholds at FCC (top) and BEL (bottom) for the 100 year return period data. The blue and red colourings correspond to 12 hour and 24 hour de-clustering windows, respectively. The dark blue/red lines illustrate the size of the confidence interval. The ‘occurrences per year’ (right hand axis) indicate the average number of occurrences in each year for values above each threshold (light blue and brown line colour).

EVS and Nagycenk geoelectric observatory data

The Nagycenk observatory (Hungary) now has more than 50 years of geoelectric (or telluric) field data. In a poster presentation at the 10th European Space Weather Week, in Belgium, in 2013, Baillie et al. (2013) applied the technique of extreme value statistics to these data to determine the 1 in 100 and 1 in 200 year peak values of the surface electric field. These results are very preliminary and are subject to further work, as described below.

At Nagycenk (IAGA code NCK) potential differences are measured using 2 m deep, low polarizing, lead electrodes in the North-South (Ex) and East-West (Ey) directions, with an electrode spacing of 500 m and data recorded at 1 sec and 10 sec sampling intervals. The data resolution is about 6.1 μV/km.

The following steps were used to process the NCK telluric data:

• Approximately 19 years (Feb 1994 to Aug 2013) of digital data were analysed
• These data were 10-second NCK geoelectric data in Ex and Ey, the 3-hour NCK K index (geomagnetic activity measure) and the three-hour T-index (geoelectric activity measure)
• The data were pre-screened using the local K index >7, producing 106 days to analyse
• For each day a least squares fit, representing the daily baseline, was removed to leave the variations
• The variations were analysed using the eXtremes software toolkit (Gilleland and Katz, 2005) through the R statistical package (R Development Core Team, 2008)
• The maximum 10-second values per 3 hour time block were used as the basic data set, providing a manageable reduction in data size, from around 1 million samples to approximately 53 000 data
• A second data set of the maximum 10-second values per day were also analysed for comparison and consistency with the first data set
• The maxima for both Ex and Ey were determined for the time-span of data. The projected GPD distribution for periods of 100 and 200 years were computed and the 95% confidence levels in the extremes were also determined

In Figure 2.9 we show comparison plots of Ex and Ey for Eskdalemuir and Nagycenk for the storm of 17th March 2013. There are similarities in the data, particularly at longer wavelengths, which is reassuring, but the overall scales are quite different. This suggests, if true also for other storms, that EVS statistics will produce larger return levels for Eskdalemuir, compared to Nagycenk, in line with estimates produced elsewhere by EURISGIC project team members.

In Figure 2.10 we then show two examples of geoelectric data for NCK. Plots of one month of data for a ‘good’ month (top) and a ‘less good’ month (bottom) are shown, with ‘good’ and ‘less good’ defined in terms of the homogeneity and overall data quality. Steps and spikes are clearly visible in the ‘less good’ data and show the need for careful quality control. We note in particular that the E-field can be relatively large even when the observatory K is only moderately disturbed.

In Figure 2.11 the return level plots for the Nagycenk Ex and Ey data, as estimated by EVS are given. Probably the most striking result is the apparent asymptotic upper limit to Ex, with this being less evident in Ey. Whether this is physical and a robust result requires further study.

Figure 2.9    Electric field in the north (Ex) and east (Ey), for Eskdalemuir (left) and Nagycenk (right), during the storm of 17th March 2013.

Figure 2.10    Electric field in the north (Ex) and east (Ey) directions, NCK T-index and observatory K index for two months of data, to demonstrate quality control issues when pre-processing the data for EVS analysis.

Figure 2.11    Return level plots in the Ex (left) and Ey (right) for Nagycenk. Units for the return level are mV/km. 95% confidence limits are shown in blue. Circles denote measured data.

In summary, what our initial look at Nagycenk data has shown is that:

• Both the 100 and 200 year return periods for NCK geoelectric data are about 45 mV/km and 54 mV/km for Ex and Ey respectively. These small levels of geoelectric field, compared to Eskdalemuir, may reflect the sedimentary nature of the region and its more conductive nature
• The return level curves tend to ‘saturate’ meaning that such ‘extremes’ are already in the data set. This is unexpected and indicates that further work is required
• Using daily maxima of 10-second data produces similar results (the second data set referred to earlier and not shown here)
• There is a sensitivity of the results to the choice of threshold (currently set at ~33% of the peak level, which is substantially less than for similar treatments of geomagnetic data, e.g. Thomson et al., 2011)

Future work will include:

• Repeating the analysis by including days when the local 3 hour geomagnetic K index is 6 (adding ~500 days more for processing and analysis), or through better automated correction of baselines, jumps and spikes
• Repeating the analysis with recently digitised (quality checked) telluric data for more of the 50 years in the NCK record
• Re-computing the NCK T index with an unbounded upper limit and analyzing this data set via EVS
• In the long term, comparing results with similar analyses of Eskdalemuir and other UK station data

In 2012, geoelectric monitoring sites were installed at each of the three UK geomagnetic observatories. Data are now being recorded at these sites with the intention, over the long term, of comparing with the Nagycenk data, to provide a wider European scale view of surface geoelectric fields.

For an initial look, in Figure 2.12 we show the distribution of electric field values from the Eskdalemuir observatory, for 1st January to 1st August 2013. It can be seen that there is a long tail in the distribution, in both the north (Ex) and east (Ey) components. Moreover there are clear difference between the Ex and Ey channels. This may be physical (local effects?) or may reflect issues with the instrumentation. We note that the local maxima during this time interval (of ~0.2V/km) compare with an estimated 5V/km for central Scotland during the October 2003 storm (Thomson et al., 2005) and even estimates of ~20V/km for genuinely extreme events in resistive geological terranes (e.g. Pulkkinen et al., 2012).

Over time, probably many years, we hope to construct a data set from which we can analyse the tail of the UK distributions using EVS, with a view to confirming, or otherwise, estimated extremes in modeled electric fields, such as given by Pulkkinen et al. (2012), or the results obtained from Nagycenk.

Figure 2.12    Electric field distributions in the north (Ex) and east (Ey) directions at Eskdalemuir, for 1st January to 1st August 2013.

Extreme GIC in the UK grid

Using the EVS results of Thomson et al. (2011), summarised in Section 2.1, Beggan et al. (2013) investigated the impact of extreme geomagnetic variations in terms of worst case GIC in the UK transmission system. Beggan et al. firstly created an updated model of the UK surface conductivity by combining a spatial database of the UK geological properties (i.e., rock type) with an estimate of the conductivity for specific formations. Secondly, they developed and implemented a sophisticated and up-to-date model for the 400 kV and 275 kV electrical networks across the whole of Great Britain and, in addition, for the 132 kV network in Scotland.

They were then able to deduce the expected GIC at each transformer node in the system based on the network topology and an input ‘extreme’ surface electric field. Beggan et al. applied these developments to study the theoretical response of the present-day UK high-voltage power grid to modeled extreme 100 year and 200 year space weather scenarios and to a scaled version of the October 2003 geomagnetic storm, approximating a 1 in 200 year event. The analysis of Beggan et al derived much from Thomson et al. (2005), Beamish et al. (2002) and McKay (2003).

As a first hypothetical model, Beggan et al. developed two synthetic electrojet model profiles: an electrojet model with an amplitude profile akin to a ‘top-hat’ function, extending from 53N to 63N in geomagnetic latitude, with a second model that had a ‘tapered cosine’ profile extending between 48N and 68N in geomagnetic latitude. They used the two different models to examine if the amplitude gradient (slope) of the magnetic field strongly affects the estimated GIC. The Top-Hat model gave a very strong gradient across its edges while the Tapered-Cosine model had a gentler gradient. Two orientations of the auroral electrojet were computed for each synthetic 100 and 200 year electrojet: (a) geomagnetically east-west aligned across the UK and (b) a geomagnetic north-south alignment (approximately following the central axis of the UK).

To scale each electrojet model to the correct amplitude for an extreme event, the results from the Thomson et al. [2011] study on the statistical predictions of extreme values in European magnetic observatory data were applied. Beggan et al. used 1000 nT/min, 3000 nT/min, and 5000 nT/min in their analysis to approximate the expected maximum in dH/dt for 30, 100, and 200 years. For an assumed driving sinusoid of period of 2 min, this led to magnetic field input strengths H of approximately 450 nT, 1350 nT, and 2250 nT, corresponding to dH/dt = 1000 nT/min, 3000 nT/min, and 5000 nT/min.

As a final, more ‘realistic’ model, a model of the magnetic field during the October 2003 event was constructed based upon the measurements from nine observatories and variometers around the United Kingdom and North Sea region, scaled by a factor of five, again according to the results suggested by Thomson et al. (2011). This provided an event with peak magnitude of 5000 nT/min at one instant (21:20 UT).

The estimated electric field for the scaled October 2003 model is shown in Figure 2.13. In Figures 2.14 and 2.15 we show the GIC results for Beggan et al, for, respectively, the synthetic electrojet models and the scaled October 2003 model.

It is clear that substantial GIC, of more than a few hundred Amps, are possible, according to these scenarios. The largest GIC typically occur in the north of the UK, and in the ‘corner’ nodes of system (e.g., southwest Wales and England) or in isolated regions (Scottish Borders). Where nodes lie close together, especially in the southern UK, there is a tendency for smaller GIC (e.g., London/southeast England), though this is not necessarily the case with other clusters of transformers (e.g., northeast England). It is also found that any North-South component to the electrojet over the UK (e.g., during a westward travelling surge) will on average increase the GIC flowing in transformer earths. In places, this GIC can be an order of magnitude greater than that for the more commonly East-West oriented electrojet.

Figure 2.13    Electric field induced in the surface for period of 120 s due to magnetic fields from an extreme version (x5) of the 30 October 2003 geomagnetic storm. The columns show the (left) Y component and the (right) X component. Nominal times (in UT) are illustrative, taken from the time profile of the October 2003 storm.

Figure 2.14    GIC in the National Grid GB high-voltage network due to a 100 year extreme scenario (120 s period) from an auroral electrojet with the following configurations: (a) Tapered Cosine East-West aligned, (b) Tapered Cosine North-South aligned, (c) Top Hat East-West aligned, and (d) Top Hat North-South aligned. Blue indicates GIC directed into the grid, red indicates GIC into the ground. Circle size represents size (relative to scale). Note that many sites have multiple transformers present.

Figure 2.15    Snapshots of GIC in the National Grid GB high-voltage network due to an extreme storm scenario (approximately a factor of 5) of the 30 October 2003 geomagnetic storm (due to an electric field with a period of 120 s). (a) Time: 19.30 h, (b) Time: 20.50 h, (c) Time: 21:20 h, and (d) Time: 22:50 h (see Figure 2.5). Blue indicated GIC directed into the grid, red indicates GIC into the ground. Circle size represents size (relative to scale). Note that many sites host multiple transformers.