## 2 Methodology
### Uncertainty Analysis  
The deterministic analysis procedure adopted in the previous section assumes precise knowledge of the system parameters. However, all problem variables are characterized by a certain degree of uncertainty, including the material strength and the seismic loading. To quantify the effect of these uncertainties on the pipeline performance, this study develops a robust fragility analysis framework based on the Monte Carlo Simulation (MCS). The latter considers a large number of soil-pipeline system samples, generated by random sampling of the input system parameters, based on their probability distribution. The seismic demands are compared to the structural capacity values associated with each performance limit state, allowing us to estimate the associated probability of exceedance conditioned to the level of ground motion intensity measure, IM, represented by the vector (\(\delta,L_{b}\)). Table 5 summarizes the input variables and their probability distribution. To investigate the effect of the soil cover depth \(H_{c}\), three different values between 1.0 m and 2.0 m were considered, representative of typical onshore pipeline installations. The soil strength parameters, i.e., the longitudinal soil reaction per unit length of pipeline \(f_{s}\), and relative soil-pipe displacement at friction sliding, \(u_{0}\) are assumed to follow a normal distribution, with the mean values derived according to [31], and a coefficient of variation COV = 30% [26; 35; 36], as indicated in Table 5. The uncertainty of the X42 steel grade pipe material was modeled considering that the yield strength follows a normal distribution, with a mean value \(\sigma_{y}\) = 290 MPa and a small COV = 3.5% (Table 5), representing the low variability of the steel properties. Figure 8: The peak axial strain magnitude in the unpressurized pipeline (\(P_{i}/P_{max}\) = 0, \(\Delta T\) = 0 \({}^{\circ}\)C) as a function of the PGD length \(L_{b}\) and displacement \(\delta\) for (**a**) tension and (**b**) compression. The dashed horizontal curves represent the strain isolines corresponding to the NOL and PIL performance limit states. The limit state functions for normal operability (NOL) and pressure integrity (PIL) for tension and compression are expressed in terms of the system demand (\(\delta\), \(L_{b}\)) and capacity (\(\delta_{cri}\), \(L_{cri}\)) corresponding to each performance criteria:  
\[g_{t1}=\max(\delta_{t1}-\delta,L_{b,t1}-L_{b}) \tag{6}\]  
\[g_{c1}=\max(\delta_{c1}-\delta,L_{b,c1}-L_{b}) \tag{7}\]  
\[g_{t2}=\max(\delta_{t2}-\delta,L_{b,t2}-L_{b}) \tag{8}\]  
\[g_{c2}=\max(\delta_{c2}-\delta,L_{b,c2}-L_{b}) \tag{9}\]  
The probability of exceeding the normal operability (NOL) and the pressure integrity limit (PIL) state, conditioned to the PGD intensity measure level (\(\delta\), \(L_{b}\)), is given by the joint union of the associated damage states for tension and compression:  
\[\text{P}[\text{NOL}\,|\,(\delta,L_{b})]=\text{P}[\max(g_{t1},g_{c1})\leq 0 \,|\,(\delta,L_{b})] \tag{10}\]  
\[\text{P}[\text{PIL}\,|\,(\delta,L_{b})]=\text{P}[\max(g_{t2},g_{c2})\leq 0 \,|\,(\delta,L_{b})] \tag{11}\]  
These probabilities can be effectively calculated using the Monte Carlo Simulation (MCS) method, as summarized in Table 6.  
The described algorithm can be easily implemented within most programming languages, like Python [37], to evaluate the fragility surfaces representing the conditional probability of the system reaching a performance limit state, as a function of the seismic demand.