A bayesian approach for estimating the post-earthquake recovery trajectories of electric power systems in Japan

Post-disaster recovery modelling of engineering systems has become an important facet of catastrophe risk modelling and management for natural hazards. The post-disaster recovery trajectory of a civil infrastructure system can be quantified using (a) the initial post-disaster functionality level, Q o ; (b) rapidity, h (i.e., the rate of functionality restoration); and (c) recovery time, R t . This study uses a Bayesian estimation approach to derive a set of probabilistic models to estimate Q o , R t , and h of electric power networks (EPNs) using post-earthquake recovery data from 16 large earthquakes in Japan between 2003 and 2022. The considered predictor (explanatory) variables include earthquake magnitude, year of occurrence, seismic intensity, and exposed population (PEX). Apart from being a simple and efficient stand-alone tool, the proposed data-driven models can be a useful benchmarking tool for simulation-based approaches for EPN recovery modelling.


Introduction
Field evidence has shown that disaster-induced damage and disruptions to critical infrastructure systems (e.g., electric power networks; EPNs) can cause significant direct and indirect socioeconomic losses, including casualties.Apart from economic losses associated with business interruptions, power outages can result in human losses, increased waste from perishable food items, failure of security systems, increased disease spread, significant direct repair and restoration costs, and several other forms of nuisance to the general public (e.g., Chang, McDaniels, Mikawoz, et al., 2007;Dugan, Byles, & Mohagheghi, 2023).
Accurately estimating power outages and recovery can help EPN management authorities and other private/public organizations reliant on electric power define effective short-term and long-term recovery strategies to improve community-level resilience.For example, business owners and homeowners can set up appropriate measures (e.g., backup power sources and power disruption insurance) if potential post-disaster recovery trajectories of EPNs at different hazard intensities are known.
The post-disaster recovery trajectory of civil infrastructure systems can be quantified using (a) their initial post-disaster functionality level (i.e., the ratio of the number of serviced customers/end users post-disaster to that pre-disaster); (b) rapidity (i.e., the rate of functionality restoration); and (c) recovery time (i.e., time to restore full functionality to the total number of serviced customers/end users).Currently, two main approaches for infrastructure recovery modelling are widely adoptedempirical and simulation-based modelling.A more detailed discussion of the relevant state-of-the-art studies is provided subsequently in this paper.Empirical models (e.g., Guikema & Quiring, 2012;Liu, Davidson, Rosowsky, et al., 2005;Nojima & Sugito, 2003) for EPN assessments are based on correlations from data analysis of historical spatial outages and recovery times.Given the data resolution used in developing these models, most empirical models do not explicitly account for the network topology and fragility (probability of a given damage level vs hazard intensity) of individual components and subcomponents in the EPNs.However, empirical models are considered helpful in forecasting regionlevel disaster-induced outages.It is noted that most existing empirical models have been developed for hurricane and storm hazards.Furthermore, available empirical models on post-earthquake recovery of EPNs do not adequately link earthquake features (e.g., magnitude, local intensities, etc.), exposed population, post-earthquake functionality level, and recovery time.Hence, developing more efficient yet simple empirical models for post-earthquake recovery of EPNs is crucial.
Simulation-based methods (e.g., Çağnan, Davidson, & Guikema, 2006;Guidotti, Chmielewski, Unnikrishnan, et al., 2016;Ouyang & Dueñas-Osorio, 2014) use network analysis to simulate the post-disaster functionality level and recovery trajectory of EPNs, while explicitly accounting for component and system fragilities.One of the fundamental challenges of simulation-based models is that, due to security reasons, EPN topologies are not publicly available data.Also, reliable information on component and system fragilities and repair/replacement times and sequences are needed to develop reliable estimates of post-disaster recovery trajectories.The adequacy of simulation-based methods can be improved through appropriate validation exercises using real-life events.Empirical models may be helpful benchmarking tools whenever such validation exercises are not feasible.This study's main objective is to enhance the disaster risk modelling and management of EPNs, filling the discussed gaps above.Specifically, this study seeks to develop probabilistic models that adequately link earthquake features (e.g., magnitude, local intensity, year of occurrence), exposed population, initial post-earthquake functionality level, rapidity, and recovery time.
To this end, this study (a) collects and aggregates data on the performance and restoration of EPNs from past earthquake events in Japan; and (b) develops simple probabilistic models for post-disaster recovery trajectory estimation using relatively easy-to-obtain information such as seismic intensity, event magnitude, exposure data (in terms of number of households exposed to each local seismic intensity level), number of serviced households, and EPN restoration trajectory data.
The paper is organized as follows.Firstly, a brief overview of existing research studies is presented, and the research gaps (and corresponding research objectives) are defined.Subsequently, the data collection and aggregation process is described.Finally, Bayesian regression analyses are carried out to develop probabilistic models to estimate the initial post-earthquake functionality level, rapidity, and recovery time of EPNs as a function of the considered predictors.

Existing research studies
The concept of post-disaster recovery of critical infrastructure systems and its analytical modelling has gained widespread attention in the last two decades.As highlighted by the global disaster risk reduction sector (e.g., UNISDR, 2015), the need for safer and more resilient communities has catalyzed significant research efforts into post-disaster recovery modelling.As briefly discussed in the Introduction above, postdisaster recovery modelling frameworks in the literature are either empirical or simulation-based.Given the scope of the current study, this section explicitly focuses on EPNs.
Before discussing state-of-the-art studies on postdisaster recovery modelling of EPNs, it is essential to highlight the critical output metrics of any recovery modelling framework.Figure 1 presents a post-disaster recovery trajectory for a component/system (or even a community) with a pre-disaster functionality level Q pre , reduced to Q o (referred to as the initial post-disaster functionality level) after the disaster at time t o .Repairing or replacing damaged system components will result in the gradual restoration of functionality.The time required to restore the system's functionality to a desired level Q R (which can be similar, close to, or better than the pre-disaster functionality level) is referred to as the recovery time (R t ).The recovery trajectory of the damaged system (i.e., from {t o ,Q o } to {t R , Q R }) can be quantified by a recovery function, Q(t).The recovery rapidity h quantifies the recovery rate.An ideal recovery modelling approach should be able to estimate Q o , R t , h, and Q(t).Furthermore, any sources of uncertainties in these metrics need to be accounted for adequately.
Empirical post-disaster recovery models for EPNs adopt field data from past events to develop statistical models to estimate the power outage level and restoration time.Proposed tools range from simple linear regression models to sophisticated machine learning tools.Empirical models are either set up for post-disaster recovery modelling at a (regular) grid or community level.For instance, Liu, Davidson, Rosowsky, et al. (2005) adopted data from three hurricanes in North and South Carolina, U.S.A., to develop negative binomial regression models that relate the expected number of outages over one square kilometre grid cells to the number of transformers, maximum gust windspeed, and affected power company.Liu, Davidson, and Apanasovich (2008) used data from six hurricanes and eight ice storms to develop a spatial generalized linear mixed modelling for estimating power outages over nine square kilometre (3 km × 3 km) grid cells.The authors concluded that maximum gust wind and ice thickness are directly related to higher outage counts in hurricanes and ice storms, respectively.The number of protective devices was also found to influence hurricane and ice storm-induced outage levels.Other studies that have developed empirical power outage models include (but are not limited to) Guikema and Quiring (2012) -for hurricanes; Cerrai, Koukoula, Watson, et al. (2020) -for snow and ice storms; and He, Wanik, Hartman, et al. (2017) -for storms.It is noted that the majority of existing empirical models were developed for hurricane and storm hazards.
It is noted that most empirical studies focus on estimating power outages.Only a few studies have developed empirical models for recovery time estimation.These studies use limited data (e.g., a single event), a small number of explanatory variables and consider other hazards (instead of earthquakes).Nojima and Sugito (2003) developed a probabilistic model for evaluating the outage level and recovery time of EPNs using data from the 1995 Hyogoken-Nanbu (Kobe), Japan, earthquake.Variables considered in the model include exposed population and seismic intensity.This model was subsequently modified using the 2011 Tohoku, Japan, earthquake (Nojima & Kato, 2014).Liu, Davidson, and Apanasovich (2007) adopted accelerated failure models, developed using data from six hurricanes and eight ice storms, to predict the post-disaster outage duration.Variables considered in the Liu, Davidson, and Apanasovich (2007) study include maximum gust wind speed, duration of strong winds, 7-day rainfall, ice thickness, and population density.Kammouh, Cimellaro, and Mahin (2018) developed probabilistic recovery curves for estimating downtime for power, water, gas, and telecommunications systems using data from significant earthquakes worldwide from 1960 to 2015.Using such data, the authors proposed recovery curves as a function of earthquake magnitude and development level of the affected country.
Over the last three decades, significant research has been conducted to develop and adopt simulation-based methodologies for post-earthquake resilience assessment of EPNs.Simulation-based methods use network analysis, whereby the power grid is modelled as a network with nodes and edges.The nodes represent the generators, substations, and load points.The edges represent the transmission and distribution lines in the network.Each component or subcomponent in the network is modelled with its own fragility models.Ang, Pires, and Villaverde (1996) developed a probabilistic model for estimating the seismic reliability of electric power transmission systems under earthquake ground shaking.Such a model adopts Monte Carlo simulation to assess the power outage probability by combining the fragility models of critical equipment and substations with ground shaking intensities at the local sites through power flow and network analyses.Cagnan and Davidson (2007) developed a simulation-based model of the post-earthquake restoration process of EPNs.The model captures geographic variability of risk and restoration resource constraints (e.g., material and workmanship) in estimating the restoration trajectory of EPNs.Ouyang and Dueñas-Osorio ( 2014) developed a restoration model that incorporates resource availability, resource mobilization and restoration sequence in estimating the recovery trajectory of EPNs.The restoration sequence algorithm in their study prioritizes EPN systems and components vital to public safety, health and welfare.Mensah and Dueñas-Osorio (2016) proposed a Bayesian Network-based framework to estimate power outages in electrical grids under hurricane hazards.Their proposed framework is combined with the Ouyang and Dueñas-Osorio (2014) restoration model to develop a resilience enhancement framework.
Based on information available in existing simulation-based studies, it is noted that most existing studies do not provide a holistic tool for estimating all relevant parameters of the recovery (or restoration) curve.Also, most current models cannot effectively evaluate the recovery function and rapidity of damaged EPNs.Furthermore, most existing models were developed for post-hurricane and post-storm recovery modelling.Some models calibrated for earthquake-induced damages (e.g., Nojima and Kato, 2014) are based on just one or two events.Therefore, a more reliable estimation of post-earthquake recovery trajectories of EPNs may be achievable by considering more earthquake events.
While simulation-based methodologies may be desirable in some instances, the fact that topology and typologies of EPNs are not publicly available information means that most studies rely on synthetic test beds.Therefore, it is essential to develop practical empirical tools that can be used to benchmark simulation-based tools.
This study introduces simple yet effective empirical models to fill the observed gaps in the recovery modelling of EPNs following earthquakes.The novelty of these models lies in their ability to: 1) be specifically tailored to earthquake events occurring in Japan; 2) incorporate a wider range of past events and explanatory variables than those considered in previous studies; and 3) estimate the initial post-earthquake functionality, rapidity, and recovery time of EPNs, including associated uncertainties, through a Bayesian estimation approach.These aspects mark a shift from earlier studies that primarily focused on estimating power outages from other natural hazards using limited data and explanatory variables.The proposed models may be helpful benchmarking tools for simulation-based approaches to estimate the postearthquake recovery of EPNs.

Data collection and aggregation
The first step is establishing a database of EPN performances/damages and recovery from past earthquakes in Japan.The data collection exercise for this study focused on major damaging earthquakes that caused significant damage to EPNs.Another criterion was the need to have sufficient information on the initial post-disaster functionality level, recovery trajectory and recovery time.Based on these criteria, 16 events between 2003 and 2022 were considered.Information collected includes the earthquake magnitude, time of occurrence, exposed population data, number of serviced households, and recovery time for the EPNs.The data obtained from the above are summarised in the table in Appendix A and are readily available for use by interested researchers.The subsequent subsections provide details on the data collection process.

Earthquake magnitude
This study adopts the Japan Meteorological Agency (JMA) magnitude (M JMA ) to quantify the magnitude of each event.More details on the JMA magnitude scale can be found in JMA ( 2014).The M JMA for each event was extracted from Japan Real-time Information for earthQuake (J-RISQ) reports (J-RISQ, 2015).The J-RISQ report database contains information for events from 1996 to date.The J-RISQ reports are typically published immediately after an earthquake.

Population data (number of households and population)
The number of households and populations in the affected regions are taken from publicly available census data (SBJ, 2020).For the number of households for the entire region, if the municipality where the affected households are located is known, the total number of households within that area is used; if not, the total number of households within the jurisdiction of the electricity company's sales office or prefecture is used.In some instances, information on the total number of households was unavailable.However, it was inferred from events with information on both the number of households and the population that the average ratio of population to the number of households equals 2.8 (i.e., 2.8 persons per household).Hence, if the number of households was unavailable, it was inferred by dividing the reported population by 2.8.

Number of serviced customers (households)
This is the number of households with continued access to power immediately after the event.The number of serviced customers and total customers were extracted from published literature.Priority was given to data published by municipal electricity companies.In cases where incomplete information is available from electricity companies, data from newspapers and journal articles were used.Q o is defined as the ratio of the number of households with continued access to power immediately after the disaster to the total number of households in the affected region.

PEX (population exposed to the earthquake) data
Data on the exposed population at each JMA (local) seismic intensity scale level were extracted from J-RISQ reports (J-RISQ, 2015).JMA intensity scale intensities range from 0 to 7.More information on the JMA intensity scales can be found in JMA (2015).
JMA ( 2015) identifies seismic intensities of '5 Lower' or more to influence electric power disruptions.Hence PEX data collection focused on the population exposed to seismic intensities '5 Lower' (P 5 L ), '5 Upper' (P 5 U ), ' 6Lower' (P 6 L ), '6 Upper' (P 6 U ), and 7 (P 7 ).Figure 2 shows an example of the extracted J-RISQ reports.The PEX data for each event was normalized by the total population in the high-intensity zones (as defined in section 2.2) (i.e., p 5 u = P 5 U /(P 5 L +P 5 U + P 6 L + P 6 U + P 7 )).

Recovery time
The recovery time was defined as the time taken to restore power to all affected households in the affected region.The recovery time was extracted from publicly available information.Priority was given to information published by electricity companies, local authorities and relevant ministries and agencies.In some instances, newspaper data (typically in the form of government announcements) were used.

Rapidity
As mentioned earlier, rapidity characterizes the recovery rate of the EPN.Opabola and Galasso (2023) proposed a recovery trajectory function Q(t) that is dependent on Q o , mobilization time t 1 (i.e., the time after which restoration work starts), R t , and rapidity coefficients g and h (See Equation ( 1)).
The values of g and h express the level of preparedness, resource availability, technical know-how and societal conditions of a community.For known values of h and Q(R t ), g can be computed from Equation (1) as: g It is noted that g must be greater than zero.
A concave recovery curve corresponds to h << 0 and represents the recovery trajectory of a community with a poor level of preparedness, a high level of resource unavailability, and unfavourable societal conditions.In contrast, h >> 0 represents a community with good preparedness and resource availability and favourable societal conditions.Figure 3 shows how the value of h influences the recovery trajectory of a recovering system.The rapidity coefficient h from each event was estimated by fitting Equation ( 1) to the field data (see Figure 4 for an example) so that the sum of squared estimate of errors(SSE) between the observed and predicted functionality trajectory during the recovery phase is minimized.

Correlation analysis
A correlation analysis for the aggregated database is conducted first.The correlation matrix is presented in Figure 5.As shown in the figure, Q 0 is negatively correlated with the PEX (more significant correlation with p 7 and p 6 u ).The trend indicates that as the proportion of a population exposed to higher seismic intensity increases, the initial post-disaster functionality level decreases.
R t is shown to be highly correlated with Q o (with correlation coefficient r = −0.84),year of occurrence and earthquake magnitude.Low Q o corresponds to increased damage to EPN components, resulting in a longer time to recover.The significant correlation of PEX with R t may be attributed to the fact that PEX and Q o are also highly correlated.
Furthermore, the rapidity coefficient h is shown to be correlated with Q o and R t with r equals 0.49 and −0.42, respectively.As earlier described, a higher level of h depicts a system with good preparedness, good resource availability, and favourable societal conditions.Hence, the positive correlation between h and Q o can be attributed to the fact that a higher Q o is also associated with good preparedness and favourable societal conditions.

Bayesian linear regression approach
Unlike frequentist approaches, a Bayesian parameter estimation yields a probability distribution of model parameters instead of a single value.There are some advantages of Bayesian regression models over standard frequentist regression, for instance: 1) they tend to perform better than standard frequentist regression models when working with a relatively small sample size (as in this study).In fact, they provide inferences that are conditional on the data and are exact, without reliance on asymptotic approximation; 2) credible intervals (or regions), which are the Bayesian version of confidence intervals, have a much more straightforward interpretation than the confidence intervals from standard linear regression models (e.g., 'the true parameter has a probability of 0.95 of falling in a 95% credible interval'); 3) they are generally more flexible than frequentist regression models.For instance, they provide a natural  and principled way of combining prior information with data, within a solid decision theoretical framework; they also provide a convenient setting for a wide range of models, such as hierarchical models and missing data problems.Various literature (e.g., Gelman, Carlin, Stern, et al., 1995;Reich & Ghosh, 2019) provides a detailed introduction to Bayesian methods, which is outside the scope of this study.This paper just provides a brief summary of Bayesian parameter estimation.Consider a standard linear regression model given as: where y n is a scalar response, β is a vector of model parameters (or coefficients) for the vector of regressors x n , and error ε n is a zero mean Gaussian noise term with a non-zero variance σ 2 .The probabilistic interpretation of Equation ( 3) can be written as: In Equation ( 4), y n is a normally-distributed random variable with mean of β T x n and a corresponding nonzero variance of σ 2 .
The Bayesian treatment of a linear regression model introduces a prior probability distribution over the model parameters β and σ.In this study, we consider the semiconjugate prior for which the prior and posterior share the same parametric family.The natural semiconjugate prior of parameters β and σ 2 is a normalinverse-gamma distribution, with the form: 5) where N is the normal distribution, and IG is the inverse-gamma distribution.V is the conditional covariance matrix of the Gaussian prior of β ; and μ is the mean hyperparameter of Gaussian prior on β.A and B are the shape and scale hyperparameters of the inverse gamma prior on σ 2 , respectively.
The prior probability distribution of a parameter can then be combined with the likelihood of the observed data to obtain the posterior distribution of the parameter (see Equation ( 6)). 6) This study adopts the Markov Chain Monte Carlo (MCMC) sampling algorithm, precisely Gibbs sampling (e.g., Gelfand, 2000) for the posterior estimation.The Bayesian regression analyses presented in this study drew 10,000 samples from the posterior distribution generated using a burn-in period of 10,000 and a thinning level of 10.It is also noted that a stepwise removal process was carried out via the Bayesian estimation method to identify the governing predictors for Q o , h, and R t .

Initial post-disaster functionality level (Q o )
The Bayesian approach outlined in section 4.2 defines the probabilistic formulation for Q o .Two data points were removed from the analysis.The 2011 Tohoku event was excluded because the functionality loss and recovery trajectory data are assumed to be significantly influenced by the earthquake and tsunami sequence (rather than earthquake only, as for the other considered events).The 2018 Hokkaido event was excluded because the total blackout was a precautionary measure due to a fire event at one power plant.The stepwise removal process shows that the main predictors for Q o are the normalized PEX data.Among the PEX data, p 5 l is determined to be the least significant in predicting Q o ; hence, it is not considered in the Bayesian analysis.
Table 1 presents the posterior summary statistics for the proposed model.The model assumes that the higher the proportion of a population exposed to a more significant seismic intensity, the lower the initial post-disaster functionality level, which is intuitive.Table 1 shows that the Gelman -Rubin convergence diagnostic (Gelman & Rubin, 1992) declares convergence (i.e., R < 1:1).
It is noted that the lower-and upper-bound values of Q o are zero and unity, respectively.Also, it is worthwhile to constrain the model parameters to negative values (i.e., to capture the reduction in Q o ).In line with the mean values, another constraint may be to ensure 6 shows the relationship between the estimated mean and measured Q o .As shown in the figure, the model provides good estimates of Q o with an R 2 of 0.81 and normalized root mean square error (NRMSE) of 0.1.

Recovery time
The stepwise removal process shows that the three main predictors for R t are M JMA , Q o , and the occurrence year (Y).As discussed earlier, Q o accounts for PEX.It is presumed that the sensitivity of recovery time to event magnitude is because the earthquake severity on other infrastructure systems and the entire community determines the resources dedicated to repairing damaged EPNs.For example, larger events could have severe impacts on other lifelines which are interdependent with EPNs.
Based on the estimated RMSE, the predictors relate better with R t in the natural log space.Hence, the Bayesian analysis was carried out in natural log space (Equation ( 7)).ln Table 2 presents the posterior summary statistics for the proposed model.The proposed model provides adequate estimates of R t with an R 2 of 0.8 and NRSME of 0.2.The model captures the increase in recovery time with an increase in magnitude and proportion of customers without service (i.e., 1 -Q o ).Similarly, the sensitivity of recovery time to occurrence year can be attributed to disaster management agencies learning from past events to reduce the post-earthquake recovery time.It is noted that a stepwise removal of occurrence year as a predictor of R t drops the R 2 from 0.8 to 0.65.This highlights the strong correlation between occurrence year and R t in the data used for this study.Hence, it was decided to retain occurrence year as a predictor.However, the availability of more EPN recovery data will be helpful for future studies to validate this relationship and refine the recovery time model.
Table 2 shows that the Gelman -Rubin convergence diagnostic declares convergence (i.e., R < 1:1).Figure 7 shows the good relationship between the estimated mean and observed R t .

Rapidity coefficient h
As shown in Figure 5, h is positively correlated with the initial post-earthquake functionality level and negatively correlated with R t .There are also significant correlations between spatial seismic intensities and h.The stepwise removal process suggests that h = f(1-Q o ).R t was excluded during the removal process due to a high correlation with Q o (r = −0.84).Furthermore, as shown in Figure 8, there is a nonlinear relationship between shape constant h and 1-Q o .Hence, the Bayesian analysis was carried out in natural log space also in this case (see Equation ( 7)).Table 3 presents the posterior summary statistics for the proposed model.The adequacy of the mean estimate    is represented in Figure 8.As shown in the figure, the proposed model provides a good estimate of h with an R 2 of 0.84 and NRMSE of 0.17.Table 3 shows that the Gelman -Rubin convergence diagnostic declares convergence (i.e., R < 1:1).

Illustrative example
As an illustrative example, this section demonstrates how the proposed methodology can generate probabilistic realizations of recovery trajectories for an EPN in a region following an earthquake.For simplicity's sake, rather than carry out a hazard analysis to generate seismic intensity on a synthetic or real testbed, the seismic intensity map and PEX data from the 2016 Kumamoto earthquake are considered for this illustrative example (see Figure 2 and Table A1).The order of applying the proposed methodology for predicting the recovery trajectory of EPNs, using posterior summary statistics presented in Tables 1, 2, and 3, is presented in Figure 9.In this illustrative example, we generate 1000 realizations of the recovery trajectories of the EPN in the region affected by the 2016 Kumamoto earthquake.
Based on the PEX data for the 2016 Kumamoto event (from Table A1), using Table 1, the mean Q o is estimated as 0.6.Figure 10a compares the actual (i.e., field data) and probabilistic distribution of the 1000 simulated Q o .Figure 10a shows that the actual Q o falls within the simulated distribution of Q o .Each simulated Q o is combined with the occurrence year and earthquake magnitude, using Table 2, to simulate 1000 realizations of R t .Figure 10b shows the actual R t falls within the simulated distribution of R t .Finally, each simulated Q o is used to simulate 1000 realizations of h. Figure 10b shows that the actual h falls within the simulated distribution of h.
For each realization of Q o , R t , h, 1000 realizations of the recovery trajectory Q(t) are simulated.To do this, Equation ( 2) is used to simulate 1000 realizations of g (where g > 0).Subsequently, for a time range [0 R t ], Q(t) is simulated using Equation (1).It is assumed that t 1 (i.e., the time after restoration work starts) equals zero for the simulations.Figure 11 compares the actual and the 1000 simulated realizations of recovery trajectories of the EPN.
As shown in Figures 10 and 11, using a frequentist approach may not help risk analysts and decisionmakers capture the randomness (i.e., both aleatory and epistemic uncertainties) present in post-disaster recovery modelling.Furthermore, Figures 10 and 11 show that the proposed approach is an efficient post-earthquake recovery modelling technique that can serve as a stand-alone tool and/or be a useful benchmarking tool for simulation-based models.  1 Simulate recovery time R t using Table 2 Simulate rapidity h using Table 3 Simulate g using Equation ( 2) for n th simulation using Equation ( 1) No Yes Discard simulation Figure 9. Analysis flowchart using the proposed models (where N is the total number of required realizations).

Discussions and conclusions
Prediction of power outages and recovery following earthquakes can help electric power network (EPN) management authorities and other private/public organizations cultivate preparedness plans and devise shortterm and long-term recovery strategies to improve community-level resilience.Simulation-based recovery modelling approaches can capture the network topology and seismic fragility of the generation, transmission, and distribution systems.However, as a result of the unavailability of EPN topology for several earthquakeprone countries, simulation-based recovery modelling studies often adopt synthetic testbeds.The adequacy of simulation-based methods can be improved through appropriate validation exercises using realistic data.Empirical models may also be helpful benchmarking tools whenever such validation exercises are not feasible.Furthermore, such empirical models can serve as simple stand-alone tools for integrating community-level and building-level recovery modelling frameworks.Most existing empirical post-disaster recovery models were developed using hurricane and storm hazard events.Hence, it is important to develop models that can be used for post-earthquake recovery modelling.To this aim, this study developed a set of probabilistic models to predict the initial post-earthquake functionality level (i.e., the ratio of the number of serviced households/end users post-disaster to that pre-disaster), recovery time (i.e., the total time to restore full functionality to the total number of serviced households/end users), and recovery rapidity of earthquake-damaged EPNs (i.e., rate of functionality restoration), using Bayesian parameter estimation to capture uncertainties.The initial post-earthquake functionality level was found to be dependent on the seismic intensity and exposure characteristics -i.e., the higher the proportion of a population exposed to a more significant seismic intensity, the lower the initial postdisaster functionality level.The recovery time was found to be dependent on the initial post-earthquake functionality level, event magnitude, and year of occurrence.The sensitivity of recovery time to occurrence year can be attributed to disaster management agencies learning from past events to reduce the post-earthquake recovery time.The recovery rapidity is dependent on the initial functionality level -i.e., communities with low EPN functionality loss recover more rapidly.
One of the key limitations of the current study is the data size.Additional data collection exercises are needed to update the probability distribution functions of the parameters with new observations.Future studies will also look at collecting data from other countries to improve the applicability range of the proposed formulations.Furthermore, future studies can extend the presented simulation-based approach to other lifelines (e.g., water and gas networks).

Figure 2 .
Figure 2. Example of reported data for the 2016 M JMA 7.3 Kumamoto earthquake on seismic intensity and exposed population (source: J-RISQ, 2015).

Figure 3 .
Figure 3. Graphical representation of different values for the the h factor for a hypothetical system with different recovery trajectories but similar initial post-disaster functional level and recovery time.

Figure 4 .
Figure 4. Example of post-recovery trajectory following the 2016 M JMA 7.3 Kumamoto earthquake showing calibrated rapidity coefficient h.

Figure 5 .
Figure 5. Correlation matrix for extracted data.

Figure 6 .
Figure 6.Relationship between observed and estimated mean initial post-disaster functionality level (Q o ).

Figure 7 .
Figure 7. Relationship between the measured and estimated recovery time.

Figure
Figure Relationship between proportion of unserved households (1-Q o ) and rapidity coefficient h.

Figure 11 .Figure 10 .
Figure 11.A comparison of the actual and simulated realizations of recovery trajectories for the EPN following the 2016 Kumamoto earthquake.

Table 1 .
Posterior summary statistics for Q o model.

Table 2 .
Posterior summary statistics for R t model.

Table 3 .
summary statistics for h model.