Nonlinear earthquake response of integral abutment bridges on vertical and batter piles: A practical approach

This paper presents a new finite element framework for seismic analysis of IABs (and other similar structures) on vertical and batter piles that takes into account nonlinear soil‐structure interaction. A previously developed macro‐element and linear impedance matrix are implemented in a time and modal domain solution, respectively. Detailed pseudo‐codes describing the implementation strategy for the macro‐element are provided. The solutions are validated using rigorous numerical models constructed in OpenSees MP. The software is further demonstrated by performing an extensive set of incremental dynamic analyses, where the effect of linear and nonlinear SSI, batter angle, pile spacing and asymmetric configuration is addressed and discussed.


INTRODUCTION
Bridges are often designed as statically determined systems which usually requires the use of bearings and expansion joints.During the life-span of a bridge, these mechanical elements require expensive maintenance routines.Integral abutment bridges (IABs) are monolithic structures without such elements.They are unique in the sense that substantially less maintenance is required compared to structures with bearings and expansion joints.In addition, the construction costs are generally lower, the durability is higher and the seismic risks are lower due to the absence of weak components.IABs are therefore popular among consultants, contractors and governmental road maintenance departments.However, monolithic structures are sensitive to restraint forces such as thermal effects, creep and shrinkage, and the monolithic connections increase the overall soil-structure interaction (SSI) during seismic excitation.In many cases, these structures are supported by deep foundations due to overlaying soft soil.Particularly poor seismic performance has been observed for deep foundation with batter piles.Until the early 1990s, batter piles were commonly used in seismic design of bridges and other large structures to improve the lateral capacity.In the following years however, batter piles became generally discouraged due to several earthquakes where batter piles experienced serious damage.In the Loma Prieta-earthquake, harbor ports supported by batter piles suffered severe damage.Several of the squared, pre-stressed concrete piles supporting the Public Container Wharf in the Port of Oakland failed in tension near the connection to the deck.Due to liquefaction, the batter piles settled and attracted large moments.Similar failure was observed at The Matson Terminal Wharf and the Oakland Outer Harbor Pier 7. The prestressed concrete batter piles supporting the Ferry Plaza Pier in San Francisco failed in tension, and some of the piles also punched through the deck.In April 1991, Costa Rica was struck by an earthquake with   = 7.2.The Rio Banano bridge was severely damaged due to liquefaction.At one of the abutments, it was observed that the battered piles in front suffered substantial bending and shear damage, whereas the vertical piles at the rear were less damaged.Today, several governing codes (e.g., Eurocode 8 1 ) advise against the use of batter piles in the seismic design of deep foundations.3][4][5][6][7][8][9] It has also been suggested that the unsatisfactory performance of batter piles is likely due to flawed design rather than inherent features of the pile. 10n recent decades, performance-based earthquake engineering (PBEE) has become a well-accepted practice in earthquake engineering.Traditionally, design codes and guidelines prescribe a fixed set of requirements that must be met in seismic design of structures.Such requirements ensure that the structure has adequate strength and ductility to safely resist the seismic forces.This is known as force-based design (FBD).In contrast, PBEE focuses on design methods that predict actual structural behavior during and after an earthquake.This approach enables solutions that consider how the structure responds to site-specific earthquake excitation based on factors such as structural configuration, materials, seismic data, and more.As PBEE involves a more realistic evaluation of structural response, SSI plays a crucial role in such design approaches.Although sophisticated software packages enable precise evaluations of most structural and geotechnical problems, detailed numerical solutions of dynamic SSI-response are generally more complex.The anticipated timeframe and the associated economic expenses related to the design of relatively simple structures such as IABs, do not facilitate for rigorous and time-consuming procedures.Consequently, practitioners must rely on simpler methods that are robust, user-friendly and easily verifiable.
A macro-element may be regarded as a nonlinear continuum model described by advanced constitutive laws.The main advantage of the macro-element approach is that the complex soil-foundation response is condensed into a single element.However, such elements require pre-defined parameters that must be calibrated, and they are often restricted to a specific foundation configuration, soil profile and soil type.Also, analyzing soil-structure problems using macro-elements inherently introduces an approximation since the analysis is fundamentally nonlinear, while kinematic and inertial effects are not assessed simultaneously.
The macro-element approach has increasingly garnered attention as a practical tool for evaluating SSI-problems for deep foundations.Correia 11 and Correia and Pecker 12 formulated a macro-element for monoshaft foundations where the formulation accounted for non-linear behavior of pile, soil and separation effects.In addition, the authors presented a rigorous calibration procedure for typical soil profiles.Liu et al. 13,14 first developed a macro-element for single piles embedded in homogeneous sand, which they later extended to single batter piles.This formulation was inspired by the work of Salciarini and Tamagnini. 15Pérez 16 presented a macro-element for vertical pile groups based on the work of Liu. 17 Page et al. 18 presented a macro-element model for mono-pile foundations based on multi-surface plasticity.The model was verified against large-scale model tests. 19More recently, Cemalovic et al. 20 developed a nonlinear macro-element for vertical and batter pile groups that considers the inelastic behavior of both pile and soil.The formulation is based on single-pile response with two separate load-displacement formulations (axial and transverse).Rotation is considered implicitly.The main advantage with aforementioned macro-element is that it does not require pre-defined failure surfaces or other parameters, and is therefore not restricted to a specific foundation configuration, soil profile or soil type.The macro-element may be calibrated using any type of non-linear pile-soil model, and the calibration procedure is rather straight-forward.Furthermore, Cemalovic et al. 21developed a simple, linear method for dynamic impedance of vertical and batter pile groups.The approach is based on a closed-form solution of a BWF-problem that takes into account pilesoil-pile interaction.The method is primarily intended for low-exaction seismic problems, vibration problems or estimates in the early-stage design process.
The main objective of this paper is to aid the industry with practical computational tools for IABs (and other similar structures) that are supported by deep foundations using both vertical and batter piles.These methods aim to capture the essence of SSI while also being straight-forward to understand, implement and apply.This paper supplements the research conducted by Cemalovic et al. 20,21 by demonstrating the developed solutions.More precisely, it involves the implementation of (1) the macro-element 20 within a time domain solution and (2) the linear impedance matrix 21 within a F I G U R E 1 Overview of the analyzed bridge.modal domain solution.Subsequently, these solutions are utilized to perform a set of incremental dynamic analyses (IDA) of an IAB where the effect of linear and nonlinear SSI, batter angle and pile spacing is evaluated.
The investigated system (Figure 1) is a monolithic, one-spanned IAB with relatively stub frame abutments.The abutments are founded on a 3-by-3 pile-group with both vertical and batter piles in soft clay.The selected pile cross-section is not based on a conventional standard-based design, but it represents a realistic case for the purposes of the research.
It should be noted that the abutment-backfill interaction, which is highly significant for the overall seismic response of an IAB, is not considered here.Addressing this interaction rigorously requires an advanced procedure, where the results are expected to depend highly on the chosen modeling strategy.One solution might be use of dedicated nonlinear macroelements for the soil-wall interaction which could be implemented in the model presented here.Alternatively, one could employ simplified calibrated springs for this purpose.The present study focuses on performance of macro-elements in a complete seismic SSI model; other refinements, as described, can be explored in future research.
This paper is divided in four sections.Section 2 describes the theory, implementation and validation of the numerical solution.Section 3 presents the IDA-analyses and Section 4 presents the conclusion.

Macro-element implementation
This section covers the theory and numerical methods used to solve nonlinear SSI-problems in the time domain using the macro-element.In essence, the macro-element serves as a function that provides the global finite element code the foundation stiffness at a certain time instance.The input of that function is a displacement increment, and the output is a coupled tangent stiffness matrix with three degrees of freedom.The macro-element is implemented into a newly developed (in-house) finite element package using the programming language Python 3. 22

Global time domain solution
The time domain solution is based on Newmark's method, 23 d−1 end while end procedure must be assumed.Newmark's method proposes the generalized equations where  and  are parameters that control the time step variation of acceleration.Combining Equations 1 and 3, the state variables at time instance   may be determined from known properties and state variables defined at the known time instance  −1 .Note that this approach assumes that the external load vector, together with the mass, stiffness and damping matrices, are not functions of the state variables.That assumption is only valid for linear systems.In most structural engineering problems, significant nonlinear effects stem from nonlinear material behavior (nonlinear constitutive models).The resisting force vector is then expressed as where the stiffness matrix is a function of the state variables.In that case, the nonlinear system must solved through an iterative procedure.The developed software uses the well-known Newton-Rhapson method [24][25][26] to restore the resisting force vector.The numerical scheme for the nonlinear time domain solution is shown in Algorithm 1.Since the frame model contains penalty constrains (imposed earthquake motion), it is recommended that relative incremental displacement is used to check convergence, that is, Here, || ⋅ || denotes the  inf -norm of the vector and   is the predefined tolerance.
It is often desirable to damp out high-frequency response, which may arise from spurious oscillations associated with the discretization of the problem.In order to maintain numerically stability and simultaneously add numerical damping to high-frequency response, it is recommended that  = 1 4 where  = 1∕2 equals no numerical damping. 27ocal macro-element solution The local macro-element solution is based on the bounding surface plasticity theory with radial mapping. 28,29A detailed description on bounding surface (or point) plasticity is provided in the literature 12,20 The essential features are repeated here for convenience.Displacement and rotation rates are decomposed into elastic and plastic components, that is, The rate of generalized forces is expressed as where  is the force and   is the elastic stiffness.The bounding load is restricted through where  is the load parameter.The evolution of plastic displacements is given by the plastic flow rule where the plastic multiplier  is zero during elastic response and greater than zero during plastic response.The hardening rule can be defined as where  is the hardening parameter.The consistency condition guarantees that the image point aligns with the bounding load at all times.It then follows that Solving for γ gives The macro-element algorithm assumes that the initial step in the creation of any loading (virgin loading, unloading or reloading) is purely elastic and that the subsequent steps within the respective loading always contain plastic deformations.
1][32] In short, this implies using Equations 7, 8, 10, 11, and 14 to obtain the elastic displacement, plastic displacement, load and load multiplier using values from the previous iteration or the last converged step.Convergence is achieved when the displacement residual and the bounding load equation both drop below a pre-defined tolerance value.Algorithm 2 presents the load predictor step.
Here,    is the single pile tangent stiffness in local coordinates,   is the trial load,   is trial load parameter,   is the minimum load parameter obtained during the previous loading steps,   (and   ) is the load (and load parameter) associated with the last initial unloading state,   (and   ) is the load (and load parameter) associated with the last initial reloading state and TANGENT_RETURN is the return mapping numerical scheme (which differs for transverse and axial response).For convenience, the subscript  that implies normalized values has been omitted from the presented algorithms.Nevertheless, normalized values are assumed in all cases.This part of the macro-element code determines if the elastic trial load corresponds to a load reversal point or not, that is, if a new loading is created or not.If the trial load is a load reversal point, then the trial load is correct.The code then proceeds to determine if the trial load is a virgin loading point, unloading point or reloading point, and the corresponding loads and load parameters are updated.If the trial load is a virgin loading point,   is updated.If the trial load is a unloading point,   and   are updated.Similarly, if the trial load is a reloading point,   and   are updated.If the trial load is not a load reversal point, the return mapping scheme is initiated.Although Algorithm 2 applies for both transverse and axial response, it should be mentioned that the predictor step must be assessed individually for transverse and axial response.
Algorithm 3 shows the return mapping scheme for transverse response.
Here, subscript  represents the iteration counter,   is the elastic displacement,   is the plastic displacement,  is the plastic multiplier,  is the displacement residual,  is the bounding load equation,   is the transverse plastic modulus,   0, is the transverse plastic modulus constant for fixed-head conditions and  , is a parameter that controls the transverse unloading/reloading curve shape.The return mapping algorithm is essentially an iteration scheme that aims to determine the nonlinear relationship between the elastic displacement, plastic displacement and the corresponding load for a given displacement increment.In the initial iteration step ( = 0), the state variables, load, load parameters are set to values obtained at the previous (last converged) time step .The plastic multiplier increment is set equal to zero.During the subsequent iteration steps , the trial load parameter, which was computed by the load predictor algorithm (Algorithm 2), determines if the load is a virgin loading point, unloading point or reloading point.The plastic modulus is computed, and the state variables, load parameter, displacement residual, bounding load equation, plastic multiplier increment and load are updated.If the displacement residual and the bounding load equation are within the prescribed tolerance, the transverse tangent stiffness is computed.Similarly, Algorithm 4 shows the return mapping for axial response.The only difference is that axial response requires a distinguished formulation for tension and compression.In both cases, V, U, R, T, and C denote virgin loading, unloading reloading, tension and compression.Once the transverse and axial tangent stiffness is computed, the foundation tangent stiffness matrix is assembled and passed to the global solution.Then, the global stiffness matrix   and internal force vector    are computed, and the global state variables are updated (Algorithm 1).A detailed description of the plastic modulus evolution for both transverse and axial response in the different loading stages, as well as the calibration procedures and its parameters, is provided by Cemalovic et al. 20 For convenience, the presented algorithms only show selected parts of the total numerical scheme.Assessment of the implicit transverse-rotational coupling, pile group assembly, and the interaction with the global code are also important aspects of the macro-element algorithm in terms of stability and efficiency.Furthermore, the pseudo-codes provided herein are simply illustrations of how the code can be implemented, and it is important to note that both the predictor step and return mapping schemes can generally be constructed in a manner that is most convenient for the task at hand.

Linear impedance matrix implementation
This section provides a procedure for establishing a diagonal impedance matrix for vertical and batter pile groups.The single pile impedance matrix may be expressed as where   is the coordinate transformation matrix and is the impedance matrix of pile  in local coordinates.The matrix entries of  *  , may be retrieved from the literature (e.g., Gazetas 33 ).We assume that there is no contact between the pile cap and the soil, that is, all forces and moments from the pile cap are transferred through the piles.In addition to the force distribution from the pile cap, pile-soil-pile interaction effects must be considered.Due to the rigidity of the pile cap, the displacement of the single pile  is approximately equal to the pile group displacement.With reference to Figure 2 ) ur, ) ( expressed as where  is the total number piles,    is the displacement of the source pile  due to its own load and is the interaction factor matrix between receiver pile  and source pile  in global coordinates.The interaction factor matrix represents the interaction between receiver pile  and source pile  in local coordinates.Closed-form expressions for the matrix entries of    are presented by Cemalovic et al. 21he horizontal displacement of pile  may be expressed as The vertical complex impedance is determined as Finally, the rotational impedance is assembled on the basis of individual vertical pile forces and moments.The vertical displacement of pile  is given as where   is the distance between the global node and pile .The rotational impedance is determined as and the diagonal impedance matrix is expressed as The modal approach is not strictly suitable for assessing linear SSI-problems since the stiffness expressions associated with foundation and soil are generally frequency-dependent.In that case, frequency-domain analysis provide the mathematically sound and accurate solution.Although there are especially developed software packages for SSI-problems in the frequency domain, 34 such solutions are not particularly convenient in practical engineering.In fact, solving structural engineering problems in the modal domain rather than the frequency domain provides several advantages.First the modal approach provides useful information in terms of modal shapes and the corresponding natural frequencies (dynamic properties of the structure), which are concealed in time-or frequency domain solutions.Second, modal analysis are the standard option in most commercial software packages, and the modal superposition concept is familiar for most practicing engineers.Third, modal analysis are an inherent part of the commonly used response spectrum analysis.6][37][38] Perhaps most practically, and especially in cases where estimates are sought, the frequency-dependent foundation stiffness may be approximated as the constant value corresponding to the fundamental frequency of the soil-structure system.The constant value may be computed by iterating on the eigenvalue-problem until the assumed frequency matches the computed eigenvalue.With reference to Equation 26, the constant foundation stiffness is then determined as and the foundation damping ratio in each mode may be estimated as Iterating on all modes (and not only the fundamental mode) provides the accurate mode shapes and natural frequencies for each individual mode.However, the eigenvectors are then not orthogonal and the modal responses cannot be superimposed.Although the stiffness term is assessed for the fundamental mode only, the damping ratios are determined separately for each mode.The damping ratio    is bound to maximum 30% of the computed value, as suggested by Ref. [38].
The presented impedance matrix is a closed-form solution limited to homogeneous deposits, while the considered profile varies linearly from a certain depth.However, as is the case in many situations, the upper part of the soil is homogeneous.Since the upper part is considered most important in terms of pile behavior (at least for transverse response), the linear impedance matrix may be used as an approximation.The authors recognize that there are better models for pile group impedances in general. 39,40However, given that the primary emphasis of this paper is on practical methods suited for the industry, the intention is to showcase and assess the very simple model formulated in Cemalovic et al. 21Furthermore, although the rotational-translational coupling is generally important for batter pile groups, the diagonal matrices are often used in practical engineering, especially for simple structures and/or in the preliminary design stages.

Damping
Rayleigh damping is commonly used for linear analysis.For nonlinear analysis however, adding damping is not a straight-forward task since (1) material damping is an inherent part of the material model and ( 2) the stiffness matrix is continuously changing throughout the system response.Moreover, the issue of damping in SSI-problems is further complicated due to the fact that the superstructure and foundation-soil region usually demand distinct damping ratios.To take into account the aforementioned issues, the damping matrix is determined by updating both the Rayleigh-coefficients and the tangent stiffness matrix at the last converged step.Furthermore, the mass and stiffness matrices are split into subdomains, enabling the utilization of different damping ratios for the foundation and superstructure.The damping matrix is thus expressed as where subscripts  and  indicate the superstructure and foundation-soil region, respectively, and   refers to the last converged tangent stiffness.[43]

Validation
The nonlinear time domain solution is validated using a fully coupled, three-dimensional finite element model constructed in OpenSees MP 44,45 and STKO 46 as described by Cemalovic et al. 20,47 It should be noted that the contact elements available in OpenSees MP produce in some cases rather noisy results for transient analyses.This may occur even for integration schemes with numerical damping such as the Hilber-Hughes-Taylor method. 48,491][52] The macroelement is thus calibrated using a model without contact elements (for the validation case only).Due to symmetry, only half of the model is considered in both models.The nodes corresponding to the plane cutting the mid-span are fixed in the vertical direction, but are free to rotate and translate in the horizontal direction.The three-dimensional finite element model is validated in Cemalovic et al. 47 Furthermore, the superstructure in both models is represented using linear beam elements, which confines nonlinear behavior to the foundation and soil.Figures 3A-D compare the macro-element solution against OpenSees MP in terms of deck displacements and accelerations when the deck is subjected to a horizontal, harmonic load for  equal to 1 and 3 Hz.The results show that the macro-element solution performs quite well for both frequencies.Figures 3E and F compare the macro-element solution against OpenSees MP in terms of deck displacement and acceleration for seismic loading.The input motion is the Kocaeli Gebze time history record with  ,30 = 790 m/s (soft rock).The input motion is scaled to  = 1 g.For the fully coupled OpenSees MP model, the displacement history is applied at the base (rock).For the macro-element solution, a nonlinear free-field analysis is first performed using the same modeling strategies as the validation model.The free-field displacement history is then applied at the base node of the macroelement, which inherently neglects kinematic interaction.The results indicate that macro-element is able to simulate the seismic response with reasonable accuracy, admittedly with some overestimation of deck acceleration.It is worth mentioning that the OpenSees model uses the Hilber-Hughes-Taylor method, while the macro-element solution uses the Newmark's method, where the former is more effective in damping out high-frequency response.
As previously mentioned, analyzing soil-structure problems using macro-elements inherently introduces an approximation since the analysis is fundamentally nonlinear, while kinematic and inertial effects are not assessed simultaneously.Consequently, the results presented in Figures 3A-D demonstrate the presented macro-element performance, while Figures 3E and F illustrate, at least suggestively, the accuracy of the macro-element approach in general.

INCREMENTAL DYNAMIC ANALYSIS
This section present a series of IDAs using the in-house code in order to evaluate the effect of linear and nonlinear SSI, batter angle and pile-spacing.IDA is a parametric analysis method that evaluates the structural performance under seismic loads, where the computational model is subjected to one or more ground motion records scaled to different levels of intensity. 54The results (IDA-curves) provide a relationship between a response parameter and intensity level, which may further be analyzed in a statistical sense.In the present context, IDA-curves provide an excellent demonstration of the developed solutions and how they might influence seismic design.It should be emphasized that the IDA analyses are not performed in a traditional sense (which puts the focus on the structure), but rather utilizing the concept for a practical presentation of the results.Furthermore, the batter angle and pile spacing is changed without making any other changes to the design of the system in order to clearly show the effect of batter angle and pile spacing in a general sense.
An important part of the SSI-analysis is the assessment of foundation input motion (FIM).Cemalovic et al. 47 explored the feasibility of using nonlinear kinematic interactions factors to estimate FIM, but this approach was based on the results of comprehensive finite element analyses.If this step were to be included in a simplified method that emphasizes practicality, it would make the overall solution strategy more complicated, potentially causing it to lose its intended purpose.The analyses presented herein are therefore performed by enforcing the ground motion displacement histories at the base, which inherently neglects the kinematic response.Furthermore, since the analyses only involve models where the input motions are directly enforced at the foundation level, the ground motion records are selected such that they represent soft soil.This is achieved by limiting  ,30 to maximum 200 m/s.Table 1 shows the selected ground motions records and their key characteristics.Here,  is the database number,  is the non-scaled peak ground acceleration,   is the earthquake magnitude and  ,30 is average shear wave velocity for the top 30 m of the soil.Figure 4 shows the corresponding unscaled response spectra.Each ground motion record is analyzed for eight -values ranging from 0.1 to 2.2 g.While it is acknowledged that such high  values are not realistic for soft soil, the analyses are carried out up till 2.2  to effectively illustrate the effect of inelastic, nonlinear behavior.The results are presented for nine different response parameters, namely, horizontal displacement residual, rotation, and acceleration for both foundation and deck, in addition to base level vertical force, shear force and moment.The displacement residual refers to the residual between deck/foundation and input motion.The IDA-curves are presented as the median for all considered models.The individual results and the 16 − 18% fractiles are shown for the nonlinear (macro-element) model only.
F I G U R E 5 IDA-curves showing the effect of SSI.IDA, incremental dynamic analyses; SSI, soil-structure interaction.

Effect of SSI
The IDA-curves are obtained for three different analyses; (1) nonlinear analysis using the macro-element (nonlinear SSI), linear analysis using the impedance matrix (linear SSI) and linear analysis with fixed base condition (neglected SSI).It is recognized that IDA is by definition nonlinear.However, in the context of this study, linear analyses were preformed to clearly demonstrate the effect of linear and nonlinear SSI in comparison to fixed-base conditions.Figures 5A-C show the IDA-curves for horizontal displacement residual, rotation, and acceleration of the foundation.Clearly, there is no residual displacement nor rotation at the foundation level for fixed conditions.The results show that nonlinear SSI substantially It is particularly interesting to note that linear and nonlinear SSI converge for low .This observations strengthens the validity and usefulness of the linear impedance matrix in combination with the modal approach as an efficient and easy-to-use tool for estimation of low-strain seismic response.However, it should be noted that when assessing low-strain vibration problems where accuracy is of greater importance, frequency domain solutions should be used.It is also interesting to note that both linear and nonlinear SSI increase several of the response parameters for low-to-mid range .The considered system is a relatively stiff structure founded on soft soil.The fixed-base fundamental period of the structure is   = 0.21 s, while the flexible-base fundamental period (using the linear SSI-approach) is T = 0.36 s.This yields which is considered to be a large period lengthening.The detrimental effect of SSI is perhaps expected when considering the response spectra shown in Figure 4.It is clear that several of the ground motions attain larger spectral accelerations for periods longer than the fundamental fixed-base period.Introducing SSI, it is thus very likely that the fundamental mode (and consequently any higher mode), contributes to an increase in spectral acceleration.Even for the nonlinear case, which may cause even larger period lengthening, but simultaneously increase the overall damping, SSI increases several of the response parameters for low-to-mid range .
To further clarify the discussion above, the IDA-curves for deck acceleration, base shear and base moment are normalized by the corresponding fixed-base response and plotted against .The results are shown in Figure 6.The lines showing fixed-base and linear SSI are of course constant with respect to the normalized response, where the former is equal to unity.It observed that linear SSI increases deck accelerations and shear forces by more than 50%.Base moments are increased by approximately 35%.Nonlinear SSI converges with linear SSI as  tends towards zero, that is, when the macro-element responds relatively linear.As  increases, nonlinear SSI yields larger period lengthening, but also additional system damping.The normalized response decays, but it is still greater than unity until  reaches approximately 1 g.Past this point, nonlinear SSI yields lower values compared to fixed-base conditions.Figure 6 gives a rather clear overview of how linear and nonlinear SSI affect the respective response parameters in terms of ground motion intensity.In essence, these plots are merely an alternative representation of the typical IDA -curves.

Effect of batter angle and pile spacing
The effect of batter angle is evaluated by comparing the results (macro-element approach) using four different combination of batter angles  1 and  2 with constant pile spacing  0 = 5  .Figures 7A-C show the IDA-curves for horizontal displacement residual, rotation, and acceleration of the foundation.Increasing batter angle decreases the displacement residual and increases rotation for high .It is interesting that increasing batter angle decreases rotation for low  and increases foundation acceleration for moderate and high .The latter may be attributed to the fact that batter piles increase the elastic domain in the horizontal direction, and hence also the accelerations as  increases.The effect of pile spacing is evaluated by comparing the results (macro-element approach) using three different pile spacing's (3  , 5  , and 10  ) with vertical piles only.Figures 8A-C show the IDA-curves for horizontal displacement residual, rotation, and acceleration of the foundation.Increasing pile spacing decreases the foundation rotation and has a negligible effect on foundation acceleration.There are no clear trends with respect to displacement residual.Figures 8D-F shows the IDA-curves for horizontal displacement residual, rotation, and acceleration at the deck level.Similar trends are observed as for foundation response.Figures 8G-I shows the IDA-curves for vertical force, shear force and moment at the base level (top foundation).Increasing pile spacing decreases vertical forces, slightly increases shear forces (at least for high ) and increases moments.

Perspectives
The macro-element formulation neglects the frequency-dependent radiation damping.For larger displacements, hysteric damping is expected to dominate the overall damping, and the macro-element may yield reasonable result without further modification.It is possible to approximate radiation damping using the expressions available in the literature. 33However, these expressions are strictly not applicable to nonlinear time domain solutions because they are dependent on both frequency and soil stiffness.As an approximation, radiation damping may be added by appropriately selecting frequency and soil stiffness for the task at hand.Since radiation damping increases with frequency, the fundamental frequency may be used in order to minimize potentially excessive damping.The constant soil stiffness may chosen to match the expected shear strain level.It should be noted that radiation damping vanishes for frequencies lower than the fundamental frequency of the soil deposit, or is negligible in the case of a rigid bedrock at relatively shallow depth.Therefore, whether or not to include radiation damping should depend on the case considered.A dedicated study on the implementation, as well as the importance, of radiation damping in the context of the macro-element solution would be valuable.The macro-element formulation also neglects pile-soil-pile interaction, which may decrease the accuracy of the macroelement when the piles are closely spaced.The limited set of analyses presented by Cemalovic et al. 20 showed that the macro-element performed well for pile spacing's down to 3  , but the analyses were limited to soft soil where it is expected that pile-soil-pile interaction is relatively low. 55,56In the context of nonlinear macro-elements, including pile-soil-pile interaction should consider the effect of yielding, preferably through the use of simplified, closed-form solutions.
The macro-element formulations assumes that the axial and transverse response are de-coupled, which inherently introduces an error related to the bounding load.Considering that the implicit transverse-rotational coupling was implemented successfully, similar strategies may be explored for transverse-axial coupling.
The interaction between the bridge abutment and the backfill soil is an important part of the overall seismic response.First, the backfill soil causes additional input motion along the abutment wall.Second, the backfill-interaction introduces additional inertia loads on the bridge, which in some cases may be of considerable magnitude, depending on the size of the backfill soil. 57Third, the backfill-soil may yield additional hysteric and radiation damping.A practical and simplified approach accounting for the above-mentioned effects would be particularly useful in the assessment of embedded (or partly embedded) structures such IABs and culverts.
Both the linear impedance matrix and the macro-element are restricted to planar analysis.It would add great value if the formulations, and particularly the macro-element, were extended to three dimensions.Furthermore, the analysis presented herein were restricted to nonlinear foundation response.In reality, nonlinear behavior extends to all parts of the structure.The in-house software may easily be extended to include inelastic behavior of the superstructure using either lumped plasticity models or fiber sections.

CONCLUSION
A new finite element framework for seismic analysis of IABs (and other simple structures) accounting for linear and nonlinear SSI was presented.Although the software was implemented and demonstrated for relatively simple bridges frequently encountered in everyday engineering, the solution is valid for any type of structure that may be represented by planar frames.The previously developed macro-element and linear impedance matrix were implemented in a time and modal domain solution, respectively.The solutions were validated using rigorous numerical models constructed in OpenSees MP.Detailed pseudo-codes describing the implementation strategy for the macro-element were provided.The software was demonstrated by performing an extensive set of IDAs, where the effect of linear and nonlinear SSI, batter angle, pile spacing and asymmetric configuration were addressed and discussed.The results revealed that the macroelement and linear impedance matrix converged as  tended to zero, which strengthened the validity and usefulness of the linear impedance matrix in combination with the modal approach as an efficient and easy-to-use tool for estimation of low-strain seismic response.Furthermore, linear and nonlinear SSI increased several of the response parameters for low-to-mid range .This was attributed to the fact that the stiff structure founded on soft soil produced large period lengthening when accounting for SSI, while several of the ground motions attained larger spectral accelerations for periods longer than the fundamental fixed-base period.Increasing batter angle decreased rotation for low  and increased acceleration for mid-to-high  (for both foundation and deck).The latter was attributed to the fact that batter piles increase the elastic domain in the horizontal direction, and hence also the accelerations as  increases.
The ability to efficiently perform IDA analyses for a number of key variables is a powerful tool in the early design stage of any structure prone to seismic loading.In the context of macro-elements, or any simplified approach for that matter, such analyses allow the engineers to efficiently form a clear overview of how a given change in the design will affect the seismic performance of the structure.This is particularly useful when considering that the complete design of structures reaches far beyond seismic design, and involves several other demands that must be fulfilled.

A C K N O W L E D G M E N T S
The first author would like to express gratitude to Sweco Norway AS and The Research Council of Norway for financial support.

C O N F L I C T O F I N T E R E S T S TAT E M E N T
The authors declare no potential conflict of interests.

D ATA AVA I L A B I L I T Y S TAT E M E N T
Data sharing is not applicable to this article as no new data were created or analyzed in this study.

1 if 1 ,
, the displacement of the single pile  in global coordinates is A L G O R I T H M 4 Return mapping scheme for axial tangent stiffness   , procedure TANGENT_RETURN_V while |  | <   and   <   do  =  +   =  −1 , Δ = 0,   =   else if   ≤   then if   > 0 then if (  >   ) and (  >  −1 ) and (   −1 > 0) then if   > 0 then    (  −     ) ) , ] ln(  ) else if   = −1 then RC →  pl  =  pl 0, [( ln(  ) + ln ( 2      (  −    ) is the horizontal displacement of pile  due to horizontal loading of pile ,    is the horizontal force in pile ,   , is the interaction factor representing the additional horizontal displacement of pile  due to a horizontal displacement of pile  and   is the horizontal pile group displacement.The  equations are solved for  horizontal forces    as functions of   , and the horizontal complex impedance is determined as can express the vertical displacement as a function of vertical forces and corresponding stiffness components and interaction factors, that is, =

F I G U R E 3
Comparison of the macro-element solution against OpenSees MP.

F I G U R E 6
IDA-curves normalized by fixed-base response.IDA, incremental dynamic analyses.increases the foundation displacement residual and rotation.SSI also increases foundation accelerations for most values.Figures 5D-F show the IDA-curves for horizontal displacement residual, rotation, and acceleration at the deck level.Linear and nonlinear SSI yield approximately the same deck displacement residual up to a certain value.As  increases, nonlinear SSI yields substantially larger values.Nonlinear SSI reduces deck rotations and accelerations compared to linear SSI and fixed conditions as  increases.Figures 5G-I show the IDA-curves for vertical force, shear force and moment at the base level (top foundation).In all cases, nonlinear SSI reduces forces and moments as  increases.

F I G U R E 7
IDA-curves showing the effect of batter angle.IDA, incremental dynamic analyses.

Figures
Figures 7D-F show the IDA-curves for horizontal displacement residual, rotation, and acceleration at the deck level.Similar trends are observed as for foundation response, but with somewhat less difference between vertical and batter pile groups.Figures7G-Ishow the IDA-curves for vertical force, shear force and moment at the base level (top foundation).In all cases, increasing batter angle increases forces and moments as  increases.It is also observed that the asymmetric configuration ( 1 = 15 0 and  2 = 0 0 ) is very close to the symmetric configuration ( 1 =  1 = 7.5 0 ) for several response parameters.The effect of pile spacing is evaluated by comparing the results (macro-element approach) using three different pile spacing's (3  , 5  , and 10  ) with vertical piles only.Figures8A-Cshow the IDA-curves for horizontal displacement residual, rotation, and acceleration of the foundation.Increasing pile spacing decreases the foundation rotation and has a Return mapping scheme for horizontal tangent stiffness   , procedure TANGENT_RETURN_H while |  | <   and   <   do ,   =  −1 , Δ = 0,   =   else if   ≤   then   >   ) and (  >  −1 ) and (   −1 > 0) then A L G O R I T H M 3 53lected earthquake records from PEER Ground Motion Database53.
TA B L E 1