Resilience of gas interchangeability in hydrogen-blended integrated electricity and gas systems: A transient approach with dynamic gas composition tracking

Green hydrogen can be produced by consuming surplus renewable generations. It can be injected into the natural gas networks, accelerating the decarbonization of energy systems. However, with the fluctuation of renewable energies, the gas composition in the gas network may change dramatically as the hydrogen injection fluctuates. The gas interchangeability may be adversely affected. To investigate the ability to defend the fluctuated hydrogen injection, this paper proposes a gas interchangeability resilience evaluation method for hydrogen-blended integrated electricity and gas systems (H-IEGS). First, gas interchangeability resilience is defined by proposing several novel metrics. Then, A two-stage gas interchangeability management scheme is proposed to accommodate the hydrogen injections. The steady-state optimal electricity and hydrogen-gas energy flow technique is performed first to obtain the desired operating state of the H-IEGS. Then, the dynamic gas composition tracking is implemented to calculate the real-time traveling of hydrogen contents in the gas network, and evaluate the time-varying gas interchangeability metrics. Moreover, to improve the computation efficiency, a self-adaptive linearization technique is proposed and embedded in the solution process of discretized partial derivative equations. Finally, an IEEE 24 bus reliability test system and Belgium natural gas system are used to validate the proposed method.

G reen hydrogen is one of the most appealing solutions to the decarbonization of energy systems. It is usually produced by PTGs by consuming surplus renewable generations [1] . The green hydrogen can be then injected into the pipelines, and transported to other locations for further use. Recently, many small trials have been implemented in China, the UK, the US, etc., which demonstrate the feasibility of blending hydrogen into natural gas pipelines (https://new.qq.com/rain/a/20221123A088JI00; https://hydeploy.co.uk/app/uploads/2022/06/HyDeploy-Close-Down-Report_Final.pdf; https://news.bjx.com.cn/html/20220927/ 1257566.shtml). For example, the hydrogen blending demonstration project in Zhangjiakou, China, 2020, is estimated to provide more than 4 million hydrogen to residential users and vehicles per year, which can reduce approximately 3000 t carbon emissions (https://news.bjx.com.cn/html/20200917/1105156.shtml).
However, excessive injections of hydrogen can affect the gas compositions, and consequently jeopardize gas interchangeability. Interchangeability is used to describe whether two gases are interchangeable. If the new gas mixtures can be used to substitute the original natural gas without affecting the operation of gas appliances, then the interchangeability of the new gas mixture can be considered qualified [2] . Due to the lower heat value of hydrogen, if too much hydrogen is blended, the new gas mixture cannot produce the same amount of heat energy compared with the original gas. Then, the performance of gas appliances, such as gas water heaters, can be affected. This issue is more severe in the case of green hydrogen. Because renewable generation is stochastic, the gas interchangeability of gas will become more unpredictable [3,4] . Therefore, it is essential to closely monitor and regulate gas interchangeability in the presence of green hydrogen.
Since gas interchangeability is dominated by gas composition, some researchers are dedicated to the gas composition tracking problem. For example, the steady-state gas composition simulation model with the injection of alternative gas is studied in Ref. [5]. The gas composition tracking problem is combined with the electricity system operation in Ref. [6] to reach a global optimum. The impacts of hydrogen produced by distributed photovoltaic generation on the gas system are investigated in Ref. [7]. The probabilistic gas flow with multiple gas types is established in Ref. [8] considering uncertainties. Though the gas composition can be simulated or optimized in these studies, the impacts of hydrogen on gas interchangeability are not quantified. As a result, we still do not know if the gas composition is qualified in terms of interchangeability, or keep the gas composition in an interchangeable range. Recently, some works start to pay attention to maintaining gas interchangeability. For example, the security of the gas mixtures is considered in the optimal operation of H-IEGS with alternative gas in Ref. [9]. It has also been integrated into the robust operation of H-IEGS with hydrogen and renewable energies in Ref. [10].
Nevertheless, the previous studies usually adopt static gas security criteria, such as the Wobbe index, combustion potential, etc., straightforwardly to calculate gas interchangeability. Though these static criteria can measure the values of gas interchangeability quantitatively under a specific given condition at a certain time point, they cannot fully reflect the dynamic abilities of the H-IEGS to maintain gas interchangeability during time-varying operating conditions. In another word, these studies cannot fully utilize the flexibility of the gas system itself to defend the gas interchangeability violations. For example, when the hydrogen injection is increased, the hydrogen fraction can increase suddenly near the injection point, and the gas interchangeability near that point may be inferior. However, the gas system may increase the gas flow rate near that injection point gradually, pushing the gas interchangeability back to the acceptable range. This ability of the gas system to defend the gas interchangeability variations, and to recover gas interchangeability to the normal level, can be defined as the resilience of gas interchangeability.
Though the resilience of the electricity system [11][12][13] or gas system [14][15][16] regarding the energy supply capability is widely studied, the resilience of the H-IEGS with respect to the gas interchangeability against the fluctuating renewable generations has not been investigated yet [17] . The evaluation of gas interchangeability resilience requires near real-time captures of the dynamics of gas flow, and more importantly, gas composition variations due to the fluctuation of hydrogen injections. These two kinds of dynamics are governed by PDEs, which are very challenging to incorporate in the optimal operation of H-IEGS efficiently. Especially for the gas composition dynamics, though recently a few studies start to incorporate them into the coordinate operation of H-IEGS [18] , the solution efficiency still needs improvement. Moreover, its impact on gas interchangeability has not been investigated either.
To address the research gaps, this paper proposes a novel resilience evaluation method for gas interchangeability in the H-IEGS. Other contributions are summarized as follows: (1) Novel metrics are defined from multiple dimensions to evaluate the resilience of gas interchangeability. Compared with traditional static gas security metrics, the proposed metrics can: (a) better reflect the capability of the gas system to defend gas interchangeability violations caused by the hydrogen injections; (b) calculate the accumulated gas interchangeability loss, which is more in line with the real gas safety regulations [19] .
(2) A two-stage gas interchangeability management scheme is proposed. It entails solving the SOEF problem first to determine the desired operating state of the H-IEGS. In the second stage, the DGCT problem is solved to determine the optimal path of reaching the desired operating state, as well as using the gas flow dynamics to mitigate the gas interchangeability loss during this process.
(3) A self-adaptive linearization technique is used to solve the discretized PDEs. The reference points of the gas flow, and gas property coefficient (i.e., specific gravity, compressibility factors, etc.) are updated based on a gap criterion, so that the computation efficiency can be improved while the accuracy can still be guaranteed.

Gas interchangeability
Gas interchangeability can be defined as the ability to substitute one gaseous fuel for another in a combustion application without materially changing operational safety, efficiency, performance or materially increasing air pollutant emissions [20] . Gas interchangeability is tightly dependent on the gas composition. With the injection of hydrogen, the gas composition may change, which may further jeopardize gas interchangeability. The traditional method, such as the Dutton method, which is widely used in the UK, Australia, etc., usually adopts three metrics, i.e., WI, ICF, and SI to measure the gas interchangeability from different aspects [19] . However, these regulations are established without considering the participation of hydrogen [21] . Since hydrogen has a higher flame speed, here we introduce another metric FS to further guarantee gas interchangeability [22] . They are introduced as follows: (1) WI: WI is designed to guarantee that gas appliances can produce the same amount of heat energy when consuming the same amount of gas under prespecified conditions. Because hydrogen has a lower heat value, WI can reflect the decrease of heat energy contained in the gas mixture brought by hydrogen. Thus, it is used as a criterion for gas interchangeability, as calculated in Eq. (1). If the WI is too high or too low, the heat energy from the combustion may exceed or be beneath the design value of the gas appliance, and are thus unacceptable. Therefore, we have upper and lower bounds for WI, as shown in Eq. (2).
The GCV and specific gravity in Eq. (1) are calculated by: (2) ICF: ICF characterizes the ability of a gas to be combusted completely without residues. It can be calculated by the empirical equations in Eq. (5). A higher ICF indicates that the gas mixture contains more propane, butane (with a higher carbon element), or inert gases, while containing less methane, ethane (with a lower carbon element), or hydrogen. In this case, the gas mixture is less likely to combust completely. Therefore, an upper bound for ICF is set, as shown in Eq. (6).
(3) SI: SI measures the soot formation in the gas appliance during the combustion, which is calculated by Eq. (7) [23] . Soot is a mass of impure carbon particles resulting from the incomplete combustion of hydrocarbons [24] . It can cause various types of cancer and lung disease. A higher SI means there will be more soot during the combustion. Therefore, an upper bound is set, as shown in Eq. (8).
(4) FS: Recently, more and more studies have reported that the integration of hydrogen can increase the flame speed of the gas mixture, which may cause serious issues like overheating, flashback, corrosion, etc. [25] FS can be used to measure this phenomenon. It describes the approximate maximum velocity with which a flame can travel in any gas-air mixture, as calculated in Eq. (9) [26] . Thus, keeping the FS within a certain range is effective in regulating the combustion dynamics, as shown in Eq. (10).

Resilience of gas interchangeability in hydrogen-blended IEGS
144 iEnergy | VOL 2 | June 2023 | 143-154 Among these gas interchangeability matrics, WI and FS are most likely to be affected by hydrogen injections, and further violate their constraints. For example, the violation of the lower bound of WI is usually caused by injecting excessive hydrogen by PTGs. A typical countermeasure is to convert more hydrogen into methane first, and then inject the methane into the gas network instead of injecting hydrogen directly. It is also feasible to just reduce the amount of hydrogen injection. For the FS lower bound violations, it does not happen constantly. It is usually caused by unqualified gas sources. For example, the gas sources from biogas have more impurities (e.g., carbon dioxide, nitrogen, etc.) than natural gas, which may cause a lower FS. To overcome this issue, we usually strictly control the source of the gas supply, preventing unqualified gas from injecting to the gas system. Besides, injecting hydrogen can also improve this situation.

Resilience of gas interchangeability
Gas interchangeability introduced in the last subsection is a static concept. For example, at a given time point, the gas composition is fixed, and then the interchangeability metrics are fixed. It focuses more on whether the gas composition violates the gas interchangeability constraints. However, the hydrogen injection from renewable generation is time-continuously stochastic. Because the gas system is also a dynamic and continuous system, it is important to measure how the gas system adapts to the hydrogen injections during continuous operation. With this idea, the concept of gas interchangeability resilience is proposed in this subsection.
The resilience of gas interchangeability is a dynamic concept, which is similar to the resilience in the power system. In power systems, resilience can be defined as the ability to defend, endure, and recover from the impacts of a triggering event [27] . In the presence of uncertain and fluctuating hydrogen injections, the resilience of gas interchangeability can be defined as the ability of the H-IEGS to maintain the gas interchangeability within the acceptable range. Or, in cases where the violation of gas interchangeability constraints is inevitable, it also refers to the ability to quickly recover to acceptable interchangeability. Compared with gas interchangeability, the resilience of gas interchangeability focuses more on the gas system itself. It focuses on evaluating the ability of the gas system to defend against or recover from potential gas interchangeability violations during the operation. Apart from the violation of gas interchangeability at a certain time point, it also focuses on the maximum violation of the gas interchangeability during the operational period (e.g., the worst scenario), the duration of the disturbance, recovering time, accumulated performance loss, etc.
Following this idea, here we define several metrics that measure the resilience of gas interchangeability. An illustrative example is shown in Figure 1. Suppose there is a spike in wind generation, as shown in Figure 1(a). Then, the hydrogen production of PTG also increases correspondingly. Figure 1(b) shows the travel of the hydrogen content in the gas network. At different times (e.g., , , and ), the distribution of the molar fraction in the pipeline is different. As the peak of the hydrogen molar fraction curve gradually moves from the injection point to distant locations, the gas interchangeability metrics at a specific location, e.g., , is shown in Figure 1(d).
Several metrics are proposed to characterize the resilience of gas interchangeability from multiple dimensions (we take WI as an example. The resilience of other gas interchangeability metrics can be derived similarly): (1) Maximum gas interchangeability loss . This metric can give information on the worst gas interchangeability due to the stochastic hydrogen injection during the operation.
(2) Settling time . This metric shows the "inertial" of the gas system in defending the hydrogen injections. If the settling time is long, it means the gas interchangeability in this gas network is difficult to be affected by hydrogen injection. It will take a longer time to make the gas interchangeability violate the constraints. If the settling time is short, it means that the gas interchangeability in this gas system is prone to follow the fluctuation of wind power, which is more likely to violate the constraints quickly. Settling time is calculated by the time length between and , as shown in Figure 1(d). is the time point when the triggering event happens. The triggering event is the emergence of the wind speed which will lead to the gas interchangeability violation based on SOEF. is the time point when the gas interchangeability violation happens. Then, the settling time can be calculated by: Length of the enduring phase . The length of the enduring phase reflects the ability of the gas system to recover from a gas interchangeability violation event. If is long, it means that the gas system will take more time to recover from a gas interchangeability violation event. It is usually associated with the settling time (or the "inertial" of the gas system). A longer settling time usually suggests a longer time length of the enduring phase. is calculated by: (4) Accumulated gas interchangeability violation . This metric indicates the accumulated gas interchangeability loss during the hydrogen injection. It is calculated by: 2 Two-stage gas interchangeability management scheme During the operation, due to the fluctuation of the hydrogen injection, gas interchangeability may be affected. The H-IEGS will take active measures to defend against gas interchangeability loss. For example, if the hydrogen injection is high, the system operator  may choose to increase the gas flow near the hydrogen injection point to dilute the hydrogen fraction, although it is not the most economic way to do so in normal states.
Considering that renewable generation is difficult to predict, here we propose a two-stage gas interchangeability management scheme, as shown in Figure 2. In the first stage, the SOEF is solved first to determine the desired operating state of H-IEGS with minimum gas interchangeability loss. However, due to the slower gas flow and gas composition dynamics, the system may not able to reach the desired state immediately. Therefore, in the second stage, the DGCT problem is performed to optimize the path toward the desired state (i.e., the real-time operating condition of the H-IEGS with less gas interchangeability loss and more hydrogen injection). During the operation, the solution of the first-stage SOEF optimization problem is passed to the second-stage DGCT to set the boundary conditions, as shown in Figure 2. The solution of the second stage does not pass information back to the first stage. Owing to this one-directional information exchange, the error will not be accumulated.

First stage: steady-state optimal electricity and hydrogengas flow problem k
In this stage, the SOEF is performed for each renewable generation capacity level to determine the desired operating condition of the H-IEGS, as illustrated in Figure 2. The objective is to minimize the total cost, including electricity generation and gas purchasing costs, as shown in Eq. (15) (the notation of is omitted for conciseness). Theoretically in our framework, the time resolution in SOEF depends on the state transition of wind, and thus does not have to be the same as the time step in the PDEs of DGCT. The SOEF only has to be performed when the generating capacity of the wind farm changes. Practically in our study, because we have wind data with high time resolution, and the computation time of SOEF is relatively small, we conduct the SOEF at the same time resolution as in DGCT [28] .
The optimization model is subject to: (1) Gas system constraints: where Eq. (17) indicates that the total heat energy of the gas demands of natural gas and hydrogen should equal the total heat energy of the original gas demand measured in natural gas; Eq. (18) indicates that the gas composition of the gas consumption of gas demand should equal the gas composition at the exact gas bus; Eq. (19) is the Weymouth equation for gas mixtures [29] ; Eqs. (20) and (21) are the nodal gas flow balances of natural gas and hydro-

Second stage
First stage Steady-state optimal electricity and hydrogen-gas ow Gas (25) and (26) are the calculation methods for the gas constants of the gas mixtures at bus and in the pipeline , respectively; Eq. (27) is the calculation method for the gas composition in pipelines; Eqs. (28)−(33) are the bounds for optimization variables.
(2) Electricity system constraints: The first stage defines the desired operating state of H-IEGS. Then, in the second stage, the DGCT is implemented to optimize the path in which the H-IEGS reaches the desired operating state. Therefore, the optimization problem is formulated on each time step, with the objective of minimizing the deviations to the desired operating state, as well as the gas interchangeability loss, as shown in Eq. (44). It is worth noting that because the gas system is a dynamic system governed by PDEs, there may exist an overshoot in the response to the hydrogen injection change. To better balance the feasibility and optimality, the relaxation factor and corresponding penalty coefficient are introduced. By this means, slight violations of the gas interchangeability constraints are allowed at the beginning. With the increase of time step , by increasing the penalty factor , the slack variable will gradually approach zero, and the gas interchangeability constraints can finally be guaranteed.
It is subject to: (1) Gas dynamic constraints: The dynamics of gas flow, as well as the traveling of specific gas content, are governed by three PDEs, namely, continuity equation, motion equation, and advective equation. Their discrete form in an isothermal and horizontal pipeline can be written as [18] : It is worth noting that the discretization forms of the variables should be chosen carefully in order to keep a good balance between accuracy, feasibility, and stability. Normally, the central difference scheme has better accuracy, but is more likely to cause oscillation issues. Therefore, the implicit differentiation scheme is applied as well.
(2) Initial and boundary conditions: The above PDEs are formulated for each pipeline in the gas system. To derive the operating state of a dynamic system, initial and boundary conditions are required. The initial condition is determined by the SOEF (for k = 1) or the operating state at the last time step (for k >1), which is obtained by solving the DGCT problem at the last time step, as shown in Figure 2: The initial conditions of gas density, gas composition, and gas flow, can be given similarly.
The boundary condition is given by the pipelines it is connected with. For example, the gas pressures at the connecting point should be equal. The gas composition at the beginning of the pipeline should be equal to the upstream gas bus. Therefore, we have: The boundary condition for the gas flow is given by the nodal gas flow balance equation, similar to Eqs. (20) and (21).
Due to the variation of hydrogen injection, the gas constant becomes variable subject to the gas composition. Thus, the gas pressure no longer has a linear relation with the gas density. Thus, the motion equation and advective transport equation become more nonlinear than the case with constant gas composition, which greatly affects the gas system dynamics.

Solution method
The optimization problems in the first and second stages are both nonlinear and nonconvex. To solve the problem more tractably, an adaptive linearization method is proposed. The idea of the adaptive linearization method is to approximate the nonlinear terms using Taylor expansion or other kinds of linearization methods, and the reference point of the Taylor expansion will be updated if the state of the system changes dramatically. With this idea, the optimization problem can become more tractable without increasing the computation burden, and the accuracy can also be improved compared with fixed reference points.

Adaptive linearizations method
In the first stage SOEF problem, the nonlinearities exist in the: (1) the nodal consumption of gas demand Eq. (18); (2) the Weymouth equation (19); (3) the nodal gas mixing equation (24). They are reformulated to: In the second stage problem, the nonlinearities exist in the: (1) motion equation (47); (2) advective transport equation (48); (3) gas state equation (55). They are reformulated as [30] : The reference point is initially selected according to the solution in the first stage SOEF problem (e.g., ) when . When , the following criterion is set to determine whether the reference point should be updated: If Eq. (61) is satisfied, then update the reference point of as .

Solution procedure
The whole solution procedure is elaborated as follows: K ∆x ∆t Step 1: Initialize the system data, including the wind speed, system structure, energy demands, physical data of components, etc [31,32] . Set the total time steps , and the time and space resolutions for PDE problems and .
Step 2: Solve the SOEF problem by replacing the gas production of PTG with natural gas according to Ref. [9]. Set the solutions as the reference points (e.g., ) for .
Step 3: For each wind level, solve the first stage SOEF problem based on Eqs.

Case studies
Mm 3 /day ε In this section, the IEEE 24 bus reliability test system and Belgium gas system are used to validate the proposed method. The topological structure of the two coupling systems is presented in Figure  3. The detailed data of the two systems can be found in Refs. [33,34], respectively. Three PTGs are installed at electricity bus #10, 17, 18, respectively, which are also connected with gas bus # 1, 5, 9, respectively. The capacities of the PTGs are 2 . The wind farm is located at electricity bus #18, with a generation capacity of 800 MW. The time and space resolutions are set to 1800 s and 10000 m, respectively. is set to 0.01.

Validation of proposed method
First, we validate the proposed method in terms of both computation efficiency and credibility. We assume only the PTG at electricity bus #18 is available. The wind generation is assumed to increase from 100 MW to 800 MW at . As a result, the hydrogen production of the PTG increases from 1.55 to 2 .
We compare our method with a general nonlinear solver (e.g., IPOPT). The total computation time of our method with is 1.53 s, which is 98.94 % faster than the IPOPT solver (143.69 s). The solutions of the nodal molar fraction of hydrogen right after the hydrogen injection, at which time the system state changes dramatically and is most likely to cause large errors, are compared in Figure 4. As we can see, the relative error between the proposed method and the IPOPT is tiny (the calculation method of the relative error can be found in the Appendix). The errors at all buses do not exceed 1%. The largest errors occur at gas bus #9 (the injection point) and #14 (where the pure natural gas and hydrogen-natural gas mixtures are mixed). At other gas buses, the errors are almost neglectable.
Moreover, we conduct a sensitivity analysis on the penalty factor to validate the effectiveness of our proposed DGCT framework. Four scenarios are set. In S1, ; in S2, ; in S3, ; in S4, . The optimization results of the molar fraction of hydrogen at the injection point gas bus #8, and the value of , are shown in Figure 5. Mm 3 /day μ As we can see, different values of lead to different results of the molar fraction of hydrogen and . Comparing S1, S2, and S3, we find that the increase in will lead to the decrease of . As a result, the molar fraction of hydrogen will be lower with a lower value of . The relative differences of the molar fractions of hydrogen when can be up to 7.85%. However, we should note that selecting a large penalty factor at the beginning is not always the best option. Although it will strictly guarantee gas interchangeability, it will increase the computation time, and lead to a local optimum and loss of optimality. For example, although the in S3 and S4 both end up in zero after , the hydrogen productions are 1.83 and 1.85 , respectively. The hydrogen production of PTG in S4 is 1.09% higher than in S3. Therefore, our selection method of can better balance the feasibility and optimality.

Gas interchangeability resilience of different gas buses
In this subsection, we investigate the gas interchangeability k = 2 resilience of different buses. We assume only the PTG at electricity bus #18 is available. The wind generation is assumed to increase from 100 MW to 800 MW at .
First, some critical gas interchangeability metrics (e.g., the Wobbe index and flame speed factor) when the hydrogen is injected into gas bus #9, are presented as an example in Figure 6. We can find that, with the injection of hydrogen, due to its lower GCV and higher flame speed, the Wobbe index and flame speed factor gradually decrease and increase with time, respectively. However, the spikes of gas interchangeability could be generated in some gas buses (gas bus #9 and 14 in this case) near the hydrogen injection time point, as shown in Figures 6(a) and 6(b). This is because the suddenly injected hydrogen can be stacked near the injection points. These spikes are not caused by numeral simulations, and cannot be neglected because they will cause unexpected temporal violations of gas interchangeability metrics. For example, the maximum and accumulated Wobbe index losses at gas bus #9 and #14 are 0.  To defend against the gas interchangeability loss, the gas system raises the gas flow rate of natural gas near the injection point to dilute the molar fraction of hydrogen. For example, as shown in Figure 6(e), the nodal natural gas injections at gas buses #10, 11, 12, and 13 increase dramatically near the injection time point. Figure 6(f) shows the settling times of gas buses in response to the hydrogen injection. It can be seen that the settling time is longer with the longer distance between the gas bus and the injection point. For example, gas buses #16 and 20 are at the ends of two pipeline routes, and their settling times can reach up to 10 and 11 hours. It shows that the transient process of gas composition cannot be neglected in the gas interchangeability evaluation. The longer settling time means those gas buses can have more time to prepare and respond to the hydrogen injection, which increases the resilience of these buses against the hydrogen injection.
Following this idea, we can study the resilience of the gas system with the hydrogen injected into other different gas buses, i.e., gas buses #1, 5, and 8. The gas interchangeability metrics are shown in Table. 1. We can observe that different gas buses have different resilience against hydrogen injection. Comparing S3 and S4, we can find that the upper stream locations are more suitable for hydrogen injections. Compared with S4, the total settling time is longer in S3, which means the gas system will have more time to respond to mitigate the gas interchangeability loss. The maximum Wobbe index and flame speed factor losses in S3 are 23.81% and 9.75% lower than that in S4, which means the gas interchangeability losses are also lower. The accumulated Wobbe index and flame speed losses in S3 are also 36.68% and 36.38% lower than in S4, which means the gas interchangeability loss will recover faster in S3.
Gas buses #1, 5, and 8 are all at the beginning of pipeline routes.
Observing these scenarios, we can find that the gas interchangeability shows various patterns when the hydrogen is blended in different gas buses. Comparing S2 with S1 and S3, we find that although the gas interchangeability loss is lower, the gas production of PTG is also significantly lower. This is because the gas flow rate near gas bus #5 is relatively low, which cannot support the large volume hydrogen injection without causing gas interchangeability loss. Comparing S1 and S3, we find that the settling time of S1 is shorter than S3. This is because when hydrogen is injected in gas bus #1, it only propagates to limited gas buses (i.e., #2-4, 6-7, and 14-16). These gas buses are relatively closer to gas bus #1, and thus the gas system will have less time window to regulate gas interchangeability. In contrast, other gas interchangeability resilience metrics in S1, such as maximum and accumulated Wobbe index and flame speed factor loss, are lower than in S3. This means the gas system will have a stronger ability to adapt if the hydrogen is injected into the gas bus #1.

Gas interchangeability resilience during daily operation
In this case, the proposed gas interchangeability resilience evaluation method is further applied to daily operations. The wind speed data are obtained from National Oceanic and Atmospheric Administration (https://www.noaa.gov/). All the PTGs are committed in this case. The operating conditions of the H-IEGS are presented in Figure 7. From Figure 7(a), we find that the hydrogen production of PTGs follows the wind speed generally. The highest wind speed appears around 3:00−4:00 and 9:00−14:00, during which the hydrogen production also reaches its maximum values. PTG #2 is less affected by the wind, which contributes to hydrogen production significantly over the operational horizon. Since PTG #2 is connected with gas bus #8, we further present the flame speed factors along the critical pipeline route which starts from gas bus #8 to #16, as shown in Figure 7(b). The distribution pattern of the flame speed factor along that route verifies the theory that we developed

Resilience of gas interchangeability in hydrogen-blended IEGS
in the last subsection. The gas buses that are closer to the hydrogen injection point are more likely to violate the upper bound of the gas interchangeability metrics (e.g., gas bus #9). For the downstream gas buses, such as #12-16, they will have more time windows to increase the gas flow rate of the natural gas to dilute the molar fraction of hydrogen, and thus the gas interchangeability is more resilient. Observing from the whole H-IEGS, not only the gas buses near #8, but also the gas buses near #5, are likely to violate the gas interchangeability constraints, as shown in Figure 7(c). This is because #5 also has hydrogen injections. Unlike gas bus #8 where the gas flow rate is large, the gas source in the gas bus #5 is relatively small. Therefore, the gas system will have less flexibility to dilute the molar fraction of hydrogen in gas bus #5, and the gas interchangeability resilience is also inferior.

Validation of proposed metrics
To validate the effectiveness of the proposed metrics, a more general case is set.  Fig. 8 shows the operating conditions of the H-IEGS with multiple wind farms. From Figure 8(a), we find that the gas production of the PTGs also follows the wind speed generally. However, compared to the last case, this case shows the spatial differences in the gas production of PTGs. PTG #1 almost does not produce any hydrogen. PTG #3 only produces hydrogen during the wind peak hours, while PTG #2 produces a large amount of hydrogen throughout the operation. This is owing to the locations of the PTGs. Although these PTGs are all at the electricity buses that are different from those of wind farms, the PTG #2 is very close to the wind farm at electricity bus #18. In contrast, the location of PTG #1 is very distant from any of these three wind farms. This phenomenon proves the importance of location for PTGs.
From Figures 8(b) and 8(c), we can see the resilience of gas interchangeability of gas buses. In Figure 8(b), almost all the gas buses violate the FS upper bound. However, it is obvious that the resilience of gas interchangeability increase, as the distance to the hydrogen injection point increase. For example, gas buses #12 and 13 are more resilient than others. With the same hydrogen injection pattern, gas bus #8 violates the FS constraints immediately at 9:00  Resilience of gas interchangeability in hydrogen-blended IEGS ARTICLE iEnergy | VOL 2 | June 2023 | 143-154 and 11:00 just after the hydrogen injection increase. In contrast, the gas bus #12 successfully defended the hydrogen injection at 9: 00, and its FS does not violate the constraint. Even facing the hydrogen injection starting from 11:00, the FS of gas bus #12 is still controlled within the normal range before 20:00, and the violations are very slight. For the gas buses #14-16, their gas interchangeabilities are controlled within the normal range throughout the operation. Moreover, to demonstrate the advantages of the proposed metrics, we compare the operating condition of H-IEGS with/without the proposed metrics. Two scenarios are established. In S1, the gas interchangeability resilience is not considered. In S2, these constraints are considered.
The operating conditions of the H-IEGS are compared in Figure  9. We find that with the consideration of gas interchangeability resilience, the hydrogen productions of PTGs are reduced, especially in peak hours. The gas productions of PTG #2 and #3 are reduced by 2.75% and 38.0%, respectively. Consequently, the FSs during the whole operation are then controlled within the upper bound, as shown in Figures 9(c) and 9(d). Especially at gas bus #12, due to a large amount of hydrogen injection, the FS constraint is violated by up to 0.52% in S1. While in S2, the FS is reduced significantly.

Conclusions
This paper proposes a gas interchangeability resilience evaluation method in integrated electricity and gas systems with the injection of hydrogen. Novel gas interchangeability resilience metrics are proposed. A two-stage gas interchangeability resilience management scheme is proposed to improve the interchangeability considering the traveling of hydrogen content across the gas network. A self-adaptive linearization technique is proposed to improve the computation efficiency while maintaining a satisfying accuracy. The numerical results show that our solution method can improve the computation efficiency by 98.94%, while the relative error is controlled with 1%. We also find that the upper stream gas bus with a relatively large gas flow rate is more suitable for hydrogen injections in terms of improving the gas interchangeability resilience. For example, the gas interchangeability losses in gas bus #8 under hydrogen injection are lower than in gas bus #1. We also demonstrate the effectiveness of our method in a multi-period daily operation of H-IEGS. With the urging requirement for a lowcarbon energy system, the gas interchangeability resilience evaluation method proposed in this paper can help the energy system operators to optimize the operating strategy of the H-IEGS in the future.