Development of Three-Dimensional Soil-Amplification Analysis Method for Screening for Seismic Damage to Buried Water-Distribution Pipeline Networks

A soil-amplification analysis method is developed that uses high-resolution ground data and a three-dimensional nonlinear dynamic finite-element method to screen for possible areas of seismic damage to buried water-distribution pipeline networks. The method is applied to a cut-and-fill developed area in Japan, whose water-distribution pipeline network was severely damaged in the 2011 off the Pacific Coast of Tohoku Earthquake. The obtained soil amplification is compared with known points of pipeline damage to check the validity of the analysis. A sensitivity test is also conducted to account for uncertainties in the properties of the ground material. From the results, it is expected that the developed soil-amplification method could be used to screen for possible damage to buried pipelines in a given area, and used to support methods for estimating damage to buried pipelines based on observations and seismic indices.


Introduction
When a large earthquake hits a city, it is important that the water distribution systems either continue being usable or are at least unavailable for as short a time as possible.Because water distribution systems work under hydraulic pressure, their seismic performance is usually assessed by excluding any damaged pipelines and conducting hydraulic simulation over the remaining network [1][2][3].Conducting many such simulations under anticipated earthquake scenarios helps identify the critical parts of the system and decide on effective countermeasures.Because the hydraulic analysis is conducted based on the distribution of pipeline damage, improving the assessment of seismic damage to buried pipelines is essential for estimating the seismic performance of water distribution systems reliably.
Dynamic ground strain caused by local soil amplification is a major cause of seismic damage to pipelines [4], so the estimation of this amplification is important for assessing the damage to buried pipelines.However, measuring local soil amplification requires many sensors and is therefore expensive.Hence, many methods have been proposed to estimate earthquake hazard empirically or statistically (e.g., [5][6][7][8]), and use them as inputs for repair rate curves for estimating pipeline damage.These methods are often used to estimate the seismic damage done to pipeline networks for water distribution [1][2][3].However, such methods have limited estimation accuracy when applied to ground with complex geometry, and it is not straightforward to consider three-dimensional (3D) ground effects [9].With the improvement in sensor technology, computational methods, and computer technology, ground data information accuracy and resolution has improved, and physics-based simulation based on such data is becoming feasible.In addition to physics-based earthquake wave propagation simulations targeting wave propagation from source to sites (e.g., [10][11][12]), methods have been developed recently to estimate local soil amplification by analyzing dynamic nonlinear 3D ground motion [13,14].In these methods, a ground model consisting of ground layers is generated from boring logs, and soil amplification is estimated by inputting acceleration waveforms into high-resolution 3D analyses.Because the dynamic nonlinear wave equation is solved, it is expected that the time-history characteristics of the input waves, as well as the material properties and geometrical features of the ground, will be reflected in the estimations of local soil amplification.Although this method cannot obtain accurate results where detailed ground information is not available, and thus statistical or empirical methods may be sufficient in such cases, this method has possibility to improve the reliability of pipeline damage estimations where detailed ground information is available.
In previous studies, 3D soil amplification analyses were applied to sites with single high-pressure gas pipelines [13,14].Although such main gas pipelines are usually placed in sites that have relatively good seismic characteristics, some parts of an entire water-distribution pipeline network could be located where there are distinct changes in the ground properties.Cut-and-fill land is such an example, where larger soil amplification may occur locally.Indeed, many pipelines in these areas in Japan were damaged during the 2011 off the Pacific Coast of Tohoku Earthquake, hereinafter referred to as the 2011 Tohoku Earthquake (see [15] for damage to buried pipelines and [16,17] for damage to other facilities).It may be possible to improve the accuracy to which seismic damage to pipeline networks is estimated by analyzing the soil amplification of the entire network area while reflecting the 3D ground information.In the present study, we developed a modeling and analysis method that is capable of estimating local soil amplification for pipeline networks over large areas with complex geometry that changes on the scales of the order of 10 m.We apply this method to a district whose water-distribution pipeline network was severely damaged in the 2011 Tohoku Earthquake.Herein, a ground model with detailed 3D geometry is generated from digital ground data, and numerical analysis is conducted using 3D solid finite elements.Because pipeline damage is correlated with the strain of the surrounding ground, the validity of the 3D soil-amplification analysis is assessed by comparing the computed strain distribution with known locations of pipeline damage.Although the estimation cannot be perfect because of uncertainties in the ground model and input waves, these simulations are expected to be useful for screening for possible damage in a given area.
The remainder of this paper is organized as follows.In Section 2, the methodology of soil-amplification analysis is explained.Here, database-generation of ground models, generation of finite-element models, and large-scale soil-amplification analysis method is explained.In Section 3, the developed method is applied to a district that incurred pipeline damage during the 2011 Tohoku Earthquake, and the computed response distribution is compared to known damage locations.In Section 4, the paper is summarized.

Overview
Because local soil amplification is strongly correlated with local geometry and material properties, the analysis should explicitly reflect ground geometry and material heterogeneity.An unstructured 3D solid finite-element method is suitable for modeling complex geometry and material heterogeneity.However, small elements are needed to attain verifiable results with numerically convergent solutions, making finite-element model generation and high dynamic nonlinear analysis costly.It is also expensive to generate high-resolution 3D models of ground geometry when targeting large areas.Thus, we developed an automated method that generates ground geometry models from digital ground data, and applied a robust mesh generator and fast 3D finite-element analysis method.Because of the uncertainties associated with pipe aging, ground models, and input waves, we targeted their soil amplification analysis at the screening of areas with possible damage to a buried pipeline network.This involved comparing the locally amplified ground strain at the position of the buried pipeline with the design strain levels.

Method of Generating High-Resolution Ground Models
Point-wise boring data or two-dimensional cross sections are often used to estimate 3D ground structures.However, a target area of a few square kilometers has a ground structure that change every few tens of meters, and it would thus be expensive to use such methods to estimate ground structure.Thus, we developed a method for estimating ground structure automatically by means of elevation maps before and after cut-and-fill land development.Figure 1 shows the process of modeling the ground geometry.The input data for the target area before land development are in the form of a 2 m contour map, whereas the data for the developed land are in the form of a digital elevation map (DEM) on 5 m structured grid [18].The original 2 m contour map was firstly traced to generate a 5-m-grid DEM.The DEM difference before and after land development was then used to obtain the cut-and-fill depth.Because a 2-m contour map was used, the maximum error involved in converting the data from a contour map to a DEM was 2 m.Thus, this method could be used to generate ground geometry models for large cut-and-fill land areas to an accuracy of ±2 m.Because the target area is large (a few square kilometers) and elements on the order of 1 m are required to attain numerical convergence of the strain response, the finite-element model has the order of 10 8 degrees of freedom (DOFs).For such a large problem, it is not trivial to generate finite-element models that reflect the complex geometry while having high-quality elements; using general-purposed mesh generators could lead to elements with squashed shapes that would cause numerical errors in the obtained solution, as well as poor convergence of numerical solvers.Often, it is necessary to change the mesh-generation parameters on a trial-and-error basis to circumvent this problem, so mesh generation is a bottleneck for conducting finite-element simulations at such a scale.Thus, we used an octree-based finite-element modeling method [19] that decomposes the domain with an octree and generates a mesh locally for each octree-leaf.Because the meshing in each leaf is independent of the others, a complex geometry can be decomposed into elements robustly.The independence between meshing of each leaf leads to a straightforward parallelization on shared-memory computers, leading to the fast generation of a large mesh.Furthermore, the octree data structure was used to enlarge the elements for homogeneous areas, which led to fewer elements and thus a lower computational cost for soil-amplification analysis.Using this method, we could model regions of a few square kilometers robustly with second-ordered tetrahedral elements of the order of 1 m.

Large-Scale Soil-Amplification Analysis Method
Following previous studies [13,14], we used a solid 3D finite-element method that can model complex ground geometry appropriately while satisfying the traction free surface boundary condition analytically.Modeling ground shaking as a wave propagating in a nonlinear material and discretizing in the temporal domain using the Newmark-β method (β = 1/4, δ = 1/2) leads to the following: Here, δu, u, v, a, and F are the vectors of incremental displacement, displacement, velocity, acceleration, and external force, respectively.M, C, and K are the mass, damping, and stiffness matrices, respectively, dt is the time step, and n is the time step number.Rayleigh damping is used; the element damping matrices are computed as using an element mass matrix M e and an element stiffness matrix K n e .Here, α and β are determined by solving the least squares problem minimize Here, f max , f min , and h are the maximum target frequency, minimum target frequency, and damping ratio, respectively.The modified Ramberg-Osgood model [20] and Masing rule [21] are used for the nonlinear constitutive model.Half-infinite absorbing boundary conditions are applied to the side and bottom of the domain, and the wave is input from the bottom of domain when analyzing the soil amplification.
As the computing cost increased considerably, we used a fast nonlinear dynamic 3D finite-element method parallelized by MPI and OpenMP.Here, the solver developed in [22,23] was used, which is an accelerated solver based on the method used in [13].This finite-element solver was verified through numerical convergence of acceleration waves and SI values [24] with errors less than 1% [22,25].
It was validated by confirming that computed waveforms match observed waveforms recorded during the 2011 Tohoku Earthquake [13] at a site with complex ground geometry.By using the robust finite-element mesh-generation method with this solver, we could analyze local soil amplification for regions with complex geometry and high material heterogeneity.

Problem Settings
Soil amplification analysis was conducted in an area whose water-distribution pipeline network was severely damaged during the 2011 Tohoku Earthquake, and the computed response distribution was compared with actual pipeline damage.The target region was developed in the late 1960s and the early 1970s by cutting and filling a hilly area in Sendai City, Miyagi Prefecture, Japan.Figure 2 shows the interface between the cut-and-fill layer as generated from ground data [18].The water-distribution pipeline is buried 1.2 m below the black lines indicated on the surface; the pipes are buried in a region whose cut-and-fill area has complex, branch-like shapes.The pipeline was laid at the time of land development, was 40-200 mm in diameter, and consisted mainly of VP (vinyl chloride pipes with adhesive joints) or DP (ductile iron pipes with type-A fittings) types.Following [26], we set the material properties of the ground as in Table 1.Because the cut-and-fill ground has high contrast in its material properties, with the complex interface changing on the order of 10 m, the soil amplification could be locally concentrated.
Figure 3 shows the generated 3D finite-element model.Here, 150 m buffer regions were added to the four sides of the 1320 × 1020 m (EW × NS) target area to avoid reflection waves that cannot be completely damped by the absorbing boundary conditions.As the non-uniformity of the ground is only near the surface, we modeled the region from a depth of −127 m to avoid reflection waves that cannot be completely damped by the absorbing boundary conditions at the bottom of the model.Modeling the domain with 2.5 m second-order tetrahedral finite elements leads to a finite-element model with 85,143,240 DOFs, 20,666,257 elements, and 28,381,080 nodes.The wave observed in the target area during the 2011 Tohoku earthquake (observation point S14, SAKR of Small Titan [27]) was pulled back to the bottom of the model (i.e., the wave at the bottom of the model is estimated based on 1D equivalent linear analysis) and input into the finite-element model to conduct nonlinear dynamic analysis.The wave was input as velocity and stress at the bottom and side boundaries of the model, and the waves that do not match with the 1D analysis results were damped by the absorbing boundary conditions.Here, we used a target frequency of f max = 2.5 Hz, f min = 0.1 Hz to set the Rayleigh damping parameters.Figure 4 shows the input wave.Here, the velocity response spectrum S v of the level-2 design wave used for the seismic design of water-distribution pipelines is shown together.It can be seen that the merged horizontal component (i.e., S 2 vx + S 2 vy ) of the input wave has similar magnitude to that of the design spectrum.Computation of 300 s of analysis (dt = 0.001 s × 300,000 time steps) took 46,993 s (13 h 03 min) using 512 compute nodes of the K computer [28].

Comparison of Computed Strain and Design Strain
Pipeline damage can be estimated from the ground strain of the surrounding soil, which is caused by soil amplification.Because the pipeline is buried near the surface, the strain response is seen initially at the surface.Figure 5 shows the maximum displacement response and maximum principal-strain response at the surface.The displacement is relatively small in the cut areas but larger in the areas with deep fill layers.This matches the intuitive expectations based on the one-dimensional ground structure directly beneath each point.The maximum principal strain becomes large at places slightly within the fill area from the cut-and-fill interface.Thus, places with V-shaped cross sections comprising a soft layer surrounded by cut ground have two lines with large strain responses (e.g., Figure 5(b-1)).A similar strain response was observed at V-shaped cross sections comprising a soft layer between hard layers during the 1978 earthquake off the coast of Miyagi Prefecture [29].Additionally, previous studies using 3D soil-amplification analyses reported similar strain responses at V-shaped cross sections [13], and there is thus consistency among observations and analyses for such V-shaped ground.By using 3D analysis, it is also possible to obtain the strain distribution at areas with more complex shapes, such as Y-shaped cut/fill interfaces (e.g., Figure 5(b-2)).The strain distributions at such areas are not easily estimated using simpler analysis methods, so an advantage of the present 3D analysis is its ability to estimate quantitatively for these complex-shaped areas.Consequently, the maximum strain response in the whole target domain was obtained in this Y-shaped area.In addition, the strain response is small for some filled areas (e.g., Figure 5(b-3)) .Thus, the strain response is not necessarily large for all areas of filled ground, but depends on the 3D shape of the surrounding ground.
Next for consideration is the relationship between the computed strain distribution and strain estimated using seismic indices.An example of seismic estimation based on seismic indices is the estimation method used in the seismic design of water-distribution pipelines in Japan [30].As a seismic index, this method uses the velocity response spectrum S v (T G ) for strain estimation; the ground is estimated as a one-dimensional stratified structure, and the displacement of an S wave vibrating at its natural period T G is obtained with Here, t is time, z is depth from the surface, S v (T G ) is the velocity response at normal period T G = 4H/V s2 , H is the depth of the surface soil layer, and V s1 , V s2 are the shear-wave velocities for the cut and fill layers, respectively.The maximum principal strain is approximated by U(t, z) propagating horizontally at a speed of V = 2V V s2 /(V s1 + V s2 ).Thus, the strain at the surface is estimated as Equation ( 6) is multiplied by coefficient µ to reflect the complexity of the ground at the target site, and the pipelines are designed to withstand this strain level (µ is set to 1.2 for complexly shaped ground sites and 1.0 elsewhere).Because the input wave in this study corresponds to a level-2 wave, the difference between the design strain and the actual strain that can happen for a large earthquake can be estimated by comparing the computed strain-response distribution.Figure 6 shows the strain estimated by the design method (µ = 1.2 is used for the whole domain).When compared with the computed strain in Figure 5b, similarities such as the strain becoming larger in the fill regions can be seen.Also apparent is the fact that the maximum strain-response amplitude is similar between the two methods.However, the distribution inside the fill region differs locally in some places.For example, the strain tends to be large at the center of V-shaped fill areas for the design strain, which differs from that observed in the 1978 earthquake and that computed by the 3D analysis.The reason for such a difference can be found by seeing the time history displacement and strain distribution of the 3D analysis as shown in Figure 7.As seen in the top row of Figure 7, the displacement is amplified at the fill regions, and the displacement response thus differs between the cut regions and the fill regions.This difference in displacement response causes a large strain in regions slightly within the fill region from the cut/fill interface (bottom row of Figure 7).The large strain occurs in regions slightly within the fill region from the cut/fill interface for all of the four time steps shown.On the other hand, as the velocity response spectrum S v of the design wave increases monotonically with the natural period (Figure 4), the estimated strain also increases monotonically with the fill-layer depth.Thus, it can be seen that the design method is sufficient for estimating the maximum magnitude of the strain response in a given area, but it may not be as accurate when predicting the position of the maximum strain response.Thus, depending on the setting of µ, the centers of the V-shaped fill domains may have excessive design strength or the edges of such domain may lack design strength.Indeed, much of the pipeline damage indicated with circles in Figure 6 occurred near the edges of the fill regions, which corresponds better with the 3D analysis results.

Screening for Pipeline Damage
Although peak ground velocity or peak ground acceleration are typically used for indices in repair rate curves of pipelines [31][32][33][34], pipeline damage is pointed out to be due to the deformation and strain of pipes [4].By using the developed 3D ground analysis method, we can directly compute ground strain; thus, we use strain as an index for damage assessment of pipelines.Because 12 out of the 17 instances of damage in the target area were due to dislocation at pipe joints, we consider the axial strain of the pipes, which is correlated with pipe-joint dislocation.Figure 8a shows the axial strain at the underground pipeline position.Here, axial strain is computed as the ground strain component in the axial direction of the pipe at the pipeline position.Compared with Figure 5b, the distribution of axial strain at the pipeline network is similar to the maximum principal strain at the surface.Thus, it can be seen that the strain distribution at the surface can be used to estimate the subsurface ground strain at the pipeline network.Because water-distribution pipelines are designed to withstand a strain of 0.004 in this area [30], the computed maximum axial strain (0.00357) is large enough for damage to occur, given the effects of long-term aging.When comparing the pipeline strain distribution with that of the damage in Figure 8a, the locations of maximum strain do not necessarily match those of the damage.However, there appears to be no damage in the cut areas.In addition, the damage is not uniform inside the fill areas; parts of the fill areas with small computed strain were not damaged (e.g., area 1 in Figure 8a).Although we can tune the ground geometry or material properties such that the estimated strain distribution becomes closer to the actual pipeline damage distribution, we can see that the results even without tuning of ground geometry and material properties could be used to screen for areas with possible damage (i.e., identify areas that have low probability of damage).
Compared to the relatively accurate ground geometry, the ground material properties may not be as accurate because they are not based on data measured at the target site.Thus, discrepancies in the location of any damage may be due to the material parameters.Here, consideration is given to the sensitivity of the ground strain response to changes in the ground material properties.Figure 8b,c show the strain response for two additional sets of material parameters (Cases 1 and 2 summarized in Table 2).Compared to the standard settings used in the previous computations, the nonlinear material parameter of the fill soil is changed to 0.0005 in Case 1.In Case 2, the nonlinear material parameter and S-wave velocity of the fill soil is changed.We chose these two parameters because they are expected to make large differences in ground motion.The values of these parameters were set based on typical measurement errors during soil material property tests.The maximum axial strain changed from 0.00357 in the standard case to 0.00345 for Case 1 and 0.00488 for Case 2. However, the positions with large strain are similar regardless of the material properties.Thus, it can be seen that the strain distribution is affected more by the ground geometry than by the material properties.In reality, backfill is used when burying pipelines, and the damage uncertainties may be due to the inhomogeneity of the pipeline joints and to aging.Further improvement in the estimation of pipeline damage can be expected by analyzing pipeline deformation by means of soil springs that account for slip conditions between backfill soil and the pipe.As the three-dimensional ground motion distribution is computed, curvature strain, which was not considered for damage assessment in this paper, can also be accounted for by using pipeline deformation analysis with soil springs.

Closing Remarks
A seismic soil-amplification analysis method has been developed for estimating seismic damage to a buried water-distribution pipeline network using high-resolution ground data and large-scale 3D finite-element analysis.Here, a method was developed to generate ground models from digital data, a robust meshing method was applied to generate finite-element models of large areas with complex geometry, and a fast dynamic nonlinear finite-element method was used to conduct seismic soil-amplification analysis of region of a few square kilometers under the 2011 Tohoku Earthquake.
The computed strain obtained from soil-amplification analysis was in accordance with that observed at V-shaped fill areas in the 1978 earthquake off the coast of Miyagi Prefecture, while a complex strain distribution was obtained in complexly shaped ground domains such as Y-shaped domains.Although the axial strain distribution, which is expected to be highly correlated with pipe-joint damage, did not exactly match the observed damage points, damage was not observed in areas with small computed axial strain.Thus, the developed method could be used to screen for areas in which damage was more likely.The obtained maximum ground strain was similar to the maximum design strain over the whole area, but the local strain distribution was different among the methods inside the fill region.The seismic design of pipelines could be improved by using such insights regarding the strain distribution.For example, the method could be used to choose the order in which old pipelines are replaced in a given area.Based on the results of this study, the probability of pipe failure can be linked with hydraulic simulation results of water distribution networks.In the future, we plan to compute the response to several input waves to estimate the sensitivity of the strain distribution to the input-wave properties.We also plan to conduct pipeline deformation analysis using soil springs that account for slip conditions between backfill soil and the pipe.With improvements in computer technology, we expect to use such high-resolution 3D analysis practically to support damage-estimation methods based on measurements and ground-motion indices as currently used in practice.

FillFigure 1 .
Figure 1.Generation of ground geometry model.Old terrain data before cut/fill of ground given as a contour map that was traced to make a digital elevation map (DEM).This was subtracted from a DEM of the current terrain to make a ground geometry model.Red in the generated ground geometry model indicate fill regions, while blue in the ground geometry model indicate cut regions.

Figure 2 .Figure 3 .
Figure 2. Target area.(a) Dark gray areas indicate cut ground (hard), while light gray areas indicate fill ground (soft).Black lines indicate position of pipeline.Here, pipeline is buried 1.2 m below the surface.Circles indicate the position of pipeline damage during the 2011 Tohoku Earthquake.(b) Elevation of cut/fill-layer interface.White lines indicate the cut/fill-layer interface.

Figure 4 .
Figure 4. Input ground motion.Here the wave observed in the target area during the 2011 Tohoku earthquake (observation point S14, SAKR of Small Titan [27]) is pulled back to the bottom of the model and used as the input wave.Velocity response spectrum of design wave of pipeline is also indicated.(a) Acceleration time history.(b) Velocity response spectrum (damping h = 0.15).

Figure 5 .
Figure 5. Distribution of the maximum displacement and maximum principal strain.White lines indicate boundary between cut-and-fill layers.(a) Maximum displacement.(b) Maximum principal strain.

Figure 6 .Figure 7 .
Figure 6.Design strain for level-2 ground motion (factor µ = 1.2 for all regions).White circles indicate the positions of pipeline damage.

Table 1 .
Material properties (standard case).Here, V p , V s , ρ, h max , and γ r , indicate primary wave velocity, secondary wave velocity, density, maximum damping, and reference strain, respectively.

Table 2 .
Material properties for Cases 1 and 2. Material properties of cut layer and bedrock are those of standard case.Here, V p , V s , ρ, h max , and γ r indicate primary wave velocity, secondary wave velocity, density, maximum damping, and reference strain, respectively.
V p m/s V s m/s ρ kg/m