Quantifying the seismic risk for electric power distribution systems

Abstract Electric power distribution systems are generally more prone to disruption from natural hazards than transmission systems due to their often less redundant circuit structures. However, seismic risk analysis for distribution systems is rare compared to the rich body of literature focussing on transmission systems. This paper proposes a seismic risk assessment framework for electric power distribution systems considering both the network topology and the functional vulnerability of distribution substations. Implicit Z-bus method is applied to solve distribution system power flow and evaluate system serviceability. Monte Carlo simulation is applied to obtain probabilities of the scale of unserved loads resulting from disconnection and abnormal voltage condition. The seismic risk is jointly quantified using multiple risk metrics, and importance measures are used to determine criticality of substation components for prioritisation of seismic retrofit. The seismic risk assessment framework is applied to the CIGRE medium voltage distribution test network and two ground motion intensity scenarios – one for peak ground acceleration values based on a scenario earthquake and the other for uniformly distributed peak ground acceleration across the network. The framework allows the quantification of different network topologies and substation configurations. This enables network owners and operators to evaluate the seismic vulnerability of their substation configuration and network topology, identify potential bottlenecks of the systems and thus inform effective planning and risk-reduction investments.


Introduction
Electric power distribution systems deliver electricity from the transmission system to individual consumers. Transmission systems, as the backbone of national grids, are built to be strongly meshed for reliability purposes. In contrast, distribution systems are typically built as weakly meshed, radial networks, mainly due to economic constraints. Therefore distribution systems are intrinsically more vulnerable to earthquakes than transmission systems due to lower redundancy in both network topologies and substation circuits. The performance of distribution substations are crucial to the entire distribution network performance, and a single component failure could lead to complete disconnection of the down-stream system, which is not always the case in transmission systems. The functionality of electric power networks can be severely affected by extreme natural hazard events. Continuing power supply after an extreme event is crucial to emergency response and recovery as multiple utilities, such as telecommunication and water supply, rely on power supply for continued function.
Therefore utility managers are interested in estimation and minimization of the disruption on service level. According to a survey on causes of large blackouts in the U.S., earthquakes are shown to have the least frequency of occurrence but the most significant impact in terms of lost power (Hines et al., 2008). However, detailed electric power infrastructure impact assessment following such events are often excluded from current industrial reliability assessment practice. This paper presents a quantitative assessment framework for estimating the impact of earthquakes on electric power distribution systems, which can assist network owners and government officials to identify system safety bottlenecks and guide investments.
The electric power network of New Zealand has been exposed to significant earthquakes in recent years. The 22nd February 2011 Christchurch earthquake caused severe damage in four substations and extensive damage to buried cables in areas with severe soil liquefaction. The earthquake resulted in power outage to over 80% of the Christchurch area and 629 million customer minutes were lost (Kwasinski et al., 2014;Massie and Watson, 2011). In the 14th November 2016 Kaikōura earthquake, although there was no significant structural damage to the electric power network, a substation bus conductor dropped off as it reached its seismic movement limit leading to disconnection of the downstream distribution systems from the grid for approximately 8 hours (Liu et al., 2017). It has been observed that earthquake damage to transmission systems in New Zealand is often minor compared with that of distribution systems.
Under most circumstances, post-event transmission level repairs are completed within hours to days whereas complete distribution system repairs can range from weeks to months (Kwasinski, et al., 2014;Liu, et al., 2017;Massie and Watson, 2011).
Studies of natural hazard impacts on electric power networks can be performed either at component level or up to system level. Substation component level damage has been evaluated by a wide variety of methods (Anagnos, 1999;Paolacci et al., 2014;Zareei, Hosseini, & Ghafory-Ashtiany, 2016. However, power outages can still happen even if there is no extensive structural damage to the system. The component structural damage assessments are necessary but alone they are not sufficient for full power system functionality assessment. A seismic reliability evaluation method based on Monte Carlo Simulation (MCS) was proposed by Pires et al.(1996) for electric power transmission networks. Ground motion intensity was mapped to component availability through fragility functions, and transmission power flow was assessed to evaluate the overall network performance. Related work was conducted by Vanzi (1996), where substation internal circuits were specifically modelled. Shinozuka, et al. (2007) presents an MCS on earthquake damage of a transmission system as well as the recovery process, where simulation results were in good agreement with observed data.
A comparative study of five seismic performance evaluation models with varying complexities are reported by Cavalieri et al. (2014). The models were compared based on both system-level and component-level measures, suggesting that simple graph models can be used for retrofit prioritization, whereas the complex flow-based models are more appropriate when actual system performance is needed. Panteli, et al. (2017) considered temporal impacts of extreme weather events on a national transmission grid, where sequential MCS was applied to simulate the system functionality over time. Resilience metrics factoring in both functional and structural aspects were proposed to better characterize performance degradation due to increasing intensities of the extreme weather event. Despite the literature on individual component and transmission network performance, the performance of distribution networks following extreme hazard events have rarely been considered in their entirety. Winkler, et al. (2010) developed a framework to assess the performance of power transmission and distribution system under hurricane events. The framework embeds component fragility functions and connectivity based network study to capture the impact of component and network topology on system reliability. However, a graph model based on node connectivity can significantly underestimate power loss, as connected consumers can still experience power outage due to voltage issues. Han et al. (2009) presents a statistical power outage risk estimation model, where principle component analysis is used to deal with collinearity and generalised linear models are employed for regression analysis. However, this statistical approach does not consider any internal mechanisms of the distribution system. This paper presents a seismic risk assessment framework designated for electric power distribution systems. The framework considers component fragility, substation configuration and network characteristics including topology and physical power flow.
Probabilistic estimates of the damage state of network components are accounted for by considering their individual fragility function and ground motion intensity. A distribution load flow that can handle both weakly meshed and radial network characteristics is conducted to evaluate post-event lost load due to abnormal voltage profile and disconnection. Monte Carlo simulation is used to obtain probability distribution of the lost load and associated risk metrics. The method is demonstrated using a modified CIGRE medium voltage distribution network for both a case study earthquake approach and a simplified representation where ground motion intensity is consistent across the network. The effect of changes in network and substation topology on a range of risk metrics and importance measures are demonstrated using these networks.

Seismic risk assessment framework for electric power distribution systems
The seismic risk assessment framework for electric power distribution systems presented in this paper specifically models both the seismic behaviour of critical substation components as well as substation internal logic due to the vital role substations play in the distribution system. In addition to substations, distribution network functionality is essential to ensure power delivery to end users. In the immediate aftermath of an earthquake, the consumers may lose mains supply due to complete loss of connections or degraded voltage levels. A distribution load flow algorithm that can effectively handle both radial and weakly meshed distribution networks is adopted to assess the scale of power loss following an earthquake. MCS is applied to estimate the probability distribution of power loss, from which several risk metrics can be defined for the seismic risk evaluation for distribution systems. The details for each part of the framework are presented in the following sections.
Typical power system reliability analysis not only involves estimation of the scale of power loss but also considers duration of the loss as the primary concern is the value of the service that the network can provide. This framework only considers immediate power loss caused by earthquake related damage, as the time of repair after an earthquake is partly an operational decision, and also dependent on road access which would require integrated power-transportation network interdependency modelling. Although the operational decision and transportation availability can vary considerably between different operators and situations, electric power system fragility is dependent on the seismic and co-seismic hazard, the design and specification of the asset and its maintenance program. This framework could be combined with approaches to estimate repair time and regeneration of the network, but this is outside the scope of the current study.

Component fragility functions
The relationship between ground motion intensity and component failure is quantified using fragility functions, which determines the probability that a component will reach or exceed a damage state for a certain ground motion intensity measure. This study focuses on the characterization of operational impact, which is the reduction in load serving capability in the aftermath of an earthquake. In this case the damage state of interest is when a component ceases to normally function and needs to be isolated for repair and replacement. Characterization of further damage states, representing more severe damage, are more useful in recovery modelling and the relationship between damage and the time and cost required to repair or replace components.
One of the most significant components of power systems that may experience damage is power transformers. The power transformer can slip and fall off from its foundation due to high intensity ground motions. Circuit breaker failure caused by ground shaking is usually due to porcelain column fracture. The most vulnerable part of a disconnector is the bottom joint of the porcelain column which is often susceptible to ground shaking. Despite the various types of damage to power systems components, these need to be classified in relation to whether the component needs to be isolated for repair and replacement. The differences between short-circuit, open circuit and lock-in faults are not specifically distinguished in the assessment framework. The performance of eight classes of substation equipment in twelve California earthquakes was studied by Anagnos (1999), with the majority of damage to circuit breakers and disconnectors. Seismic response of a vertical disconnector is investigated in Paolacci, et al. (2014) based on experiment results, and corresponding fragility curves evaluated. Zareei, et al. (2017) proposed multivariate fragility analysis to ascertain more reliable assessment of the seismic vulnerability, using a circuit breaker as an example.
Analytical fragility functions were developed for a power transformer in Zareei, et al. (2016) to characterize seismic vulnerability. HAZUS (FEMA, 1999) has a multiinfrastructure seismic fragility function library, including the fragility functions for substations obtained from a comprehensive earthquake damage database. However, the damage states for the substation fragility function are defined with respect to the percentage of components being damaged within the entire substation, which ignores the fact that different substation components can have different responses to ground motion.
The HAZUS method is thus more suitable for characterizing physical rather than functional damage. In order to characterize the impact of earthquake intensity on substation component functionality, this paper takes fragility functions for a range of components from Vanzi (1996), where the substation components fragility functions are obtained from shaking table testing. The parameters of the fragility functions for substation components are summarised in Table 1.

Substation modelling
Distribution substations are the bridges between transmission and distribution networks.  The energizing status of outgoing feeder 1 in Fig.1 is dependent on the status of B3 and B4, as well as that of the bridging components. The functional statuses of B3 and B4 depend on those of B1 and B2, and the transformer circuit group components.
The functional statuses of B1 and B2 depends on whether power can be fed from the incoming lines through the circuits. The logic governing the functionality of outgoing lines is divided into three hierarchical levels: a top level fault tree governing the availability of outgoing feeders (Fig.2) , medium level governing the availability of B3 and B4 (Fig.3) and bottom level determining the status of B1 and B2 (Fig.4).
In order for the feeder to be active, CB4 and D12 must remain connected.
Additionally either D10 or D11 needs to be active, in order to be connected with one of the bus bars, B4 or B3.  The substation will not be able to supply any load if both B1 and B2 cannot be energised via incoming lines.

Figure 4 Fault tree of bottom level events
Note that due to the "N-1" security requirement, the two transformers are designed to be 50% loaded under normal operation condition so that even a single transformer has sufficient capacity to support the entire network under emergency or maintenance conditions on the other. Therefore the availability of a single transformer as well as the associated switch gears is required to ensure the availability of outgoing feeders.
Minimal cut sets up to triple components for the availability of outgoing feeder 1 are listed in Appendix.

Distribution network power flow
Power flow is widely used for power system operation and planning, and is one of the most important analyses for determining whether the power system can be normally operated (Saadat, 2010). The power flow determines the magnitude and phase angle of voltages at each bus, as well as the real and reactive power flow in each line. A normally functional distribution system must be able to hold the voltage at each bus within a specific range, out of which the consumers will experience power outage. For each bus i, the real and reactive power balance equations are as follows: is solved within each iteration and is used to update the solutions until convergence criteria is met: The Newton-Raphson method is conventionally applied to solve power flow of transmission systems, and has been employed to solve network flow based models in Cavalieri, et al. (2014), Pires et al. (1996), Shinozuka et al. (2007) and Vanzi (1996).
However, it is known to have convergence problems for solving distribution power flow due to the special features of distribution systems. Transmission networks are highly reactive (i.e. low resistance to reactance ratio (R/X)), which results in strong interactions between P ,  and Q , V . This makes the Jacobian matrix diagonal dominant, facilitating computation of its inverse. Distribution networks are much more resistive compared with transmission networks, meaning the R/X ratio is higher. This can lead to singularity of the Jacobian matrix, forbidding the application of Newton-Raphson method. In order to overcome the difficulty, the forward and backward sweep method (Kersting, 2012) was proposed to solve distribution power flow. The forward and backward sweep method proves particularly effective for radial distribution networks, but can be tedious for the contingent power flow. Teng (2003) proposes a method that directly applies Kirchihoff's laws to obtain the load flow results without using Newton-Raphson or forward backward sweep techniques. However, the BIBC, BCBV matrices used by the method containing network topological and electrical characteristics need to be formed by extraction of local topological information making it difficult for computer implementation and extension. The implicit Z-bus method (Chen et al. 1991) only requires the formation of traditional Y-bus matrix for which it is convenient to represent removal of network components therefore facilitating computation of the contingent distribution power flow. Furthermore, the convergence of implicit Z-bus method is independent of the R/X ratio and thus is adopted as the solution algorithm for the distribution power flow.
Consider an n-bus distribution network, the bus current injections can be expressed as: ( 1) ( ) ( 1) n n n n Inj bus The current and voltage relations can be separated by substation swing bus and load P-Q buses:

Loss of load
Three metrics are defined to characterise the risk of loss of load: • Expected load not served (ELNS) -reflects average loss of load. Representation of the loss in an average sense does not account for that which results from extreme hazard events.
• Conditional load at risk (CLaR) -represents the expected loss sustained beyond load at risk accounting for pessimistic possibilities. ELNS is defined as the expectation of and is estimated by the average of multiple independent samples of ( ) determined by independent samples of . The sample size is adaptively determined against the coefficient of variation ̂√ , in order to determine an appropriate sample size and allow for trade-off between accuracy and computation time.
LaR at confidence level is defined as the value at risk (VaR) of at confidence level , which represents the minimum loss of load that will not be exceeded with at least probability . In order to account for the risk beyond LaR, CLaR at confidence level is defined as the conditional value at risk (CVaR) of at confidence level , which is the expectation for loss greater than or equal to load at risk at confidence level .

( ) { | } loss loss loss
CVaR P E P P VaR  = It can be shown that CVaR  is the optimal value to where VaR  is a solution (Pflug, 2000;Rockafellar and Uryasev, 2002 Computation of importance measures is a relevant outcome of the probabilistic seismic risk assessment of distribution systems. Two metrics are defined to characterise the criticality of different substation components: • Risk Achievement Worth (RAW) -the ratio of the risk when a component is assumed to be always failed to the base case model risk. The RAW represents maximum relative possible increase in total risk due to failure of a certain component. It can be used to identify components that need to be seismically retrofitted in order to reduce loss of functionality.

Case Study Applications
The proposed seismic risk assessment framework is demonstrated using a modified CIGRE medium voltage test network (Rudion et al. 2006), details of which are summarised in Fig.6. The network has a rural character and supplies a small town and the surrounding area. The rated voltage level of the network is 20 kV, and it is supplied from a 110 kV transformer substation. The network can be operated either as weakly meshed by closing switches S1, S2 and S3 in Figure    Two approaches for the representation of the hazard are used to assess this case study network using the seismic risk framework. In the first approach, a scenario earthquake with associated geospatial variability of ground motion intensities is used.
Peak ground acceleration (PGA) is used as the ground motion intensity measure for the study. This approach evaluates the impact of a specific earthquake scenario on the network, and could be a useful approach if there is a fault that dominates the hazard in a certain region. This could be expanded to assess a range of earthquake scenarios, however this is beyond the scope of the current paper. The second approach applies a range of uniformly distributed PGA values across the entire network and the influence of various PGA levels on the seismic risk is assessed. Seismic risk is quantified with the two hazard representation approaches in terms of network risk metrics and substation component importance measures. This includes the effect of bus bar configuration in substation design, and the effect of a meshed and radial network topology.

Benchmark scenario
For the earthquake scenario approach, the network is assumed to be situated in a region in the South Island of New Zealand that was affected by the 2016 Mw7.8 Kaikōura earthquake. The geospatial PGA distribution, interpolated based on recorded ground motions, is applied across the test network as shown in Figure 7. This is not meant to represent any real-world network, with the PGA distributed used to demonstrate the application of the framework. Substation circuit configuration in this benchmark scenario is the typical double bus-bar system as described in Section 2.2. The network is operated as a weakly meshed network where the switches are closed. The numbering in this figure matches the numbering of each bus in Figure 6. Figure 7 Geospatial representation of the case study network topology and variation of peak ground acceleration across the area.
Probability distribution of the unserved load is estimated using the assessment framework and is summarised in Fig.8. It can be seen that the loss of load distributes itself into three clusters. Under the majority (70%) of the simulated scenarios the unserved demand is between 0 and 5.5 MW (referred to as minor loss of load). The unserved demand is also distributed between 21.6 and 25.9 MW for 13% of the simulated scenarios (referred to as medium loss of load), and between 44.6 and 44.7MW for the remaining simulated scenarios (referred to as major loss of load). The simulations show that minor loss of load is most likely to happen resulting from network disconnections due to bus failures and abnormally low voltage levels due to changes in topologies, leading to unexpected long distribution branches. The major loss of load has the lowest frequency of occurrence (17%), but these are high impact scenarios with loss of power to almost the entire network. This is mainly because of the loss of substation feeders. The expected load not served (ELNS) for the entire distribution is 14.1MW, which can be over-optimistic and does not show the losses from the tail of the distribution. The 80%-LaR of the system subjected to the simulated earthquake is 24.9MW. This implies that with 80% confidence the earthquake will not result in loss of load at a scale greater than 24.9MW. In pessimistic cases when the loss of load exceeds 24.9MW, the expected loss of load is 40.0 MW (80%-CLaR).

Assessment of substation design
The effect of substation design is assessed through the comparison of a single bus bar design and a double bus bar design. The single bus-bar layout is the same as the double bus-bar configuration shown in Fig. 1, except that one of the bus and the associated circuit breakers/disconnectors are removed. The loss of load histogram of both single and double bus bar layout with meshed networks are shown in Fig.9. Fig.9 (b) presents probability distributions of unserved demand when a single bus-bar is used. In contrast to the double bus bar design, where the scenarios of minor loss of load have the highest probability of occurrence, the single bus bar system undergoes major loss of load for 48% of the simulated scenarios. The expected load not served for the single bus bar configuration is 25.5MW, which at 80% is significantly higher than that for the double bus bar configuration. The 80%-CLaR for the single bus bar is 44.7MW, indicating that in a pessimistic case the substation could lose all its load supply capability, which is worse than the double bus case where the 80%-CLaR is 40.0MW (89.5% of the full load supply) Figure 9 Histogram of unserved load for the single and double bus bar substation configurations as part of a weakly meshed network.

Assessment of network topology
The seismic risk of different network topologies are assessed through the structure of the network as weakly meshed or radial. This is dependent on the status of three switches in the network as described in Section 3. Fig.10 shows histograms for unserved load when the components of the system are exposed to ground motions in the earthquake scenario (Fig. 7). For single bus substation design, the expected load not served for meshed and radial network are 25.5and 25.8MW respectively. For double bus substation design, the ELNS for the two network layouts are 14.10 and 14.07MW respectively. This shows that the effect of network topology is not as significant as that of substation configuration for this test network and scenario.
Figure 10 Histogram of unserved load for different combinations of single/double bus substation layouts and radial/meshed network topologies.

Criticality assessment of substation components
The RAW and RRW indices for the double bus bar substation configuration and weakly meshed network topology are summarised in Fig. 11 and 12 respectively.. It can be seen from Fig.11 that loss of components CB4, D12, CB8 and D24 has a significantly higher contribution to increasing the risk than the rest of the substation components. This is primarily because failure of components CB4 and D12 will directly lead to failure of feeder 1, with CB8 and D24 leading to the failure of feeder 2. Seismic strengthening of these components should be prioritised as most other circuits have greater redundancy.
The second highest RAW group of components constitutes the four bus bars (B1-4), loss of which leads to almost 1.5 times the risk of other remaining components. This is followed by the transformer circuits group (T1 with D6 D7, T2 with D18 D19). Their loss of functionality will result in approximately 1.3 times the risk of other remaining components in this scenario. It is therefore not sufficient to assess component criticality solely based on RAW or RRW index. Rather, both of them should be considered simultaneously as well as the component status.

Uniform intensity approach
Instead of applying an earthquake scenario with a specific spatially distributed ground motion intensity, the sensitivity of power disruption to ground motion intensity is assessed by applying a range of uniform PGAs across the network.

Assessment of Substation Design
The result obtained from the previous comparison of single and double bus bar system for the scenario earthquake approach is also demonstrated here where the system is subjected to uniformly distributed PGA. Fig. 14 shows that for a specific PGA, the expected load not served for double bus-bar is lower than that of single bus-bar design. Fig. 15 shows a lower cut-off PGA for the single bus configuration (0.19g) than the double bus configuration (0.29g), suggesting a lower risk with double bus bar layouts.
The double bus bar arrangement is assessed to be more reliable in this case as one bus bar can provide an alternative route when the other is damaged. This finding, in conjunction with the results associated with RAW and RRW (Section 3.1.4), suggests that double bus bar designs should be given preference over single bus bars when installation occurs following new substation set-up in areas with high seismic hazard.
By quantifying the network risk metrics, the effect of different substation designs on potential loss of load can be clearly demonstrated. However, given that bus bars are particularly vulnerable components to ground motion damage, and that repair time for damaged double bus bars can be longer than single bus bars (Li, 2014), newly installed double bus bars in such areas should offer high seismic resilience. Seismic retrofitting of existing bus bars should be considered in substation maintenance and investment programmes, as should potential alterations to substation design to incorporate double bus bars in place of existing single bus bars. Such approaches will reduce the seismic risk to electric power distribution systems, and estimates of loss of load can inform the justification of these decisions.

Assessment of Network Topology
The scenario earthquake suggested that the effect of network topology is not as significant as that of substation configuration. The result is also highlighted by inspecting a range of uniformly distributed PGA intensities across the network as shown in Fig.16, where the two network fragility curves for both single and double bus bar substation configuration overlap one another, despite the different network topologies. Fig.17 shows that the cut-off PGA, where 95% LaR/CLaR reaches full load serving capacity, is relatively independent of network topology, rather it is more closely related to the substation configuration. As discussed previously, distribution networks are not typically built to be meshed due to economic reasons, in contrast to that typically seen in transmission systems. The network assessed in this paper is a typical medium voltage distribution network and it is constructed to be weakly meshed in order to provide alternative conducting paths during normal operations. It is clear in this case that change in network topology would not contribute to seismic risk reduction, however in other cases this may not be the case, particularly in situations where parts of the network or particular substations are exposed to higher seismic hazards, whether ground motion or others such as landslides, liquefaction or surface fault rupture. This falls outside the scope of this paper but this framework allows for the assessment of these different topologies and range of scenarios. If assume all lost load is of the same importance, could identify an acceptable level of loss of loadbased on range of factors. This could also be linked to repair time and availability estimates. Across these examples we have effectively assumed loss of load is the same regardless of the type of facility across the network. If critical facilities are connected to the network, this approach could be modified by changing the metrics used to account for the variation in criticality/importance.  The proposed risk evaluation method can be applied to assess the seismic risk of various substation configurations and distribution network topologies, providing distribution utility managers with crucial insights for network design and informing rational investment decisions.