Assessment of the design spectrum with aggravation factors by 2D nonlinear numerical analyses: a case study in the Gemlik Basin, Turkey

The response spectra of multidimensional analyses are compared with one-dimensional (1D) local models to couple the irregular soil stratification effect at a site. In recent studies, the surface motion 2D/1D or 3D/1D spectra ratios are defined as the spectral aggravation factors for each region of a site. Particularly in alluvial basins, where the soil media is typically formed by fault ruptures or topographic depressions filled with sediments, the inclination of the rock outcrop at the edge of the basin has a considerable effect on the site response, and such an effect has not yet been taken into consideration in recent seismic building codes and general engineering applications. In this study, the natural alluvial basin near the North Anatolian Fault in Gemlik, Marmara region, Turkey, was investigated by 40 seismic site tests and 4 validation borings. The 2D and 1D nonlinear response history analyses in the north–south and east–west directions in the Gemlik basin were performed by numerical models on a finite difference scheme considering nonlinear elastoplastic material behaviors and geometric discontinuities. Twenty-two strong ground motions recorded at the rock site were excited vertically as SH waves. The numerical results showed that narrow basin effects were derived not only by reflection, refraction, and shifting behavior but also by the focusing and superposition of the seismic waves propagating from opposite basin edges. As a result, the site-specific spectral aggravation factors SAF2D/1D defined by the ratio between the 2D and 1D acceleration response spectra for each period and any location on the site were proposed for the Gemlik basin. The variations in the aggravation factors were observed as values increasing to 1.2–2.2 near the edge and at the basin center.


Introduction
The design of earthquake-resistant structures for living quarters is one of the significant objectives of civil engineering. The prediction of hazards to all living facilities to protect lives by reducing the destructive effect of earthquakes has always been crucial in terms of engineering. Estimating a design earthquake is the main phenomenon in building design and has a crucial role in examining the existing building performance to minimize the loss of life and property. Site response investigations performed in recent decades have pointed to four main aspects that shape strong ground motions. The first aspect is the amplification of displacement, which exists when a seismic wave progresses through an interface to lower rigidity layers from higher rigidity layers. The second is the resonance of the flat layers developed mechanically at specific frequencies. The third originates from the nonlinearity of soil stress-strain behavior and the nature of inhomogeneity and anisotropy in the material. The last effect derives from the wave propagation variation in soil half-space, which has multilayered site conditions with stratigraphic heterogeneities.
In early studies, to reveal the effect of plane-incident waves in soft two-dimensional basins, Aki and Larner (1970) and Wong and Trifunac (1974) estimated the surface motion of valleys with perfectly elastic material for incident plane SH rays by a semianalytical method. Smith (1975) performed analyses using finite difference and finite element techniques to study the effects of irregular layer interfaces systematically. Bard and Bouchon (1980) examined the formation of surface waves and edge effects. King and Tucker (1984) noted significant differences in surface peak horizontal accelerations at the center and near the edges of the valleys. Yamanaka (1989) performed an in situ investigation and numerical analysis to study the propagation of seismic waves within the deep sedimentary layers of the southwestern Kanto district, Japan. Papageorgiou and Kim (1991) investigated the effect of bedrock slopes on amplifications by establishing a 2D model in the Caracas basin in Venezuela during the Caracas earthquake on 29 July 1967. Zhang and Papageorgiou (1996) worked in the Marina basin in San Francisco, California, to predict the impact of the Loma Prieta earthquake on 18 October 1989. Pei and Papageorgiou (1996) studied the surface waves created by motions traveling up from the basin bottom based on records from Gilroy seismograph arrays placed on the surface along the Santa Clara basin in California. Kawase (1996) used the 2D finite element model to analyze the basin edge effect along the damaged zone in the Kobe basin observed in the Hyogoken Nanbu earthquake on 17 January 1995. Graves et al. (1998) suggested that the increase in amplification was directly caused by surface waves generated from the basin edge. Bielak et al. (1999) evaluated soil amplification and structural damage together in a small valley in Kirovakan depending on basin conditions. The performed 1D wave propagation analyses could not provide sufficient results for the existence and spatial distribution of the 1988 Armenia (Spitak) earthquake damage over a wide area in Kirovakan. In the progressive stages of multidimensional response studies, the concept of the aggravation factor was proposed by Chávez-García and Faccioli (2000). Bakir et al. (2002) analyzed the strong ground motion in the Dinar district, where the 2D model was established on the edge of an alluvial basin in southeastern Anatolia during the Dinar earthquake. Somerville and Graves (2003) asserted that commonly used empirical approaches do not reflect the additional effects in sedimentary basins. In relevant studies, namely, Semblat et al. (2005), Iyisan and Hasal (2011), Iyisan and Khanbabazadeh (2013), Abraham et al. (2016), Khanbabazadeh et al. (2016), Riga et al. (2016), Makra and Chávez-García (2016), Chávez-García et al. (2018), Stambouli et al. (2018), Cipta et al. (2018), Moczo et al. (2018), Zhu et al. (2018), the results found in the two-dimensional numerical analyses support this situation. In the current researches, Hasal et al. (2018), Khanbabazadeh et al. (2016, 2019), and Ozaslan et al. (2020 studied time-domain dynamic analyses in idealized 2D models of the Dinar and Düzce basins.
This study presents the effects of the heterogeneities in both vertical and lateral directions on the local seismic response by concentrating on the earthquake response analysis of the basin, which is laterally confined and filled with sediment. In this aspect, nonlinear soil behavior was examined under risk-targeted levels consisting of design earthquakes and maximum considered earthquakes. For this purpose, the proposed procedure is based on the correlation between detailed numerical modeling solutions to extract the contributions due to 2D effects in addition to 1D soil behavior. The impacts of a multiaxial stress state in the soil, formed by a 2D-3D decomposed shallow stratigraphy, need to be investigated since they can play a remarkable role on the resulting nonlinear strains. Thus, in this study, comprehensive in situ investigation and numerical analyses were performed to study the site response and propagation of seismic waves within the shallow sedimentary basin in the Gemlik district of the western Marmara Sea, Turkey. In the scope of the study, the effects of lateral discontinuities of the basin geometry and bedrock slope on design earthquake accelerations were investigated by nonlinear numerical analyses on natural basin conditions.

Seismic array tests and borings
The Gemlik basin in Bursa, located between latitudes of 40°26′20″-40°25′15″ and longitudes of 29°09′06″-29°11′10″, has been determined as the research area. The basin, which is approximately 1.7 km wide in the north-south direction and 2.5 km long in the east-west direction, is enclosed by two fault segments on the southwestern branch of the North Anatolian Fault. Major faults in the Anatolian Plate and the fault segment passing from the southern border of the research area in the east-west direction are given in Fig. 1. The alluvial deposit and bedrock in the basin have been investigated by a large number of microtremor array measurements and validation borings. Site experiments also includes one active-source test Multichannel Analysis of Surface Waves (MASW), which has been performed due to the lack of enough plain area.
Microtremor array measurements were made at a total of 40 different research points with two instrument sets. The recorded microtremors were analyzed by performing the SPAC method using Fortran codes to create a detailed underground model of the basin (Okada 2003;Yamanaka 2005). In addition, confirmatory drillings and Standard Penetration Tests were performed at 4 sites. The coordinates of all site investigations are presented with emphasis on the methodologies and tools employed in Figs. 1 and 2, respectively. In the seismic experiments, a circular array consisting of two equilateral triangles having three corner stations and three mid-point stations located around the central station was preferred because of its feasibility in highly intensive residential areas. In this method, the phase velocities of Rayleigh waves are computed on cross-correlations between each station pair on the circular array by analyzing the vertical components of microtremors. The observed phase velocities are employed to estimate the shear wave velocities and layer thicknesses by performing inversion with the Genetic Algorithm and Simulated Annealing method, which generates optimization in the mutations, crossover, and selection of individuals in a population (Yamanaka and Ishida 1996a, b).
The applied sizes of the array circle ranged from 20 to 50 m in radius. The phase velocities of Rayleigh waves were estimated by 20-30 min-long microtremor data, and the shear wave velocities of the soil layers were determined by the inversion method. The soil structure consists of two layers lying on the bedrock with mean S-wave velocities are approximately 180 m/s and 500 m/s, which increase with depth as shown in Fig. 3. Microtremor   Fig. 3, the phases of the SPAC method performed at point MA-BBB9 are shown. The colored curves on the center of the figure give phase velocities, and the circles on the red curve represent data points forming the dispersion curve. Since dispersion is a function of phase velocities that varies depending on the frequency, it provides the shear wave velocity of the soil structure by the inversion phase. In Fig. 3, given Vs-profile was obtained from the inversions with one standard deviation indicated by the dotted lines. However, there is doubt about whether a circular array is sufficient for the application of the SPAC method. Many previous researchers using the SPAC method have employed a circular array with one station at the center and the other three at the circumference (Kudo et al. 2002;Matsuoka et al. 1996;Okada 2003). Yamamoto et al. (1997) and Okada (2006) concluded that a minimum 3-station array provided phase velocities in the frequency range from 5 to 15 Hz with an error of approximately 5% relative to those calculated from a known shear wave velocity profile at an experimental site.
The accuracy of the soil section models created in the study was controlled by the consistency of the data collected in both seismic and penetration tests. For this purpose, BH1, BH2, BH3, and BH4 investigation borings and penetration tests were performed at the same points as the microtremor measurements, MA BBB-5, MA BBB-6, MA BBB-7, and MA BBB-10, which were made in the basin center and near the edges. In the borings, BH1 and BH2 medium-to high-plasticity clay and silty clay units are observed within depths of 8-12 m. The unit in question is generally very weathered, very weak, and fragmented colluvial soil because it is in the fault zone. It is a fully weathered residual soil (sandy silty clay), and the weathering degree decreases toward the end of the borings. In the third borehole, a silty clay unit with a fine sand band that is grayish in color with medium to high plasticity passed between 2 and 35 m. These units are Quaternary-aged alluvium consisting of gravelly coarse metasedimentary soil and metavolcanic rocks. Under the alluvial layers, there is coarse-grained soil (sandy gravelly clay) deeper than 35 m until the end of the borings, and it is predicted to have formed after the complete weathering of the Triassic Abadiye Formation that is part of the Karakaya Group which consists of meta-claystone, metasiltstone, metasandstone, and metalimestone. In the 4th borehole, shells belonging to sea creatures are observed in the boring sample after 20 m in the formation (Okay and Göncüoğlu 2004). The alluvial deposits in the main lithological units detected in boring logs BH1-BH4 in the Gemlik basin are given in Fig. 4. The correlations suggested by Imai (1977), JRA (1980), and Iyisan (1996) were used to calculate the variation in shear wave velocity of the Standard Penetration Test data. The average values of the results calculated by three different SPT correlations V s = 102N 0.292 , V s = 100N 0.33 , and V s = 51.5N 0.516 (where N is the number of blows) were compared with the results of the seismic tests in Fig. 5. In borings BH3 and BH4, which are closer to the basin center, the standard penetration test blow counts (SPT-N) have lower values that range from 5 to 15 in the upper layer within depths of 10-40 m, while these values increase to 20 to 40 after 35 m of depth. In the borings located near the edges of the basin, the mentioned increase starts at a depth of 10-20 m.

The soil condition and basin geometry
A geographical information system tool (GIS) was used to integrate all collected data on storing, evaluating, and presenting the field investigations. It provided the visualization of the layers created with the interpolation analysis of the spatial locations of the whole field study in terms of 2D and 3D maps and helped the preparation of sections by establishing relations between data. The variation in the bedrock depth along the bottom of the basin detected in simultaneous microtremor array measurements is presented in Fig. 6. In The Gemlik basin is a relatively shallow sedimentary structure with gentle slopes in the east-west direction, while the north-south direction has highly inclined rock outcrops on the edges. The fault segment passing from the southern border of the research area shapes the slope of the bedrock/sediment interface, thereby increasing the α 3 and α 4 angles at the southern edge of the basin illustrated in Fig. 7. The basin has a complicated geometry, where the depth of the soil deposit is estimated to be 140-180 m at the center, and the basin overlies rigid bedrock and is surrounded by an inclined rock outcrop. The 2D models of the shear wave velocity exhibited bedrock inclinations of the basin edges that extend α 1 = 18° and α 2 = 10° in the north direction and α 3 = 16° and α 4 = 32° in the south direction, and width of the flat interface between the bedrock and the soil deposit is about 400 m. In the east-west direction, there is a broader flat base having 1400 m width and edge slopes of α 5 = 9° and α 6 = 10° in Fig. 7.

Description of the finite difference-based nonlinear method
In the nonlinear 1D and 2D response history analyses of the basin, the explicit finite difference method was performed by the Fast Lagrangian Analysis of Continua 3D (FLAC3D code). Contrary to previous studies, this method provides elastoplastic soil nonlinearity under shear and compressional wave propagation by considering the straindependent nonlinear constitutive rule and yielding criteria at any time during dynamic In this way, the averaged strain rate calculations over all subzones in the soil space mesh are performed before any computation steps by regarding constitutive law functions of the soil materials without any other damping requirements.
In the numerical analysis, the soil constitutive model needs to reflect nonlinear elastoplastic material properties under cyclic loading in the time domain. In the applied constitutive model, energy dissipation is provided by the hysteretic model with the degradation of the shear modulus (G/G max ) and cyclic damping (D) at small strain levels. Furthermore, the plastic deformations of soil materials at high strain levels are determined by the Mohr-Coulomb model. The combination of the strain-dependent damping ratio and secant modulus functions are derived by the given equations and illustrated in Fig. 8a. A loop tracked on an initial cycle of unloading/reloading is illustrated in Fig. 8b.
Subsequently, the yielding level of the hyperbolic rule must be higher than the Mohr-Coulomb yield stress. The state is provided with the condition given in the equation below (Cundall 2001).
where γ ref is the ultimate value of τ m /G maks in the hysteretic function, G max is the maximum shear modulus, and τ m and γ m are the constant yield stress and shear strain, respectively.
In the elastic range γ c < γ m , the modulus reduction factor, is defined by Eq. 3: In the plastic range, γ c ≥ γ m, The energy dissipation in a cycle, ΔW, is expressed as the total contributions from the elastic ΔW H and plastic ΔW MC ranges.
The maximum stored energy, W, and the damping ratio, D, in a cycle are given by the equations below.

Boundary conditions and damping
The entire range of problems potentially encountered in geotechnical engineering exists in semi-infinite soil space, and the solutions are performed on simulated finite media by discretization with numerical methods. Thus, the boundary conditions for the solution of the problem are crucial and need to be provided suitably with the surveyed wave propagation for infinite site conditions. In the dynamic analysis, advanced boundary conditions, which are improved to prevent the waves from being trapped inside the model without reflections that exist at the finite model borders, are executed and given in Fig. 9. The effectiveness of this type of energy-absorbing boundary has been demonstrated in both finite-difference and finite-element models (Cundall 2001). At the bottoms of the models, the quiet boundaries that prevent the reflection of outward-propagating waves back into the model were assigned. On the lateral boundaries of the models, the nonreflecting free-field boundaries were settled in a continuum finite difference scheme by coupling the main grid to the free-field grid by viscous dashpots. The assigned dashpots produce viscous normal and shear stress tractions along the model boundaries by Eqs. 11-12.
where t n and t s are the normal and shear stresses tractions, ρ is the mass density, C p and C s are the pressure (P) and shear (S) wave velocities, and v n and v s are the normal and shear components of the velocity at the quiet boundary.
where F x , F y, and F z, are gridpoint tractions of the free-field boundary; v x m , v y m and v z m are x, y, and z velocities of gridpoints in the main grid at the side boundary; v x ff , v y ff and v z ff are x and y velocities of the gridpoint inside the free field; A is the area of influence of the freefield gridpoint; F x ff , F y ff and F z ff are free-field gridpoint forces with contribution stresses of the free-field zones around the gridpoint (Itasca 2017).

2D and 1D models of the Gemlik basin
The soil material in the basin has a low shear wave velocity (V s ) ranging from 160 to 220 m/s in the near surface, and the underlying layers become more rigid with increasing V s from 300 to 550 m/s due to the influence of the effective vertical stress until the bedrock with V s greater than 750 m/s at the base. The shear wave velocities of the sublayers defined in the models are the average of the data collected along the sections. The basin consists of alluvial deposits and old river and sea sediments composed of sandy-silty clay layers with medium plasticity index (PI) values are in the range of 15-25%. The significant V s contrast between deposits and bedrock is not very common in real basin conditions, even if the sediments are poorly consolidated soils overlying bedrock, because of the increasing overburden pressure on the deeper layers (Bard and Bouchon 1980;Zhu and Thambiratnam 2016). The final models were generated by constituting sections in the north-south (A-A) and east-west (B-B) directions that  Table 1. Following the research topic, the preferred software allows investigating the wave propagation produced by multiple physical phenomena, such as refraction, reflection, and resonance, in infinite soil media subjected to nonlinear soil behavior. In 2D and 1D analysis, when performing this type of analysis, the sizes of the discriminated elements in the soil media need to permit the transition of the applied seismic waves. Similarly, the smallest time step required for each calculation needs to ensure the propagation of the highest frequency component of input motion in the model. Only under these conditions can the motion be transferred accurately to the surface throughout the defined finite media. Figure 10 presents free-field boundary zones and the finite difference scheme of the Gemlik basins in 2D plane strain models, which are built and investigated by numerical analyses. In this way, plane waves propagate upward and sustain no distortion at the boundary because the free-field zones supply identical information to those in an infinite model.   Kawase and Aki (1989)

Verification of the wave propagation
In this study, a trapezoidal model that is identical to that indicated by Kawase and Aki (1989) was produced to validate wave propagation. The same acceleration pattern was captured across the model surface and is given in Fig. 11. The result of the Kawase and Aki (1989) analysis was also tested by Iyisan and Khanbabazadeh (2014) and Gil-Zepeda et al. (2003). The properties of the materials in the model have been defined by 1000-2500 m/s shear wave velocities (V s ), and the unit weight of materials have been set as the same to provide a constant impedance value concerning the compared model. The shear and bulk moduli were considered by taking the Poisson ratio as 1/3. The propagation and distortion of the wave, the change in the wave characteristics in the material environment, and the spectrum distribution have been investigated using the wavelet. As the input motion, the Ricker pulse has provided the opportunity to definitely examine wave propagation for different frequency (f c ) ranges and defined amplitudes u(t). Consequently, it was determined that the numerical method and boundary conditions properly provide wave propagation, distraction, and reflection that are the same as the verification model. Two-dimensional plane strain and 1D soil column models of the Gemlik basin have been built and investigated by the same numerical method. When generating the finite difference scheme, the maximum defined zone size (l max ) has been restricted to be equal to or less than 1/10 or 1/8 of the lowest wavelength (λ min ) defined in the model. In other words, it depends on the wavelength of the wave with a maximum frequency (f max ) transmitted in the softest layer of the media. The zone dimensions are determined by the equations given below (Cundall 2001;Itasca 2017).
In the models, considering the assumed relations, the maximum zone sizes (Δl max ) are regulated as 1-5 m by fitting to the shear wave velocity of the soil layers by seeing that the earthquake records have the highest frequency components, up to 15 Hz.

Input motions
Computing response spectra for several varying strong ground motions and averaging them lead to a smoother target spectrum. Producing such smoothed spectra is a significant step in improving a design spectrum. The ground motions are defined probabilistically as risk-targeted spectra by current seismic code provisions through seismic hazard analyses. The results of hazard analyses are used to determine the site-specific Maximum Considered Earthquake (MCE R ) and the Design Earthquake (DE R ) spectrum and the sitespecific design acceleration parameters of the short-period SDS and 1-s period SD1. The underlying methods of site-specific ground motion analysis are necessarily highly technical and require a unique combination of geotechnical, earth science, and probabilistic expertise (FEMA 2020). Across the basin models, both 1D and 2D site response analyses were performed under the sets of ground motion data that were selected by matching to the level of (16) min ≤ V smin f max (17) Δl max ≤ min 10 the MCE R spectrum and the level of the DE R spectrum defined by exceedance probability of 2% in 50 years and 10% in 50 years, respectively, in Fig. 12. The input motion selection also considers near-fault and transitional regions with given distances in the specification of the earthquakes in Table 2. In total, 22 earthquakes were selected by matching two levels of target spectra. Each strong ground motion set contains 11 earthquakes filtered by a 25 Hz low-pass filter and corrected baseline (Bommer and Acevedo 2004;Katsanos et al. 2010). To extract the effect of soil layers from selected accelerograms, the motions recorded on stiff soil and rock sites during real earthquakes were chosen. Strong ground motions were selected from the PEER ground motion database, COSMOS Virtual Data Center, and AFAD earthquake catalog.

Comparison of 2D and 1D basin response
Two seismic code-based levels of excitations have been used in models to obtain the seismic response of the Gemlik basin in two directions. Thus, not only 2D and 1D models but also the maximum accelerations in input motions have been compared. In the assessment stage of the processed data, SH wave propagation in one-dimensional soil columns created for each point that was selected with 50 m intervals throughout the basin surface was investigated. In 2D basin models, the effects of the stratigraphic two-dimensional discontinuity produced by refraction and reflection of SV waves were analyzed. Furthermore, the surface motions on 2D models were recorded by synthetic seismographs located at equal intervals on the surface, similar to 1D analyses at 50 m intervals. Consequently, the spectral aggravation factors SAF 2D/1D = Sae(T) 2D /Sae(T) 1D have been defined as the ratios between the response spectra of 2D and 1D models by considering locations and periods. A total of almost 1000 dynamic time history analyses were performed, and all results are presented in detail.
On the other hand, the effects of the distance between opposite basin edges to each other on surface motion for the narrow section have also been examined. In this study, the ratio of the depth (H) of the basin to its width (L) is smaller than H/L = 1/10, in contrast to the prior recent studies carried out by Riga et al. (2016), Khanbabazadeh et al. (2019), and Zhu et al. (2018), in which the edge effect was examined. It is significant to investigate the effects that are formed by refracted, reflected, and shifted waves from   inclined bedrock in narrow basins because the aggravation is higher than that in wider basins Riga et al. (2018).
The response spectra ratios calculated on the surface for different periods across the basin in the 2D and 1D models are given in Figs. 1 and 14. The surface acceleration spectra obtained in the artificial recorders s10, s16, s20, s24, and s30, which were placed on the surface in the point projections where the changes in the edge slope at the bottom and the basin center have been compared. In the A-A direction, where the opposite sides of the basin are closer to each other for the E7 earthquake at the DE R level, the aggravation factor increases to 2 in the basin center, especially at the T = 1 s period, while in the larger period of T = 1.2 s, it reaches 1.5 above the edge regions. Similarly, in earthquake E20 at the MCE R level, the aggravation factor reaches the maximum value at the center in the T = 1 s period and the DE R level shifts to the edge region in a larger period. In lower periods, at T = 0.2 s and T = 0.4 s, multiple peak values of the highest aggravation factors occur close to the edges of the basin, as shown in Fig. 13.
On the other hand, in the site response analysis conducted in the east-west direction, where the basin is wider than the other directions and the edge slopes are approximately In Fig. 14, artificial recorders s11, s21, s31, s38, and s45 are used to compare the results of earthquakes E7 and E20. In section B-B, where the inclination of the outcrop is relatively low, the aggravation at high frequencies remains under 1.2, and the 1D and 2D response history analysis results in the center of the basin are almost the same, contrary to section A-A. The highest aggravation factor values were calculated as between 1.25 and 1.5 in the region close to the eastern edge. This situation confirms that in accordance with wave propagation phenomena previously described in the Semblat et al. (2005), Riga et al. (2016), andZhu et al. (2016) studies, the increase in the width of the basin reduces the interference of the waves dispersed into the basin independently of the earthquake intensity levels. In addition, the decreasing edge slope reduces the effects of the refraction, reflection, and shifting produced at high frequencies in the regions close to the edge.

Maximum spectral aggravation factors for DER and MCER
It is obvious that the basin effects are 3D in nature. However, numerical analyses of 2D models are preferred due to the computation time, cost of the analysis and software. Most of the previous studies confirm that the two-dimensional (2D) basin effects are mainly caused by irregular interfaces between the soft soil layer and underlying bedrock. Although the mechanism of the basin effect is clear, there are still some uncertainties about quantitative descriptions of the influential area of basin effects, the effects of the inclined input motion, and 3D basin geometry.
Recent studies assert that it will be possible to bring the findings obtained in these 2D models, which are examined according to the variables reflecting numerous situations, into practice only with generalized findings. Otherwise, this type of specific research, which requires considerable time, cost and expertise, should be performed specifically in each basin to build earthquake-safe structures. The increase in the number of sitespecific analyses, as in this study, will direct future studies to bring the general results into practice.
For this reason, in addition to the seismic code provisions, it is assumed that aggravation factor charts grouped according to the primary variables defining the soil classes and the geometric structure of a basin by making some assumptions with the procedure followed in several engineering calculations will help the application. Thus, 1D analyses that are currently used can be included in basin response analyses by developing them with additional aggravation factors. Figures 15, 16, 17 and 18 present the highest values of spectral aggravation by grouping the strong ground motion levels consisting of 11 real earthquakes. Thus, both the effects of the differences in earthquake levels and lateral discontinuities that cause a change in the aggravation factor depending on the location in the model are explained. Figures 15 and 16 illustrate that the narrow section of the Gemlik basin in the north-south direction increases the amplifications under the input motions at both levels that cannot be neglected. Considering the results given in Figs. 17 and 18 in the east-west direction, this difference is interpreted as the superposition of the refracted and reflected seismic waves due to the interface, which has a higher slope depending on location and frequency.  Figure 15 clearly shows that the greatest spectral aggravation reaches 2 in near-edge regions x = 350-450 m and x = 1200-1450 m in the period of 1.2 s. The aggravation for the center of the basin takes values between 1.25 and c1.5, which is illustrated by the green region, at values higher than 1.5 s period. In Fig. 16, the same model under the DE R earthquake level, which has lower PGA values, produces higher aggravation values that are greater than 2 at the center of the basin. The results dramatically reveal that the levels of strong ground motion and frequency components form the aggravation deviation. The findings also put forward the crucial role of the time domain analysis of the nonlinear elastoplastic soil model.
On the other hand, the maximum aggravation factors are within 1.25-1.5 on the eastern edge of the B-B section. It is clear that the higher slope reflects the edge effects dominantly on the surface motions, and regardless of the earthquake level the behavior is similar to the 1D analysis results obtained where the bedrock outcrop does not reach the surface.
As the spectral aggravation factor charts are unique for the Gemlik basin, the maximum values of the response spectra ratios are calculated by 1D and 2D time history analyses for 22 earthquakes in both directions at all points and for each period. This spatial distribution of the maximum aggravation across the narrow direction is visible in Fig. 19.
For a 1 s period, a restricted region in the center has aggravation factors greater than 2-2.25, and the second peak is observed as 1.50-1.75 at a 1.2 s period, while on the basin edges, the aggravations are noticed as 1.75-2.00 between 1 s and 1.2 s with narrowing basin widths. In Fig. 20, which is totally different from the A-A section, it can be noticed that in section B-B, the maximum aggravations increase to 1.25-1.5 in periods within 1-1.5 s only in the region close to the eastern edge. As is clearly visible in the figures of maximum spectral aggravation, the distribution of the ground motion amplitudes in the 2D model is dominantly shaped by the interaction of scattered waves in the soil half space and basin geometry.

Conclusion
In the Gemlik basin models created in the north-south and east-west directions by seismic tests and boring investigations, the aggravation factors have been defined by site-specific response analyses performed considering the DE R and MCE R earthquake levels. The progression of the surface waves derived from both edges in the narrow direction of the Gemlik basin into the center of the basin creates higher amplification, particularly at lower frequencies, in 2D plane strain analyses with respect to the results of the 1D soil column method. In contrast, in the wider direction where the bedrock slope is lower, it is recognized that the aggravations take lower values only near the edge region, and the 2D and 1D analysis results become similar as they move away from the edge.
Consequently, the results of the research can be used to include the basin effect by aggravation factored site-specific response spectrum in addition to the values provided by current seismic codes. As the second important result, 1D site response analyses are not sufficient. In 1D estimations, dynamic properties of the soil deposit and earthquake characteristics have a remarkable effect on strong ground motion. On the other hand, in addition to mentioned effect for 1D soil column model, the discontinuities of soil layers and distracted-interacted propagation of seismic waves have a more significant role in distributing the peak ground accelerations of basins. Particularly at alluvial basins, different regions along the site surface are affected to various degrees. Therefore, the aggravation coefficients proposed for each specified region in the basin could be used more feasibly by calibrating to the 1D design spectrum. Multidimensional response analysis methods necessitate a combination of technical knowledge of geotechnical, earth science and physical wave propagation phenomena. Hence, the study method can develop basin-specific aggravating factors, and the suggested charts can be used with seismic code provisions. Finally, it is asserted that further studies that numerically define the regions where the basin affects mainly emerge and which periods are critical will significantly contribute to revealing the uncertainties about the subject.