Financial allocation and network recovery for interdependent wastewater treatment infrastructure: development of resilience metrics

1

ABSTRACT

The disrupted performance of wastewater treatment plants (WWTPs) when facing floods can cause complications in urban systems’ functioning. Therefore, utilities need to improve the resilience of their facilities and establish asset restoration plans. The present study aims at investigating two metrics to improve the WWTP’s resilience. An attribute-based resilience metric is devised using a fuzzy multi-criteria decision-making (FMCDM) framework to attain a budget allocation proposition. A network model representing the interdependent infrastructure is also formulated to develop a performance-based metric. The network’s post-disaster behavior is simulated to obtain an asset restoration sequence (RS). The results indicate that considering imprecisions in experts’ views lowers the possibility of underestimation or overestimation of resilience improvement when allocating funds. Furthermore, prioritizing facilities’ restoration will increase WWTPs’ resilience. The proposed methodology, which is tested on Brooklyn, provides a resilient solution to some challenges in asset management and can be applicable to other coastal settings.

KEYWORDS:

• Resilience assessment
• flood resilience improvement
• wastewater treatment plants
• asset management

1. Introduction

Urbanization and increased expected quality of service have created complex interdependencies in urban infrastructures. Disruptions in the functionality of these systems caused by natural hazards such as floods can profoundly affect urban societies. Furthermore, exposed infrastructures such as wastewater treatment plants (WWTPs) often placed near water in coastal cities are more vulnerable and need additional provisions for flood resilience improvement (Karamouz et al., 2016). The flood resilience of a system is its ability to withstand floods, mitigate the related adverse impacts, and return to normal conditions (Simonovic & Peck, 2013). Therefore, a WWTP’s resilience depends on its intrinsic strength and resources needed to withstand flooding before and during the hazard, as well as the rate and speed of its recovery pertaining to interdependencies after the event. As a result, the flood resilience of a WWTP can be improved by increasing the flood resistance aspects of these systems and prioritizing the restoration of interdependent assets to minimize the interruption of their serviceability.

Moreover, a pressing concern currently addressed by infrastructure asset managers is the lack of well-defined, resilience-informed asset operational planning in response to cascading interdependent failures (Imani & Hajializadeh, 2020). Therefore, the proposed resilience improvement measures need to be aligned with the current asset management programs (AMP) of WWTPs to initiate and promote new and adaptive pathways toward sustainable infrastructure management. An AMP entails all the events occurring over the life-cycle of an infrastructure’s assets (Malano et al., 1999). Although extreme events and natural disasters have long-lasting effects on the performance of an infrastructure’s assets, the role of resilience improvement measures in the AMP of urban infrastructures is among the less emphasized topics in the literature. To improve the resilience of a system, metrics need to be developed to quantify this system’s characteristic. Vugrin et al. (2017) indicated that resilience assessment methodologies could be classified as attribute-based or performance-based.

The attribute-based resilience metric of a system is often developed by distinguishing the characteristics of a system that can influence its resilience (Duijnhoven et al., 2017). A limitation of an attribute-based assessment is that it does not provide any estimation or confidence in how well the system will operate in the event of a disruption (Vugrin et al., 2017). Therefore, this metric can only provide insight into the preparedness of the system itself, not its response after the failure. However, attribute-based assessments are powerful techniques for resource allocation and quantifying a system’s resilience. This is often made possible by including subjective expert opinions on the pairwise comparison of different resilience elements. In this regard, multi-criteria decision-making (MCDM) methods have proven effective in identifying and analyzing the attributes contributing to a system’s resilience (Karamouz & Movahhed, 2021).

In several fields, MCDM techniques have been widely used in prioritizing resources and alternatives in the overall evaluation of systems (Anaokar et al., 2018; Tahmasebi Birgani & Yazdandoost, 2018; Zamani et al., 2020). Based on the input information, MCDM approaches are categorized into crisp MCDM and fuzzy MCDM (FMCDM; Sitorus et al., 2019). The outcome of Crisp MCDM methods usually relies significantly on decision-makers’ subjective judgments and preferences (Kabir & Ahsan Akhtar Hasin, 2011). Applying the fuzzy set theory in MCDM enables the incorporation of imprecisions associated with the vague judgments of experts (Chen et al., 2011).

The fuzzy analytical hierarchy process (FAHP) technique, as a part of FMCDM methods, can be viewed as an advanced analytical method extended from the traditional AHP developed by Saaty (1987). Generally, it is impossible to reflect the decision makers’ imprecise preferences through crisp values. Chang (1996) introduced an FAHP approach that employs triangular fuzzy numbers (TFNs) for the pairwise comparison scale instead of crisp values. FMCDM methods have been frequently used in water and wastewater infrastructure management.

Noori et al. (2021) used a grouped FMCDM method for water supply choice optimization. Comparing their proposed method with conventional MCDM methods, their results indicated that utilizing FMCDM methods stands more consistent with the local experts’ opinions. Dewalkar and Shastri (2022) evaluated different wastewater treatment systems by integrating life cycle cost parameters in an FMCDM framework. Results from their study showed that the experts mostly prefer natural treatment systems. Thus, FMCDM methods have proven more effective in the context of prioritizing investments and improvement measures for urban water and wastewater infrastructures.

One widely used framework for measuring urban communities and infrastructures’ attribute-based resilience has been proposed by Bruneau et al. (2003), which is based on quantifying the system’s resilience dimensions, including rapidity, robustness, resourcefulness, and redundancy (four Rs). Rapidity is the system’s ability to achieve the expected level of performance in the shortest possible time, while robustness is the system’s resistance to maintaining its function against stressors. Furthermore, resourcefulness is defined as the ability to employ material and human resources to meet the expected level of service, and redundancy is the units available for replacement in different parts of the system to ensure its continued operation (Bruneau & Reinhorn, 2007; Schoen et al., 2015). By identifying the flood resilience attributes of WWTPs, a resilience quantification framework built on the preference ranking organization method for enrichment evaluation (PROMETHEE) was designed by Karamouz et al. (2018). Their findings showed how resilience could be increased by improving the values of factors that affect resilience, mainly through financial allocation.

Conversely, the performance-based method is generally applied by monitoring the performance of a simulated system after a disruption with the help of a defined resilience metric (Hosseini et al., 2016; Sun et al., 2020). However, developing a performance-based metric can be complex and require a significant amount of data (Vugrin et al., 2017). Thus, utilizing this method for simulating the preparedness and response of a complex system such as a WWTP to assess its resilience can create complications and could be time-consuming. Sweetapple et al. (2022) suggested that performance-based resilience assessment methods can aid in shifting the emphasis from the hazard itself to the system process by providing an insight into system characteristics that increase its resilience. In this regard, interdependencies can bring about further disruptions in urban systems’ functioning because of the cascading failure effects (Guidotti et al., 2016).

Rinaldi (2004) examined the critical interdependencies between urban infrastructures by citing power outages in early 2001. He concluded that identifying and analyzing interdependencies between infrastructures to manage them effectively is essential as these infrastructures become more centralized and complex. Approaches developed to date for modeling the performance of interdependent infrastructure systems following disruptive events are generally placed into five categories of agent-based, system dynamics-based, network-based, economic theory-based, and others (Ouyang, 2014). Employing Ouyang’s classification, the topology-based simulation network, one of the implementations of the network-based method, has been selected for the presented study to derive a performance-based resilience metric. Network modeling, which follows the graph theory, is a suitable method for simulating an interdependent infrastructure system since most urban systems are composed of networks of physical assets (Sinha, 2013).

Water supply systems, power grids, communication infrastructures, and their patterns were modeled using network theory and numerical analysis by Kong and Simonovic (2016). Based on three network resilience features, a new resilience metric was represented that could improve the decision-making processes. He and Cha (2020) classified the interdependencies among urban infrastructures into three levels of system-to-system, system-to-facility, and facility-to-facility, depending on the resolution. Using a dynamic network-based model, they concluded that recovery periods would be underestimated if interdependencies between systems were ignored. Based on the reviewed literature, the role of interdependencies in the performance of WWTPs during and after natural hazards is not widely investigated.

According to Henry and Emmanuel Ramirez-Marquez (2012), a generic resilience assessment must integrate both methods mentioned above to consider measurements for avoiding disruptions and facilitating the recovery from disruptions. In this regard, and based on the concept of resilience, an FMCDM process is proposed to develop an attribute-based metric to increase the WWTPs’ ability to withstand the additional load of a storm through financial allocation. Also, a network model is utilized to devise a performance-based metric to improve the recovery rate of WWTPs’ serviceability after the disaster by prioritizing the restoration of interdependent assets. The proposed resilience improvement measures are chosen by considering the relevant features of an AMP (allocation of financial resources and operation and maintenance planning). It should be noted that other resilience improvement measures, such as the application of nature-based solutions (Lafortezza et al., 2018), monitoring asset performance against predefined service delivery targets (OECD, 2021), and increasing existing components’ interconnections (Shin et al., 2018), are not considered in this study. Imprecisions in linguistic terms are handled in the first method, and uncertainties in the state of damaged facilities are addressed in the latter. Consequently, different aspects of infrastructure characteristics and interdependencies are considered as both parts of a financial allocation scheme and a prioritized asset restoration plan. In the following sections, the methodology is discussed thoroughly. Then, different specifications of the study area are presented, followed by discussions of the results from applying the methodology to the case study. As the last step, a summary is provided, and conclusions are drawn.

2. Methodology

The methodology of this study consists of identifying different systems characteristics and data preparation at first. An attribute-based resilience metric is then obtained by combining the fuzzy set theory with an already established MCDM framework by Karamouz and Hojjat-Ansari (2020) to determine each WWTPs share when allocating funds. In order to simulate the process of spreading interference in the serviceability of network components and the subsequent recovery of failed components after a flood event, the Interdependency Network Model (INM) is developed as a performance-based assessment method. Eventually, the restoration of INM components is prioritized based on maximizing the resilience of WWTPs. The proposed framework provides a practical solution to achieve a more resilient AMP for WWTPs at the risk of flooding.

The flowchart of the proposed framework for financial allocation based on improving the interdependent WWTP’s resilience before the disaster and the prioritized asset restoration for post-disaster recovery is shown in Figure 1. The flowchart consists of a data preparation section covering topological, physical, and systems characteristics, two resilience quantification methods, and a final segment for improving each resilience metric.

Figure 1. Proposed framework for WWTP’s flood resilience improvement by financial allocation and prioritized asset restoration.

2.1. Attribute-based resilience assessment

For quantifying an attribute-based resilience metric, the WWTPs’ attributes affecting their resilience already identified by Karamouz and Hojjat-Ansari (2020) will be utilized. Therefore, the identified attributes (hereafter called sub-criteria) should be categorized into the four Rs criteria. The corresponding sub-criteria associated with the WWTP’s resilience are presented in Figure 2. Additional descriptions of the sub-criteria are provided in Table A1 of Appendix A. Rapidity and robustness criteria each consist of eight sub-criteria. Resourcefulness has three sub-criteria while four sub-criteria are assigned to redundancy. As mentioned before, the attribute-based metric is developed to improve a WWTP’s resilience prior to and during the flood. Therefore, attributes (sub-criteria) only affecting its performance in resisting the flood load are addressed by five additional sub-criteria representing WWTPs’ interdependencies to power infrastructure.

Figure 2. Sub-criteria defined for quantifying WWTPs’ flood resilience (Adapted from Karamouz & Hojjat-Ansari, 2020).

2.1.1. Imprecise judgments handling

The AHP method has been extensively used for assigning weights for the sub-criteria based on expert knowledge. However, due to the subjectivity inherent in human judgment, imprecision inevitably enters this process (Xiao, 2019). Therefore, the conventional AHP method cannot handle the ambiguity associated with mapping one’s judgment to a number (Liu et al., 2020). Incorporating fuzzy set theory into the framework of an MCDM technique will assist the decision-maker when transforming vague linguistic expressions made based on subjective information into quantitative assessments (Kaya et al., 2019).

In this regard, Chang (1996) introduced a new AHP approach by combining it with the fuzzy set theory. This method considers the imprecisions associated with the assessment of linguistic terms by quantifying the pairwise comparison scales with the use of TFNs. First, the scales of relative importance used in the pairwise comparison matrix will be defined. Nine TFNs are defined with respect to the classical scale system. Each TFN can be represented with $M=\left(l,m,u\right)$, where $l\le m\le u$, $l$ and $u$ stand for the lower and upper value of the support of, respectively, and $m$ for the modal value. The linguistic terms for describing the levels of importance for scales 1 through 9 are defined in Table 1. Corresponding TFNs of each linguistic value for pairwise comparison are also listed in Table 1.

Table 1. 9-level fuzzy scale for the pairwise relative importance.

Crisp scale	Linguistic scale	Triangular fuzzy scale
1	No difference	(1,1,1)
2	Intermediate	(1,2,3)
3	Moderate importance	(2,3,4)
4	Intermediate	(3,4,5)
5	Strong importance	(4,5,6)
6	Intermediate	(5,6,7)
7	Extreme importance	(6,7,8)
8	Intermediate	(7,8,9)
9	Outmost importance	(9,9,9)

In the second step, pairwise comparison matrices are generated concerning the experts’ opinions. The geometric mean method is used to integrate the fuzzy pairwise comparison matrices from different experts (Buckley, 1985). The fuzzy geometric mean value of $n$ criteria can be calculated as

(1) $\bar{M}=\left(\bar{l},\bar{m},\bar{u}\right)=\left({\left(\prod _{i=1}^n{l}_i\right)}^{1/n},{\left(\prod _{i=1}^n{m}_i\right)}^{1/n},{\left(\prod _{i=1}^n{u}_i\right)}^{1/n}\right).$(1)

Consider the TFNs $M_1$ and $M_2$, where $M_1=\left(l_1,m_1,u_1\right)$ and $M_2=\left(l_2,m_2,u_2\right)$. The operational laws for adding, multiplying, and inverting TFNs are done using Equations 2(2) $M_1\oplus M_2=\left(l_1+l_2,m_1+m_2,u_1+u_2\right),$(2) to 4 as

(2) $M_1\oplus M_2=\left(l_1+l_2,m_1+m_2,u_1+u_2\right),$(2)

(3) $M_1\otimes M_2=\left(l_1{l}_2,m_1{m}_2,u_1{u}_2\right),$(3)

(4) $M^{-1}\approx \left(\frac{1}{u},\frac{1}{m},\frac{1}{l}\right).$(4)

The fuzzy weight of the $i^{th}$ criterion can then be calculated based on Equations 1(1) $\bar{M}=\left(\bar{l},\bar{m},\bar{u}\right)=\left({\left(\prod _{i=1}^n{l}_i\right)}^{1/n},{\left(\prod _{i=1}^n{m}_i\right)}^{1/n},{\left(\prod _{i=1}^n{u}_i\right)}^{1/n}\right).$(1) through 4(4) $M^{-1}\approx \left(\frac{1}{u},\frac{1}{m},\frac{1}{l}\right).$(4) as

(5) ${\tilde{w}}_i={\bar{M}}_i\otimes {\left({\bar{M}}_1\oplus {\bar{M}}_2\oplus \dots \oplus {\bar{M}}_n\right)}^{-1}.$(5)

In order to investigate the effect of imprecisions in the quantification of resilience, the crisp weights of the identified sub-criteria are also calculated. For this purpose, experts’ linguistic opinions are transformed to the crisp scales mentioned in Table 1 and combined according to the conventional AHP method (Karamouz & Hojjat-Ansari, 2020). Comparing the crisp resilience values with the fuzzy results can highlight the practicality of the FAHP method and assist the decision-maker in selecting a range of variations as far as the actual state of resilience is concerned.

2.1.2. Fuzzy resilience metric

For quantifying resilience, the PROMETHEE method, which is an efficient decision-making method for a limited set of alternative measures, can be used (Behzadian et al., 2010). PROMETHEE uses geometrical analysis to combine subjective weights and values in order to compare the alternatives with each other. A detailed description of the applied approach to quantify resilience in a deterministic manner is explained in Karamouz et al. (2016).

In this study, the PROMETHEE method is extended to the fuzzy environment adapted from Tong et al. (2020).

Let $A=\left\{A_1,A_2,\dots ,A_m\right\}$ be a set of alternatives, $D=\left\{D_1,D_2,\dots ,D_l\right\}$ be a set of the decision-makers and $C=\left\{C_1,C_2,\dots ,C_n\right\}$ be a set of criteria concerning their fuzzy weights $\tilde{W}=\left\{{\tilde{W}}_1,{\tilde{W}}_2,\dots ,{\tilde{W}}_n\right\}$. Here, only the weights of the FMCDM problem are treated as TFNs. Therefore, the proposed fuzzy PROMETHEE (FPROMETHEE) method has no differences from the conventional method until calculating the aggregated preference function. The fuzzy aggregated preference index, $\tilde{\pi }\left(x_i,x_k\right)$ can be calculated as

(6) $\tilde{\pi }\left(x_i,x_k\right)=\left[\sum _{j=1}^n{\tilde{w}}_j{P}_j\left(x_{ij},x_{kj}\right)\right]/\sum _{j=1}^n{\tilde{w}}_j,$(6)

where $i\ne k$, $i,k=1,2,\dots ,m$, and $x_i$ and $x_k$ are given preferences by the decision-makers for $i^{th}$ and $k^{th}$ alternatives, respectively. $P_j\left(x_i,x_k\right)$ is the indifference interval linear priority rule preference function over the $j^{th}$ criterion. Then, the fuzzy leaving (positive) flow, ${\tilde{\phi }}^{+}\left(x_i\right)$ and entering (negative) flow, ${\tilde{\phi }}^{-}\left(x_i\right)$ for the $i^{th}$ alternative can be obtained respectively as

(7) ${\tilde{\phi }}^{+}\left(x_i\right)=\frac{1}{m-1}\left[\sum _{k=1}^m\tilde{\pi }\left(x_i,x_k\right)\right],$(7)

(8) ${\tilde{\phi }}^{-}\left(x_i\right)=\frac{1}{m-1}\left[\sum _{k=1}^m\tilde{\pi }\left(x_k,x_i\right)\right],$(8)

where ${\tilde{\phi }}^{+}\left(x_i\right)$ represents the degree to which alternative $x_i$ dominates all other alternatives. Contrastingly, ${\tilde{\phi }}^{-}\left(x_i\right)$ represents the degree to which alternative $x_i$ is dominated by all the other alternatives (Abdullah et al., 2019). The net outranking flow of the $i$ th alternative, $\tilde{\phi }\left(x_i\right)$, will be

(9) $\tilde{\phi }\left(x_i\right)={\tilde{\phi }}^{+}\left(x_i\right)-{\tilde{\phi }}^{-}\left(x_i\right)=\left(l_i^{\phi },m_i^{\phi },u_i^{\phi }\right).$(9)

The fuzzy resilience metric for each alternative (WWTP) is calculated by transferring the values of the net flow to a range of $\left[0,100\right]$ to express the resilience metric in percentage. Lower and upper values of the fuzzy resilience metric can be interpreted as the minimum and maximum values for the resilience of $i^{th}$ WWTP. Rather than aiming to be analytically precise, this method provides relative and indicative measures for a resilience comparison between similar systems. Eventually, crisp resilience metrics are also developed using crisp weights obtained in the previous step and the conventional PROMETHEE II method explained by Karamouz and Hojjat-Ansari (2020). This is done to see how effective the inclusion of imprecision in the experts’ opinions is when quantifying an attribute-based resilience metric.

2.1.3. Resilience improvement

One way to increase this resilience metric is by improving sub-criteria values. For this purpose, economic sub-criteria (ESC), the sub-criteria that can potentially be improved by investment, are identified and listed in Table 2, along with their corresponding marginal cost in improvement.

Table 2. Economic sub-criteria with their corresponding marginal costs in USD. (Adapted from Karamouz & Hojjat-Ansari, 2020).

Economic sub-criteria	Unit	Marginal cost
Plants design capacity (Ra3)	MLD	$2.64 million
Generator capacity (Ra8)	1000 kW	$0. 8 million
Percent of not-at-risk equipment (Ro4)	%	$1 million for sandbagging, $1.75 million for elevating critical equipment, and $1.25 million for waterproofing the critical equipment
Number of technical staff (Rs1)	-	$0.043 million for each annual staff wage ($21 hourly wage)
Availability of dewatering facilities (Rs2)	-	$7 million
On-site storage (Rd3)	m^3	$159 million
Additional generators (Rd4)	2000 kW	$1.6 million

ESCs are then targeted for maximizing the resilience improvement of the $i^{th}$ WWTP for a total budget, $TB,$ according to Equation 10(10) $max\mathrm{\Delta }{R}_i=R_i\left(ES{C}_1,ES{C}_2,\dots ,ES{C}_{NE}\right)-R_i\left({\widehat{ESC}}_1,{\widehat{ESC}}_2,\dots ,{\widehat{ESC}}_{NE}\right)$(10) . Allocation constraints for each ESC are also conditioned in the optimization problem.

(10) $max\mathrm{\Delta }{R}_i=R_i\left(ES{C}_1,ES{C}_2,\dots ,ES{C}_{NE}\right)-R_i\left({\widehat{ESC}}_1,{\widehat{ESC}}_2,\dots ,{\widehat{ESC}}_{NE}\right)$(10)

$Subjectto$

(11) $ES{C}_j=\left\{\begin{array}{c}{\widehat{ESC}}_j+\left[\frac{Allo{c}_2}{M_j}\right]ifj=7\\ {\widehat{ESC}}_j+Allo{c}_j/M_{j}otherwise\end{array}\right.$(11)

(12) $\sum _{j=1}^{NE}Allo{c}_j\le TB$(12)

(13) $Allo{c}_j\ge 0\hspace{0.17em}j=1,2,\dots ,NE$(13)

(14) $Allo{c}_1\le 0.35TB$(14)

(15) $Allo{c}_2\le M_2\left(R{a}_7-R{a}_8\right)$(15)

(16) $Allo{c}_3\le 0.6TB$(16)

(17) $Allo{c}_4\le 0.1TB$(17)

(18) $Allo{c}_5=\left\{\begin{array}{c}{M}_{5}ifR{s}_1=0\\ 0otherwise\end{array}\right.$(18)

(19) $Allo{c}_6\le 0.35TB$(19)

where $\mathrm{\Delta }{R}_i$, is the increase in resilience after any investment in the $i^{th}$ plant. $R_i$, is determined from FPROMETHEE. $Allo{c}_j$, is the part of $TB$ that is dedicated to $ES{C}_j$, and ${\widehat{ESC}}_j$, is the value of $ES{C}_j$ before the investment. Also, according to Table 2, $M_j$, is the cost corresponding to the unit increase of the ESC, and $NE$, is the number of distinguished ESCs, which is equal to six. The incremental improvement of sub-criteria values after financial allocation is shown in Equation. 11(11) $ES{C}_j=\left\{\begin{array}{c}{\widehat{ESC}}_j+\left[\frac{Allo{c}_2}{M_j}\right]ifj=7\\ {\widehat{ESC}}_j+Allo{c}_j/M_{j}otherwise\end{array}\right.$(11) . With Equations 14(14) $Allo{c}_1\le 0.35TB$(14) and 19(19) $Allo{c}_6\le 0.35TB$(19) , the total allocation to a plant to increase its treatment capacity and on-site storage is limited to 35%. In order to compensate for the power shortage, the on-site generators’ capacity will be increased according to the difference between the power consumption of the plant and its capacity (Equation 15(15) $Allo{c}_2\le M_2\left(R{a}_7-R{a}_8\right)$(15) ). Because WWTPs typically have 1,000–2,000 KW generators, according to Table 2, it is assumed that for every 2,000 KW increase in capacity, one generator is added, which costs $1.6 million. This sum has already been allocated to the capacity of the generators. Equation 16(16) $Allo{c}_3\le 0.6TB$(16) determines the allocation limit for sandbagging, waterproofing, and elevating critical equipment. These three items are the measures that need to be taken in order to increase the percentage of not-at-risk equipment during flooding. Equation 17(17) $Allo{c}_4\le 0.1TB$(17) is used to determine the maximum increase of the personnel. Equation 18(18) $Allo{c}_5=\left\{\begin{array}{c}{M}_{5}ifR{s}_1=0\\ 0otherwise\end{array}\right.$(18) indicates that the investment needed to equip a plant with dewatering facilities is $7 million, assuming the plant does not already have these facilities. It should be noted that Equations 14(14) $Allo{c}_1\le 0.35TB$(14) through 19(19) $Allo{c}_6\le 0.35TB$(19) are defined using preliminary information. For the actual implementation of the proposed methodology, specialists from each system can do a more thorough analysis to properly define these limits for specific ESCs.

Initially, the resilience improvement of each plant is evaluated for different figures through optimization. Total budgets (TB) of $5, $10, $50, $100, and $200 million are allocated to each plant to obtain their corresponding resilience increase separately in order to form their utility functions (UFs). Here, the UF is defined as the rate of increased resilience resulting from financial allocation. Three UFs representing the fuzzy improved resilience metric’s lower, middle, and upper values are fitted to the incremental allocation increase against the percentage of resilience improvement. In order to attain a single UF for each WWTP, Equation 20(20) $R_i=\frac{l_i+4\times m_i+u_i}{6}$(20) is utilized to defuzzify the TFNs (Wu et al., 2017).

(20) $R_i=\frac{l_i+4\times m_i+u_i}{6}$(20)

2.2. Financial allocation

In this step, the derived UFs for each WWTP are employed to determine the final budget allocation by maximizing the Nash product (MNP). When reaching an agreement is challenging in allocation problems, the solution approaches the Nash bargaining solution with utilities that reflect the incentive to settle and have a determined disagreement vector (Binmore et al., 1986). The unique solution for the presented allocation problem based on maximizing Nash product can be achieved by

(21) $maxZ=\prod _{i=1}^N\left(U_i-t_i\right)$(21)

(22) $Subjectto$(22)

(23) $U_i>t_i$(23)

(24) $A{B}_i\le A{B}_i^{\ast }$(24)

(25) $A{B}_i\ge 0$(25)

where N is the number of WWTPs, $U=\left(U_1.U_2.\dots .U_N\right)$ is the vector of systems’ UF, which is a function of $A{B}_i$, the allocated money for the $i^{th}$ WWTP, and $A{B}_i^{\ast }$, is the allocation at which the rate of resilience improvement of the $i^{th}$ WWTP becomes insignificant. $t=\left(t_1.t_2.\dots .t_N\right)$ is the disagreement vector. Here, the disagreement factor is assigned to be half of the funding recommended by the WWTP’s management department. In practice, each plant manager can determine the minimum amount of funding it needs as the point of disagreement in the allocation process.

2.3. Performance-based resilience assessment

In this assessment method, a functionality metric that conforms with the study’s objective must be defined, and its behavior monitored through the simulation process of the system’s performance after a hazard occurrence.

2.3.1. Flood inundation

Floodplain delineation and routing are essential in analyzing the vulnerabilities of an infrastructure network. The Gridded Surface Subsurface Hydrologic Analysis (GSSHA) model, developed by the US Army Corps of Engineers (USACE), was used for this purpose. GSSHA is a 2D hydrologic model that uses digital elevation model (DEM) and geographic information system (GIS) tools to illustrate flooded areas and inundation depth. GSSHA is capable of combining the impact of rainfall and storm surge for overland flow simulation during coastal flooding events (Karamouz & Fereshtehpour, 2019). In this model, spatial analysis with a raster format representing data in DEM grid cells with different land-use and soil types is employed to produce the inundation map of the watershed (Downer & Ogden, 2004; Karamouz et al., 2014). The inundation depths at different infrastructures’ components over the storm’s duration will be obtained from the GSSHA model.

2.3.2. Interdependency network modeling

While there have been numerous studies on modeling the aftermath operation of interdependent water infrastructure systems following a disruption, the majority of them have employed separate physical and deterministic damage models working jointly to simulate the recovery process. In order to refine and extend the previous studies, the present study proposes the Interdependency Network Model (INM). The INM further advances interdependent infrastructure modeling under a storm surge condition in a single model that uses a facility-to-facility resolution of interdependencies. Henceforth, this study refers to the network components in the INM facilities in the wastewater and power infrastructures as assets in the infrastructure network.

The proposed model is a topology-based network that simulates the interactions of network components based only on their topologies and physical connections. Typically, the INM is comprised of nodes and links. Nodes indicate source/sink facilities (assets) connected by links that allow services, i.e., power or wastewater, to flow between them. The links also define the dependency and connection between the ends of the two nodes (Almoghathawi et al., 2021).

Nodes. Each node in the INM has two indices; operability and serviceability. The operability of a network node is its own independent functionality, which is determined based only on the node’s condition in a binary manner. If node $i$ is considered failed at time step $t$, its operability, $O_{i,t}^{Node}$, will be equal to zero, and if it is considered functional, it will equal one. The serviceability of node $i$ at time step $t$, $S_{i,t}^{Node}$, is defined as its actual output ratio relative to its desired service level (when it performs at its routine level). Thus, the serviceability of a sink node depends on its source node(s) supplying the desired service levels and its own operability. The damage and recovery of network nodes over time are measured using operability and serviceability indices.

Links. Each network link in the INM only has an operability index. The operability of a link between node $i$ and node $j$ at time step $t$, $O_{i,j,t}^{Link}$, represents the functioning physical connection that allows services to be sent from node $i$ to node $j$. Each link will have a significance factor for its sink node in the form of an interdependency weight. This factor represents the dependence of the serviceability of the sink node on its source node’s serviceability. Thus, the sum of the link weights reaching a node should equal one. The adjacency matrix summarizes the interdependency weights of all the network links.

(25) $A=\left[\begin{array}{ccc}\begin{array}{c}\begin{array}{cc}0& A_{1,2}\end{array}\\ \begin{array}{cc}{A}_{2,1}& 0\end{array}\end{array}& \cdots & \begin{array}{c}{A}_{1,n}\\ A_{2,n}\end{array}\\ ⋮& \ddots & ⋮\\ \begin{array}{cc}{A}_{n,1}& A_{n,2}\end{array}& \cdots & A_{n,n}\end{array}\right]$(25)

where $n$ is the number of all the nodes in the INM. Each element of the adjacency matrix, A_ij, determines the relative importance of the link between source node i and sink node j in providing the desired service level of node j, relative to other links connected to node j. Figure 3 illustrates an example network comprised of two substations and one pumping station that provide the required power for wastewater delivery and WWTP operation. When the second substation fails due to the impact of a hazard ($O_{2,1}^{Node}=0$), serviceabilities of the dashed nodes are affected.

Figure 3. An example network before and after the failure of one of its nodes.

2.3.3. System’s post-disaster behavior

A system’s post-disaster behavior can be simplified and illustrated in Figure 4. Prior to a hazard, a system’s serviceability is at a routine level (100%). When a hazard occurs at time $t_0$, as illustrated in Figure 4, the system’s serviceability drops due to the damaged components in its network. This decrease in the system’s serviceability happens during the resistance phase. In this phase, the system is still in shock, cascading failures due to interdependencies have affected the system, and the recovery process is not yet started. Then, the system enters the reaction phase at time, during which the essential damaged components of the system (crucial nodes and links) should be restored as soon as possible to bring the overall serviceability up to the minimum operational capability (MOC) at time $t_2$. This is the absolute minimum level of operation, and it is not sustainable. MOC in a WWTP can be defined as when the wastewater is only primarily treated, and further treatment will require restoration of its assets.

Figure 4. Post-disaster behavior and the system’s recovery curve. (Adapted from Bristow, 2015).

The response phase begins once MOC is achieved. This phase includes a series of system reconstruction and restoration measures required to achieve a minimum sustainable capability (MSC) at time $t_3$. This is the level of serviceability at which the system starts to be compliant and stable (Hay, 2016). For a WWTP, MSC can be defined as when the system behaves sustainably by treating the wastewater to the partial secondary level. Defining MOC and MSC based on the treatment level of WWTPs is used as a standard for evaluating the accuracy of the simulated recovery process.

The rehabilitation phase is the system’s return to the routine serviceability level. At the end of this phase ($t_4$), the system’s serviceability is restored to the pre-disaster state, and each component’s serviceability is recovered to 100%.

A system’s performance behavior after a hazard consists of the abovementioned phases. This behavior can be monitored by the network’s serviceability versus the time curve (hereafter called the recovery curve). The recovery curve of an asset shows the infrastructure’s initial resistance and its subsequent service improvement process for returning to its pre-disaster condition. Thus, it can be a measure of the asset’s resilience (He & Cha, 2018). According to Reed et al. (2009), the time average of the area under the recovery curve of node $i$, can be used as a metric of its resilience, $Re{s}_i^{Node}$ as expressed in Equation 26(26) $Re{s}_i^{Node}=\frac{{\sum }_{t=t_0}^{t_4}{S}_{i,t}^{Node}}{t_4-t_0}.$(26) .

(26) $Re{s}_i^{Node}=\frac{{\sum }_{t=t_0}^{t_4}{S}_{i,t}^{Node}}{t_4-t_0}.$(26)

A system’s recovery curve also provides a general overview of its four Rs properties based on the findings of Kong et al. (2019).

2.3.4. Network modeling of resilience

Before the hazard’s occurrence, the INM functions at the routine level, and the operability and serviceability of all network components are equal to one and 100%, respectively. When the storm occurs, a certain number of damaged assets in the INM will be considered failed (their initial operability is changed to zero). The initial operability of a specific component will be determined by considering the uncertainties both in the random nature of the hazard and the design and operation of the assets (discussed in section 2.3.3). Considered failed components are then removed from the network. This procedure will result in the network’s structure to alter accordingly, resulting in a modified adjacency matrix. The serviceability of an operating network node can be lessened in the following time steps due to the cascading failure effect.

After the calamity, a failed component would be recovered and its operability restored to one when a restoration task is performed on the asset for the required time. In order to prioritize the restoration of INM components, only one restoration task was assigned to each infrastructure type at a particular time (one task on wastewater and one for power infrastructure). The recovery coefficients, $R_i^{Node}$ and $R_{i,j}^{Link}$, determine the number of time steps needed for node $i$, and link $i,j$ to be restored completely, respectively. Therefore,

(27) $if\hspace{0.17em}{O}_{i,t}^{Node}=0,\hspace{0.17em}then\hspace{0.17em}{O}_{i,t+R_i^{Node}}^{Node}=1,$(27)

(28) $if\hspace{0.17em}{O}_{i,j,t}^{Link}=0,\hspace{0.17em}then\hspace{0.17em}{O}_{i,j,t+R_{i,j}^{Link}}^{Link}=1.$(28)

When a network component is restored, it will be reinstated in the network structure. Then, the remodified adjacency matrix will be used to obtain the updated serviceability of all INM nodes based on the new network structure. Hence, the serviceability of a node $j$ in time step $t$, $S_{j,t}^{Node}$, is defined as

(29) $S_{j,t}^{Node}=\left[\sum _{i=1}^n\left(A_{i,j,t}\times S_{i,t}^{Node}\right)\right]\times O_{j,t}^{Node},$(29)

where $n$ is the number of all the nodes in the INM, and $A_{i,j,t}$ is the weight of the link between nodes $i$ and $j$ obtained from the remodified adjacency matrix at time $t$.

Since the resilience of the WWTPs is this study’s focus, the serviceability of these nodes (treating the wastewater and stormwater) will be considered to derive the recovery curve and quantify the resilience of each. It is worth mentioning that at this stage, the utilized sequence for the restoration of failed INM components is determined randomly to resemble the recovery process of interdependent assets when there are no prioritization plans by the utilities. The INM will simulate the recovery process until the return of serviceability of all WWTPs to the normal condition (routine level). Recovery curves are then produced in this step to calculate the resilience values of WWTPs as a baseline when there are now specific restoration strategies.

The genetic algorithm optimization technique is then utilized to determine the optimized sequence of restoration tasks on the failed assets to improve the developed performance-based resilience metrics. In this regard, after quantifying the baseline resilience, the INM will be reset to its condition at the start of the simulation process and commence once more with the same initial operability values for its components. This time, the restoration sequence (RS) will be optimized by maximizing the resilience of all the WWTPs. Consequently, the objective function of this optimization (Equation 30(30) $maxRes=\prod _{i=1}^{N}Re{s}_i,$(30) ) is intended to maximize the product of all the resilience values, $Res$, for the WWTPs.

(30) $maxRes=\prod _{i=1}^{N}Re{s}_i,$(30)

where $Re{s}_i$ is the resilience value of the $i_{th}$ WWTP obtained from Equation 26(26) $Re{s}_i^{Node}=\frac{{\sum }_{t=t_0}^{t_4}{S}_{i,t}^{Node}}{t_4-t_0}.$(26) , and $N$ is the total number of WWTPs. After the optimization, corresponding recovery curves for the WWTPs are developed again, and the improved resilience values are calculated accordingly. Also, the optimized RS will be used to classify the restoration priority of different assets in later steps.

2.3.5. Uncertainty in the state of network components

When the hazard occurs, the initial operability of damaged INM components will be changed to zero according to their failure probabilities. Since many factors affect an asset’s performance, the resulting uncertainties need to be modeled when determining the failure probability of an asset in a performance-based method (Olyaei et al., 2018). In the case of INM nodes (except for WWTPs), their failure probabilities are determined based on the maximum inundation depth at their location during the storm’s period, and existing fragility curves (FC) utilized from previous studies. The FC of an asset relates its probability of failure to the severity of the hazard while considering uncertainties in its capacity (vulnerability of its components) and the demand (inundation depth or high-speed winds; Karamlou & Bocchini, 2015; Ma et al., 2020). Because a WWTP is a complex system, it is necessary either to quantify its failure probability using a more precise approach or to assume that it is either zero or 100% based on accounts of its failure during similar hazards in the past.

Links, unlike nodes, have lengths, and they would fail if they sustain physical damage at any point along their length (He & Cha, 2018). Consequently, the failure probability of the INM transmission lines is determined relative to their length. Overhead power transmission lines are created by connecting electric conductors through distribution poles. A damaged distribution pole or conductor will result in the failure of the entire transmission line. For the sake of simplicity, the failure probability of conductors is not considered in this study. Therefore, transmission lines are divided into segments equal to the number of distribution poles, and their failure probability can be calculated based on Equation 31(31) $P{r}_{tr}=1-{\left(1-P{r}_{pole}\right)}^p,$(31) proposed by He and Cha (2018),

(31) $P{r}_{tr}=1-{\left(1-P{r}_{pole}\right)}^p,$(31)

where $P{r}_{tr}$ is the failure probability of a transmission line, $P{r}_{pole}$ is the failure probability of a distribution pole obtained from its corresponding FC, and $p$ is the number of poles contained in the transmission line.

After determining the failure probability of all network components, a uniformly distributed random number in the interval $\left(0,1\right)$ is generated for each INM component and compared against its failure probability. Suppose the generated random number is less than or equal to the failure probability of the component. In that case, it will be considered failed, its initial operability will drop to zero, and it will be removed from the INM. Otherwise, it will keep operating after the hazard. If a transmission line is considered failed, the INM will calculate the number of failed poles by comparing the $P{r}_{pole}$ to random numbers separately. According to Ouyang and Dueñas-Osorio (2014), the average distance between two wooden poles is 42 meters. Hence, transmission lines in the INM are divided into segments by poles 42 meters apart to establish the number of poles.

After determining the initial operability of INM components, the recovery process discussed in the previous section is simulated, and the serviceability values of WWTPs over time, in both baseline and optimized states, are saved. Then, the INM is reinitialized to its status before the hazard’s occurrence (all components performing at the routine level). Again, failure probabilities of network components are compared to generated random numbers, and new initial operability values for network components are determined. Newly failed components are then removed from the network, and the recovery process will initiate again until the network is fully restored. This process will be repeated for several realizations with different initial states. Each realization of the INM can result in different network structures and recovery paths. Finally, serviceability values of the WWTPs in each time step are averaged over all the realizations to attain final recovery curves and corresponding resilience values for the WWTPs in both baseline and planned states.

2.4. Prioritized asset restoration

Prioritizing the restoration and maintenance of the damaged assets after a hazard plays a critical role in the resilience of urban systems due to resource and first responder limitation issues. In this step, the optimized RSs attained from each realization are utilized in classifying the assets’ priority in the restoration process. This prioritization will help achieve/keep the desired level of service after the flood event. When each INM realization ends, all failed nodes of an infrastructure type are graded based on their position in the RS, according to Equation 32(32) $\left\{\begin{array}{c}{G}_{i,n,m}^{Node}=I_n^{Node}-P_{i,n,m}^{Node}+1,\left(O_{i,0}^{Node}=0\right).\\ G_{i,n,m}^{Node}=1,\left(O_{i,0}^{Node}=1\right).\end{array}\right.$(32) . Nodes with their operability not affected will receive a grade equal to one.

(32) $\left\{\begin{array}{c}{G}_{i,n,m}^{Node}=I_n^{Node}-P_{i,n,m}^{Node}+1,\left(O_{i,0}^{Node}=0\right).\\ G_{i,n,m}^{Node}=1,\left(O_{i,0}^{Node}=1\right).\end{array}\right.$(32)

where $G_{i,n,m}^{Node}$ is the grade for the $i^{th}$ node in the infrastructure, at the $m^{th}$ realization, $I_n$ is the total number of nodes in the $n^{th}$ infrastructure (wastewater or power), and $P_{i,n,m}^{Node}$, is the position of the $n^{th}$ node belonging to the $n^{th}$ infrastructure, in the optimized sequence of the $n^{th}$ realization.

A similar grading technique shown in Equation 33(33) $\left\{\begin{array}{c}{G}_{i,j,n,}^{Link}=I_n^{Link}-P_{i,j,n,m}^{Link}+1,\left(O_{i,j,0}^{Link}=0\right),\\ G_{i,j,n,m}^{Link}=1,\left(O_{i,j,0}^{Link}=1\right),\end{array}\right.$(33) is also applied to the links of the INM.

(33) $\left\{\begin{array}{c}{G}_{i,j,n,}^{Link}=I_n^{Link}-P_{i,j,n,m}^{Link}+1,\left(O_{i,j,0}^{Link}=0\right),\\ G_{i,j,n,m}^{Link}=1,\left(O_{i,j,0}^{Link}=1\right),\end{array}\right.$(33)

where ${\mathrm{G}}_{i,j,n,m}$ is the grade for the link between the $i^{th}$ and the $j^{th}$ node in the $n^{th}$ infrastructure, at the $m^{th}$ realization, $I_n^{Link}$ is the total number of links belonging to the $n^{th}$ connection type (wastewater pipeline or power transmission line), and $P_{i,j,n,m}^{Link}$ is the position of the link type $n^{th}$, connecting $i^{th}$ and $j^{th}$ nodes together in the optimized RS of the $m^{th}$ realization. Then, based on Equations 34(34) $S{c}_{i,n}^{Node}=\frac{{\sum }_{m=1}^M{G}_{i,n,m}^{Node}}{I_n^{Node}\times M}\times 100,$(34) and 35(35) $S{c}_{i,j,n}^{Link}=\frac{{\sum }_{m=1}^M{G}_{i,j,n,m}^{Link}}{I_n^{Link}\times M}\times 100,$(35) , a final score will be assigned to each INM component.

(34) $S{c}_{i,n}^{Node}=\frac{{\sum }_{m=1}^M{G}_{i,n,m}^{Node}}{I_n^{Node}\times M}\times 100,$(34)

(35) $S{c}_{i,j,n}^{Link}=\frac{{\sum }_{m=1}^M{G}_{i,j,n,m}^{Link}}{I_n^{Link}\times M}\times 100,$(35)

where $S{c}_{i,n}^{Node}$, and $S{c}_{i,j,n}^{Link}$ are the final scores of the $i^{th}$ node and the link connecting $i^{th}$ and $j^{th}$ node together, respectively, and $M$ is the total number of INM realizations. The priority of each INM component in the restoration process is determined based on its score and the classification scheme depicted in Figure 5. Restoration crews can use the derived prioritization after the hazard as a guide to help WWTPs recover their lost serviceability, hence increasing their resilience. This prioritization is configured in a way that factors such as the failure probability of an asset and its role in the interdependent infrastructure can contribute to its priority. The proposed framework can potentially improve the flood resilience of WWTPs by providing better flood resistance attributes of the plant during the flood and facilitating the recovery and serviceability of that after the event (see the elements of the recovery curve in Figure 4).

Figure 5. Classification of INM components according to their score.

3. Case study

The study area selected for testing the proposed methodology should have specific features and available data. These include wastewater and power infrastructure and interdependencies characteristics and technical information, as well as hydrological data and accurate observations of previous storms. Also, the selected area should be at severe risk from coastal storms and have potential resilience improvement for its infrastructures. Therefore, three sewersheds in New York City that were most affected by Superstorm Sandy are selected as the case study of this paper. These are the Coney Island, Owls Head, and Rockaway sewersheds and their corresponding WWTPs. Table 3. lists the WWTPs’ names, capacities, drainage area, receiving waterbody, and the population they serve (Balci & Cohn, 2014).

Table 3. Specifications of the WWTPs in the study area (Retrieved from NYCDEP (New York City Department of Environmental Protection), 2013).

Plant’s Name	Design Capacity (MLD)	Population Served	Receiving Waterbody	Drainage Area (km^2)
Owls Head	454.25	758,007	Upper New York Bay	52.39
Coney Island	416.39	596,236	Jamaica Bay	61.05
Rockaway	170.34	90,474	Jamaica Bay	37.47

The New York City Department of Environmental Protection (NYCDEP) is responsible for operation and maintenance of the mentioned WWTPs as well as the allocation of financial resources to implement resilience improvement measures. Values of the sub-criteria needed for calculating the attribute-based resilience metric are gathered from the reports issued by the New York City’s Office of the Mayor (NYCDEP (New York City Department of Environmental Protection), 2013) and those provided in the study by Karamouz and Hojjat-Ansari (2020). Table 4. lists the values assigned to the identified sub-criteria for the three WWTPs in the study area. After Superstorm Sandy, NYCDEP allocated 41.6 million dollars to the specified WWTP to improve their resilience. It is reexamined in this paper how those funds could have been reallocated to the 3 WWTPs.

Table 4. Values of sub-criteria for WWTPs in the study area.

Sub-criteria	ID	Owl’s Head	Coney Island	Rockaway
Hurricane flood elevation	$\mathrm{R}{\mathrm{a}}_1$ (m)	2.82	3.07	3.47
level of serviceability	$\mathrm{R}{\mathrm{a}}_2$ (%)	0	100	100
Plant design capacity	$\mathrm{R}{\mathrm{a}}_3$ (MLD)	454.25	416.39	170.34
Post-stress recovery	$\mathrm{R}{\mathrm{a}}_4$ (hr)	16	112	180
Population served	$\mathrm{R}{\mathrm{a}}_5$ (#)	758	596	90
Untreated effluent discharge	$\mathrm{R}{\mathrm{a}}_6$ (MCM)	0.144	0.134	44.8
Average electricity consumption	$\mathrm{R}{\mathrm{a}}_7$ (kW hr)	2645	3639	2600
Back-up generator capacity	$\mathrm{R}{\mathrm{a}}_8$ (kW hr)	2750	6733	2000
Additional load	$\mathrm{R}{\mathrm{o}}_1$ (MLD)	454.25	416.39	170.34
Critical flood elevation	$\mathrm{R}{\mathrm{o}}_2$ (m)	4.42	4.72	4.42
Maximum inundation depth	$\mathrm{R}{\mathrm{o}}_3$ (m)	1.22	0.91	2.13
Critical assets not at risk	$\mathrm{R}{\mathrm{o}}_4$ (%)	71	73	62
DMR violations	$\mathrm{R}{\mathrm{o}}_5$ (%)	1.4	2.1	1.9
Damage from hurricane	$\mathrm{R}{\mathrm{o}}_6$ ($)	48.4	84.9	49.3
Flood elevation at unit substation	$\mathrm{R}{\mathrm{o}}_7$ (m)	0.73	1.63	1.74
Power loss period	$\mathrm{R}{\mathrm{o}}_8$ (hr)	8	28	72
Asset operators	$\mathrm{R}{\mathrm{s}}_1$ (#)	68	68	41
Availability of dewatering facilities	$\mathrm{R}{\mathrm{s}}_2$ (-)	0	0	0
Risk avoided for every dollar spent	$\mathrm{R}{\mathrm{s}}_3$ ($)	13.8	22.6	13.1
Underground tunnel system	$\mathrm{R}{\mathrm{d}}_1$ (-)	1	1	1
Availability of neighboring WWTPs	$\mathrm{R}{\mathrm{d}}_2$ (km)	7.9	9	9
On-site storage	$\mathrm{R}{\mathrm{d}}_3$ (m^3)	0.5	0	1.51
Number of back-up generators	$\mathrm{R}{\mathrm{d}}_4$ (#)	4	5	2

The selected study area is a distinct region that still uses an overhead radial power supply system with wooden poles connected to the national grid. During Superstorm Sandy, the study area experienced widespread outages primarily due to downed trees that caused interruptions in the wastewater infrastructure’s performance (Kurtz et al., 2013). Table 5 lists the data sources used in this study to determine the location and technical information of the assets. WWTPs’ locations and their interdependencies with electric power are illustrated in Figure 6. The study area depicted in this figure consists of 12 electric substations, five pumping stations, three WWTPs, and 37 links, including transmission lines and wastewater pipelines between these nodes.

Figure 6. The study area with its wastewater infrastructure and its related power facilities.

Table 5. Hydrological and technical data sources used in the study.

Data	References
Rainfall hyetograph (JFK airport)	https://www.ncdc.noaa.gov/
Soil type, land use, Boundary	https://opendata.cityofnewyork.us
Dem (30 m)	https://earthexplorer.usgs.gov
Wastewater infrastructure	https://openseweratlas.tumblr.com
Power infrastructure	https://hifld-geoplatform.opendata.arcgis.com

During Superstorm Sandy, the Owls Head sewershed was the least affected, while the other two sewersheds in the study area experienced higher levels of inundation due to their location. Hydrological data required for flood inundation modeling by the GSSHA model are also tabulated in Table 5. When simulating a flood hazard, due to the spatial and temporal characteristics and changes in the state of the waves, utilizing gauged storm hydrographs from the nearest station (Battery Park) is not quite applicable due to the spatial variation of storm from south of Manhattan (where Battery station is) to the southern Brooklyn shoreline where the case study is. To solve this issue, four ungauged storm hydrographs around the case study boundaries, developed by Karamouz et al. (2022), are reexamined and utilized.

Based on the results of the GSSHA model and the FCs obtained from Sánchez-Muñoz et al. (2020) and Yu et al. (2021), the failure probabilities of assets within the network infrastructure are calculated. The sources used in this study for obtaining failure probabilities of INM components (with the exception of wastewater pipelines), along with their recovery coefficients for simulating the recovery process of the INM, are summarized in Table 6.

Table 6. Sources used for determining the failure probability and recovery coefficients of different INM components.

	Failure Probability	Recovery Coefficient
WWTP	NYCDEP (New York City Department of Environmental Protection; 2013)	NYCDEP (New York City Department of Environmental Protection; 2013)
Substation	FC (Sánchez-Muñoz et al., 2020)	Depth-damage, and damage-time curves (Sánchez-Muñoz et al., 2020)
Pumping stations	FC (Yu et al., 2021)	Depth-damage functions (FEMA, 2013) Damage-time functions (FEMA, 2003)
Transmission line	Averaged FCs (Salman & Li, 2016)	~N (3.5,2.5; Ouyang & Dueñas-Osorio, 2014)

As mentioned earlier, the failure probabilities of WWTPs are assumed to be zero or 100% due to their complexity. Because the serviceabilities of Coney Island and Rockaway WWTPs were interrupted during Superstorm Sandy, their failure probability was considered 100% (NYCDEP (New York City Department of Environmental Protection), 2013). In contrast, the failure probability of Owls Head plant was assumed to be zero. A typical flood FC for an electrical substation is illustrated in Figure 7. Utilized FCs from literature are derived from HAZUS’s FEMA (2013) flood model and are applicable to the study area. The failure probability of electric transmission lines in the INM was estimated based on Equation 31(31) $P{r}_{tr}=1-{\left(1-P{r}_{pole}\right)}^p,$(31) , the wind gust of 38 m/s in the study area during Superstorm Sandy (Blake et al., 2013), and the findings of Salman and Li (2016). They calculated the fragility curves of various aged wooden poles under extreme wind speeds. These curves were averaged for the same wind speed values over different poles’ ages to provide a generic fragility curve for the presented study.

Figure 7. For a typical electrical substation: a) Fragility curve, b) water depth versus damage, and c) repair time versus damage (From Sánchez-Muñoz et al., 2020).

The recovery coefficients of electrical substations in the INM were determined according to the inundation depth versus damage and damage versus repair time curves of this facility type, according to Figure 7. The recovery coefficients of WWTPs were obtained based on the actual recovery time of these plants after Superstorm Sandy. The pumping stations’ recovery coefficients were estimated, similar to substations, using FEMA (2013), FEMA (2003), and FEMA (2003) earthquake model was utilized since, to date, no studies have been identified that contain damage versus repair time functions of a pumping station in the event of a flood. The recovery coefficients of distribution poles were considered to follow a normal distribution with a mean of 5 hours and a variance of 2.5 hours (Ouyang & Dueñas-Osorio, 2014).

4. Results

The proposed methodology provides a practical solution to increase the resilience of WWTPs in two separate steps. First, by allocating funds proportional to the resilience improvement potential of WWTPs prior to the flood. Second, by prioritizing the restoration of the damaged assets in the interdependent infrastructure after the flood event. In the following sections, the results of the proposed methodology are presented.

4.1. Attribute-based resilience assessment

The FAHP technique was utilized along with the conventional AHP to calculate the fuzzy and crisp weights of the identified sub-criteria based on experts’ opinions. The criteria and sub-criteria weights derived by both methods are tabulated in Table 7. As listed, the weights of the four Rs are not the same. Rapidity and robustness were assigned higher weights. Also, crisp weights are close to their modal values in the fuzzy state.

Table 7. Estimated fuzzy and crisp weights calculated for the sub-criteria.

Rapidity (Ra)				Robustness (Ro)				Resourcefulness (Rs)				Redundancy (Rd)
Fuzzy			Crisp	Fuzzy			Crisp	Fuzzy			Crisp	Fuzzy			Crisp
0.211	0.334	0.576	0.347	0.167	0.298	0.487	0.302	0.113	0.172	0.263	0.173	0.124	0.196	0.300	0.178
Ra1				Ro1				Rs1				Rd1
0.082	0.148	0.263	0.148	0.057	0.102	0.176	0.098	0.183	0.298	0.519	0.307	0.183	0.306	0.537	0.315
Ra2				Ro2				Rs2				Rd2
0.045	0.078	0.128	0.072	0.091	0.158	0.292	0.161	0.245	0.447	0.732	0.452	0.141	0.236	0.393	0.238
Ra3				Ro3				Rs3				Rd3
0.087	0.148	0.260	0.148	0.049	0.087	0.148	0.082	0.161	0.255	0.449	0.262	0.120	0.215	0.372	0.216
Ra4				Ro4								Rd4
0.071	0.120	0.214	0.119	0.054	0.095	0.167	0.091					0.144	0.244	0.399	0.244
Ra5				Ro5
0.083	0.136	0.223	0.133	0.099	0.177	0.310	0.178
Ra6				Ro6
0.099	0.175	0.309	0.176	0.108	0.176	0.296	0.176
Ra7				Ro7
0.057	0.099	0.163	0.094	0.057	0.101	0.177	0.098
Ra8				Ro8
0.056	0.095	0.161	0.091	0.060	0.104	0.176	0.100

Utilizing the PROMETHEE framework in fuzzy and crisp states, values and weights for the sub-criteria listed in Tables 4 and 7, respectively, were utilized in Equations (6(6) $\tilde{\pi }\left(x_i,x_k\right)=\left[\sum _{j=1}^n{\tilde{w}}_j{P}_j\left(x_{ij},x_{kj}\right)\right]/\sum _{j=1}^n{\tilde{w}}_j,$(6) –9). Figure 8 illustrates the attribute-based resilience metrics of the WWTPs in both states. The lower and upper values are discerned as their imprecise boundaries of the resilience estimation. Lower, modal, and upper values of the fuzzy resilience metric are shown on the left of each WWTP, while the values on the right represent the crisp resilience metrics (shown with red lines).

Figure 8. Fuzzy and crisp resilience metrics obtained for the WWTPs in the study area.

According to Figure 8, the Owls Head WWTP has the highest resilience, followed by the Coney Island WWTP. The Rockaway WWTP has the lowest resilience due to its flood-prone location and lack of flood resistance equipment. The modal value for the Owls Head WWTP is about 58% which is less than its value in the crisp form (61%), showing experts’ tendency to overestimate the resilience of this plant when not considering the imprecisions. Instead, the modal value of the fuzzy resilience metric for the Rockaway plant is 43%, which is four percent higher than the obtained crisp value. This indicates the potential for underestimating the resilience of this WWTP when using a deterministic method. The fuzzy metric calculated for the Coney Island WWTP has the widest range between its lower and upper values. For this reason, the quantification of its resilience is more imprecise than the other two plants. However, the modal value of its fuzzy resilience (47%) is very close to its crisp value.

Based on identified ESCs (Table 2) and the resilience improvement optimization formulated in Equations 10(10) $max\mathrm{\Delta }{R}_i=R_i\left(ES{C}_1,ES{C}_2,\dots ,ES{C}_{NE}\right)-R_i\left({\widehat{ESC}}_1,{\widehat{ESC}}_2,\dots ,{\widehat{ESC}}_{NE}\right)$(10) through 19(19) $Allo{c}_6\le 0.35TB$(19) , the allocation for the maximum increase in the fuzzy resilience metric of each WWTP was obtained. Optimization was carried out for different total budgets ranging from 5 to $200 million. The resilience improvement values obtained from these optimizations were utilized to plot the utility functions (UFs) for the WWTPs illustrated in Figure 9. A logarithmic UF in the form of $a\ast ln\left(x\right)+b$ was fitted for each of the improved resilience metrics’ lower ($l$), modal ($m$), and upper ($u$) values since a logarithmic expression had the highest correlation with the data points, as is evident from the values of their corresponding coefficients of determination, $R_l^2$, $R_m^2$, and $R_u^2$, respectively. Another logarithmic curve was also fitted to the points corresponding to the crisp improved resilience values ($y_c$) with the coefficients of determination of $R_c^2$. The area between the upper and lower fitted curves represents the imprecision in the outcome of the financial allocation process. Next, by utilizing Equation 20(20) $R_i=\frac{l_i+4\times m_i+u_i}{6}$(20) , the three UFs were defuzzified and transformed into single crisp functions ($y_d$ in Figure 9) to be utilized in the final financial allocation step.

Figure 9. Fuzzy and crisp UFs for the three WWTPs of Owls Head (a), Coney Island (b), and Rockaway (c).

Based on Figure 9a, the defuzzified UF of Owls Head plant is lower than its crisp UF. This indicates the possibility of overestimating the effects of financial resource allocation on improving the resilience of this WWTP if the imprecisions are not taken into account. According to the values of $y_d$ and $y_c$ functions in Figure 9b, Coney Island WWTP has the least potential for resilience improvement. In addition, due to having a wider gap between the lower and upper UFs of this plant compared to the other two, its allocation process is more prone to deal with imprecision. According to Figure 9c, Rockaway WWTP has the highest potential for resilience improvement. Also, since of this plant is close to its lower UF in the fuzzy form, there is a possibility of underestimating the impact of financial resources on the resilience improvement of this facility by the experts. In general, it can be inferred that by increasing the amount allocated to each of these systems, the effect of the allocated financial resources on the improvement of resilience decreases. Simultaneously, imprecisions intensify as the allocated budget increases, resulting in a more ambiguous resilience improvement estimation.

4.2. Financial allocation

The allocated funds were distributed among the three WWTPs in the study area based on the MNP technique mentioned in Equations 21 to 24. The final allocation and corresponding growth in the resilience values obtained from MNP compared to the actual allocation recommended by the NYCDEP (New York City Department of Environmental Protection; 2013), given the $41.6 million of the total funding, are illustrated in Figure 10.

Figure 10. Final financial allocations and their corresponding potential resilience improvement.

The final financial allocation numbers based on the obtained fuzzy and crisp UFs are in good agreement compared to the amount that NYCDEP actually allocated. The allocated amount for the Owls Head WWTP based on its defuzzified UF is lower than its amount from using a crisp UF. The resilience of the Owls Head WWTP in the fuzzy and crisp mode is increased by 6.6% and 7%, respectively, when allocating the values of 12.15 and 10.43 million dollars. The Coney Island plant, which had the lowest potential for resilience enhancement, receives high allocation in both states and has a minor improvement (around 5%). This is because, in the MNP approach, the tendency is to minimize dissatisfaction among alternatives. With its low resilience and high growth potential, the Rockaway WWTP receives the most considerable allocation ($17.04 million) in the crisp state, which increases its resilience only by 9.2%. Meanwhile, the amount of allocated budget based on its defuzzified UF ($14.61 million) is lower than the crisp value and the amount suggested by NYCDEP, which increases this WWTP’s resilience by 11.4%. Due to the underestimation of the Rockaway plant’s resilience when using crisp weights, more funds are allocated to increase its resilience less than when the FMCDM is utilized.

4.3. Performance-based resilience assessment

As the first step in developing a performance metric, the flood inundation map of the study area during Superstorm Sandy was obtained. Ungauged storm hydrographs for different shoreline sections of the study area, along with DEM of 30-meter resolution and other hydrological data mentioned in Table 5, were used in the GSSHA hydrologic simulation model. The GSSHA model was already calibrated for the study area in previous studies. Figure 11 demonstrates the simulated maximum inundation depth during the storm.

Figure 11. The simulated inundation map of the study area obtained from hydrological modeling.

The Rockaway sewershed experienced the most severe flooding, while the Owls Head sewershed was the least affected. Values of the corresponding inundation depths of different nodes of the INM are summarized in Table 8. All of the assets located in the Rockaway sewershed witnessed severe inundations, with the Rockaway WWTP itself undergoing 1.74 meters of water. Meanwhile, the storm did not affect two substations in the Owls Head sewershed and one in Coney Island. With 1.63 meters of water depth at the Coney Island WWTP, it had the highest inundation depth among the interdependent assets in its sewershed. The Owls Head WWTP experienced 0.73 meters of inundation when substations S1 and S2, which provide the power for the whole district, were more affected by the storm.

Table 8. Maximum inundation depth at network nodes under simulated Sandy Superstorm.

Node	Water depth (m)	Node	Water depth (m)
S^a 1	1.32	S 11	1.74
S 2	1.05	S 12	1.1
S 3	0	P^b 1	1.48
S 4	0	P 2	0
S 5	1.2	P 3	1.05
S 6	0.2	P 4	1.92
S 7	0	P 5	1.16
S 8	1.08	WWTP 1	0.73
S 9	1.42	WWTP 2	1.63
S 10	1.22	WWTP 3	1.74

^aelectric substation, ^b pumping station

The study area’s INM, depicted in Figure 12, was developed using MATLAB software. Gathered FCs, along with the estimated inundation depths at network nodes and the maximum wind gust during Sandy, were utilized to determine the failure probability of electric substations, pumping stations, and transmission lines. At the start of each realization, the initial operability of network components was determined based on their corresponding failure probabilities and generated random numbers, with the exception of wastewater pipelines. The failure of wastewater pipelines during a storm will result in combined sewer overflow. Hence, for a wastewater pipeline link to function without failure, both of its endpoint nodes (the pumping station and the WWTP) need to be operational. Therefore, if any node connected to a wastewater pipeline was in the failure state, the initial operability of this link was also considered to be zero.

Figure 12. The INM of the study area developed in MATLAB.

In the next step, network components with zero operability underwent the recovery process according to their corresponding recovery coefficients. Eventually, wastewater pipelines with zero operability were reinstated in the INM as soon as their connected non-operating nodes were restored. Prior to any optimization, the system’s post-disaster recovery process resulting from a random RS was simulated. Serviceability values in each time step were updated according to Equation 29(29) $S_{j,t}^{Node}=\left[\sum _{i=1}^n\left(A_{i,j,t}\times S_{i,t}^{Node}\right)\right]\times O_{j,t}^{Node},$(29) . After the recovery process’s completion, baseline serviceability values of Owls Head, Coney Island, and Rockaway WWTPs over time were determined. Then, the operabilities of INM components were reset to their initial states, and the genetic algorithm was utilized to obtain the optimized RS and the planned serviceability values based on Equation 30(30) $maxRes=\prod _{i=1}^{N}Re{s}_i,$(30) . Finally, the INM was reinitialized to its pre-disaster status.

To account for uncertainties in the state of the affected assets, the recovery process discussed above was simulated for 100 realizations. Based on the size of the INM and the computational effort and time for each simulation, 100 realizations were deemed enough. For actual implementation, 500 realizations are suggested. Each realization commences with comparing the failure probabilities of network components against newly generated random numbers to determine initial operability states. This allows considering various disruption scenarios that result in different recovery paths for the network infrastructure. Eventually, baseline and planned recovery curves for the three WWTP nodes in the INM were obtained by finding the expected value of their serviceability over all realizations in each time step. Figure 13 illustrates the recovery curves and corresponding optimized resilience values for the WWTPs.

Figure 13. The recovery curves and their corresponded resilience metrics obtained from the INM before and after the optimization.

According to Figure 13, the Owls Head WWTP’s serviceability does not drop below 50% at any point in the simulation period, denoting that the primary treatment would not be stopped. In this plant, the MSC is achieved in a shorter time compared to the other two WWTPs. After the optimization, the resilience of this WWTP increased from 66.3% to 89.1%. The optimization improves the resilience of Coney Island WWTP the highest, experiencing a 27.7% increase in its performance-based resilience metric. This plant reaches its MOC and MSC approximately two days sooner when the optimized RS is implemented.

The Rockaway WWTP is entirely offline in the first 53 hours of the simulation period. This can be due to more prolonged and higher levels of inundation in the Rockaway sewershed and the vulnerability of the power infrastructure in this area. Both substations serving this plant were disrupted and damaged severely during the storm (Balci & Cohn, 2014). However, the asset restoration optimization that least affected this plant improved its resilience from 26.7% to 45.5%. The Rockaway WWTP has limited interdependencies, so prioritizing its network components leads to a relatively insignificant outcome for the resilience improvement process.

Overall, MOC and MSC levels show good agreement with reports on the serviceability interruptions of WWTPs during Superstorm Sandy. After the storm, the Owls Head plant’s secondary treatment was only interrupted for approximately five hours. Meanwhile, the Coney Island WWTP was completely offline for 12 hours and regained its routine serviceability after four days. The Rockaway WWTP was severely inundated during Superstorm Sandy. Despite the many preventive measures taken by the plant’s staff and their efforts to recover the equipment, it took two full weeks to return to the routine service level (NYCDEP (New York City Department of Environmental Protection), 2013).

4.4. Prioritized asset restoration

By prioritizing the restoration of failed assets, the WWTPs’ service levels after flooding will be improved noticeably. After each realization, assets within the infrastructure network were graded based on obtained optimized sequences and Equations 32(32) $\left\{\begin{array}{c}{G}_{i,n,m}^{Node}=I_n^{Node}-P_{i,n,m}^{Node}+1,\left(O_{i,0}^{Node}=0\right).\\ G_{i,n,m}^{Node}=1,\left(O_{i,0}^{Node}=1\right).\end{array}\right.$(32) and 33(33) $\left\{\begin{array}{c}{G}_{i,j,n,}^{Link}=I_n^{Link}-P_{i,j,n,m}^{Link}+1,\left(O_{i,j,0}^{Link}=0\right),\\ G_{i,j,n,m}^{Link}=1,\left(O_{i,j,0}^{Link}=1\right),\end{array}\right.$(33) . The final scores of each INM component were calculated according to Equations 34(34) $S{c}_{i,n}^{Node}=\frac{{\sum }_{m=1}^M{G}_{i,n,m}^{Node}}{I_n^{Node}\times M}\times 100,$(34) and 35(35) $S{c}_{i,j,n}^{Link}=\frac{{\sum }_{m=1}^M{G}_{i,j,n,m}^{Link}}{I_n^{Link}\times M}\times 100,$(35) . The restoration priorities of network components were then classified into low, moderate, high, and crucial classes based on the classification scheme in Figure 5. The restoration priority of each network component in the study area is illustrated in Figure 14.

Figure 14. Prioritized network components for restoration in the study area.

Based on Figure 14, the restoration priority of most of the assets in the Rockaway sewershed is classified as crucial. The Owls Head sewershed experiences the lowest inundation. However, due to the high dependency of other components in the INM on the performance of the assets located in this sewershed, two components have been given a relatively crucial priority while the majority are classified as low or medium. Substations S5 and S8 are classified as crucial because of their flood-prone location, while other substations placed at the Coney Island sewershed have been given a low priority. Moreover, Coney Island and Rockaway WWTPs are both classified as crucial, while the Owls Head plant’s restoration priority is low.

In order to validate the proposed scoring scheme, the performance-based resilience metrics are calculated once more based on the final RS. Resilience values for WWTPs Owls Head, Coney Island, and Rockaway were calculated as 83.3%, 61.8%, and 42.1%, respectively. These values show about 90% acceptance with the resilience values obtained from averaging the optimization results.

The overall framework of this study provides a template for improving WWTPs’ resilience during and after flood hazards. Attribute-based and performance-based resilience metrics are developed, and improvement measures are proposed to assist in the planning of resilient AMPs. The Coney Island WWTP had the lowest attribute-based metric resilience improvement, while its post-storm recovery resilience (performance-based metric) could be significantly enhanced. The Rockaway WWTP showed the exact opposite effect. In contrast, the Owls Head WWTP’s resilience improvement was moderately affected by both approaches. These observations demonstrated both methods’ range and applicability in composing resilience enhancement strategies. The results have also benefited from incorporating uncertainties and imprecisions in the modeling process. The findings of this study revealed how WWTP’s resilience could be enhanced by improving the resource allocation schemes of these plants and prioritizing the restoration of the infrastructure network. It provides new insight into the decision-making process by providing practical tools and resilience measures.

5. Summary and conclusion

Treatment failure of coastal WWTPs and interruptions in maintaining their service levels due to floods will result in grave environmental consequences. Therefore, the resilience of these systems has become of great concern. The present study proposes an attribute-based resilience metric that addresses the imprecisions associated with the subjective inputs of experts by utilizing an FMCDM method. Next, a framework has been designed to allocate limited financial resources to improve their resilience. This part of the study looks at the overall resistance of the system in facing an extreme flood event such as Superstorm Sandy.

By considering the interdependencies among wastewater and power infrastructures, a performance-based resilience metric methodology is also calculated to prioritize the restoration of assets within the infrastructure network. This part of the study concentrates on recovery efforts after the storm. Interconnected assets in the study area are classified into four states of urgency to restore. Restoration crews can plan their restoration efforts according to this classification in order to improve the WWTPs’ resilience. The study area selected for implementing this methodology is the flood-prone districts of Southern Brooklyn and Rockaway in New York City.

The fuzzy attribute-based resilience assessment placed the Owls Head WWTP as the plant with the highest resilience metric and the Rockaway as the least resilient plant. $41.6 million that the city had planned to allocate to these plants after Sandy was redistributed among the WWTPs and compared with their actual share. According to this study, the reallocation of $14.6, $14.8, and $12.2 million to WWTPs of Owls Head, Coney Island, and Rockaway, respectively, resulted in increases in their attribute-based metric by 6.6%, 5.1%, and 11.4%. Moreover, by comparing the fuzzy attribute-based method with its corresponding crisp method, it can be stated that utilizing FMCDM methods can show the range of overestimation or underestimation of resilience. Therefore, having the range of resilience variation, decision makers can consider a value in that range depending upon their risk-taking attitude.

The simulated post-disaster behavior of the infrastructure network showed the performance-based resilience metrics of the three WWTPs to be in the same order as was obtained from the attribute-based method. This seems relevant even though the two approaches are not entirely compatible. Recovery curves obtained from 100 INM realizations illustrated the improved resilience values when the restoration of damaged assets was prioritized. The Coney Island WWTP demonstrated the most considerable improvement in its resilience since it is well-connected to other assets in the infrastructure network. By contrast, the Rockaway plant’s resilience was least improved due to the limited number of interdependencies and its flood-prone location, even though network components in this sewershed are classified as either crucial or high to be repaired.

The results of this study show how well the food hazard could be simulated, and the WWTPs’ resilience could be quantified in line with the consequences that NYCDEP documented. In addition, it is demonstrated how the resilience metrics can be used for resource allocation and asset restoration planning. The proposed methodology can be applied to other coastal regions. One limitation of the proposed attribute-based method is the sensitivity of its results to the type of the chosen fuzzy membership function. Therefore, a sensitivity analysis should be performed by considering various membership functions and fuzzy scales for actual implementation of the proposed mythology. The suggested performance-based approach can be enhanced to reflect the recovery process in practice more accurately. One suggestion for enhancement is to take into account other interdependent infrastructures, such as communication and transportation. Results from the INM can also benefit by creating unique FCs for each WWTP. Furthermore, the accuracy of recovery curves can be further enhanced by considering partial operability states in the interval $\left[0,1\right]$ for the assets (especially WWTPs).

Abbreviation

AHP, Analytic Hierarchy Process; AMP, Asset management programs; DEM, Digital Elevation Model; ESC, Economic Sub-Criteria; FC, Fragility Curves; FAHP, Fuzzy Analytic Hierarchy Process; FMCDM, Fuzzy multi-criteria decision-making; GIS, Geographic information system; INM, Interdependency Network Model; MNP, Maximizing the Nash product; MOC, Minimum Operational Capability; MSC, Minimum Sustainable Capability; RS, Restoration Sequence; GSSHA, The Gridded Surface Subsurface Hydrologic Analysis; NYCDEP, The New York City Department of Environmental Protection; TNF, Triangular Fuzzy Number; UF, Utility Function; WWTP, Wastewater Treatment Plants

Disclosure statement

The Coalition for Disaster Resilient Infrastructure (CDRI) reviewed the anonymised abstract of the article, but had no role in the peer review process nor the final editorial decision.

Additional information

Funding

The Article Publishing Charge (APC) for this article is funded by the Coalition for Disaster Resilient Infrastructure (CDRI).

Appendix A

The detailed description of the sub-criteria utilized in this study are tabulated in Table A1.

Table A1. Descriptions of sub-criteria.

Sub-criteria	Description
$\mathrm{R}{\mathrm{a}}_1$	Flood elevation based on North American Vertical Datum of 1988 (NAVD88)
$\mathrm{R}{\mathrm{a}}_2$	Adverse environmental impacts on the surrounding area due to treatment failure
$\mathrm{R}{\mathrm{a}}_3$	Capacity as function of the plant users
$\mathrm{R}{\mathrm{a}}_5$	Number of residents that are served by the plant
$\mathrm{R}{\mathrm{o}}_1$	The difference between WWTP capacity for the total maximum wet and dry weather flow (Maximum wet weather flow is the maximum flow received during any 24 hours. Maximum dry weather flow is the maximum daily flow during periods without rainfall.)
$\mathrm{R}{\mathrm{o}}_2$	100-year flood elevation +30 inches for expected sea-level rise by the 2050s (This is determined based on the FEMA’s 100-year flood elevation maps)
$\mathrm{R}{\mathrm{o}}_3$	Due to the flat terrain of the plant, several areas may be flooded by up to this value of water during the critical flood event
	The percentage of discharge monitoring reports that resulted in effluent violations (The Discharge Monitoring Report (DMR) violation percentages indicate how well each treatment plant can cope with daily operational stresses.)
$\mathrm{R}{\mathrm{o}}_6$	Damage cost from the most severe historical hurricane without flood protection for the plant
$\mathrm{R}{\mathrm{s}}_1$	Number of plant’s technical staff
	Availability of facilities to drain sludge to decrease 90% of its liquid volume
$\mathrm{R}{\mathrm{s}}_3$	Total risk avoided for every single dollar spent on improving operation and maintenance over 50 years
$\mathrm{R}{\mathrm{d}}_2$	Distance from the closest WWTP
$\mathrm{R}{\mathrm{d}}_3$	Volume of lakes in the WWTP’s zone

Gogus and Boucher (1998) showed how the consistency ratio in the FAHP method can be estimated. They also generated random indices using a sample size of 400. The same procedure was used in this paper.