Seismic site-effects assessment in a fluvial sedimentary environment: case of Oued Martil floodplain, Northern Morocco

In Northern Morocco, seismic site-effects in general and liquefaction hazard in particular can occur in the event of a major earthquake due to the thick sedimentary cover characterizing the peripheral Neogene basins of the Alboran Sea. An example is Martil plain, which was the subject of important economic development during the last two decades. In this regard, we present in this study an assessment of seismic site-effect hazard through the horizontal-to-vertical spectral ratio (HVSR) method and the vulnerability index (Kg). Multiple Analysis of Surface Waves (MASW) method and core-drilling data are also used to complete our analysis and interpret the spatial distribution of Kg maps. Our findings suggest more vulnerability to liquefaction in the Southern segment of the basin, which can be explained by the asymmetrical geometry of Quaternary sedimentation, due to tectonic uplift that influences also the surface and subsurface hydrology processes. However, in this study, the main difficulties and the challenges that we faced were related to data acquisition during windy days and in the areas near industrial zones or near high traffic roads, which have significant effects on the accuracy of the HVSR analysis. In this study, HVSR microtremor measurements are recorded at 197 stations in the Martil floodplain to generate a seismic vulnerability index (Kg) map. According to the analysis results, the predominant frequency (F0) values range from about 0.2 to 10.6 Hz and the peak amplitudes (A0) values are in the range of 0.39–10.4. As a result, some districts, especially those classified as economically disadvantaged, are found to be the most exposed to this hazard (with a Kg > 10), which must be taken into consideration in future risk reduction and mitigation plans. We conclude the existence of significant seismic effects potential despite the moderate seismicity of the area. Therefore, our research needs to be completed by scenario-based seismic hazard modeling to investigate the capacity of seismic events at the region to produce the above-suggested amplifications.


Introduction
Over the previous decades, many earthquakes have manifested worldwide, leaving behind a large number of injuries, deaths, and significant material damage (Jaiswal and Wald 2010;Daniell et al. 2012). Consequently, the scientific community largely investigated and recorded seismic events since the nineteenth century thanks to the professionalization of science (Agnew et al. 2002). However, the instrumental recording of seismic events in Morocco only dates back to the early 1960's following the devastating event of Agadir 1960, which left behind an excess of 12,000 casualties and completely destroyed 75% of the city (Cherkaoui et al. 1989). Previous earthquake occurrences are the subject of historic and literary descriptions suggesting that Morocco, just like other countries of the West Mediterranean region, was exposed to and affected by historic earthquakes such as the infamous 1755 occurrence (El Mrabet 2005;Kaabouben et al. 2009). In Northern Morocco, recent seismic activity is mainly concentrated in the vicinity of Al Hoceima city, which has experienced in the space of 10 years two violent earthquakes, one in 1994 and the other in 2004. The latter (Mw = 6.3) was more pronounced and caused the death of 629; injured 926 and left behind 15,230 homeless and 2539 destroyed buildings (El Alami et al. 1998;Jabour et al. 2004;Stich et al. 2005;Akoglu et al. 2006). A more recent event took place in 2016 with less severe consequences (Kariche et al. 2018).
According to the historic seismic catalog provided by the Spanish National Geographic Institute (IGN), the Tangier peninsula is characterized by less frequent seismic activity with moderate intensities and less expressed consequences. In fact, most of the seismic hypocenters are concentrated in the Alboran Sea domain (Comas et al. 1999;Vernant et al. 2010;DeMets et al. 2015;Grevemeyer et al. 2015) or in the Gulf of Cadiz region Gutscher 2004;Gutscher et al. 2006). This means that significant earthquake events are more likely to be produced 100's of kilometers away from the major cities of the said peninsula (35° 15′ 46.642″ N 5° 33′ 42.221″ W). Despite this, the unconsolidated Neogene sediments in the peripheral basins of the Alboran Sea (Guerra-Merchán et al. 2014) could induce significant site-effects even in the case of such far-away occurrences.
In fact, the seismic risk for a region depends on the probability of occurrence of an earthquake causing a certain level of ground shaking (seismic hazard) and on the structural vulnerability of the built structures. It varies widely depending on population density and on the economic potential of the investigated area. Judging from the ongoing economic development in the Northern Rif in general and Tetouan city in particular, many issues at stake could be identified which include but are not limited to important infrastructures and dam projects.
Despite this, no extensive scientific research was conducted in the area in order to assess the seismic behavior of the above-described peripheral basins of the Alboran Sea where most anthropogenic activities are concentrated. In effect, seismic hazard microzonation studies use deterministic techniques such as the horizontal-to-vertical (H/V) ratio method (Kanai and Tanaka 1961;Nakamura 1989Nakamura , 2008Bonnefoy-Claudet 2004) to obtain such results. In Southern Spain, the large alluvial basins of the Andalusia region were intensively researched with results suggesting strong local site-effects manifested by larger predominant seismic periods for the unconsolidated alluvial and floodplain sediments (> 0.5 s), thus suggesting a greater risk of liquefaction (Navarro et al. 2002(Navarro et al. , 2007(Navarro et al. , 2013García-Jerez et al. 2008;Gaspar-Escribano et al. 2010). Nevertheless, the southern part of the Alboran Sea is poorly characterized in terms of local site-effect assessment and seismic microzoning. A recent publication by Arab et al. (2021) attempts to address this problem but fails to define the nature and spatial extension of the investigated material. Therefore, their findings are deemed less useful for understanding the local site-effects in their regional context. Also, the intra-basin frequency variance is not taken into consideration, which largely reduces the significance of the results. A second attempt was carried out in Nador city (Chaaraoui et al. 2021), the results of which cannot be extrapolated on the Tangier peninsula region because of their contrasted geological settings. To our knowledge, no other research provides further information to fill this knowledge gap between the Northern and Southern regions of the Betico-Rifain arc.
In this regard, we attempt to characterize seismic site-effects at the Martil floodplain. This large basin with thick sedimentary fill deposits constitutes a good example for investigating seismic site-effects in general and liquefaction hazard in particular. In addition, the increasing urbanization and population growth at the area with little to no regard to the possibility of a future earthquake may cause large material and human losses due to liquefaction in this water-rich environment. In general, fundamental frequencies obtained by microtremor measurements cannot be related to the total thickness of Quaternary sediments, but are highly influenced by layers of conglomerate or bedrock (Gosar 2010). This represents a challenge for a site-effects study, because no simple relation between sediments thickness and the related seismic response can be expected. In the other hand, the main difficulties arise from traffic and industrial noise. In addition, in the populated areas, the free-field space between houses was also very limited. We avoided taking measurements on windy days, because the noise introduced by wind can severely affect the reliability of HVSR analysis.
Subsequently, an assessment of soil liquefaction potential in this earthquake-prone region is a crucial step both for understanding such phenomena in their regional context as well as the identificating hazardous zones. In the present research, we measure the site response at different locations of the Martil plain through the HVSR technique. The MASW method along with geotechnical data from boreholes is also used to complete the analysis. As a result, some districts are found to be the most exposed to the seismic hazard, explained by the existence of significant seismic effects potential despite the moderate seismicity of the area, which must be taken into consideration in future risk reduction and mitigation plans.

Geological and tectonic setting
The Northern provinces of Morocco are among the most active in terms of seismic and tectonic activity, with moderate rates of deformation estimated around 2 to 4 mm/yr (Vernant et al. 2010;Chalouan et al. 2014;Poujol et al. 2014;Agharroud et al. 2021). They are part of the collision zone between the Eurasian and African plates where compressive and extensional deformation exist simultaneously (Platt and Vissers 1989;Gutscher et al. 2012;Fernàndez et al. 2019). In the Rif (Northern Morocco), based on stratigraphic and tectonic criteria, three paleogeographic domains are distinguished: Internal Domain, the Flysch Domain, and the External Domain (Fallot et al. 1937;Durand-Delga 1960;Didon et al. 1973). Its structure is in the form of tectonic nappes with west verging subhorizontal thrust faults, the distribution of which follows the arched morphology of the chain.
The latest tectonic phase that sculpted the morphology of the northern Rif is dated to the upper Burdigalian since no later sedimentation was reported in the area during the rest of the Miocene age (Feinberg et al. 1990). Neogene sedimentation resumed in the Pliocene era with shallow marine and then deep marine sedimentation in the peripheral paleoria basins of the Alboran Sea (Guerrera et al. 1993;El Ouahabi and Fagel 2009;El Kadiri et al. 2010b;Hlila et al. 2014;Merzeraud et al. 2019). The Quaternary cover of these basins is exclusively continental with thick alluvial deposits (Fig. 1).
The plain of Martil, an example of such peripheral basins and the focus of our study is located in the internal domain of the Rif chain, which is characterized by its tectonic and stratigraphic complexities. This domain mainly formed by Paleozoic material of the Sebtides and Ghomarides units (Durand Delga and Kornprobst 1963;Didon et al. 1973) was subject to extensional deformation due to the triggering of pull-apart systems in the late Miocene by a N-S oriented compressive phase (Benmakhlouf and Chalouan 1994;Romagny 2014). Since then, the Martil plain was the subject of thick sedimentation that constitutes the Quaternary cover of the Paleozoic substratum. Seismic profiles and field data suggest that the sedimentary deposits cover a graben and horst structure that is related to two extensional deformation phases, one in the late Oligocene (Ouazani Touhami et al. 1994) and the other in the late Miocene (Romagny 2014). In fact, the latter phase is related to a N-S compressive phase that triggered the collapse of the Oued Martil basin by a pullapart system.
The lithology of the said cover is essentially formed by fluvial detrital deposits in the form of unconsolidated gravel and pebbles constituting the basis of the sedimentary Fig. 1 Geological map of the study area sequence, clayey deposits forming its stratigraphic top, and travertine and tuffaceous layers deposited near carbonate rocks of the Dorsale calcaire (Benmakhlouf 1990;Benmakhlouf and Chalouan 1994).
In topographic terms, the Oued Martil basin is essentially a flat terrain covering an area of about 87 km 2 with little to no topographic variations (Himi et al. 2017) (Fig. 1). It is bounded to the North and South by two horsts manifested by cape-like landforms denominated, respectively, Cabo Negro and Ras Mazari (Romagny 2014). In Koudia Es-skirej hill (Fig. 1), the substratum outcrops near Oued Allila stream with low to moderate slopes and little topographic manifestations. Uphill, the basin is thinner with Paleozoic formations delimiting its geometry (Fig. 1). Near Tetouan city no more Paleozoic material is observed in the cluse of Tetouan, with Flysch nappes and Dorsale calcaire units outcropping in the peripheries of the basin (Benmakhlouf 1990).

Socio-economic setting
The study area is part of the North Moroccan coastline, which was the subject of significant socio-economic development since the late 1990's. According to data provided by the High Commission for Planning through the general census of population and housing of 2014 (RGPH 2014), Martil city houses currently 63 853 inhabitants with a projected population of 119 128 in 2030. This rapid increase in population will result in a higher demographic density especially knowing that medium-rise multi-apartment buildings are replacing the individual housing units to respond to the increasing demand on real estate in the area (RGPH 2014). Consequently, the management of natural risks becomes more difficult with more issues at stake being concentrated in increasingly smaller spaces.
The economy of the city depends on the tourism sector with very small industrial units located 2 km away from the city. Therefore, the population of the city varies depending on the seasons, with very high densities during the summer (i.e., the active season) and relatively low densities during winter.

The horizontal-to-vertical spectral ratio (HVSR) technique
HVSR technique was invented in the seventies by Nogoshi (1971) based on the initial studies of Kanai and Tanaka (1961) and was later developed and popularized by Nakamura (1989). These authors demonstrated that the ratio between horizontal and vertical ambient noise records depends on the fundamental frequency of the soil strata and its amplification factor.
In fact, this approach compares the ratio of the H/V spectrums at a given site to that of a reference site, usually a rocky substratum outcrop. In this study, we chose the HVSR technique because of its ease of use and low cost as well as its ability to directly estimate the frequency of sedimentary strata without having to know the geological structure and the shear wave velocity of each subsoil level (Herak 2008;Triwulan et al. 2010;El-Hady et al. 2012). Accordingly, ambient noise measurements were carried out in 2018 in 197 sites of Oued Martil basin ( Fig. 1), using a portable Tromino seismograph (Moho 2018), with three orthogonal velocimetric sensors (two horizontal components and one vertical) and a laptop computer to digitize and store the data. A few reference measurements were also 1 3 performed in the phyllite and flysch substratums. Each acquisition station was equipped with an internal GPS which allowed us to know the precise geographical position of the recording (in UTM coordinates). The data processing was realized using Grilla software in accordance with international guidelines of Sesame European Research Project (SESAME 2004).
The Tromino devices use the Fourier spectra calculated for selected windows using the fast Fourier Transform (FFT) algorithm (Nakamura 1989). The H/V spectral ratio of the two horizontal Fourier and one vertical Fourier spectra can be computed using Eq. (1).
where S(ω)NS, S(ω)EW and S(ω)V are the ratio of the horizontal and vertical spectra of the ambient seismic noise, respectively, and ω is the angular frequency.

Multi-channel analysis of surface waves (MASW) technique
The MASW method was first initiated in geophysics by Park et al. (1999Park et al. ( , 2005 and is widely used in recent years to measure the shear-wave velocity (Vs) of geological materials in shallow subsurface studies. This method is not destructive and was shown to be more cost-efficient than destructive techniques (seismic refraction, for example), which makes it suitable for a wide range of applications. It provides a good characterization of the subsoil and a very useful method for the assessment of site-effects in a given area (Mohamed et al. 2016;Roy & Jakka 2017;Pamuk, et al. 2018).
The material used in this study consists of an impulse source (hammer and metal plate), a seismic acquisition unit (Tromino), and a series of geophones spaced 3 m from each other deployed along a linear profile ( Fig. 1). At each station, series of shots are generated at the edges of the device, which allows us to obtain seismic profiles (Fig. 2). From the seismic section ( Fig. 2), the arrival times of the surface waves are determined to obtain a dispersion diagram from which we estimate the dispersion curve of the analyzed waves. Muting was also performed to filter out noise and to enrich the surface wave signal. Once all the dispersion curves were generated by analyzing the Rayleigh-type surface waves on a multichannel record, they were further inverted to generate a 1-D Vs (average shear wave velocity) profile using the inversion algorithm given by Xia and Park (1999). These 1-D Vs profiles were interpolated to generate a 2-D Vs profile of each site (see example in Fig. 2). The latter data processing steps were performed using Grilla (SESAME 2004).
(1) . 2 Interpreting a MASW profile to generate the S-wave velocity profile obtained by inversion

Topographic data
Topography is essential to delimit and understand the geometry and spatial extent of sedimentary basins. Although Shuttle Radar Topographic Mission (SRTM) data provide reasonably accurate Digital Elevation Models (DEM), they tend to perform poorly in flat terrains because of the low spatial resolution of the data, which does not reveal landforms in gentle morphologies. Hence, the topography used in this study is derived from digitizing the elevation contours of the 1:25,000 topographic maps of the area, which were analyzed using a geographic information system (GIS) platform to generate a medium resolution (5 m) DEM. These data were exploited to construct the 2-D topographic profiles presented in this study as well as the calculation of the depth of the stratigraphic boundary separating the soil from its Paleozoic substratum.

Boreholes
In this study, 34 boreholes were exploited which cover different parts of the study area in order to reconstruct the geometry and geology of the Oued Martil basin. The same database was partially or totally used for other purposes such as groundwater vulnerability assessment in coastal aquifers of the Martil plain (Himi et al. 2017;Benabdelouahab et al. 2019) and geotechnical characterizations (Ahniche 1997). These boreholes were drilled by the Loukkos Hydraulic Basin Agency (ABHL) between the years 1980 and 2005 and were stopped in either Quaternary or Pliocene strata without reaching the Paleozoic substratum, with depths varying between 30 and 80 m. An exception is the borehole that reached the Paleozoic substratum with a depth exceeding 100, which can be found in the North of the study area (X: 510,400, Y: 558,940). In order to make use of information provided by these core drilling, microtremor measurement is realized near borehole sites so as to analyze the relationship between sediment thickness and fundamental frequency in the study area.

Methodology
Since Nakamura (1989, several studies (Alfaro et al. 2002;Garcia et al. 2006;El-Hady et al. 2012;Navarro et al. 2013;López Casado et al. 2018) have used the microtremor measurements for site-effect evaluation as well as seismic micro-zoning studies. Based on our literature review, the HVSR technique, which is developed by Nakamura (1989Nakamura ( , 1997Nakamura ( , 2008, is one of the most popular methodologies that is frequently used worldwide to identify and delimit areas of high liquefaction hazard.
The Nakamura technique has been adopted to compute the fundamental frequency response with respect to the peak amplification with the help of source tool Grilla software for different microtremor data recorded in the Martil floodplain. Thus, the obtained results from the HVSR analysis have been checked for the reliability and clarity of the peak by using the SESAME guidelines condition. Further, the Tromino devices use the Fourier spectra (Eq. 1) calculated for selected windows using the FFT algorithm (Nakamura 1989; Pandey et al. 2018Pandey et al. , 2022; with respect to azimuth have been performed for all recorded microtremor data in the study region. The range of rotation is kept from 0 to 180 degrees for the horizontal components of motion and also the azimuth direction is selected clockwise in the north direction. Dynamic soil response during earthquake ground motion is the object of the present work, which is assessed through HVSR method (Fig. 3). Consequently, the liquefaction phenomena are a process in which soil loses (Meshram et al. 2022), its strength and stiffness due to ground shaking during an earthquake. This can result in the ground behaving like a liquid, causing structures and buildings to sink, tilt or even collapse. In our case study, microtremor measurements were executed at 197 stations in the Martil plain. For all stations, microtremor was measured and H/V was calculated according to the procedure used in the Grilla software prepared by the group SESAME (2004). This software generates H/V curves based on the Nakamura (1989) approach. From the H/V interpretation, the predominant frequency (F0) and the amplification factor (A0) are easily obtained of each site (Fig. 3). These parameters were then used to compute the vulnerability index (Kg) (Eq. 2), which quantitatively assesses liquefaction hazard intensity (Nakamura 1996;Beroya et al. 2009). Subsequently, areas with Kg > 10 are considered more vulnerable to liquefaction. Although Kg threshold values are arbitrary, a value of 10 is deemed reasonable by many authors (e.g., Livao et al. 2017; Akkaya 2020). Huang and Tseng (2002) concluded that the higher vulnerability index is observed in the locations where the chance of liquefaction is certain.
where Kg is the vulnerability index, A 0 is the fundamental amplitude, and F 0 is the fundamental frequency.
The useful of Kg values is to detecting weak point and for determining the dangerous damage possibility could occur before an destructive earthquake in the areas under study (Beroya et al. 2009;Akkaya 2020).
(2) Fig. 3 Schematic model of the methodology used in this study

Basin geometry from core-drilling logs
The core-drilling logs were exploited in this study to build 2-D longitudinal and transverse geological cross sections of Martil basin (Fig. 4). Overall, the geometry of the basin is defined by thick quaternary alluvial deposits resting directly over Pliocene marly material. Across the study area, a said basal layer mainly formed by semi-consolidated gravel pebble and/or sand deposits separates these two layers (Stitou 1995;Himi et al. 2017;Benabdelouahab et al. 2018) (Fig. 4). The geological substratum in the form of Paleozoic phyllite and Greywacke material (Fig. 4D,E and C) is never reached by the core-drillings and is believed to exist in a much deeper position according to nearshore seismological profiles provided by ENIEPSA exploration company in 1981. The same can be said about the Flysch substratum in profiles A and B (Fig. 4) which is not revealed by any borehole at the area.
A general thickening of Quaternary material toward the shoreline marks the geometry of the basin ( Fig. 4D and E). In some areas, a neotectonic vertical motion believed to be caused by normal faulting translates into significant and sudden changes in the thickness of Quaternary layers (Fig. 4D). Such changes can have important implications with regard to seismic site-effects in Oued Martil floodplain.
The transverse cross sections reveal on the other hand the asymmetrical sedimentation processes that characterize the basin, with Quaternary material being thicker in one side of the basin or the other depending on the location (Fig. 4A and B).

Vs30 measurements across the study area
With regard to the shear wave velocity characteristic of each geological material present in the area, our results reveal relatively low Vs30 values (around 280 m/s) for the Quaternary material compared to those of the paleozoic substratum (in excess of 500 m/s) (Fig. 5).

Fig. 4
Geological cross sections of Martil basin Interpreted from core-drilling data. The location of each section is shown in Fig. 1 Also, the Vs30 variance for the Quaternary material is found to be small despite its abovedemonstrated lithological heterogeneity (Fig. 5). Based on this, one can conclude that such material can induce significant amplification effects due to its low shear wave velocity; a fact that was previously demonstrated in other areas of the world (Castellaro et al. 2008;Graizer 2009;Laurendeau et al. 2013).
Pliocene deposits are also marked by similar Vs30 values (around 360 m/s) with an even smaller variance (Fig. 5), which exposes their vulnerability to seismic site-effects given the right circumstances (Fig. 5).

HVSR tests
The microtremor measurements yield various H/V spectral ratio peaks that range from low to high values in Oued Martil basin. Typical examples of H/V curves are illustrated in Fig. 6.
In some H/V measurement sites, the analysis reveals no clear H/V peaks (flat curves) (station T3P9 in Fig. 6) which can be due to a multitude of factors such as the proximity of the seismic bedrock from the surface (Machane et al. 2014), or the fact that some lithological facies are known to produce such flat curves due to their low impedance contrast (Guillier et al. 2005;Vella et al. 2013). According to borehole (111-2) drilled near such material, the marly material outcrops near the surface with the absence of clear H/V peaks in its corresponding H/V curve. Such curves are indicative of little to no amplification of seismic waves (A0 < 2).
The rest of the H/V stations are marked by a clear F0 peak (e.g., "T6P25" in Fig. 6) suggesting a strong impedance contrast between the soft sedimentary cover and its consolidated stratigraphic substratum (SESAME 2004). In fact, a large peak with a high amplification is indicative of a strong impedance contrasts (SESAME 2004;Gosar 2010;Oubaiche et al. 2012) which is susceptible to amplify ground vibrations.
In a few particular cases, the HVSR graphs show a plateau or bulging (station 'T5P46' in Fig. 6), suggesting low to medium contrasts (Guillier et al. 2005). These plateau curves usually have an amplitude factor inferior to 3 and are probably due to the lithological heterogeneity of the Quaternary cover.
Some H/V curves do not show one but two peaks as illustrated in Fig. 6 (stations T3P21, T5P35), with an amplitude varying from 2 to 6, sometimes even surpassing 10. They are indicative of the presence of a consolidated layer in the stratigraphic sequence of the sedimentary cover, which translates into a second seismic interface (SESAME 2004;Guillier et al. 2005;Fnais et al. 2010;Gosar 2010). Such a layer can be attributed to the consolidated Villafranchien conglomerate layer described in Sect. 5.1. (Fig. 4). Among the more complex forms, a few measurements exhibit several peaks that are not well separated (very broad peak) (station T6P31 in Fig. 6), which reflect a complex structural setting as shown by Gosar (2010). For such cases (two clear peaks), we used the value of the peak with the highest amplitude.
Subsequently, all of the results shown above are exploited to prepare iso-frequency and amplitude maps of Martil plain.

Microzonation Maps of F0 and A0
Based on our literature review, the H/V method plays a significant role in seismic hazard evaluation studies, especially in urban areas. In our case studies, the interpolated F0 values and their corresponding amplification factor (A0) are shown to be ranging from 0.41 to 10 Hz and 0.45 to 10, respectively (Fig. 7).
The F0 map (Fig. 7a) also shows that high F0 values dominate the central part of the said basin, reaching their maximum near Koudiat Es-Skirej where the stratigraphic substratum outcrops. To the South, F0 values decrease gradually and come close to zero near the main water streams at the area (F0 < 2.80 Hz). These low values are present in the Mallalyienne basin as well along Oued Allila stream.
To interpret these results, one can point to the fact that resonance frequency becomes higher in areas where the depth of the sedimentary cover is greater (Parolai et al. 2001). According to Teves-Costa et al. (1996); Parolai et al. (2001); Gosar (2017) and Rezaei & Choobbasti (2017) this situation is characteristic of bedrock-soil environments because of the presence of bedrock strata far away from the surface.
In the iso-amplitude map (Fig. 7b), it can be seen that several areas with lower amplification values (0.39 and 2.62) were found in the center of the city (Martil basin) as well as in the northwestern periphery with amplitude factors around 3. Conversely, the second range indicates an important amplification effect (3 to 7), mostly concentrated in the southern part of the plain, and high amplitudes (> 7) found mostly in the peripheries of the basin. In addition, a strong amplification is observed near the coastline area at Koudiat ES-Skirej. Such results can simply be due to the influence of marine waves, which can induce errors that lead to high amplification readings (noisy measurements). Also, we suspect the presence of artifacts due to interpolation in areas with low H/V station density. Fig.7 a Iso-frequency map of sediments resonance frequency and b maps of H/V peak of the study area 1 3

Microzonation maps of the vulnerability index
From the F0 and A0 data, the ground vulnerability index (Kg) is calculated for all measurement stations conducted in the Martil floodplain (Fig. 8). In fact, this parameter (Kg) can be interpreted quantitatively, which allows estimating the soil weaknesses in seismic hazards studies, and subsequently forecast the effects of a potential destructive earthquake (Nakamura 1997).
In the study area, most of the computed Kg values are inferior to 5 (Fig. 8). These values are representative of the Mallalyienne sub-basin, with low to moderate vulnerability index values. In other areas, Kg values are ranging from 5 to 10 (Fig. 8). For terrains with high vulnerability index values (kg > 10), their spatial distribution follows that of Martil river. They are concentrated in the Southern segment of the basin, where the geology is dominated by asymmetrical Quaternary alluvial deposits (Fig. 8). Possible reasons explaining such a spatial distribution are discussed below.

Socio-economic implications
The optimal goal of seismological research is mitigating the impact of earthquakes on human societies. In this regard, this work constitutes a preliminary characterization of liquefaction vulnerability, which is not only important to geologists and geomorphologists alike, but can also save lives and reduce economic damage to the community. In addition, the Kg map produced in this study can be considered as an important indicator in any future risk management and mitigation plan especially in coastal areas that are densely populated such as Martil plain.
In effect, the densification of human settlements alongside coastal areas is a worldwide phenomenon (Small et al. 2000;Small and Nicholls 2003;Barragán and De Andrés 2015). In Morocco, the population growth in the coastal areas has been very significant. For instance, the coastal population of Northern Morocco grew from 9.4 million inhabitants in 1982 to 14.8 million in 2000, thus marking a 60% increase in only 16 years (Mrini and Nachite 2008). However, this growth is not uniformly distributed alongside the Moroccan Mediterranean coast. It is mainly focused on coastal areas of the Tangier peninsula in general and Tetouan and Tangier cities in particular (Mrini and Nachite 2008). Thus, while the average annual growth rate from 1994 to 2004 was 1.4 percent at the national scale, Tetouan province was marked by a 2.1 percent rate, which is among the highest nationally (RGPH 2004). Also, according to the Center for Demographic Studies and Research (CERED), the coastal municipalities show the highest annual growth rates with 5.18%, a rate that can be maintained for a long time according to future projections (e.g., 5.24% is projected from 2004 to 2025).
This change in the demography and population density distributions translates into higher probability of natural hazard occurrence and introduces additional challenges for authorities responsible for risk management and mitigation. This is especially true in the Martil floodplain because of its relatively faster growth of 300% from 1971 to 2004. The population doubled again from 2004 to 2014 according to RGPH (2014), thus making it the fastest-growing town/city of the country. This creates high demand for real estate in the area, thus exposing the population to more risk. This is more true in the case of the economically disadvantaged people who are forced to adopt unregulated and sometimes ill-adapted construction practices that can cause great economic damage in future. A clear example of this situation is the Diza district, which is built on backfill material that covers old swamp deposits (Khabali et al. 2003). This is proven by our results that show high Kg values in the Diza area due to the presence of a high concentration of moisture in the subsoil layers.
Despite this, two major earthquakes that affected the Al Hoceima city in 2004 (Mw = 6.3) (El Alami et al. 1998;Jabour et al. 2004;Stich et al. 2005;Akoglu et al. 2006) and later in 2016 (Mw = 5.2) (Kariche et al. 2018), did not cause any significant site-effects in the study area. The seismically active zone, which produced these earthquakes that left behind major socio-economic damages, is far away from Martil plain, and thus, the effects of the seismic occurrences are not fully felt. Other than these two occurrences, historic events such as the 1755 tsunamigenic earthquake Gutscher et al. 2006) and the infamous Ghomara earthquake (El Mrabet 2005) may have induced significant site-effects at the study area. However, their impacts are not documented and thus may never be revealed.
Our best bet is the use of scenario-based seismic hazard modeling techniques such as the Peak Ground Acceleration (PGA) approach, which may provide important information regarding the impact of major earthquakes and the possibility of triggering destructive liquefaction processes. This would help to identify areas that are most vulnerable to liquefaction during earthquakes and inform the development of appropriate measures to mitigate these risks. Data from Southern Spain suggest high PGA values including site-effects near the coastline . Given the geological similarities between both the latter region and the study area (Didon et al. 1973), similar results are expected in Northern Morocco. However, a study is needed to prove the presence of such high levels of seismic hazard.

Seismic site-effects in Northern Morocco
The spatial distribution of the vulnerability index map shows that the largest kg values are located in the southern parts of the plain, more precisely along the Oued Martil main stream. Based on Eq. (2), the value of kg is proportional to that of A0. And knowing that the highest values of A0 reflect a strong impedance contrast between soft, water-rich sedimentary layers and their respective bedrock (Gosar 2010), such a distribution is indicative of a thick sedimentary cover at this portion of the basin. In fact, the distribution of the kg values in Martil basin is irregular, even within the same Quaternary deposits. The high kg values found in area characterized by high amplitude (areas near the coast or rivers). Differences in the values may be related to the irregularities of the lithology and the existence of the weathered soils and diverse sedimentary deposits of Oued Martil. Moreover, a satisfactory correlation between fundamental frequency and thickness of the soft sediments is observed, though with large dispersion. This may be due either to errors of the bedrock depth estimation or to possible existence of additional superficial strata (e.g., soft colluvium, man-made deposits) which were not included in the subsoil profile. Targeted geological, geophysical and geotechnical data collection is needed to shed light on this issue and refine our understanding with respect to seismic risk of the Martil basin. This can be explained by the fact that the bedrock blocs in the Northern segment of the floodplain underwent a relative uplift in consequence to neotectonic processes. This fact was evidenced by the topographic position of marine terraces which were shown to be more elevated in the N compared to their southeastern counterparts (Rampnoux et al. 1977;El Kadiri et al. 2010a;Romagny 2014). Also, at the historic Roman site Tamuda (Fig. 1), a Tensifian (Moroccan quaternary subdivision) fluvial terrace is buried under current alluvial deposits, which suggests a recent neotectonic disturbance (El Gharbaoui 1977). Being that this occurred in the Southern part of the said basin, it further exposes the above-mentioned asymmetrical deformation along the major Tetouan fault, marked by differential subsidence rates. The latter mechanisms explain the current surface and subsurface flow that favor the migration of water streams toward the South of the plain, thus increasing the potential for seismic site-effect occurrence.
With regard to the technique used in this paper, after Nakamura (1989;, many researchers have used the same method (HVSR) for the purpose of seismic site-effect characterization (Navarro et al. 2002(Navarro et al. , 2013Gaspar-Escribano et al. 2010;Gosar 2010Gosar , 2017Choobbasti et al. 2014;Akkaya 2020;EL Hilali et al. 2021). In Northern Morocco, only two research papers can be considered with regard to the matter of seismic hazard characterization using this technique.
The first is authored by Arab et al. (2021), who mainly focused on the central Rif region. However, we believe that the data points used by these researchers are not sufficient to characterize such a vast and geologically heterogeneous area. In fact, these researchers only present one H/V measurement per sedimentary basin, which does not allow establishing realistic seismic microzonation maps. Compared to our research which presents the site-effect intensity variation across the whole sedimentary body of the Oued Martil floodplain, these authors try to compare site-effects across multiple basins based on one H/V station per basin, which may induce large sampling errors that mischaracterize their seismic behavior. Therefore, we cannot compare our findings to theirs. Moreover, in the scientific literature dealing with site-effects assessment and seismic microzonation, the usual practice is to plot the spatial distribution of site amplification frequencies as contour maps (Bonnefoy-Claudet et al. 2009), which is not achieved by these authors.
The second research effort in the area is presented by Chaaraoui et al. (2021), who used a more or less similar approach to the one adopted in this study. However, the geological and geomorphological setting of the eastern Rif region is drastically different to that of our study area, and thus, our research presents the only reference in the Northern Rif that is yet to be validated by future research.
In southeastern Spain which is characterized by the same seismicity as our study area, Navarro et al. (2002Navarro et al. ( , 2013Navarro et al. ( , 2014; López Casado et al. (2018), used Nakamura's method seismic microzonation purposes in the context of urban development planning. According to these authors, the HVSR method is a very useful that provides good information regarding ground motion. The same can be said about our study area where the performance of this method was deemed sufficient to reach our research aims.

Limitations of the model
According to studies that use ambient noise measurements (Lermo and Chávez-Garcia 1993;Lachetl & Bard 1994;Parolai et al. 2004), they general associated their findings to some limitations related to the performance of the tests under various uncontrolled conditions. Among latter, the amplification functions estimated from the microtremor measurements are only acceptable for simple stratigraphic geometries characterized by horizontal material boundaries. Also, Massa et al. (2014), who were focused on the effectiveness of different spectral ratio techniques in detecting topographic effects, signaled that the H/V technique may be ineffective in identifying the fundamental frequency, especially, in the case of topographic irregularities, since they do not reflect the 1-D configuration. However, in our case study, none of the above-described technical limitations were confronted due to the subhorizontal topography and the simple stratigraphic disposition of Oued Martil basin.
However, Guéguen et al. (2007), who use the H/V method on Grenoble Valley, reported considerable errors when using this technique in the thinnest sediment layers, especially at the edges of the basin where the greatest differences occur, the soil is more heterogeneous and the bedrock topography is outstanding. They also exhibit that the geometry of the basin can have a crucial effect on the resonance frequency estimated by the H/V technique. In addition, this study shows errors of about 10% concerning the determination of the depth of the bedrock directly from the interpretation of experimental H/V frequency values. Nevertheless, in our case study, the resulting Kg value is deemed reasonable and concordant with field observations and geomorphological data. That being said, it is also worth mentioning that the heterogeneity of the Quaternary material in the study area (Fig. 4) and its anisotropic morphology confused the interpretation of some H/V curves because of the secondary H/V peaks that result from the lithological boundaries within the sedimentary cover.
Also, Bonnefoy et al. (2009) point out the performance of HVSR technique is found to be better in case where a clear H/V peak is identified. Conversely, in the case of broadband amplification of strong ground motion, this technique performs poorly. In such cases, the H/V method provides only the lower amplification frequency, as previously observed by Lebrun et al. (2002). This problem was also not evident in our case study since no broadband peaks were obtained in any of our measurement stations.
However, in our case study, the main challenges that we faced were related to data acquisition during windy days and in the areas near industrial zones or near high-traffic roads, which have significant effects on the accuracy of the H/V analysis. To limit such effects, we carefully chose the acquisition time as explained earlier. In fact, the reason for this limitation is still debated. For Konno & Ohmachi (1998);and Bonnefoy-Claudet, Cotton, et al. (2006), it could be related to the horizontal polarization of Rayleigh waves, which is involved in the origin of the H/V peaks. However, Arai & Tokimatsu (2000) and Bonnefoy-Claudet et al. (2008) link such phenomena to the influence of the Love waves, which have a significant effect on the H/V ratio in certain the frequency ranges. Either way, this poses a limitation that needs to be addressed carefully before proceeding to the acquisition step. In highly developed areas, this can significantly the acquisition time and the budget of research projects, which needs to be considered in planning for the execution schedule.

Conclusion
In this research, we use once again the HVSR technique for seismic site-effect assessment in an alluvial sedimentary environment. Our findings, although preliminary, suggest interesting connections between the basin geometry and subsoil structures with HVSR and Kg index spatial distributions. These effects need to be further investigated for they can explain not only seismic hazard in the study area but also in similar geological settings across the Alboran Sea region. The fact that hazardous areas seem to be concentrated mainly in the Southern part of the basin testifies to the importance of conducting multi-station investigations, since the hazard may not be evenly distributed across the same sedimentary body. The more HVSR measurements are made the more reliable the results are. Our findings also stress the importance of not relying in single-station measurements for characterizing vast heterogeneous areas as was done by previous authors. Eventually, we see the importance of investigating the seismic potential of the area to see whether or not seismic events can reach Martil basin and cause the above-demonstrated effects.
Finally, scenario based on seismic hazard modeling techniques such as the PGA approach can be used for understanding the seismic site-effects in Martil floodplain, and for assessing the liquefaction hazard on the built environment. Overall, additional microtremor measurements inside buildings can enhance our understanding of soil-structure resonance and help engineers design more resilient buildings, as well as to develop effective emergency response plans for earthquake events.