The evaluation of the adjacency effect on the seismic response of homogenous hills

A brief review of previous studies shows that topographic features play an important role in the amplification (de-amplification) of seismic waves. Most researchers who have focused on topographic effects have paid attention to the distribution of amplification on isolated features, whereas what is seen in nature is often more complicated. It is known that 2D or 3D geological structures further contribute to the amplification of seismic ground motion. The main objective of this study is an investigation of the influence of adjacent hills on seismic responses. First, a parametric study has been performed, considering five models (consisting of one to five 2D symmetric homogeneous semi-sinusoidal hills), which are subjected to vertical incident waves. The 2D numerical results show that the largest values of spectral amplification for a series of 2D semi-sinusoidal hills are observed at the crest stations of the central hill, and it increases when the number of hills increases symmetrically on both sides of the central hill. This means that the spectral amplification increases from 3.3 for single-hill model to 4.37 for five-hill model. In the following, to validate the parametric results, three adjacent hills located in the city of Arak are studied using the horizontal-to-vertical spectral ratio analysis of microtremor recordings. The comparison with the amplification results obtained from the distinct element (UDEC) and boundary element (HYBRID) numerical approaches indicates that the central hill is more affected by adjacency, rather than the side hills. The spectral amplification for the crest station on the central hill is higher than the two other stations.

Most researchers who have focused on topographic effects have concentrated on the distribution of amplification over single hill features, whereas what is seen in nature is often more complicated. It is known that 2D or 3D geological structures and rock weathering contribute significantly to the amplification of seismic ground motion (Buech et al. 2010;Havenith et al. 2002;Kamalian et al. 2007a, b;Sohrabi-Bidar et al. 2016;Spudich et al. 1996;Ulysse et al. 2018).
In the current paper, the terms "adjacency effect" and "the effect of the adjacency" have been used repeatedly. In general, these effects describe those of nearby geological features such as hills, valleys, and sedimentary basins on the seismic response of particular sites; the present study is focused on the effect of nearby hills on the topographic amplification of a specific hill. Accordingly, the main objective of this paper is to study the effects of adjacency on the seismic response of hills.
First, a parametric study is employed in which the seismic responses of adjacent homogeneous half-sinusoidal hills are evaluated using a boundary element method (BEM) scheme named as HYBRID (Kamalian et al. 2006). Moreover, the results are expanded to natural features, using numerical simulation and field surveys on three adjacent hills located south of Arak City (Fig. 1). Arak is a populated city in central Iran, located 260 km southern west of Tehran. Arak is surrounded by mountains in the south, west, and east. This city has been affected by numerous damaging earthquakes, such as Bou'in-Zahra (1962 and 2002), Borujerd (2006), and Kermanshah (2017).
Arak City is located at 34° 5 min 30 s north and 49° 41 min 30 s east. This city is located in an embayment related to the Zagros Mountains, on a drainage basin of the Meighan Wetland.
A section of the 1:100,000 geological map of Arak City where the features are studied is shown Fig. 1 A Location of Arak City. B Elevation contour map of 3 hills. C Location of three adjacent hills H1, H2, and H3 in Fig. 2. As it can be seen in the map, the three studied hills are composed of limestone belonging to the Cretaceous.
Based on the field surveys in the region, a thin surface layer of weathered rocks covers the mountains of the region, and in the lower part, there are relatively uniform rocks of limestone. In the numerical studies, the weathered layer is not considered due to its very small thickness.
The main purpose of field studies was to follow the results of parametric study in natural features. Numerical simulations of the hills of Arak were carried out by two numerical modeling programs named as HYBRID and UDEC (Cundall 1971;Board 1991), which are taking advantage of the boundary element and distinct element methods respectively. The frequency domain analysis results of numerical studies then have been compared with the analysis of the horizontal to vertical spectral ratio curves (HVSR) of microtremor vibrations obtained from field measurements.

Methodology
The evaluation of the adjacency effect was carried out by combining a numerical approach with experimental observations. The numerical study was carried out using a 2D simulation of the adjacent hills.
HYBRID, developed by Kamalian et al. (2003a, b), is a two-dimensional boundary element code which is employed for parametric study. The application of this code for seismic response analyses of geological features has been already examined and well documented in the previous studies (Kamalian et al. 2003a(Kamalian et al. , 2007aLovati et al. 2011;Babaadam et al. 2021).
In the next part, for the numerical modeling of the three adjacent hills of Arak City, UDEC has been employed as well. UDEC is a well-known program based on distinct element algorithm and is developed by ITASCA. This application of this software for seismic studies is also very well published before (Havenith et al. 2003;Ulysse et al. 2018;Luo et al. 2020). Fig. 2 A Section of the 1:100,000 geological map of Arak City. B The red circle show studied hills (Radfar and Kohansal 2004) Moreover, single-station microtremor surveys using the Guralp CMG-6TD seismogram have enabled us to investigate the amplification frequencies of the studied locations using horizontal to vertical spectral ratio (HVSR) of the vibrations. In the following, all these methods are described in details.

Numerical modeling
Numerical modeling is currently widely used to simulate seismic wave scattering problems in complex environments and has therefore become a common part of site effect studies, as it allows us to analyze the simultaneous effects of various parameters on the seismic response.
Some of the best-known numerical approaches to assess the seismic response of topographic features can be categorized as follows: i) Finite difference method (FDM), which is a method with computational efficiency and simple code structure, however, could present difficulties with the modeling of irregular geometry, (Oprsal 2002;Republic 2007;Blanch 2012). ii) Finite element method (FEM), which can accurately capture the effect of complex topographies, such as mountainous areas, without refining the grid resolution, but has several disadvantages such as the need for huge power of calculation and the need to apply numerical damping (Amorosi et al. 2010; Jahromi and Karkhaneh 2019). iii) Boundary element method (BEM), direct (DBEM) (Chaillat et al. 2012;Kamalian et al. 2003a, b;Mogi and Kawakami 2007;Sohrabibidar et al. 2010), and indirect (IBEM) (Griffiths and Bollinger 1979;Sánchez-Sesma et al. 1993;Gil-Zepeda et al. 2003;Lee 2013): the reduction of problem dimensionality is the most important advantage of BEM, it means that only boundary line and surface domain are discretized for twodimensional plane and 3D program, respectively.
Compared to volumetric analysis, a boundary analysis results in a significant reduction in data preparation and a much smaller system of equations to be solved. This method which is used in the present study will be discussed in detail in the following sections. iv) Spectral element method (SEM) is the combination of the flexibility of the finite element method with the precision of the pseudo-spectral method.
The important feature of this methodology is the diagonal mass matrix by construction, which considerably simplifies the performance and decreases the computational cost because one can use an obvious-time integration scheme without having to invert a linear system (Komatitsch and Vilotte 1998;Chaljub et al. 2007;van der Meijde et al. 2020). v) The discrete element method (DEM) originally developed by Cundall and Strack (1979) for the analysis of rock mechanic problems is significantly popular. Its use is widespread across multiple disciplines (e.g., Richards et al. (2004), Fleissner et al. (2007), Ketterhagen et al. (2009), Horabik andMolenda (2016)). There is a large body of existing research using DEM, as well as a number of detailed publications on the use of DEM (e.g., Barreto and O'Sullivan (2012), Havenith et al. (2013), Thornton (2015)). The usefulness of DEM is recognized; however, it is not easy for the novice user to understand in detail its capabilities, limitations, and potential pitfalls. As mentioned in the previous section, in this research, two programs HYBRID and UDEC were used for numerical modeling taking into account the material properties presented in Table 1.
In 2014, the International Institute of Earthquake Engineering and Seismology of Iran has conducted studies on risk estimation and seismic geotechnical zoning of Arak City. Prepared reports include the following: single-station and array Microtremor survey, downhole, seismic reflection, and geoelectric. All the geotechnical and geophysical parameters needed to prepare numerical models in the current research are extracted from these reports. The theoretical Ricker wavelet is used as an incident wave for numerical modelling, i.e., double frequency Ricker wavelet with predominant frequencies of 2 and 5 Hz for the UDEC simulations and Ricker wave with predominant frequencies of 1, 3, and 5 Hz for the HYBRID simulation (Fig. 3). Both of them have the same frequency content, in the HYBRID, the defined mechanism for applying Ricker is considered single-frequency by default. In order to obtain a combined range of desired frequencies, three single-frequency waves (1, 3, and 5 Hz) have been applied and then the overlapping results have been given.
All models are subjected to the vertical incident wave. This wave has a vertical propagation direction and a horizontal displacement. The vertical component received on the surface is caused by the interaction of waves in the topographic features, which have very small values, for this reason, only the responses of the horizontal component are mentioned in this research.
In the boundary element method (HYBRID), there is no need to use absorbing boundaries, but in order to reduce the impact of disturbing waves caused by the edges of the model, the lateral boundaries of the model are extended at least 10,000 m from each side. In UDEC modeling, the width of the studied area is considered to be 2600 m, and dynamic boundaries are used as absorbing boundaries on the edges and bottom of the model, which absorb energy and prevent the reflection of waves within the model environment.

Parametric study
In the first part of this research, the seismic behavior of surface irregularities was studied within the framework of a parametric study, in which 5 homogeneous semi-sinusoidal hills. To describe the surface irregularity, symmetrical hills with height = 500 m and width = 1000 m are considered. Figure 4 shows the schematic configuration subjected to the incident wave, where b denotes the half-width and h the height of the hill, respectively. Dimensionless frequencies define Ω = ꞷb/πc 2 , where ꞷ presents the angular frequency of the wave and b and c 2 denote the half-width and shear wave velocity of the hill, respectively.
The geometries used for parametric studies are semi-sin hills, which are the most popular surface irregularities in nature. The geometry of 2D semisinusoidal hills is defined as follows (Eq. 1): where b and h denote respectively the half-width and the height of the hill.
The results are presented as time and frequency domain plots in which amplitude and spectral

Fig. 3 A Double frequency
Ricker incident wave with predominant frequencies 2 and 5 Hz using for UDEC simulation. B Ricker incident wave with predominant frequencies 1, 3, and 5 Hz using for HYBRID simulation amplifications are plotted, as a function of dimensionless time (tc 2 /2b) and frequency (ꞷb/πc 2 ), respectively.

Applied numerical study
In order to verify the results obtained from the parametric study, the effect of adjacency on three adjacent hills located south of the city of Arak is studied, using both numerical modeling and field measurements. The numerical modeling of the Arak hills is carried out with the HYBRID and UDEC programs. In order to study the effect of adjacency, it was necessary to study each hill individually and also in the presence of other hills. For this purpose, 4 two-dimensional profiles, which include three profiles with one hill (model H1, model H2, and model H3) and one profile with all three hills, model H1-2-3, are used for modeling. The same cross sections are used for the both modeling techniques (Fig. 5). To prepare the cross sections using the digital elevation model with 30-m resolution (DEM, 30 m), first, a cross section including the three studied hills was extracted, then one or two of the hills were removed from the model and the remaining model has been studied. Using this method, the effect of neighbor hills on the seismic response can be studied. DEM image is downloaded from https:// earth explo rer. usgs. gov/.
Since the HYBRID code is developed with the boundary element method, it does not need internal meshing, while internal meshing (triangular meshing) is used in the preparation of UDEC models (Fig. 6). To reduce the time of analysis, the lower part of the model is meshed using a structured mesh, and the upper part, due to the complex surface topographic is filled with an unstructured mesh.

Microtremor survey
Ambient vibrations (ambient noise) at the surface of the earth are produced by various vibrational sources (road traffic, industrial works, ocean waves, wind, etc.). These vibrations are most often surface waves (Rayleigh waves or Love waves) which propagate on the surface. Low-frequency waves (below 1 Hz) are generally called microseisms, and highfrequency waves (above 1 Hz) are called microtremors (Negulescu et al. 2013). The horizontal-tovertical spectral ratio technique which is proposed by Nakamura for site effect research overcomes the difficulties of traditional use of a reference site. This technique provides new ideas and methods for many aspects of research. However, the original goal was  (H1, H2, and H3). The model shows the mesh generated for the model, as well as the location of crest stations S3-P1, S7-P1, and S11-P1 to quantify site effects based on the Fourier amplitude spectrum ratio of the horizontal and vertical components of the ground pulsation at the same surface measurement point (Nakamura 1989).
The use of HVSR analysis of ambient noises has been widely considered to study topographic effects in recent decades (De Guevara et al. 2022;Köhler et al. 2012;Maghami et al. 2021;Panzera et al. 2011;Theodulidis and Bard 1995;Ulysse et al. 2021). Figure 6 shows the location of H1, H2, and H3, three adjacent hills located south of the city, as well as a 2D cross-section, along the hills, field measurements were made on this profile. This study presents the seismic data obtained from the microtremor measurements by the Guralp CMG-6TD broadband weak motion seismometer (velocity meter) with an instrumental resonance period of 30 s to evaluate the seismic responses of the hills (Fig. 7). Ambient noise recording and Fig. 9 The effects of the number of hills on the spectral amplification (A) and wave propagation (B) of five hills, H1 to H5, subjected to the vertical incident wave corresponding HVSR analysis were carried out in accordance with SESAME guidelines (Bard et al. 2004).
Among the numerous single-station microtremor recordings in the area, the results for three stations along a cross section profile have been studied here. Three ambient noise components were recorded for at least 25 min at each station. Since the study area is not located near industrial areas, microtremor measurement was done during the day. Microtremor data was processed according to the standard HVSR procedure by Geopsy 3.3.6 (https:// www. geopsy. org/). The HVSR curve is calculated for each recording, using 74 windows with a 5% overlap and squared average for horizontal components. The Konno-Ohmachi type of smoothing with a smoothing constant of 1 Hz is considered.

Parametric study
To study the seismic response of a single homogeneous hill with the shape ratio of 1 versus in the vicinity of other hills, five models subjected to the incident SV wave are investigated by HYBRID (BEM) in the linear domain. In the case of Sv, the vertical component is zero and only the horizontal component is designed. At the beginning, the behavior of a single hill (M1) was analyzed and then similar hills were added from one side (M2-M5) to study the effect of the adjacency on the seismic behavior. Figures 8A and B depict the horizontal and vertical components of hill seismic response, respectively, when subjected to the Sv incident wave. The maximum amplification of the horizontal component occurs at the crest (indicated by the red dashed line) and gradually decreases towards the base of the hills. Conversely, the amplification of the vertical component initially increases from the crest to the base of the hills, but subsequently decreases.
A comparison of the vertical and horizontal components at stations S13 and S9 in M1 and M5 models (Figs. 8C and D) reveals that the vertical component at crest stations is minimal (i.e., close to zero) compared to the horizontal component. As previously mentioned, the Sv incident wave lacks a vertical component, and the small values observed at the surface are due to the interaction of waves with topographic features. Furthermore, the values of the vertical component are unaffected by adjacency, and even at station S9, which exhibits the highest horizontal component spectral amplification, the vertical component remains near zero.
Consequently, this paper only considers the seismic response of the horizontal component and ignores the vertical component.
The modeling results are analyzed to study the effect of adjacency on seismic behavior in the time and frequency domains (Figs. 9 and 10). Figure 9A illustrates the spectral amplification of five models in the dimensionless frequency range of 0.2 to 3. Stations S13, S11, and S9 are indicated by a dash line and M5 shows the location of each hill. The color bar represents the spectral amplification range.
The wave propagation pattern for each model in the dimensionless time range of 0 to 15 is shown in Fig. 9B. The amplitude range is displayed by the color bar.
A comparison of the seismic response of the hills suggests that, when the number of hills is odd, the middle hill displays higher spectral amplification than the side hills. On the other hand, when the number of  2 and 3). D The spectral amplification of three central stations S13, S11, and S9 ◂ Table 2 The maximum spectral amplification of three stations S9, S11, and S13 for models M1-M5 with different number of hills hills is even (especially in M4), due to the absence of a specific hill as the focal point in the center of the model, two middle hills show higher spectral amplification than the side hills due to wave scattering, as shown in Fig. 13 in the Appendix. It means that, after reflection and scattering of the waves from the hills, the waves concentrate in both middle hills. This can be clearly seen at dimensionless time of 5 in M4 (Fig. 13B in the Appendix). In this case, the waves are also scattered and propagate towards the lateral hills (H5 and H2), which makes the maximum spectral amplification of the central station of H5 (S13) in M4 greater than M3 (Table 2). Table 2 presents the maximum spectral amplification of the crest stations, S13, S11, and S9 located on H5, H4, and H3 respectively. The stations S11 and S9 show higher amplifications rather than nearby hills when located in the center of the model. For example, in M5 model, the amplification of S9 is equal to 4.3, while this value is 3.6 for S11. Moreover, the amplification of S11 is equal to 3.4 in the M3 model while for S13 this is 2.6.
Regarding S13, while the spectral amplification of the first three models decreases from 3.3 to 2.6 from single-hill models to three-hill models, the M4 shows an increase up to 3.9, followed by another decrease in the M5 (equal to 3.1). As mentioned above, the main reason for this increase is the distribution of scattered waves throughout the model, as there is no specific hill as a focal point in the center of the model to converge the waves.
As it can be deduced from Table 2, adjacency has no major effect on the main peak frequency of the S9, S11, and S13 and at all stations, a main peak in dimensionless frequency range of 0.4-0.7 can be followed. Figure 10 examines the seismic response of models M1, M3, and M5, which have an odd number of hills, in more detail. The pattern of wave propagation in all these three models is shown in Fig. 10A and in the dimensionless time range, from 0 to 15.
Stations S13, S11, and S9, which are located in the center of each model, are also marked with a dashed line, and the color bar indicates the wave amplitude. Then Fig. 10B shows the duration of each station separately.
The duration of vibration at each station shows that with the increase in the number of hills, in addition to the increase in fluctuations, new peaks are observed in the signals. It can be seen in the wave propagation pattern that when the model is subjected to the incident wave, the wave refracts, scatters, and converges towards the center of the model. Since the number of peaks is odd, the wave in each model (S9, S11, and S13) converge towards the central station, and because of this convergence, new peaks appear in the time domain results.
Another question is that, does the scattering of waves from the adjacent hills affect the peak frequency of the station under study or not? For this purpose, the frequency domain responses are examined in Figs. 10C and D. First, the overall model response in the dimensionless frequency range, from 0.2 to 3, is shown. In these diagrams, the colored bar is dedicated to spectral amplification. The increase in the spectral amplification of the central station is initially slow (M3 compared to M1), but as the number of hills increases from 3 to 5, the spectral amplification of the central station gets 1.2 times higher (S9 in model M5 compared to S11 in model M3). According to Fig. 14 in the Appendix, adjacency has no effect on the peak frequency of any station. Only S13 shows a new peak at the dimensionless frequency of 0.55 due to the adjacency, while the peak frequency for the first peak is still constant at dimensionless frequency of 0.4 for all three models. Therefore, it can be suggested that the changes in frequency peak shown in Fig. 10D are mostly related to the location of the stations and not their adjacency.
Due to the simplification of the models, most researches only studied a specific feature, and few studies have paid attention to the adjacency. In between, Kumar et al. (2018) and Afzalirad et al. (2019) have investigated the effect of multi-ridge mountain ranges on wave amplification. It can be concluded from the results presented by Afzalirad et al. (2019), that the crest of homogeneous adjacent hills with a similar shape ratio, presents larger maximum amplifications in comparison with the crest of single hills. Although this difference is increased by increasing of the shape ratio, it is negligible within the studied shape ratio. Kumar et al. (2018) suggested that neighboring topographic features affect the response of the subridge/valley due to the occurrence of interference of the reflected and diffracted waves from the adjacent topographic features with the wave field within the considered ridge/valley. Therefore, the results of the current research show a good agreement with the results of the previous research.
Due to the increasing development of numerical methods in the field of numerical modeling, the study of the effect of adjacency on seismic response should be developed using other numerical methods. In addition, since 2D studies have limitations in 3D environments, it is recommended to study the adjacency effects using 3D methods.

Applied numerical study
In the following, the seismic response of three adjacent hills located in the urban area of Arak has been investigated for a wider investigation of the adjacency effect on the seismic response. Figure 11 shows a 2D profile (P1) on three hills and the location of the crest stations S3, S7, and S11 on H1, H2, and H3, respectively. The purpose of this section is to evaluate the effect of the adjacent hills and the number of hills on the seismic responses of the mentioned stations. For this purpose, each station is studied in single mode (H1, H2, and H3) and in the vicinity of other hills (H1-2-3). In Fig. 11, the numerical results of UDEC and HYBRID in the frequency domain and in the frequency range of 0.5 to 10 Hz show that, in general, at each station, the H1-2-3 model shows higher spectral amplification compared to the response of the single hill models (H1, H2, and H3). Increasing the number of hills causes an increase in the amplification of all stations, while this amplification was not the same in all stations and the middle hill shows a greater increase than the side hills. Since the codes are produced by different methods, it can be expected that the spectral amplification obtained from HYBRID and UDEC do not perfectly match. For example, the UDEC shows the spectral amplification higher than HYBRID (Table 3) and also the content of frequency peak in UDEC shows almost 0.5 Hz more than HYBRID, but both HYBRID and UDEC graphs in the frequency domain present the same trend in the range of 0.5 to 10 Hz.
As mentioned in the parametric results section, when three hills are placed together, the crest station of the middle hill shows the maximum spectral amplification (Table 2), now the Arak hills follow this rule and Station S7-P1 (located in the middle hill) shows the highest amplification, which means that the amplification at this station increases from 1.6 to 2 (UDEC) and 1.49 to 1.8 (HYBRID) for H2 and H1-2-3 model respectively (Table 3).
In view of the spectral amplification diagrams presented in Fig. 11, in this study, the peak frequency is constant with the increase in the number of hills from single mode (H1, H2, and H3) to H1-2-3. This means that the adjacency effect only increases the spectral amplification at the peak of the graph and has no effect on the frequency. Figure 12 shows the amplification obtained from microtremor surveys using the HVSR method at stations S3, S7, and S11 to compare with numerical simulation by HYBRID and UDEC. Recording at the stations was carried out using Guralp CMG-6TD broadband weak motion seismometer. Since it is not possible to have single station data for individual hills, experimental studies have been compared with numerical studies related to multiple hills. As expected, similar to what was obtained in parametric studies (Table 2), station S7, which is located on top of the middle hill, shows higher amplification values (equal to 2.1) rather than stations S3 and S11 (equal to 1.8) that is shown in Table 3.
In general, the comparison of numerical analysis and microtremor survey in Fig. 12 shows that the results of numerical modeling are in good agreement with the field observations and in all three methods, the central station (crest station) in the middle hill has a higher spectral amplification than it shows two other stations, but it does not mean three graphs are completely matched. The microtremor survey in the Arak hills is influenced by various factors such as the impact of the sedimentary basin of Arak City, nearby mountains, and subsurface conditions, which cannot be removed from the recorded signals. In contrast, the numerical modeling is done by two-dimensional methods and the environment is studied in a simplified form (completely homogeneous media without environmental effects). Since real conditions cannot be fully applied in numerical modeling, the difference between numerical and experimental results is undeniable. And, considering the similar trends of all three graphs in the studied frequency range, it is worth saying that all three graphs confirm the results of the parametric studies presented in the previous section.
The maximum spectral amplifications presented in the previous sections are fully presented in Table 4. The main purpose of this study, which is to investigate the effect of adjacency on seismic response, is well shown in this table. The results in each section show that despite Fig. 11 Schematics of two-dimensional cross-sections including three hills together and single hills, and comparison of the amplification of top hill stations S3, S7, and S11 located on profile 1 (P1) for both single and adjacent hills using two programs, UDEC and HYBRID  . 12 HVSR from Guralp recording of top station located on Profile 1. S3-P1 located on H1, S7-P1 located on H2, and S11-P1 located on H3, and compared with UDEC and HYBRID numerical modeling the difference in the value of spectral amplifications, the middle hill always shows a higher amplification.

Conclusion
In this research, the effect of adjacency on the seismic response of homogeneous hills has been investigated using numerical modeling and microtremor measurements. Numerical modeling has been done using HYBRID and UDEC programs in a homogeneous domain and applying mathematical Ricker wave with different predominant frequencies.
The parametric study of the effect of adjacency using five models, M1 to M5, showed that when more than one hill is subjected to the seismic wave, because of the scattering, the waves converge toward the center of the model. When the number of hills is odd, significant amplification is seen in the center of the model (middle hill). For example, a spectral amplification of 3.4 and 4.3 is observed for the stations S11 and S9, respectively, which are at the center of the models M3 and M5. However, when the number of hills is even, due to the lack of a specific hill as a focal point in the center of the model, the two middle hills show higher spectral amplification than the side hills.
In addition, the results clearly show that as the number of hills increases, the magnification in the middle hill also increases (comparison of the central station in M3 and M5 models).
Numerical and experimental studies on Arak hills also show that the different numerical methods and microtremor measurements provide comparable results in a certain frequency range despite the differences in the value of spectral amplification and the frequency of the main peak, which confirm the results of parametric studies. The results obtained from both the numerical modeling and field surveys further show larger spectral amplifications for the middle hill than for the side hills.
Finally, all numerical modeling results related to the parametric and to Arak's hill study show that with an increasing number of hills, the spectral magnification may also slightly increase.
Author contribution This research was carried out as a part of Niloufer Babaadam's doctoral thesis as the main researcher. Other authors have been in charge of guiding different parts of this research as supervisors and advisors.

Data availability
The data used in this research will be available upon the request of the readers and after correspondence with the authors.

Declarations
Ethics approval and consent to participate In conducting this research, the authors consider themselves committed to complying with all the principles of scientific ethics and declare their full satisfaction with the publication of this scientific research. All the publishing rights of this research are reserved for the Journal of Seismology.