Introducing novel risk-based indicator for determining transmission line tower's backflashover performance

Abstract Transmission line (TL) backflashover (BF) performance has been traditionally ascertained using a single number, the so-called backflashover rate, which is measured in number of BF events per 100 km-years. This paper aims at presenting a novel indicator for assessing performance of high-voltage TL tower's ability for tolerating direct lightning strikes without provoking BF events. This indicator is also defined as a single number, which can be computed for any TL tower (by means of EMTP simulations) and measures, in a novel way, its tolerance against BF events. It is given in terms of the risk of the BF occurrence, which means it is statistical in nature and depends on the total sum of conditions governing the BF events. The BF risk, as an indicator, is obtained from the probability density function of the shield wire(s) incident lightning currents and the cumulative distribution function of the BF currents statistical distribution. Hence, it merges complete probabilities of obtaining lightning currents striking a tower with probabilities of those currents provoking BF events on that tower. This novel risk-based indicator can be related to the price of that risk and the associated investment costs, enabling the cost-effective optimisation solutions to the problems of TL arrester applications and station insulation coordination design. This seems appropriate, considering the fact that the investment in surge arresters and related protective measures is commonly perceived as buying insurance.


Introduction
High voltage (HV) transmission lines are exposed to lightning strikes, where only direct lightning strikes (to shield wire(s), phase conductors and tower tops) are of engineering concern. Transmission line (TL) performance in relation to direct lightning strikes is of paramount importance, for several different reasons [1,2]: (i) determining line's yearly outage rate for reliability purposes, (ii) line insulation coordination, (iii) surge arresters (TLA) application on the line, (iv) incident transformer station (or switchyard) insulation coordination design, (v) optimising TL insulators archorn gaps, (vi) optimising tower geometry for lightning shielding of the line. Of particular importance, for several of the above-mentioned instances, is the TL backflashover (BF) performance, emanating from the direct lightning strikes to the tower tops and shield wire(s), [1][2][3][4][5]. Transmission line BF performance is traditionally ascertained using the backflashover rate (BFOR), which is a single number representing an entire line, expressed as the number of expected BF events per 100 km-years, [5]. The backflashover probability, as a feature of the BFOR, is usually obtained from the repeated numerical simulations (i.e. Monte-Carlo method), e.g., [6][7][8][9][10][11] using the ElectroMagnetic Transients Program (EMTP), [12,13], or by other means (i.e. simplified analytical treatment). When this probability is combined with the number of expected BF events, traditionally determined from the application of the electrogeometric model (EGM) of lightning attachment to TLs, it yields a backflashover rate.
Analytical methods of backflashover analysis are extensively described by the IEEE WGs [5,14] and CIGRE WGs [15]. A comparison between these recommendations is given by Nucci in [16]. Nowadays, it is far-more common to treat the backflashovers on TLs in terms of the numerical simulations, carried-out by means of the EMTP, e.g., [12,13,17]. With the numerical approach to the transient analysis of TL lightning surges, detailed models of the TL components are needed, some of which exhibit non-linear behaviour or frequency-dependence [17][18][19]. Interested reader is at this point advised to consult the extensive treatment of modelling guidelines for TL lightningsurge numerical simulations provided in Ref. [20,Ch. 2]. Further important simulation details, concerning the backflashover analysis on HV transmission lines, can be found in Refs. [21][22][23][24][25].
This paper aims at presenting a novel indicator for measuring performance of HV transmission line tower's ability to tolerate direct lightning strikes without provoking a BF event. It can be determined for any single tower or given for the entire line (using a "representative" tower). This indicator is also a single number, which leverages powerful EMTP simulations in its computation, and measures, in a novel way, tower's tolerance against BF events. It is given in terms of the risk of the BF occurrence, which means it is statistical in nature and depends on the total sum of conditions governing the BF events [26]: lo-cal keraunic level, statistical depiction of lightning-current parameters (including statistical correlation between the parameters), EGM of lightning attachment, frequency dependence of TL parameters and electromagnetic coupling effects, TL span length, statistical distribution of lightning strokes along the TL span, tower geometry and surge impedance, tower grounding impulse impedance (with soil ionization if present), lightningsurge reflections from adjacent towers, non-linear behaviour of the insulator strings flashover characteristic, and power frequency pre-strike voltages.
Proposed BF indicator, in addition, utilizes the so-called curve of limiting parameters (CLP), which is derived from repeated EMTP simulations using an original algorithm developed by the authors and described in detail in Ref. [27]. The CLP itself brings into direct relationship shield wire(s) incident lightning currents with the "critical" currents for the BF occurrence. Developed algorithm minimizes the number of EMTP simulation runs using systematic simulations approach, unequal wave-front time increments, and bisection search for finding "critical" current amplitudes. Consequently, it is orders of magnitude faster then the traditional Monte-Carlo method application.
Furthermore, pseudo-random shield wire(s) incident lightning currents, necessary for the statistical treatment of the phenomenon, are generated from the appropriate bivariate statistical probability distribution, by means of the Gaussian Copula approach [28][29][30][31][32]. The Copula approach, when combined with the CLP method, provides for an extremely efficient way of obtaining pseudo-random BF currents, unlike the more traditional way of using the Monte-Carlo method (which is known to be very time consuming), e.g. see [9]. These BF currents are in-turn used to derive a probability density function (PDF) of their statistical distribution, by means of the kernel density estimation (KDE) procedure [33]. The KDE employs Gaussian kernels, with bandwidths determined using the grid search and cross-validation of the estimator performance.
The BF risk, as an indicator, is finally computed from the PDF of the shield wire(s) incident lightning currents and the cumulative distribution function (CDF) of the BF currents statistical distribution. Hence, it merges complete probabilities (instead of working with their point estimates) of obtaining lightning currents striking a tower with probabilities of those currents provoking BF events on that tower (while accounting for any tower peculiarities as such).
The paper is organised in the following manner. In Section 2, a brief outline of the TL model for the BFOR analysis is provided, which is suitable for the implementation in the EMTP software package. Section 3 provides a thorough and comprehensive statistical treatment of lightning currents. This section presents the Gaussian Copula approach, along with the EGM application, for obtaining the bivariate probability density of lightning currents incident to transmission lines. It also features BF currents probability distribution, obtained from using the CLP from the prospective tower. Section 4 introduces a novel risk-based BF indicator. A test case of the HV transmission line, along with the sensitivity analysis, is provided in Section 5, which is followed with the conclusion in Section 6.

Transmission line modelling for backflashover analysis
The EMTP model of the HV transmission line for lightning surge transient simulation in general, and backflashover analysis in particular, has been thoroughly studied and widely published, see Refs. [17-19, 34, 35]. A brief outline of the model, as employed for the purpose of this paper, will be presented in this Section. The model consists of several components: (i) TL phase conductors and shield wire(s), including spans, line terminations and power frequency voltage, (ii) TL tower, (iii) tower grounding impedance, (iv) insulator string flashover characteristic, (v) lightning current and lightning-channel impedance.

Phase conductors and shield wire(s)
High voltage transmission line phase conductors and shield wire(s) are modelled as distributed-parameters, untransposed, frequency-dependent, multiphase transmission line, by means of the so-called JMarti model [12,36]. Phase conductors and shield wire(s) positions on the tower (from the mostrepresentative tower within the TL route) are used, along with their maximum allowed sags, cross-sectional dimensions, DC resistances, ground resistivity of the ground return path, etc. Several spans of the transmission line, at each side of the tower being struck by lightning, are modelled in this way. Longer sections are then added on both sides of this chain in order to suppress further reflections. The line model is finally terminated by an ideal, grounded, power-frequency, three-phase voltage source (with fixed angles) in order to account for the pre-strike phase voltages.

Transmission line towers
The steel-lattice towers of HV transmission lines are represented as a single conductor, distributed-parameter, frequencyindependent lines. The single value of the tower surge impedance is computed from the analytical expressions, based on the theoretical background provided by Sargent in [37], which depend on the tower configuration and can be found in Ref. [20,Ch. 2]. The velocity of the surge propagation along the steellattice tower is assumed to be equal to the speed of light in free space.

Tower grounding impedance
The model of the lightning-struck TL tower grounding impedance is probably the most important single factor influencing the subsequent formation of the overvoltage on its tower top, due to subsequent reflections of travelling waves formed by the said lightning strike, where reflections from the tower base will arrive much sooner at the tower top then reflections from adjacent towers. Hence, the influence of the (apparent) TL tower footing (i.e. grounding) surge impedance on the tower top transient voltage is determined by its response time and current dependence. In cases where the tower grounding can not be regarded as being concentrated, the frequency-dependent tower grounding impedance model is needed. Its implementation can be difficult [24]. On the other hand, in case of the concentrated grounding, the tower grounding impedance exhibits only current dependence, which can be modelled in accordance with the guidelines provided in Ref. [18]. The following equation is utilised for that purpose: where R 0 is the tower grounding resistance at low frequency and low current magnitudes, (Ω), I t is the lightning current through the tower footing impedance, (kA), and I g is the lightning current level which determines the soil ionization inception process, (kA). This current is determined from the following expression: where E 0 is the soil ionization electric field gradient, provided in (kV/m) and ρ is the apparent soil resistivity in (Ωm). Grounding impedances of the adjacent towers are modelled with a simple lumped resistance, the value of which is equal to the tower grounding resistance at low frequency and low current magnitudes.

Insulator strings flashover characteristic
The insulator strings flashover characteristic is a non-linear function of the applied impulse voltage and it exhibits complicated behaviour in nature, which is rather difficult to fully reproduce inside EMTP computational framework. Interested reader is advised to consult Refs. [22,[38][39][40] for more information. It is usually modelled in the EMTP by means of the voltage-controlled switches (closure of this switch represents a flashover of the insulator strings). The flashover physics provides the "logic" behind the switch and is programmed using the MODELS language [36]. A leader progression model is used here for the insulator flashover, which is based on the solution of the following differential equation: where: d g is the insulator strings length, (t) is the leader length, v(t) is the actual (absolute value) voltage on the insulator strings, and k, E 0 are constants which are found to be dependent on the type of the insulator.

Lightning channel
A direct lightning strike to the TL tower top is modelled as a ideal current source in parallel with the resistance, the value of which represents the lightning-channel surge impedance. The current source can be of different types [36]: ramp-slope type, CIGRE type, double-exponential type, and Heidler type. The Heidler type of lightning-current source will be utilised for the purpose of this paper.

Bivariate Gaussian Copula
The bivariate Gaussian copula cumulative distribution function can be described by the following relationship [29]: where Φ 2 is the CDF of the bivariate standard normal distribution and Φ −1 is the inverse cumulative distribution function of the univariate standard normal distribution. Also, bivariate Gaussian copula density function is given by: with: where ρ is the (linear) correlation coefficient between the variates. Then, in accordance with the Sklar's theorem, an arbitrary bivariate CDF, F(x, y), can be described in terms of its marginal distributions and a copula function, [29]: where F 1 , F 2 are marginal CDFs. Furthermore, an arbitrary bivariate PDF, f (x, y), can be described by the following expression: with newly introduced functions f 1 , f 2 denoting marginal PDFs. Moreover, the usage of the copulas approach enables for easy and fast generation of correlated random variates from the arbitrary bivariate distributions with differing and unrelated marginals [28,30,31]. First, a pseudo-random sample (u, v) is generated from the bivariate Gaussian copula, which holds a desirable correlation structure. Next, using the inverse CDFs of independent marginal distributions, the random sample from the appropriate bivariate probability distribution is obtained as follows [30]: In this way, marginal distributions are independent and completely separated from the correlation structure of the final bivariate probability distribution, but the correlation is preserved nonetheless.

Ground level lightning currents
It is well-known that each individual parameters of a general downward (i.e. ground level) lightning current follows a Log-Normal distribution, with a PDF defined as follows [41]: where x µ is the median value and σ ln x is the associated standard deviation of the ln x, which holds for any of the lightningcurrent parameters seen individually. Its associated complementary CDF can be described with the expression: where erfc stands for the complementary error function from the general statistics and Table 1 features basic parameters of the individual statistical distributions of negative downward lightning currents, which will be utilised hereafter [41].

Incident lightning currents
However, transmission line tower's incident lightning currents need to be drawn from the appropriate bivariate statistical probability distribution, which accounts for the geometry of the TL tower (in accordance with the EGM) and statistics of lightning-current parameters (including correlation). This is achieved by, first, generating a pseudo-random sample (u, v) from the bivariate Gaussian Copula which holds a desirable correlation structure. In case of TL tower's incident lightning currents described here, one of the marginal distributions in (10) is the Log-Normal distribution of lightning wave-front times, obtained from (12) with the introduction of appropriate parameters (see Table 1). The other marginal distribution is the statistical distribution of TL tower incident lightning-current amplitudes, obtained from the application of the EGM of lightning attachment to the transmission line at hand.
According to [42,Ch. 7], PDF of shield wire(s) incident lightning-current amplitudes can be obtained from the following expression: with where D g (I) and D g (I) are exposure distances for the shield wire(s), while I m is the maximum shielding failure current, all of which depend on the tower geometry in accordance with the EGM theory. It is known that the EGM brings into functional relationship the lightning current amplitude with its striking distance to phase conductors, shield wire(s) and to ground surface, and is of the form r = A·I b (m), where A and b are model parameters and I is the lightning-current amplitude in (kA). Moreover, CDF of lightning-current amplitudes incident to shield wire(s), which is in-fact a second marginal distribution used in (10), can be described by the following expression [42,Ch. 7]: The associated maximum shielding failure current (I m ) can be obtained from: with [42,Ch. 7] and [20,Ch. 6]: where: r c and r g depict striking distances to phase conductors and ground surface, respectively; h is the height of the shield wire(s) on the tower; y is the height of the phase conductor(s) at the tower and a is the length of the tower arm(s) carrying the phase conductor(s), all supplied in meters. Finally, exposure distances for the shield wire(s) can be, in accordance with the EGM model, determined from the following expressions [42,Ch. 7]:

Backflashover currents
Backflashover lightning currents probability distribution is derived by further utilizing the so-called curve of limiting parameters; see Ref. [27] for more information. The CLP brings into relationship incident lightning currents with the "critical" currents for the BF occurrence and can be constructed in the coordinate space of lightning-current amplitudes and wavefront times using the systematic EMTP simulations approach. Namely, for each wave-front time there is a single value of the lightning-current amplitude (i.e. critical current) that causes a backflashover (based on the analysis using the EMTP program). Any amplitude above this critical level, for that particular wavefront time, will certainly cause a BF event (due to the determinism of the EMTP computational framework); any amplitude below this threshold will not.
The CLP, furthermore, depends on the transmission line itself, featuring all main aspects of the BF phenomenon, from the insulator strings flashover characteristic to the tower footing impulse impedance [26,27]. Hence, by superimposing the appropriate CLP (obtained from the EMTP simulations of BF events on the particular tower) on the shield wire(s) incident lightning (generated from the appropriate bivariate probability distribution), a statistical distribution of the BF current amplitudes directly follows. This will be demonstrated on a concrete example. In a sense, the CLP curve can be seen as a filter, which passes through only the BF current amplitudes, from the total pseudo-random population of shield wire(s) incident lightning currents. This "filter" also preserves the correlation structure between the lightning current parameters (wave-front times and amplitudes), which has been incorporated from the very beginning.
The obtained BF current amplitudes can then be fitted by a probability distribution function-by means of the kernel density estimation (KDE) procedure, as follows [33]: wheref (x) is the estimated PDF, with K(•) being the kernel (e.g., Tophat, Cosine, Gaussian, Epanechnikov) and h > 0 is the bandwidth. Gaussian kernels are used here. The problem of bandwidth selection is tackled numerically, by means of implementing the grid search and cross-validation procedures for evaluating estimator performance, [33]. This PDF forms the basis for obtaining (numerically) the associated CDF, as follows: The backflashover event as a stochastic phenomenon, which is governed by the (i) probability of lightning (wave-front and amplitude combination) striking the tower and (ii) probability of that same lightning subsequently evoking a BF on that tower, is now statistically fully described.

Estimating number of dangerous events
The problem of estimation of the number of direct lightning strikes to transmission lines has been tackled by means of the EGM of lightning attachment for quite some time, and the literature on the subject is exhaustive, e.g. see Refs. [42][43][44]. According to the EGM of lightning attachment to transmission lines, and considering only vertical lightning strikes, the number of lightning strikes to phase conductors and shield wire(s) depend on their exposure areas, which are determined in terms of the lightning striking distance and tower geometry. Following expression for estimating the number of direct lightning strikes to shield wire(s) can be obtained [42,Ch. 7]: where newly introduced variables have following meaning: N g = 0.04T 1.25 d in (km −2 year −1 ) is the annual average ground flash density (T d is the long-term average annual number of thunderstorm days), f (I) is the PDF of the ground level lightning current amplitudes distribution defined by (11), L is the transmission line length and S g is the distance between shield wires (zero in case of single wire).
Lightning strikes are not uniformly distributed along the TL span length, which has been shown by Mousa in Ref. [45]. If one considers only lightning strikes to TL tower tops (and their near-vicinity), which is often the case with the BFOR analysis, then the expression (23) needs to be corrected with the appropriate coefficient which takes into the account the actual statistical distribution of lightning strikes along the TL span length. For example, according to the analysis provided in Ref. [45] following expression is valid for HV transmission lines (above 230 kV): P t = 5.486 · 10 −7 2 − 7.452 · 10 −4 + 0.3587 (24) where is the length of the TL span in meters; for HV transmission lines below 230 kV expression the P t value needs to be increased by a factor of 0.02. If one considers lightning strikes collected by the TL tower top and the shield wire(s) up to the 55 % of the half-span from each side of the stricken tower, then, in accordance with [45], it yields that P t = 0.6.

Backflashover risk indicator
By using the PDF of the shield wire(s) incident lightning current amplitudes statistical distribution ( f gw ) from (14) with the CDF of the backflashover amplitudes statistical distribution (F BF ) from (22), one can establish the risk of the TL backflashover, obtained from: as the area underneath the bell-shaped curve formed by the product of distributions f gw (I) andF BF (I). Furthermore, by bringing into the relationship the risk of backflashover with a yearly expected number of direct lightning strikes to the TL tower (i.e. number of dangerous events per year), enables one to arrive at the mean time between backflashovers (MTBBF) as a measure of the tower's tolerance to BF events. The expected number of direct lightning strikes to the tower is obtained from the EGM application to the tower geometry (including keraunic level, span length, conductor sag, statistical distribution of lightning strikes along the span length), in accordance with (23). Hence, the mean time between backflashovers can be obtained as follows: This measure, which is a single number, can be used to create a rang list of towers and identify those "rogue" towers with excessive BF risk (having, statistically speaking, a small value of MTBBF years). A similar notion of the mean time between failures (MTBF) is well established in the field of station insulation coordination; see IEC 60071-2 for more information [1].
In addition, in the case of transmission line arrester (TLA) applications, following iterative procedure can be followed: (i) preselect the surge-protective measures; (ii) estimate the BF risk with the preselected design; (iii) calculate the price (C f ) associated with the total risk (which includes repair costs, lost revenue due to outage time, and any additional surcharges), and determine if R BF · C f > C p (27) where C p is the cost of the TLA design (sum of expenditures for surge-protective measures), and return to (i) as necessary.
A computer program has been developed for automating the heretofore presented BF risk computation procedure. The program is written in Python language and interacts with the EMTP computational framework-where the TL model is created and maintained-through manipulating its input and output files (these are well designed files with known structure), [13,36]. The main computational effort is spent on deriving the CLP for a particular tower (or TL itself if "representative" tower is used). This can be achieved with less then a thousand systematic EMTP simulation runs, using specially designed algorithm; see Ref. [27] for more information. For comparison, Monte-Carlo method needs in excess of 100,000 EMTP simulation runs to obtain similar information [9]. Once the CLP has been obtained, filtering out BF currents from the large population of shield wire(s) incident lightning currents (and constructing BF currents probability distribution) takes only a couple of seconds.

Test case transmission line and sensitivity analysis
Heretofore presented methodology for estimating the BFOR on HV transmission lines will be demonstrated on the typical single-circuit 110 kV transmission line with vertical conductor configuration and steel-lattice towers. Tower geometry is typical for wind pressures between 750 − 1500 N/m 2 , with individual spans of 350 m, typical for 750 N/m 2 wind pressure and 65 N/m 2 of maximum allowed conductors tensile strength [10]. Tower geometry follows: h = 27 m, y = 24 m; span between tower consoles (arms) is 2 m; top console length a = 2.5m, middle console length is 3 m and bottom console length is 3.5 m. Maximum phase conductor sag equals 3 m while that of the shield wire equals 2 m. Phase conductor DC resistance is 0.114 (Ω/km) with 10.95 mm diameter, and that of the shield wire is 0.304 (Ω/km) with 8 mm diameter. Insulator string (i.e. archorn) length equals 0.9 m.
Furthermore, different models of insulator flashover characteristic and tower grounding transient impedance have been considered, as variants provided in Table 2. Simple switch model depicts here a TACS switch which is closed (signifying flashover) when the voltage across the insulator string exceeds a fixed value of 605 d, with d being the length of the insulator strings (there is no physics-based "logic" behind this switch). The Weck's model, as already stated in Subsection 2.3, assumes concentrated tower grounding system and current dependence due to soil ionization. Simple (tower grounding) resistance model, on the other hand, assumes that the tower footing surge impedance is equal to its low-current and low-frequency value and does not change during the simulation. Obviously, some of these models are oversimplified and even inappropriate (although still employed sometimes in literature), but have been introduced here in order to provide for the sensitivity analysis. Firstly, Fig. 1 graphically presents the bivariate Gaussian Copula probability distribution (scatter and density plot), which holds the desired correlation structure (see Table 1). The associated marginal distributions are the uniform distributions, as expected.  Next, Fig. 2 depicts the bivariate (i.e. joint) PDF of lightningcurrents incident to the TL shield wire, obtained with the EGM according to Brown and Whitehead for the assumed ground flash density. This bivariate PDF has been constructed using the Gaussian Copula approach, with the individual marginal distributions defined in Subsection 3.2. It should be emphasised that this is not a bivariate Log-Normal distribution! As can be seen, Fig. 2 presents a scatter plot of lightning data with superimposed bivariate PDF, marginal distributions of the amplitudes (right-hand part) and wave-front times (top part), along with the original underlying Gaussian Copula (in the top right-hand part of the figure). It can be seen that the statistical correlation between the variates has been preserved in this non-typical bivariate distribution. Not all of these lightning-currents incident to TL shield wire will cause a BF event; only a certain portion of them will, in accordance with the curve of limiting parameters of a particular tower. Figure 3 presents curves of limiting parameters obtained for "Variant A" of TL model and several different values of TL tower footing impedance. A Blackman filter has been applied to the CLPs in order to smooth the curves and reduce the influence of unequal wave-front time increments and determinism of the EMTP computational framework [26,27]. It can be nicely seen from this figure that the critical current amplitudes increase (which means that the BF probability decreases) with the lowering of the tower footing impedance [27]. This is a known phenomenon and a reason for reducing the tower footing (impulse) impedance. It is also clear that the critical currents increase as the wave-front duration is increased, which is also expected. In fact, for very long wave-front times the associated amplitudes attain very high values, meaning that the flashover is highly improbable, regardless of the tower footing impedance. Resistance values in Fig. 3 represent tower grounding resistance values at low-frequency and low-current magnitude, while their impulse behaviour is determined during the EMTP simulation, in accordance with the Weck's model (see Table 2 and Subsection 2.3).
In addition, Figs. 4 and 5 present CLPs obtained for different TL model variants described in Table 2, with two different values of low-current and low-frequency tower footing impedances as input parameters. It can be readily discerned from these figures that the CLPs (and by association the BF probability) strongly depend on the modelling of these two crucial TL components: tower grounding and insulator flashover. This has already been established by many researchers and cor- roborates the need for detailed, sophisticated TL tower grounding models that can reproduce observed backflashover performance. It is quite clear from these figures that the simple switch model of the insulator flashover (with its stepwise CLP that flattens for longer wave-front times) fails to reproduce the flashover phenomenon at the wave-tail of the incident overvoltage impulse waveform. In other words, the up-turn of the CLP, that is so characteristic of the leader progression models, can not be reproduced by the simple switch model. Furthermore, Fig. 5 gives a clear visual demonstration of the effect of soil ionisation on the tower BF performance-which emerges from the comparison of CLPs produced from Variant A and Variant B of the TL model-and establishes the fact that inclusion of the soil ionisation process increases the critical current levels. This is a known phenomenon which stems from the fact that onset of the soil ionisation essentially dynamically lowers the impulse impedance of the tower grounding. As can be seen from Figs. 4 and 5, the effect of soil ionisation is more pronounced for larger values of the tower grounding impedance, which is again expected.
One can now easily extract BF currents from the shield wires(s) incident lightning currents, by applying any one of the presented CLPs onto the bivariate distribution from the Fig. 2. This is illustrated in Fig. 6 for the "Variant A" and a case of 50 Ω tower footing resistance, in terms of the scatter plot with superimposed histograms for each of the marginal distributions. In this figure, grey shaded marginal PDFs represent shield wire currents probability distributions, while orange shaded PDFs represent BF currents probability distributions. It ought to be emphasized once more that the BF currents obtained by this process account for the lightning statistics (including correlation), TL geometry, EGM according to Brown and Whitehead, and EMTP transmission line model specificity (such as the insulator flashover model and the model of the tower grounding impulse behaviour).
At the same time, Fig. 7 presents only BF currents (from Fig. 6), once they have been filtered out from the total shield wires incident lightning currents. On the same figure are the marginal PDFs of the BF currents, along with superimposed CLP which have been used (red line) for the purpose of extracting the BF currents. In addition, CLPs for the used "Variant A", with several different values of tower footing resistances, have been provided in this figure (gray lines), for comparison. It should be mentioned that Figs. 6 and 7 hold some of the same information, presented somewhat differently for clarification and enhanced visual exposition of the computational procedure.
The approach depicted graphically in Figs. 6 and 7 can be easily repeated for any value of the tower footing impedance and any of the TL model variants, using its appropriate CLP curve. It is straightforward and numerically very fast, even with a large initial population of shield wire incident lightning currents (e.g., a sample of 2x500,000 elements). In that regard, Fig. 8 depicts the PDFs, and the associated CDFs, of the BF current amplitudes for the "Variant A" with two different cases of tower footing resistances-which are the final products of this process (notice the usage of the two vertical (ordinate) axes in this figure). The PDF of this statistical distribution has been produced by means of the KDE procedure (see Subsection 3.4). It can be seen from Fig. 8 that the lower value of the tower grounding impedance shifts the mode of the associated BF currents statistical distribution to the right and produces a heavier tail. It can also be observed from CDFs of Fig. 8 that the increase in tower footing impedance causes, for any particular lightning current amplitude, associated increase in the BF probability. This relationship is highly non-linear and dependent on TL model specificities. It should be borne in mind that the lightning current amplitudes of Fig. 8 represent at the same time the BF currents (i.e. critical currents for BF events), with incorporated full account of the associated (i.e. statistically correlated) wave-front times.
Finally, by employing the PDF of the shield wire(s) incident lightning-current amplitudes statistical distribution ( f gw ) with the CDF of the backflashover amplitudes statistical distribution (F BF ), one can establish the risk of the TL backflashover. This has been carried-out here, for the "Variant A" with two different TL tower footing impedances, and graphically presented in Fig. 9. This figure depicts, at the same time (using two vertical (ordinate) axes), the PDF of shield-wire incident lightning amplitudes and the CDFs of the BF currents (for two different treated tower footing resistances). It also graphically presents the product of these two distributions ( f gw ·F BF ), the area below which represents the BF risk (see the inset figure), in accordance with (25). It can be seen that the BF risk increases with the increase in the tower footing impedance, when all other aspects of the transmission line are being held constant. This is, of course, expected. Risk can be seen as a measure of the tower's ability to tolerate direct lightning strikes without provoking a backflashover. The higher the risk is, the lower will be the tower's ability to tolerate direct lightning strikes without provoking a BF event, and vice-versa.
When the BF risk is associated with the yearly expected number of dangerous events, it yields the associated MTBBF, in accordance with (26). It ought to be mentioned that the BF risk, along with the MTBBF, depends on the applied EGM, and that a Brown and Whitehead model was used in this paper. The expected number of direct lightning strikes to the TL tower accounts for the applied EGM, tower geometry, span length, distribution of strikes along the span, and ground flash density. Hence, Table 3 presents the BF risk and the MTBBF for the "Variant A" of the transmission line, obtained for several dif-  It should be mentioned that the concrete values presented in Table 3 can change slightly between different runs, due to the statistical nature of the indicator. Also, further deviation is to be expected with applications of different possible EGMs. However, all these are rather small and do not influence the final conclusions. On the other hand, MTBBF will scale in the inverse proportion to the ground flash density (see Table 3). At the same time, lightning current statistical parameters, including correlation, exert important influence on the BF risk, as noticed in [27]. The "tolerable level" of the MTBBF can be determined either by examining the BF influence on the substation insulation coordination [1,2,32], or by means of the full financial risk analysis.

Conclusion
This paper proposed a novel indicator of transmission line tower's BF performance, which is risk-based and utilizes entire probability distributions-of incident lightning currents and generated BF currents-instead of working with their point estimates. The shield wire(s) incident currents probability distribution is obtained from the downward negative (ground level) lightning statistical distribution and EGM application to the particular tower geometry. The BF currents statistical distribution takes full account of the TL model particularities (featured within the EMTP computational framework) and statistical correlation between lightning current parameters.
Any TL model specificities, that are introduced within the EMTP framework, are directly manifested in the CLP curve and carried-on forward to the BF risk assessment. Also, TL tower geometry and particular EGM variant directly influence the probability distribution of shield wire(s) incident lightning currents, which in-turn implicate the BF risk estimation. Ground level lightning-current statistical parameters provide only a starting point from which the particular tower's nonstandard bivariate probability distribution is constructed-by means of the Gaussian Copula approach-and from which subsequent BF currents probability distribution is filtered-out (using the CLP from that concrete tower). However, statistical parameters used for the ground level lightning distribution, including correlation, have important influence on the subsequent BF currents distribution.
This novel indicator, being also a single number, can be used to create a rang list of TL towers and identify certain towers with a high risk of BF occurrence. These would then become prime candidates for applying dedicated measures for decreasing their BF risk, such as lowering their impulse grounding impedance or installing surge arresters (i.e. TLA application). This could be of importance for the towers emanating from the substation, in relation to the substation lightning insulation coordination design, or it could be of use for identifying the socalled "rogue" towers and/or selecting appropriate tower candidates for the TLA installation.
Proposed risk-based BF indicator can be easily related to the price of that risk and the associated investment costs, thus enabling the cost-effective optimisation solutions to the mentioned problems of insulation coordination design and TLA applications. A BF performance indicator based on the risk analysis is deemed appropriate in these circumstances, considering the fact that the investment in surge arresters and related protective measures is commonly perceived as buying insurance.
A further advantage of the proposed indicator stems from the fact that it is obtained from systematic EMTP simulations (i.e. using CLP curves) and not from the Monte-Carlo approach which has been traditionally used for computing the BF probability. This means that it uses far less computational time than would be needed for the appropriate Monte-Carlo based approach. This is particularly important when large number of scenarios are considered or when individual towers are studied one-by-one (which is the case in station insulation coordination studies and in TLA application studies).