Evaluation of the Direct-Stroke Shielding Performance of Substations Through Stochastic Modeling of Lightning Attachment

This article introduces a methodology for assessing the direct-stroke shielding performance of high-voltage substations. The proposed methodology considers both negative and positive lightning polarity and can be employed by any lightning attachment model that follows a probabilistic approach. A stochastic lightning attachment model is used, which considers the lightning discharges probabilistic progression. Stochastic model results on lightning discharge propagation, lightning attachment, and induced electric fields are in satisfactory agreement with field observations; comparisons are also made and discussed with respect to other well-established lightning attachment models. An application of the proposed methodology is performed to estimate the lightning incidence rate and shielding failure flashover rate of a 69 kV substation; the effect of lightning activity at the substation location is elucidated.

Abstract-This article introduces a methodology for assessing the direct-stroke shielding performance of high-voltage substations.The proposed methodology considers both negative and positive lightning polarity and can be employed by any lightning attachment model that follows a probabilistic approach.A stochastic lightning attachment model is used, which considers the lightning discharges probabilistic progression.Stochastic model results on lightning discharge propagation, lightning attachment, and induced electric fields are in satisfactory agreement with field observations; comparisons are also made and discussed with respect to other well-established lightning attachment models.An application of the proposed methodology is performed to estimate the lightning incidence rate and shielding failure flashover rate of a 69 kV substation; the effect of lightning activity at the substation location is elucidated.

I. INTRODUCTION
L IGHTNING strikes impose a high risk of power system outages [1], [2], [3].Substations, critical part of a power system, shall operate with high reliability and experience a practically zero outage rate [4].Thus, lightning-related studies regarding the direct-stroke shielding system design to substations and analysis of the arising surge overvoltages at substation have been a subject of extensive research [5], [6], [7], [8], [9], [10].
The complexity of shielding analysis and difficulties in formulating lightning performance assessment to substations into a generalized methodology are mainly associated with: 1) The intricacy of the actual 3-D geometry and multiple basic insulation levels and insulation types involved; 2) the lack of up-to-date and reliable field data on substation failures and lightning parameters with seasonal and geographical variation; 3) the inherent stochastic factors characterizing lightning discharge propagation.Thus, this work contributes to the development of a methodology for assessing the direct-stroke shielding performance of substations through stochastic modeling of lightning attachment.
The proposed approach considers both negative and positive lightning polarity and can be employed by any lightning attachment model that follows a probabilistic approach.In this study, a stochastic (fractal-based) lightning attachment model is employed, which considers the physical processes involved in lightning attachment.The adopted stochastic model can consider: 1) The actual 3-D geometry of substations and substation buses operating voltage; 2) the stochastic progression of lightning discharges in a 3-D domain; 3) lightning incidence to substations against direct lightning strikes of positive polarity; 4) the dynamic variation of downward to upward leader velocities supported by recent field observations.Thus, the effects of lightning polarity, object height, and statistical dispersion of lightning interception parameters are elucidated.The stochastic model validity is demonstrated against field observations on lightning discharge propagation, lightning attachment, and induced electric fields.Comparisons are also made and discussed with respect to other well-established lightning attachment models from literature.
The proposed methodology is applied to a practical 69 kV substation scheme to estimate its lightning incidence rate and shielding failure flashover rate; the effects of ratio of negative to positive lightning flashes, peak current distribution parameters, and ground flash density are quantified and discussed stressing the need for acquisition of up-to-date lightning data based on lightning location systems (LLSs) [34] and lightning strike measurements [35].

II. METHODOLOGY FOR ASSESSING THE LIGHTNING PERFORMANCE OF SUBSTATIONS
In this section, the methodology for assessing the direct-stroke shielding performance of substations is formulated; it is noted that the proposed methodology can be applied by employing any lightning attachment model that follows a probabilistic approach on lightning incidence estimation.The direct lightning incidence rate, N D (flashes/yr), and the shielding failure flashover rate, SFFOR (flashovers/yr), on an annual basis to a substation can be given as where 1) N G (flashes/km 2 /yr) is the ground flash density at the substation location.2) A (km 2 ) is the lightning collection area of the substation.3) A(I) of complex structures such as substations can be calculated by employing the following procedure in any CAD software; in this work, the AutoCAD software was employed.a) Lightning interception zones with radius equal to the interception radius R are drawn from substation equipment and air terminals for each lightning peak current I.
The interception radius is defined as the maximum lateral distance between the air terminal and the lightning interception point [36] and is commonly computed as a function of I and object height, H; R values in this work were obtained through stochastic lightning attachment model simulations (see Section V-A).b) The lightning collection area, A, is equal to the ground surface area defined by the contour of the union of the drawn circles of step 1. c) Collection area values are estimated as a function of the lightning peak current; A(I) in ( 1) and ( 2) is the fitting curve of the obtained values per I. 4) f(I) is the lightning peak current distribution with a probability density function, which is log-normally distributed as [37] f with Ī and σ ln denoting the median value and the standard deviation of the natural logarithm, respectively.5) I S (kA) is the minimum (critical) current causing flashover of substation insulation, which can be roughly estimated as [11], [38] where 1.1 is a factor to account for the reduction of stroke current terminating on a conductor as compared to zero impedance earth [11], BIL (kV) is the basic insulation level that substation buses insulators are designed to withstand, obtained based on their electrical characteristics, CFO (kV) is the negative polarity critical flashover voltage of the insulation (0.94 factor is used to represent a 6% reduction in CFO values in order to produce a withstand level roughly equivalent to the BIL rating), and Z S (Ω) is the surge impedance of substation buses (1/2 is used to account for the assumed equal split of the stroke current), which can be calculated according to the procedure described in (11, Annex C) based on the conductor's height and radius and the equivalent corona radius.Accurate estimation of I S requires electromagnetic transient (EMT) simulations due to the nonstandard waveform of lightning overvoltage surges [39] and reflections at substation equipment of high surge impedance [11].6) p SF (p.u.) is the probability of shielding failure computed from stochastic model simulations; p SF is defined as the ratio of direct lightning strikes to substation buses to the total number of direct lightning strikes to the substation's protected and protective equipment.p SF is obtained through the proposed stochastic model results (see Section V-A) or employing another probabilistic lightning attachment model.Equations ( 1) and (2) can be generalized for both lightning polarities, considering the ratio, r (p.u.), of the number of negative to total lightning flashes in a region.Thus, the total lightning incidence, N Dt , and shielding failure flashover, SFFOR t , rates are defined as where N D− , SFFOR − and N D+ , SFFOR + are the corresponding computed rates for negative and positive lightning.A typical value of r = 0.9 [40] is commonly adopted; nevertheless, r may vary depending on seasonal and geographical characteristics.N D and SFFOR values depend significantly on the lightning activity characteristics at the substation location [N G , f(I), r] as well as on the shielding effectiveness of the adopted airtermination system (p SF ).Thus, the need for acquisition of accurate lightning data based on reliable LLSs is important.A flowchart of the methodology is presented in Fig. 1.The proposed methodology, following a stochastic approach, may assist in overcoming simplifications of the methods adopted by the IEEE Std.998 [11] associated with: 1) Estimation of the substation lightning collection area employing geometrical approaches; 2) estimation of the substation lightning incidence and flashover rate neglecting possible occurrence of positive lightning flashes; 3) determination of the shielding failure rate of substations based on semiempirical purely deterministic methodologies.

III. STOCHASTIC MODELING OF LIGHTNING ATTACHMENT
A description of the employed stochastic model algorithm, developed in MATLAB, is presented in Section III-A.Model parameters for simulation of lightning discharges of negative and positive polarity are given in Section III-B.

A. Model Description
The 3-D simulation domain is initially discretized employing a uniform grid spacing.The electric potential, V, is computed at all points of the simulation domain (i, j, k), at each algorithm step, solving the Laplace equation.For that purpose, a computational electromagnetics method is employed combining the finite-difference method (FDM) solved by the successive over-relaxation (SOR) iterative procedure [41], [42] V where ω is the over-relaxation parameter (1 ≤ ω < 2, 1.85 in this study).ω can be determined with the aid of theoretical equations [41] or through dedicated studies.Appropriate boundary conditions were selected comprising the following: 1) Dirichlet type (V = const.)for the cloud, ground, airtermination system, and substation equipment; 2) Neumann type (∂V /∂x = 0, ∂V /∂y = 0) for the lateral boundaries.Electric potential calculation comprises a significant task in the employed model as it drives all the physical calculations involved; thus, accurate and fast computations are important.Accurate estimation of the electric potential requires the use of a very fine grid spacing for the domain discretization as well as a low converging tolerance in FDM and SOR; this may lead to an excessively large number of grid points, and thus, computational cost as also reported in [20] and [22].Hence, in this work, code optimization techniques, developed in [42], were employed to significantly reduce the total simulation time (approximately 1-3 hours per simulation run) retaining at the same time the accuracy in electric field computations.
At each simulation step, a single node L is added to the lightning discharge path provided that the average electric field with the parent node L, is higher than a critical electric field, E cr .The next leader discharge node is selected from the following: where η is a model parameter affecting the random progression of the simulated lightning discharges, N is the number of the lightning discharge nodes (n = 1 …N), and is the number of the potential leader nodes that satisfy E LL ≥ E cr .The denominator in (8) refers to the summation of all possible candidate points M n from all discharge nodes N and as a result normalizes the propagation probability distribution.Thus, ( 8) confers a probabilistic character on the proposed model to account for the several random factors during natural lightning discharge propagation.point is selected (for example point L 24 in Fig. 2(a), equal to the new leader node L 3 at simulation step y+1, Fig. 2(b); the same algorithm steps are then continued for N = 3 to select the next leader discharge point.It is noted that an eligible point L may be associated with two different probabilities from two neighboring leader points (as for example nodes L 11 and L 21 from L 1 and L 2 , respectively) allowing, thus, the creation of leader branches.The same procedure is also followed for the case of a 3-D simulation domain and for upward leaders' propagation, as depicted in Fig. 2(c) and (d), illustrating a downward and an upward leader at two successive simulation steps (y and y+1) and the corresponding neighboring and eligible nodes for leader propagation.
After leader node selection, the potential at the newly added node is updated considering an average discharge gradient, E ch , based on the following formula [43]: where V 0 is the leader potential at the starting point of the simulation, s is the sign of the charge carried by the leader ("−" for leaders of negative polarity and "+" for leaders of positive polarity), and d total is the total lightning discharge path computed as a sum of all the discrete steps from the starting point till the newly added node L .Then, the electric potentials at the remaining points of the domain are recalculated; thus, the competition among upward leaders and their interaction with the DL is considered.These algorithm steps apply for both the downward and each one of the multiple upward leaders.
The proposed stochastic model belongs to the single-element category, i.e., only one point is added to the lightning discharge path per algorithm step; this approach was adopted to consider the interaction of all leader discharge branches based on the dynamic recalculation of the electric potential at all points of the simulation domain after new point addition.

B. Parameter Selection
The employed stochastic model parameters are listed in Table I for both lightning polarities.These parameters are based on physical grounds of lightning discharges as well as on field measurements and observations.The DL tip potential was selected according to the formula proposed by E. M. Bazelyan [44]; as an approximation, this expression can also be used for positive lightning [45].The potential of the DL tips in the simulated tortuous and branched channel was considered by means of the average discharge gradient, E ch , as Bazelyan's formula refers to a straight leader channel.Considering the selected value of E ch , tips that deviate from the primary DL channel are characterized by a reduced leader tip potential.Thus, the differences among the charge distribution of the main channel and branches have been implicitly considered employing a voltage drop along the leader channel.The average critical electric field, E cr (kV/m), was selected as the electric field at the streamer zone of the downward and upward leader tips, respectively, required for stable streamer propagation.The average discharge gradient was chosen E ch = 10 kV/m, as this value is considered representative of a typical leader channel internal field [48].
The upward leader inception criterion proposed in [46] and [47] was selected as it can integrate the atmospheric conditions and altitude variation into the proposed expression and is validated against field observations [51].Actually, this criterion dictates the uniform grid spacing used for the 3-D domain Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.
discretization.Values of 1.4 and 2.9 m were selected for negative [29] and positive [30] lightning, respectively, in line with the critical streamer length for stable upward leader inception given by Petrov and Waters [46], [47]; it is noted that as an approximation the same spatial step was adopted for the DL and upward leaders although the descending stepped leader propagates in significantly larger discrete steps in natural lightning [47].These values are small compared to the total kilometer-range lightning discharge path, thus, the employed code acceleration techniques [42] are vital for the required multiple stochastic simulation runs.
The velocities of the downward and upward leaders were considered in the employed model based on their corresponding ratio; this ratio is a parameter that considerably affects lightning incidence results.A velocity ratio was adopted in this work and not absolute velocity values as an engineering approach; it is noted that most lightning attachment models from literature consider a constant speed ratio.However, recent high-speed video records [52], [53], [54], [55], [56] have shown a dynamic variation of downward and upward leader velocities; actually, velocity ratios even lower than unity were recorded [52], [54], [55].In this work, a dynamic velocity ratio formula is employed based on the expression proposed in [12] and modified in [49] to account for recent field observations on downward and upward leader velocity ratios (allowing also for values lower than unity) and consider the effect of atmospheric conditions.As an approximation, the same dynamic velocity ratio formula was employed for both negative and positive lightning as shown in Table I.Nevertheless, it is important to note that the adopted velocity ratio formula can consider the effect of lightning polarity.This is due to the fact that in case of negative lightning, upward leaders incept earlier (at higher altitudes of the DL) when compared to positive lightning [47], due to the larger upward leader inception field required in the latter case; hence, the velocity ratio variation from upward leader inception up to the final jump phase is higher for negative than positive lightning, since for the former ΔE (see Table I) takes generally lower values.This observation is in line with the modeling analysis of Dellera and Garbagnati [12]; however, definitely more and reliable field data on positive lightning attachment are needed emphasized on records of DL growth through advanced monitoring systems [57].

IV. MODEL RESULTS AND COMPARISON WITH FIELD DATA
In this section, the employed stochastic lightning attachment model is validated against field observations and data.In addition, results of the stochastic model regarding basic shielding design parameters (striking distance and interception radius) are compared with well-established lightning attachment models from literature.

A. Field Observations of Lightning Attachment
The stochastic model predicts predominantly lightning attachment through a tip-to-tip connection of downward and upward leaders as shown in Fig. 3(a) and (b).It must be noted that there are simulation results indicating cases, where the DL intercepts to the lateral surface of the upward connecting leader (UCL) Fig. 3. Interception cases for negative lightning.Tip-to-tip interception: (a) simulation result and (b) field observation; adapted from [24].Tip-to-side interception: (c) simulation result and (d) field observation; adapted from [24].Side flash: (e) Simulation result and (f) field observation; adapted from [32].below its highest tip [see Fig. 3(c)], as recently observed in natural lightning discharges [24], [58], [59] [see Fig. 3(d)] as well as scale-model experiments [60]; stochastic model simulations can also predict attachment cases, although rarely recorded in field observations, where the UCL intercepts to the lateral surface of the DL above its tip [58], as in the lightning attachment simulation model of [61].In addition, it can reproduce cases of side flashes [see Fig. 3(e)] as the ones recorded in recent field observations [32] [see Fig. 3(f)].

B. Fractal Dimension of Lightning Discharges
The fractal dimension is a parameter that is used to quantify the complexity of a fractal pattern characterized by selfsimilarity.The fractal dimension of the simulated lightning discharges was estimated based on the box-counting and sandbox algorithms [62], [63].The employed fractal estimation methods were used to compute the 2-D fractal dimension of the simulated 3-D lightning discharges as the mean value from different projections [63].Thus, comparison with photographic records of natural lightning is enabled.Fractal growth of the simulated lightning discharges depends strongly on the propagation parameter, η, in (8).Thus, to allow for proper parameter selection, Fig. 4 depicts the variation of fractal dimension with the propagation parameter; based on Fig. 4, η was selected equal to 8 (dotted frame) as for this value the associated mean fractal dimension and its statistical dispersion lie within the limits (1.1-1.4)[64] reported in field observations of natural lightning.

C. Lightning Interception Points
Stochastic model predicts different lightning interception points in every simulation run.Interception points for the case of a single lightning mast create an interception area above the air-terminal that increases with increasing lightning peak current and has a statistical notion [65], as shown in Fig. 5. Dots in Fig. 5 depict the various interception points for the case of a 60 m tower for positive lightning polarity; there are few cases that interception occurs lower than the air-terminal tip, in line with field observations (side flashes) [24], [32].It is important to note that field data from Nakamura et al. [66] (red crosses) derived from rocket-triggered lightning experiments to a 60 m transmission tower are located within the "cloud" of interception points (blue dots) obtained from stochastic model simulation results.

D. Striking Distance Field Data
In this work, the following two interpretations of the striking distance parameter were taken into account: 1) S inc , which represents the distance of the DL and the airterminal at the moment when the UCL begins [as depicted in Fig. 6(a)]; 2) S int , which is the distance between the point where the DL and UCL intercept and the air-terminal [as shown in Fig. 6(b)].Fig. 7 depicts stochastic model results on striking distance together with field data on natural lightning (first strokes) to 60 m towers from Eriksson [67] and reproduced by Uman [68] and Visacro et al. [69], referring to inception definition; this figure also shows striking distance data of Saba et al. [53] derived from natural lightning (subsequent strokes) to a 52 m building derived based on both striking distance definitions.It should be mentioned that although striking distance field data are associated with a huge dispersion even for similar lightning peak currents, field data lie within the range of striking distance values obtained from stochastic modeling; error bars denote the range of striking distance values based on the total number of simulations.Stochastic model results of Fig. 7 may explain the weak relationship between striking distance and lightning return-stroke peak current also found from field observations of Mazur and Ruhnke [70] indicating that striking distance cannot be described by purely deterministic expressions; this nondeterministic striking distance behavior may also be associated with the probabilistic variation of the downward leader velocity as shown in [14].

E. Electric Field at Ground Level
Stochastic model results regarding the electric field amplitude at ground level prior to return-stroke are depicted in Fig. 8 for different distances from the lightning termination point against field observations from Jerauld et al. [71] for a 51.9 kA negative CG lightning.Error bars denote the minimum and maximum values derived from stochastic model simulation runs since the electric field at the ground level depends on the stochastic variation of the downward leader path prior to lightning termination to the ground; the derived statistical dispersion of the electric field amplitude at the ground level was also shown in another simulation study [21].Thus, points on the ground level may be associated with higher electric field amplitude values if the lightning discharge propagates predominantly above these points during the final steps prior to lightning termination to the ground.The inset figure in Fig. 8 (enclosed blue dotted frame) depicts the time variation of the electric field, measured at a distance 260 m from the lightning termination point.A satisfactory agreement was found for the average values of stochastic model results (see Fig. 8) with respect to field data, as well as for the electric field time variation at 260 m (inset graph); this enhances the validity of the employed stochastic model in terms of the adopted electric field computational method and model parameter selection.

F. Comparison With Other Lightning Attachment Models
Striking distance and interception radius values (see Fig. 6) obtained through stochastic modeling are compared with other lightning attachment models from literature.Fig. 9 shows the striking distance, for a mast height of ∼20 m, based on the EGM of the IEEE Std.998 [11] as well as other well-established lightning attachment models.The EGM employing a solely  current-dependent striking distance expression may underestimate striking distance values especially for higher objects.Nevertheless, results based on the considered models lie in the statistical dispersion obtained from the stochastic model; an excellent agreement was found between stochastic model results and SLIM [73].
The interception radius, R, represents the largest lateral distance from the air terminals at which a lightning strike can occur [as illustrated in Fig. 6(b)].A comparison of R obtained through the stochastic model and other models from the literature is displayed in Fig. 10; a good agreement with Rizk's model [15] exists.From Fig. 10, it can be observed that results on interception radius obtained through stochastic modeling are lower than the widely adopted EGM of IEEE Std.998 in the low lightning peak current region (5-10 kA), and then, become higher for lightning peak currents in the region of 20-30 kA.A possible explanation for this trend is that stochastic modeling considering the inception and stochastic propagation of upward leaders, the object height, and the operating voltage of substation Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.buses predicts generally higher interception values than the EGM of IEEE Std.998 [11]; the latter is purely deterministic and neglects the aforementioned parameters.Nevertheless, results from these two models are in good agreement in the low lightning peak current region (5-10 kA) where the shielding system is designed for the majority of the substations; in this region, the effect of upward leaders' propagation is minimal and does not significantly influence the lightning attachment phenomenon.

V. LIGHTNING PERFORMANCE ASSESSMENT OF A 69 KV SUBSTATION
The lightning performance of a 69 kV substation against direct lightning strikes has been evaluated based on the introduced methodology; it comprises a typical substation of sectionalized bus arrangement, commonly found in practical substation schemes, presented in IEEE Std.998 [11].

A. Stochastic Modeling of Lightning Attachment and Results
The 69 kV substation under study comprises three buses of different heights.As a case study, an air-termination system designed based on the fixed angle method (45°protection angle) was adopted; the relative position and dimensions of the lightning masts were employed based on Annex B of [11].Modeling of the substation was performed in MATLAB as points of fixed electric potential considering the substation layout and dimensions as well as the grid spacing employed for 3-D domain discretization (see Table I).A top view of the substation layout and the MATLAB model employed for simulations are depicted in Fig. 11(a) and (b), respectively.Fig. 11(b) also shows eligible upward leader inception points, from the substation layout, associated with an elevated electric field; points at the lateral surfaces of masts and substation components were also considered allowing, thus, for side flashes [see Fig. 3(e)].
Fig. 12 shows the probability of shielding failure, p SF (p.u.) of the 69 kV substation for negative and positive lightning.Lightning discharges associated with low prospective return stroke peak currents are less likely to initiate an UCL from an air-terminal (masts and/ or catenary wires) due to the reduced electric field induced to its region; thus, the downward leader need to approach closely to the direct-stroke shielding system of air terminals in order the upward leaders to be able to satisfy the inception criterion, hence resulting in potentially larger number of shielding failures.
It is noteworthy that shielding failure cases do exist for currents larger than the maximum shielding failure current, I MSF , estimated as 7.2 kA according to the EGM of IEEE Std.998 [75].Practically all shielding failure currents would cause flashover since the critical current causing flashover of substation insulation, I S , is 2.6 kA according to (4) (BIL = 350 kV, Z S = 300 Ω [11]).p SF was found significantly lower for negative lightning since stable positive upward leaders can incept earlier from masts when an average electric field is satisfied equal to the electric field required for positive streamer propagation, 500 kV/m, within a critical positive streamer length according to [46].For the case of positive lightning, inception of stable negative upward leaders can be achieved when an average electric field is satisfied within a critical negative streamer length, which is more than double the corresponding one of positive upward leaders, with typical electric field values required for negative streamer propagation ranging between 1000 and 1600 kV/m [76].As a result, the positive downward leader needs to approach very closely to a mast to enhance sufficiently the electric field in its vicinity, leading, thus, to potentially higher probability of shielding failure.Thus, for the case of positive lightning, masts initiating late and short upward leaders exhibit shorter striking distance and interception radius when compared to negative lightning [30], [47].
In addition, Fig. 12 indicates that I MSF has a probabilistic notion and shall be treated as a statistical parameter, rather than determining a unique value corresponding to zero shielding failure probability based on a deterministic approach.Thus, the employed stochastic model considering the effect of lightning polarity, object height, and random characteristics associated with lightning attachment, could predict cases of shielding failures documented by utilities in IEEE surveys [11], [77]; these refer to recorded unexpected cases in substations with a shielding system designed based on IEEE Std.998 procedures.Simulation results of shielding failures based on stochastic modeling are shown in Fig. 13.II and a ratio of negative to total lightning flashes r between 0.5 and 1; N G = 10 flashes/km 2 /yr.

B. Shielding Performance of the 69 kV Substation
The direct lightning incidence rate, N D (flashes/yr), and shielding failure flashover rate, SFFOR (flashovers/yr), on an annual basis to the 69 kV substation were computed [see ( 1)-( 6)]; the effect of lightning activity at the substation location was also investigated.Recorded lightning peak current distributions and lightning activity parameters for different regions around the world are listed in Table II.
Fig. 14 shows N D values for the distributions of Table II and a ratio of negative to total lightning flashes, r, between 0.5 and 1. N D for negative lightning (r = 1) varies considerably, up to ∼2.7 times, among the distributions; the highest incidence rate is estimated for distribution no.4, which exhibits the highest median value, commonly found in tropical regions.For the case of higher percentage of positive lightning (1−r), thus for decreasing r, lightning incidence rate decreases due to the reduced interception radius associated with positive lightning [30], [47].II and a ratio of negative to total lightning flashes r between 0.5 and 1; N G = 10 flashes/km 2 /yr.Dotted horizontal line denotes a typical acceptable SFFOR due to direct lightning strikes (0.05% [11]).II).Dotted horizontal line denotes a typical acceptable SFFOR due to direct lightning strikes (0.05% [11]).Fig. 15 demonstrates the SFFOR of the 69 kV substation for the distributions of Table II.SFFOR for negative lightning (r = 1) varies up to ∼70 times, among the distributions; the highest flashover rate is estimated for distribution no. 1, which exhibits the lowest median value, commonly found in temperate zones.Thus, although a larger number of direct lightning flashes (N D ) is expected for substations in tropical regions, the SFFOR may be bigger in temperate zones, with the same N G , since the probability of low-intensity currents that may cause shielding failures is higher.For the case of the higher percentage of positive lightning, thus for decreasing r, SFFOR increases drastically due to the high shielding failure probability (p SF ) associated with positive lightning (see Fig. 12).Besides, the importance of accurate knowledge of the lightning parameters at the substation location is stressed by the fact that the same substation configuration for a ratio of negative to total lightning flashes r = 0.5 exhibits an acceptable SFFOR of 0.05% (1 per 2000 years [11]) Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.
for distributions no. 3 and no. 4, whereas SFFOR is higher for the other regions with lightning peak current distributions of a lower median value (see Table II).
To further demonstrate the prominent effect of lightning activity parameters on the substation lightning performance, Fig. 16 shows N D [see Fig. 16(a)] and SFFOR [see Fig. 16(b)] values for the 69 kV substation based on lightning activity data from Table II (typical values per region).It is noted that the 69 kV substation when located at a region such as Colombia with intense lightning activity (N G ≥ 20 flashes/km 2 /yr and r ≤ 0.8) would exhibit an unacceptable SFFOR (> 0.05%) in contrast to other locations, which exhibit a practically perfect direct-stroke shielding performance (< 1 flashover per 10 000 years).
Nevertheless, lightning activity parameters may vary significantly even for the same region depending on seasonal and geographical factors; thus, up-to-date lightning activity studies are of importance.It is also noteworthy that the total lightningrelated failure rate of substations is not only attributed to shielding failures due to direct lightning strikes to the substation.Incoming lightning overvoltages from the interconnected overhead lines are the main cause of lightning-related failures, and thus, may substantially affect the performance of substations against lightning [82], [83].

VI. CONCLUSION
In this work, a methodology for assessing the direct-stroke shielding performance of high-voltage substations has been introduced; the methodology can be utilized employing any lightning attachment model following a probabilistic approach.A stochastic lightning attachment model has been employed, considering the probabilistic progression of lightning discharges in a 3-D domain as well as lightning polarity.The stochastic model has been validated against field observations on lightning discharge propagation, lightning attachment, as well as induced electric fields.An application to a practical 69 kV substation configuration has been made.Results of this study have revealed the following: 1) Striking distance shall not be treated as a deterministic but as a statistical quantity depending on object height and lightning polarity.Thus, the statistical dispersion of striking distance field records can be explained.2) Shielding failure flashover events are rare but of high impact; the employed model, following a stochastic approach for lightning discharge propagation, can predict shielding failure cases (rare events) as documented in surveys conducted by the IEEE that cannot be explained by the conventional shielding analysis.Reliable field data regarding shielding failures to substations are required for further validation of lightning attachment models and shielding design methods.3) SFFOR of substations of the same configuration and shielding design exhibit significant variation when located in different regions due to the various lightning activity characteristics.Thus, the use of a universal distribution as proposed by international standards may result in reduced accuracy estimates of the SFFOR and the implementation of a fixed shielding design system to a substation does not necessarily guarantee an acceptable shielding performance at a global level.This stresses the need for the acquisition of reliable and up-to-date lightning data considering the seasonal and geographical variation of lightning activity parameters.4) Substations in tropical regions with commonly intense lightning activity and a high percentage of positive lightning need an advanced shielding system to achieve an acceptable SFFOR.Substations located in temperate zones, commonly associated with a relatively low median value, may experience a high SFFOR due to the higher probability of low-intensity currents that cause shielding failures.5) Stochastic analysis on direct-stroke shielding of substations may be adopted by international standards especially for mission-critical applications, such as substations feeding hospitals, data centers, and government and military infrastructure.This is more so considering that the available computing power will be significantly expanded within the next years allowing for more sophisticated modeling of lightning attachment through advanced computational techniques.

APPENDIX
In this Section, the methodology for the estimation of the coefficient u in the critical propagation electric field formula of Table I is presented.The methodology considers the dependence of the negative and positive leaders on altitude and atmospheric conditions.The following analysis is based on a combination of theoretical analysis and experimental or/ and simulation results from literature studies [46], [50], [76].
Based on [46] the dependence of pressure, p, temperature, T, and absolute humidity, γ, on altitude, z, can be given as where p(0) = 760 mmHg is the pressure at ground level and z 0 a constant equal to 8 km by assuming a linear temperature change where T(0) = 293 K is the temperature at ground level and k a constant equal to 6 K km −1 , where γ(0) = 11 gm -3 is the absolute air humidity at ground level and z H a constant equal to 3 km.In addition, the electric field required for stable positive streamer propagation can be given as [50] E st+ (δ, γ) = 464δ • 1 + 1.37 100 where δ is the relative air density as a function of pressure and temperature and standard sea-level conditions given as [46] δ = p 760 Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.By substituting (A1)-(A3) and (A5) into (A4), the coefficient u can be roughly estimated as u = 7.5 km and E cr0 = 464 kV/m, for streamers of positive polarity as can be deduced from the parameters of the fitting curve of Fig. 17  Following the same analysis also for streamers of negative polarity and by considering that the dependence of the average stable negative streamer propagation electric field on relative air density is given as [76] E st− (δ) = 1100 • δ (A6) the coefficient u can be roughly estimated as u = 10 km and E cr0 = 1100 kV/m, for negative streamer polarity as can be deduced from the parameters of the fitting curve of Fig. 17 (red curve).

ACKNOWLEDGMENT
Results presented in this work have been produced using the Aristotle University of Thessaloniki High-Performance Computing Infrastructure.The authors would like to thank Prof. N. Kantartzis for the fruitful discussions and assistance with the employed computational electromagnetics methods.The publication of the article in OA mode was financially supported by HEAL-Link.

Fig. 2 (
a) and (b) depicts a 2-D schematic representation of the next point selection procedure for the downward leader (DL) at two successive simulation steps y and y+1 based on the cumulative probability density function P derived from (8).Based on Fig. 2(a), N in (8) is equal to 2 (points L 1 and L 2 ), thus, the new DL discharge point is selected among the eligible points (empty dots) surrounding both L 1 and L 2 ; hence, for n = 1 (point L 1 ), eligible leader nodes are L 11 and L 12 (M 1 = 2), whereas for n = 2 (point L 2 ), eligible leader nodes are L 21 , L 22 , L 23 , L 24 , and L 25 (M 2 = 5).Then, a cumulative probability density function P is created, as shown in Fig. 2(b), where P 1 = p 11 , P 2 = p 11 + p 12 … P 7 = 1, and based on a random number m, the new leader discharge

Fig. 2 .
Fig. 2. (a) and (b) 2-D schematic representation of the next point selection procedure for the DL at two successive simulation steps y and y +1 based on the cumulative probability density function P derived from (8).(c) Corresponding part of a 3-D simulation domain.Algorithm step y; L 1d and L 1u (full dots) and L 1d and L 1u (empty dots) represent the current and possible next nodes for discharge propagation of the downward (subscript d) and upward (subscript u) leaders, respectively.(d) Corresponding part of a 3-D simulation domain.Algorithm step y+1; discharge evolution to leader nodes L 2d and L 2u .

Fig. 4 .
Fig. 4. 2-D fractal dimension, D f , as a function of the propagation parameter η; negative lightning, cloud-to-ground (CG).Solid and dashed vertical bars denote the maximum and minimum limits for the box-counting and sandbox method, respectively.Horizontal lines indicate the recorded fractal dimension values (1.1-1.4) based on field observations [64]; 100 simulations runs per η.

Fig. 5 .
Fig. 5. Lightning interception points of a 60 m tower for positive lightning.Blue dots denote stochastic model simulation results and red crosses field data from Nakamura et al. [66] derived from rocket-triggered lightning to a 60 m transmission tower (1000 runs per lightning current level).

Fig. 6 .
Fig. 6.Adopted definitions of striking distance.(a) Inception, S inc .(b) Interception, S int .R is the interception radius and H the object height.

Fig. 7 .
Fig. 7. Striking distance results based on the stochastic model and field data; negative lightning, 500 runs per current level.Empty points: Simulation results, and full points: Field data.Results with green color refer to the inception definition.

Fig. 8 .
Fig. 8. Electric field amplitude at ground level prior to return stroke; comparison of stochastic model results (100 simulation runs) with field data of Jerauld et al. [71] for different distances from the lightning termination point.Dots for the stochastic model represent average values.Inset figure (enclosed blue dotted frame) depicts the variation of the electric field with time.

Fig. 9 .
Fig. 9. Striking distance (inception definition) according to the stochastic and other lightning attachment models from literature; mast height 21.3 m, negative lightning, and 1000 runs per current level.With green solid line the fitting curve of stochastic model results (triangles); and vertical error bars denote the maximum and minimum striking distance values obtained through the stochastic model.

Fig. 10 .
Fig. 10.Maximum interception radius, R, according to the stochastic and other lightning attachment models from literature; mast height 21.3 m, negative lightning, and 1000 runs per current level.

Fig. 11 .
Fig. 11.(a) Top view of the 69 kV substation under study [11].(b) 3-D model of the substation in MATLAB; squares denoted upward leader inception points.

Fig. 14 .
Fig. 14.N D (flashes/yr) for the distributions of TableIIand a ratio of negative to total lightning flashes r between 0.5 and 1; N G = 10 flashes/km 2 /yr.

Fig. 17 .
Fig. 17.Dependence of the critical propagation electric field of negative (red curve) and positive (blue curve) streamers on the altitude above ground level.
Evaluation of the Direct-Stroke Shielding Performance of Substations Through Stochastic Modeling of Lightning Attachment Alexios I. Ioannidis , Student Member, IEEE, and Thomas E. Tsovilis , Senior Member, IEEE

TABLE I EMPLOYED
PARAMETERS FOR DOWNWARD AND UPWARD LEADERS

TABLE II LIGHTNING
ACTIVITY PARAMETERS FOR DIFFERENT REGIONS