Design Analyses For A Large-Span Tunnel in Weak Rock
Design Analyses For A Large-Span Tunnel in Weak Rock
Design Analyses For A Large-Span Tunnel in Weak Rock
Bhaskar B. Thapa Jacobs Associates, San Francisco, California Johannes Van Greunen Jacobs Associates, San Francisco, California Yiming Sun Jacobs Associates, San Francisco, California Michael T. McRae Jacobs Associates, San Francisco, California Hubert Law Earth Mechanics Incorporated, Los Angeles, California
ABSTRACT: The proposed Caldecott fourth bore will consist of a two lane highway tunnel along California State Route 24 near the City of Oakland. The proposed design and construction sequence for the 15-mdiameter tunnel are based on the New Austrian Tunneling Method (NATM). The initial support system incorporates combinations of shotcrete, rock dowels, lattice girders, spiles, and grouted steel pipe canopies. The nal lining is cast-in-place reinforced concrete. A waterproong membrane and drainage system are placed between the initial and nal linings. State Route 24 is a lifeline route, required to be open to emergency vehicles within 72 hours after a major earthquake, dened as having a return period of 1,500 years and a peak ground acceleration of 1.2 g. Although the seismic design criteria are stringent, the design of the tunnel lining system is ultimately controlled by static ground loads in the weak rock along the alignment. INTRODUCTION Project Background The existing Caldecott Tunnel complex includes three bores along State Route 24 (SR 24) through the Berkeley Hills in Oakland, California. The California Department of Transportation (Caltrans) and the Contra Costa Transportation Authority (CCTA) propose to address congestion on SR 24 near the existing Caldecott Tunnels by constructing a fourth tunnel that will provide two additional trafc lanes. The proposed horseshoe-shaped fourth bore is 1,036 m (3,399 ft) long, 15 m (50 ft) in diameter, and 9.7 m (32 ft) high. The project will include short sections of cut-and cover tunnel at each portal, seven cross-passageway tunnels between the fourth bore and the existing third bore, electrical substation buildings, and a new operations and control building. State Route 24, considered a lifeline route by Caltrans, is required to be open to emergency vehicles 72 hours after an earthquake with a return period of 1,500 years and a peak ground acceleration of 1.2 g. Construction of the fourth bore is anticipated to begin in the summer of 2009 and be completed in 2014. GEOLOGY Major Geologic Formations and Structure The geology of the alignment is characterized by northwest-striking, steeply-dipping, and locally overturned marine and non-marine sedimentary rocks of the Middle to Late Miocene age. The western end of the alignment traverses marine shale and sandstone of the Sobrante Formation. The Sobrante Formation includes the First Shale, Portal Sandstone, and Shaly Sandstone geologic units as identied by Page (1950). The middle section of the alignment traverses chert, shale, and sandstone of the Claremont Formation. The Claremont Formation includes the Preliminary Chert, Second Sandstone, and Claremont Chert and Shale geologic units (Page, 1950). The eastern end of the alignment traverses non-marine claystone, siltstone, sandstone, and conglomerate of the Orinda Formation. Major formations and geologic units within these formations are shown Figure 1. The geological structure of the project area has been characterized as part of the western, locally overturned limb of a broad northwest-trending syncline, the axis of which lies east of the project area. The fourth bore alignment will encounter four major
417
Figure 1.
Geologic formations and geologic units Seismicity The San Francisco Bay Region is considered one of the more seismically active regions of the world, based on its record of historical earthquakes and its position astride the tectonic boundary between the North American and Pacic plates. During the past 160 years, faults within this plate boundary zone have produced numerous small-magnitude (M<6) earthquakes, and more than a dozen moderate- to large-magnitude (M>6) earthquakes. Major faults that comprise the 80-km-wide plate boundary in the San Francisco Bay Region include the San Gregorio, San Andreas, Hayward, and Calaveras Faults. The active Hayward Fault, located 1.4 km (0.9 mi) west of the Caldecott Tunnel, is the closest major fault to the project site, capable of producing a magnitude 7.4 earthquake. The southern segment of the Hayward Fault produced the 1868 Hayward earthquake of estimated magnitude 6.8 that was accompanied by 30 to 35 km (19 to 22 mi) of surface faulting. INITIAL SUPPORT DESIGN The initial support system design is based on the Sequential Excavation Method (SEM), also known as the New Austrian Tunneling Method (NATM). NATM provides the required exibility to accommodate the variable ground conditions and weak, folded, and faulted rock that will be encountered along the Caldecott fourth bore alignment. The design approach involves classication of ground along the alignment into several ground classes, development of corresponding support categories,
inactive faults, which occur at the contacts between geologic units. These faults strike northwesterly and perpendicular to the tunnel alignment. In addition to the major faults, many other zones of weak ground will be encountered, such as smaller-scale faults, shears, and crushed zones. West of the fault contact between the Preliminary Chert and Shale and the Second Sandstone, the bedding encountered in the fourth bore generally dips predominantly northeast. East of this fault contact, the bedding dips southwest. Several joint sets occur within each geologic unit, and random joints occur in almost all orientations in all geologic units. Intrusive sandstone dikes and hydrothermallyaltered diabase dikes occur most frequently in the Claremont Chert and Shale, but may be encountered less frequently in other geologic units. The structure of the rock mass units along the alignment varies from blocky in the best ground to disintegrated or crushed in the poorest-quality rock. Average RQD ranges from 5 to 81. Rock Mass Ratings (Bieniawski, 1989) and Q values (Barton, 1988) at the tunnel scale vary from 20 to 65, and 0.006 and 10.5, respectively. Rock strength varies from weak to moderate along the alignment. Average values of measured unconned compressive strength in the various geologic units vary from 5.2 MPa (750 psi) to 21.6 MPa (3190 psi). Mudstone, siltstone, and shale in the Orinda and Claremont Formations are expected to exhibit swelling behavior. The fourth bore has been classied as a gassy tunnel by the California Occupational Safety and Health Administration.
418
Figure 2.
Longitudinal displacement proles for SC I is modeled using Mohr-Coulomb elastic-plastic continuum elements in FLAC3D. The hardening of the shotcrete lining is modeled as the tunnel top heading and two bench cuts advance at prescribed rates and lags to represent the early age creep effects of shotcrete described in Thapa et al. (2007). Ground relaxation factors are estimated based on a tunnel longitudinal displacement prole (LDP) and a ground reaction curve (GRC). The tunnel LDP (see Figure 2) demonstrates the development of tunnel radial displacement as a function of distance along the length of the excavation, and can be generated from FLAC3D analysis results. The GRC (see Figure 3) shows the tunnel radial displacements as a function of support pressure, and can be generated from a two-dimensional FLAC analysis. To estimate the ground relaxation factor, a FLAC3D analysis of the entire excavation sequence was performed. From this analysis, three LDPs were generated, one corresponding to each stage of excavation. From each LDP, the drift radial displacement (ur0) prior to installation of initial support was estimated. Then, a FLAC2D analysis was performed to generate the GRC for the excavation stage under consideration. Next, the radial displacement (ur0) estimated from the LDP was used to locate the corresponding support pressure on the GRC. The ground relaxation factor (GRF) for the drift under consideration was estimated as follows: GRF = ( 1 R ) 100%
and denition of criteria for application of the support categories during construction. Four major and two minor ground classes, and corresponding support categories, have been developed for construction of the fourth bore. Support category I applies to the best quality rock mass and Support Category IV applies to the poorest quality rock mass. For the fourth bore design, the team also developed additional support measures to be used if unexpected geologic conditions are encountered during construction or monitoring reveals unexpected, unfavorable ground behavior. A description of the general design approach is provided in Thapa et al. (2007) and is not repeated here. However, the sections below describe some of the specic design analyses including both two- and three-dimensional convergence-connement analyses that were performed with FLAC (Itasca, 2005) to evaluate specic design issues for the NATM initial support design. These design issues are: Stress relaxation ahead of the tunnel heading Face stability Lining loading across weak zones
Stress Relaxation Ahead of the Face FLAC3D models of the full NATM excavation and support operation were developed for each support category to estimate the amount of relaxation in the ground ahead of tunnel face. The FLAC3D models explicitly represent the sloping core used for face support and spiling presupport. The shotcrete lining
419
Figure 3.
Figure 4.
Schematic illustration of estimation of ground relaxation factor using FLAC3D results horizontal displacements near the springline of the bench drifts. The above approach used in estimating the GRF is consistent in principle with the current practice in tunnel design (Carranza-Torres and Fairhurst, 2000 and Graziani et al., 2005). Tunnel displacements (or strains dened as the radial displacements divided by the tunnel radius)
The GRFs for other stages of tunnel excavation were estimated in the same way. Figure 4 illustrates schematically the use of the LDP and GRC to estimate the ground relaxation factor ahead of the tunnel face. The radial displacements utilized in generating the LDPs and the GRCs are the vertical displacements near the crown of the top-heading drift and the
420
Figure 5.
Strength envelopes for various factors of safety in Support Category I The procedure for calculation of the face stability factor of safety begins with initialization of the model to top heading equilibrium conditions following excavation and support installation using the actual rock mass strength envelope. Factors of safety for the face region are then calculated by iteratively reducing the rock mass strength (corresponding to increasing factors of safety), agging failure zones corresponding to the iteration factor of safety and contouring zones with the same factor of safety (see Figure 6). This iteration is repeated until the model fails to reach mechanical equilibrium or a predetermined number of increments in the F value is reached. The range of F values evaluated varies from 1.0 to 5.0. During iteration, failure of a zone representing the rock mass is dened as the non-convergence of the zone velocity to a value of less than 106 m/s. Face stability is evaluated in a region that extends to four tunnel diameters longitudinally and to the model limits vertically and transversely. Predicted factors of safety against general shear failure for the top heading drift ranged from 3.2 in Support Category I without any face support to 1.3 in Support Category III with a sloping core for face support. Lining Loading Across Weak Zones A 180-m-long reach adjacent to the fault contact between the Second Sandstone and Claremont Chert and Shale geologic units is expected to have subreaches of varying ground quality ranging from very poor to fair. This reach occurs under the highest cover along the tunnel alignment. Shotcrete lining thickness requirements for this reach were evaluated using FLAC3D to account for the effect of longitudinal
calculated from FLAC3D and FLAC2D are generally in good agreement. Iterative calibration of the GRF ensures that the displacements or tunnel strain of the supported tunnel from FLAC3D at the section where a plane-strain condition is reached match those from FLAC2D. GRF estimates from the above analyses ranged from 58% to 65% for various support categories. Face Stability FLAC3D was also used to evaluate face stability by determining the factor of safety (F) against global shear failure for the top heading drift. In the FLAC3D face stability analysis, F is the factor by which the rock mass shear strength must be divided to bring the drift face to the verge of failure. The resulting factor of safety (F) is the ratio of the actual rock mass shear strength to the reduced shear strength at failure, which can be expressed as: o F = ---r where, o is the actual shear strength and r is the reduced shear strength at incipient failure. Figure 5 shows the actual and reduced rock mass strength envelopes for Support Category I, corresponding to various factors of safety. The strength envelopes shown in Figure 5 are based on four Hoek-Brown failure criterion strength parameters. The four Hoek-Brown criterion parameters for each factor of safety were determined using a cubic spline interpolation scheme built into Microsoft Excel.
421
Figure 6. Factor of safety around face of top heading in Support Category I (Figure shows longitudinal section through tunnel centerline. Top heading half-width =7.5 m, height=5.5 m.)
Figure 7. Stresses in shotcrete lining through ground with varying material properties (Tunnel width=15.0 m, height=10.5 m.) arching on lining loads. Figure 7 shows the lining loads developed across this reach from a FLAC3D analysis. Comparison of FLAC3D results to a FLAC2D analysis that does not account for longitudinal arching shows that lining loads computed using FLAC3D are about 30% lower than the loads calculated using a plane-strain FLAC2D analysis. It is noted that the FLAC3D results are in general agreement with the FLAC2D results in other reaches with uniform ground conditions. FINAL LINING AND SEISMIC DESIGN Final Lining System The Caldecott fourth bore uses a double lining system consisting of an initial support system (discussed above) and a cast-in-place reinforced concrete nal lining (Figure 8). A waterproong membrane with a geotextile backing layer for drainage will be installed between the initial support and the nal lining. The initial support system is
422
Figure 8.
Final lining and are considered temporary. In most of the tunnel, the rst 50 mm (2 in.) of shotcrete lining are applied as a ash-coat and considered sacricial. In the First Shale reach of the Sobrante Formation, the rst 100 mm (4 in.) of shotcrete is considered sacricial because of the high sulfate concentration in the groundwater in this reach. The remaining shotcrete layers are also expected to deteriorate to some degree over time. In addition, the initial shotcrete lining is assumed to have no exural capacity due to possible deterioration of any reinforcing embedded therein. Thus, as these components deteriorate over time, the nal lining will support a signicant portion of the ground load. Analyses were performed to assess the effect of the degradation of the initial support and to determine the part of the ground load that will be transferred to the nal lining. Analyses were performed using FLAC 5.0 (Itasca, 2004). Key assumptions of the analyses were: Rock dowels were completely deteriorated. Initial shotcrete lining thickness reduced by neglecting the sacricial layer as described above. The modulus of the reduced shotcrete lining was degraded to 60% of its original design value. The initial shotcrete lining has no exural capacity after degradation. Due to the presence of the waterproong membrane, the interaction between the initial and
designed to carry the ground loads that develop during construction, while the cast-in-place reinforced concrete nal lining is designed to carry long-term ground loads and any additional loads resulting from nishes or anchored equipment. The nal lining will also accommodate seismic deformations and provide a durable and sound tunnel lining. General descriptions of loads and load combinations, ground loads, and seismic demands have been described in Thapa et al. (2007) and are not repeated here. The following key aspects of the nal lining design are discussed below: Load sharing between the initial and nal linings Wave scattering analysis Pseudo-static time history analysis of seismic demands
Ground Loads from Load Sharing The initial shotcrete lining and the nal concrete lining will behave as a combined lining system. The long-term performance of the system will depend not only on the nal lining, but also on the long-term load-carrying capability and the durability of the initial shotcrete lining. During construction, the initial support will carry the ground load. However, two essential components of the initial support, the rock dowels and the shotcrete lining itself, are expected to deteriorate with time. The rock dowels proposed for the project are not protected against corrosion
423
nal linings was modeled with stiff radial springs and no tangential springs. The analysis was performed using two methods: Method A: The full tunnel excavation was modeled in one step with both the degraded initial lining and the nal lining in place. This is a conservative assumption that would be expected to transfer a somewhat higher portion of the ground load to the stiffer nal lining. Method B: The excavation and initial support installation sequence was modeled to develop the forces in the initial support. Thereafter the installation of the nal lining and deterioration of the initial support were modeled. Simply changing the properties of the shotcrete lining when the lining-ground system is in equilibrium is not a viable analytical approach using FLAC. Nodal forces are still in balance and no redistribution of forces would occur. Therefore, to force the model to perform a meaningful solution cycle, the forces (thrusts, shears, and moments) in the shotcrete lining were changed at the same time as the shotcrete lining properties. The forces in the shotcrete lining were reduced by the same factor used for the degradation of the material properties, thereby limiting the ground load supported by the initial lining and forcing redistribution of loads. The shotcrete lining forces change during cycling, and therefore have to be updated after each cycle to make sure that the resulting forces always satisfy the following criterion:
o ~ p = fp
~ represents the new parameters (thrust, where: p shear, and moment); po corresponds to the original parameters; f is the reduction factor, = (1 percent degradation ) = 0.6 here. The results of both analyses with Method A and B indicate that the nal lining will attract a maximum of 50% of the ground load supported by the initial lining. The nal lining was conservatively designed to support 2 3 of the ground load supported by the initial lining. Wave Scattering Analysis Wave scattering analyses were performed to calculate ground deformations around the tunnel opening in response to seismic wave propagation. This analysis accounts for the effect of local conditions such as tunnel geometry, adjacent tunnel cavities, geology, topography, and variation in rock quality. The time history of ground distortions around the tunnel
obtained from wave scattering analyses were used as input for the pseudo-static analysis of the tunnel nal lining seismic demand (as described later in this paper). The scattering analysis was performed using elastic material models with properties adjusted for small strain dynamic conditions. At large shear strains modulus reduction and increased damping were considered in the analysis. The tunnel lining was not included in the models as the linings are signicantly more exible than the ground and, therefore, only the properties of the ground determine the deformation of the tunnel opening. The scattering analyses were performed using QUAD4M, a nite element computer program (EMI, 2007). The nite element models include a transmitting boundary capable of minimizing seismic wave reection at the nite element boundary, which is used to model a semi-innite space outside the nite element domain (Hudson et al., 1994) (Lysmer and Kuhlemeyer, 1969). Three transverse cross sections of the mined tunnel, a cut-and-cover cross section at the west portal, and a longitudinal section of the tunnel were evaluated for wave scattering effects. The three transverse cross sections were selected to represent the critical combinations of cover and ground properties along the fourth bore alignment. Since the project is part of a lifeline route, ground motion criteria consistent with other important facilities on the same route including the Benicia-Martinez bridge, San FranciscoOakland Bay bridge (SFOBB) and Yerba Buena Island (YBI) tunnel, were selected for design. Thus, the ground motion adopted for the Safety Evaluation Earthquake (SEE) and Functional Evaluation Earthquake (FEE) are the 1,500 and 300 year return period uniform hazard spectra respectively. The performance requirements for the SEE are that the fourth bore will be open to emergency vehicle trafc within 72 hours following an SEE. Performance requirements for the FEE are that the fourth bore remains fully operational and experiences minimal, if any, damage. Three sets of earthquake time histories were developed to spectrum match the reference SEE rock spectra; Figure 9 shows three component time histories of the SEE reference rock motion for one of the three sets of earthquake time series used in the design. Results of the scattering analyses are illustrated in Figure 10 showing a comparison of the intensity of the computed motions at the tunnel perimeter and the reference motion in terms of acceleration response spectra. The spreading of the response spectra clearly indicates non-coherent wave propagation, which results in differential motion around the tunnel cavity. The computed acceleration time
424
Figure 9.
Reference rock motion (SEE Set 1) for the Caldecott fourth bore elements to model ground behavior. The two-dimensional beam-spring models were used for design to calculate strains, stresses, and forces in the fourth bore lining and cut-and-cover structures, and to ensure that the results were within acceptable stress and ductility limits. The two-dimensional beam-continuum models were used to verify the results obtained from the beam-spring models. All of the numerical models were initialized with gravity loads (rock loads and rock wedge loads) before the simulation of the seismic events. Two-dimensional beam-spring SAP2000 (CSI, 2006) models of the nal lining were developed for all support categories. The lining was represented by linear beam elements while the ground was modeled with equivalent springs, considered to be compression-only to simulate the passive support the ground will provide to the lining. The stiffness of the springs was based on the spring tributary area and the ground modulus of elasticity. The static mean modulus of elasticity was used for all analyses of static
histories at the nodes around the tunnel perimeter were integrated twice to yield displacement time histories which served as the multiple-support input motions to the tunnel lining Pseudo-Static Time History Analysis of Seismic Demand
Approach
State-of-the-art beam-spring and beam-continuum models were used to perform pseudo-static time history analyses of the tunnel nal lining using multiple-support displacement time histories from the scattering analyses described above. Two types of numerical models were used to calculate lining strains, stresses, and forces: two-dimensional SAP2000 (CSI, 2006) beam-spring models with nonlinear support springs (gap elements) to model ground behavior, and two-dimensional beam-continuum models using both FLAC (ITASCA, 2005) and ADINA (ADINA R&D Inc.) with elastic continuum
425
excursions outside the interaction envelope. However, the calculated reinforcing steel stresses and concrete strains are well within the allowable limits. CONCLUSIONS Design of initial support required several threedimensional evaluations. These evaluations were performed using FLAC3D and the results were combined with traditional two-dimensional and closedform-solution analyses. The FLAC3D evaluation of relaxation ahead of the face justied the use of high relaxation factors which resulted in lower support loading, and contributed to the selection of more realistic support requirements. FLAC3D evaluations of face stability showed that typical closedform solution evaluations can be unconservative and that three dimensional numerical analyses help assess more realistic face support requirements. The FLAC3D evaluation of lining loading across weak zones was unique and key to evaluation of support requirements in high cover reaches. The FLAC3D evaluation in weak zones showed the proposed shotcrete lining thickness was sufcient and the thicker lining required by a two-dimensional analysis was not necessary. State-of-the-art seismic design analyses were performed on this project due to the critical lifeline classication of the facility. The design considered high levels of shaking, several ground motion time histories for each design event, non-coherence due to wave scattering, and pseudo-static time history analysis of the lining response. The analyses showed that the low cover portal sections of the tunnel were subject to more severe seismic demands than interior sections with high cover. The design analyses demonstrated that a 381-mm (15-in) nal lining with 35 MPa (5,000 psi) concrete can support the ground loads and accommodate the seismic deformations. Seismic demands do not control the thickness of the nal lining, despite the close proximity of the project to a major active fault and seismic design criteria corresponding to an earthquake with a 1,500-year return period and a peak ground acceleration of 1.2g. ACKOWLEDGMENTS The authors would like to acknowledge Geomatrix Consultants for their work on the site geology, ILF Consultants for independent reviews of the initial support designs, and SC Solutions for their work on seismic demand analysis. The contents of this paper were reviewed by the State of California, Business, Transportation and Housing Agency, Department of Transportation and
Figure 10. Response spectra at fourth bore tunnel opening for Station 107+60 under SEE Set 1 motion loads. For seismic analyses the dynamic modulus of elasticity was determined by increasing the static moduli by a factor of between two and three. Gravity loads were applied through the support springs by displacing the xed ends of the springs and then iterating to achieve structural equilibrium with the required load in each of the support springs. The pseudo-static time history analyses were performed by imposing displacements, calculated at each time step through the scattering analyses described above, to the nal lining.
Results
The results of the analyses indicate that a 381-mm (15-in) nal lining with 35 MPa (5,000 psi) concrete can support the ground loads and accommodate the seismic deformations. The nal lining thickness was selected for constructability and is controlled by the thrust resulting from ground loads in the high cover section of Support Category III. In general, the analyses indicated that the maximum bending moment in the nal lining, as calculated with the different models, are not sensitive to the modulus of the ground. However, the lining thrust was generally signicantly higher for cases using the upper bound modulus of elasticity. Figure 11 summarizes the maximum lining thrust and moment for one of the critical support categories using the upper-bound ground modulus. In this support category, the thrust and bending moments due to seismic deformations result in some
426
Figure 11.
Support Category IV interaction diagramverication analyses for SEE Lysmer, J. M. and Kuhlemeyer, R. L (1969), Finite Dynamic Model for Innite Media, J. of the Engineering Mechanics Division, ASCE, Vol. 95, No. EM4, August, pp 859877. SAP2000 Version 10, Computers and Structures Inc, Berkeley, California, 2006 ADINA, ADINA R&D, Watertown, Massachusetts Earth Mechanics Incorporated, Caldecott Improvement Project, Technical Memorandum No. 16: Results of Scattering Analyses, 2007 Bieniawski, Z.T. (1989), Engineering Rock Mass Classications, Wiley, New York Barton, N., 1988, Rock Mass Classication and Tunnel Reinforcement Selection Using The Q-system, Rock Classication System for Engineering Purpose, ASTM Special Publication 984, American Society for Testing Materials, Page 5988. Thapa, B.B. et al. (2007) Preliminary Design of the Caldecott Fourth Bore, Proceedings of the Rapid Excavation and Tunneling Conference, Toronto. Carranza-Torres, C. and Fairhurst, C., 2000. Application of the Convergence-Connement Method of Tunnel Design to Rock Masses That Satisfy the Hoek-Brown Failure Criterion. Tunneling and Underground Space Technology, Vol. 15, No. 2, pp. 187213. Graziani, A., Boldini, D., and Ribacchi, R. 2005. Practical Estimate of Deformations and Stress Relief Factors for Deep Tunnels Supported by Shotcrete. Rock Mechanics and Rock Engineering, 38 (5), 345372.
the Contra Costa Transportation Authority. The contents of the paper reect the views of the authors who are responsible for the facts and accuracy of the data presented herein. The contents do not necessarily reect the ofcial views or policies of the State of California or the Contra Costa Transportation Authority. This paper does not constitute a standard, specication or regulation. REFERENCES Page, B.M. (1950), Geology of the Broadway Tunnel, Berkeley Hills, California, Economic Geology, Vo. 45, No.2 John, Max and Mattle, Bruno. Shotcrete Lining Design: Factors of Inuence. 2003 RETC Proceedings. Itasca Consulting Group Inc., 2005, Fast Lagrangian Analysis of Continua (FLAC) Version 5.0, Minneapolis. Hashash, YMA et al., Seismic Design and Analysis of Underground Structures, Tunneling and Underground Space Technology 16 (2001) pp. 247293, Elsevier. Hudson, M., Idriss, I. M. and Beikae, M. (1994), Users Manual for QUAD4M, A Computer program to Evaluate the Seismic Response of Soil Structures Using Finite Element Procedure and Incorporating a Compliant Base, Department of Civil & Environmental Engineering, University of California, Davis.
427