Characterization of the shallow subsurface structure across the Carrascoy Fault System (SE Iberian Peninsula) using P-wave tomography and Multichannel Analysis of Surface Waves

The seismicity in the SE Iberian Peninsula is distributed parallel to the coast in a well-developed strike-slip fracture system known as the Eastern Betic Shear Zone (EBSZ). This work focuses on the characterization of the shallow subsurface structure of the Algezares-Casas Nuevas Fault, within the Carrascoy Fault System of the EBSZ. The Carrascoy Fault borders the Guadalentín Depression to the south, which is a densely populated area with extensive agricultural activity. Therefore, this fault system represents a seismic hazard with significant social and economic implications. We have constructed two velocity-depth models based on P-wave tomography and Multichannel Analysis of Surface Waves (MASW) acquired from seismic reflection data. The resulting velocity models have allowed us to interpret the first ~250m depth and have revealed: i) the thickness of the critical zone; ii) the geometry of the Algezares-Casas Nuevas Fault; iii) the depth of the Messinian/Tortonian contact and iv) the presence of blind thrusts and damage zones under the Guadalentín Depression. Our results have also helped us to estimate an apparent vertical slip rate of 0.66±0.06m/ky for the Algezares-Casas Nuevas Fault since 209.1±6.2ka. Our results provide a methodological and backflow protocol to study the shallow subsurface of active faults, complementing previous geological models based on paleoseismological trenches, and can be used to improve the seismic hazard assessment of tectonically active regions around the world.


INTRODUCTION
High-resolution seismic imaging of the subsurface is a long-established technique for fault detection (Ivanov et al., 2006;Mansfield and Cartwright, 1996), which is a challenging task due to near-surface effects, strong lateral velocity contrasts across steeply dipping structures, presence of soft soils, velocity inversions and/or fault-zone heterogeneities, among others (e.g. Alcalde et al., 2013;Bruno et al., 2010;Catchings et al., 2014;Improta and Bruno, 2007). Despite those difficulties, these studies provide very valuable information about the shallow subsurface structure in active tectonic zones, allowing to characterize fault zones and their associated structures, map blind faults and damage zone areas, image colluvial wedges (Improta and Bruno, 2007;Mattson, 2004;Villani et al., 2017), and to define the depth of the critical zone (e.g. Handoyo et al., 2022;Parsekian et al., 2015).
An interesting area to study the shallow subsurface structure of active faults is the Eastern Betic Shear Zone (EBSZ), in the SE of the Iberian Peninsula (Fig. 1A, B). The EBSZ is an approximately 450km long network of NE-SW-oriented faults (Banda and Ansorge, 1980) that runs parallel to the coast of the Mediterranean Sea (De Larouzière et al., 1988;Silva et al., 1993;Silva, 1994). The EBSZ is considered the most active fault system in the Iberian Peninsula (García-Mayordomo, 2005;García-Mayordomo et al., 2007;Gómez-Novell et al., 2020). This area is characterized by an intense earthquake activity which has caused significant damage during historical times (e.g. Gómez-Novell et al., 2020) (e.g. the 1518 AD Vera (VII-IX) (e.g. Silva et al., 2003); 1522 AD Almería (IX) (e.g. Gràcia et al., 2006) and 1829 AD Torrevieja (IX-X) (e.g. Alfaro et al., 2002) earthquakes; Fig. 1A, B). One of the most recent and destructive earthquakes took place on the 11th of May 2011. It had a Mw of 5.1 and affected the city of Lorca causing severe damage and 9 fatalities. The earthquake was rapidly located and associated to the Alhama de Murcia Fault (Martínez-Díaz et al., 2012a, b;Vissers and Meijninger, 2011), an oblique reverse fault that cuts across the city of Lorca and extends for ca. 100km in a NE-SW direction (Martínez-Díaz, 2002Masana et al., 2004;Sanz de Galdeano et al., 2020) (Fig. 1A, B). However, the Alhama de Murcia Fault is not the only structure that represents a seismic hazard in that area. To the south of the Guadalentín Depression, the Carrascoy and Algezares-Casas Nuevas Faults (Fig. 1B, C) also represent a major threat for the city of Murcia and other populations nearby (García-Mayordomo et al., 2007;Gaspar-Escribano et al., 2008;Gómez-Novell et al., 2020), as well as to the agricultural industry. Specifically, for the city of Murcia the official seismic hazard map of Spain shows a peak ground acceleration on the rock of 0.23g with a 10% probability of exceedance in 50 years (475-year return period) (Instituto Geográfico Nacional-Universidad Politécnica de Madrid (IGN-UPM), 2013).
Despite the great amount of works carried out in the EBSZ, the shallow subsurface geometry of the southern border of the Guadalentín Depression, remains uncertain, a matter that is of paramount importance for producing reliable models for seismic hazard assessment. Revealing its subsurface characteristics, locating fault strands at depth and identifying blind faults, can help to understand its structure and how it accommodates the deformation. The INTERGEO project, supported in 2013 by the Spanish research program, was conceived to unravel the morphology of different active faults of the EBSZ by the acquisition of six high-resolution controlled source multichannel seismic transects. Here, we show the model obtained from one of these lines, which samples the Algezares-Casas Nuevas Fault, an oblique reverse fault formed in a strike-slip regional context (Martín- Banda, 2020;Martín-Banda et al., 2016, 2021Rodríguez-Escudero et al., 2014).
The tectonic relevance of active faults makes them singular targets for high-resolution seismic imaging approaches. Faults and related damage zones feature variations in the physical properties that can be sensed by seismic methods (Ando and Yamashita, 2007;Iacopini et al., 2016;Martí et al., 2006b). Large-scale fault zones can be studied by assessing the changes of the elasticity properties across the fractured area. If the contrast in physical properties at both fault blocks is strong enough, the volume, geometric distribution and extent of these structural elements can be traced by indirect seismic methods. Geophysicists have been using surface waves to characterize the subsurface at a variety of scales; from the entire lithosphere (e.g. Palomeras et al., 2017), to highresolution near-surface studies (e.g. Socco et al., 2010). Specifically, at our scale of interest (i.e. a few hundreds of meters), the Multichannel Analysis of Surface Waves (MASW) is perhaps one of the most extensively used schemes for subsurface characterization (Foti et al., 2018;Ivanov et al., 2006;Miller et al., 2000;Park et al., 1999;Xia et al., 1999Xia et al., , 2000Xia et al., , 2002. The MASW method has been applied to detect a wide range of subsurface features, and it has proven to be successful in other regions worldwide, reaching a resolution depth 1m to >50m (Ivanov et al., 2006;Miller et al., 2000). Complementary, seismic tomography inverts first arrival travel times to generate the P-wave velocity distribution in the subsurface. This technique is broadly used to unravel the seismic velocity structure of faults showing recent activity (e.g. Feenstra et al., 2016) as well as to image colluvial wedges associated to paleoearthquake surface ruptures, much deeper than the range reached by standard paleoseismic trenches (Gaždová et al., 2015).
To unravel the near-surface structure across the Algezares-Casas Nuevas Fault we present here a P-and S-wave seismic velocity (Vp and Vs, respectively) depth profiles based on P-wave velocity tomography and MASW. We have established a MASW processing and interpretation flow to obtain reliable Vs models that combined with the P-wave velocity tomography study provided accurate models of the subsurface down to ca. 250m. The resulted Vp and Vs models provide new constraints on the structure of the Algezares-Casas Nuevas Fault, the location of blind faults, the alluvial thickness in the Guadalentín Depression and the depth of the stratigraphic contacts. Finally, in this contribution we thoroughly examined the results of 1D S-wave velocity and 2D pseudo-S-wave velocity from the surface waves inversion that combined with the P-wave tomography method to characterize the shallow subsurface structures in an active fault system. This detailed analysis allowed us to assess the extent to which this method combination can be used to successfully characterize geological features in the subsurface, which is critical to improve seismic hazard assessments in this type of active areas.

GEOLOGICAL SETTING
The EBSZ is the longest fault system of the Betic Cordillera ( Fig. 1A), which together with the Tell Fold Belt in north Africa, accommodates the convergence between the Nubian and Eurasian plates since the late Neogene (DeMets et al., 2015;De Larouzière et al., 1988).
At the end of the Messinian-beginning of the Pliocene, the tectonic convergence under a NNW-SSE regional shortening direction (DeMets et al., 2015;Echeverria et al., 2011;Khazaradze et al., 2008;De Larouzière et al., 1988) inverted previous Neogene basins, and the Guadalentín Depression, or tectonic valley (Sánchez-Roldán et al., 2021), formed with a NE-SW orientation in an area that was previously uplifted (De Larouzière et al., 1987;Montenat et al., 1987;Vissers and Meijininger, 2011). This regional shortening direction has kept steady in the region since the late Miocene (Martínez-Díaz, 2002) to the present, as evidenced either by sea-floor magnetic anomalies (e.g. De Mets et al., 2015) or GPS data (e.g. Argus et al., 2010). The current tectonic regime explains the fault kinematics and the strike-slip and normal component of the active faults in the area (Sánchez-Roldán et al., 2021).
The study area is located in the northern part of the EBSZ (Fig. 1B). There, the Guadalentín Depression is bordered to the north by the Alhama de Murcia Fault and to the south by the complex Carrascoy Fault system, a 33km long NE-SW structure (Silva, 1994). The Carrascoy Fault System is formed by two overlapping segments, with different geometrical, structural and kinematic characteristics (Martín- Banda et al., 2016). The NE segment extends in a N50E trend for 16km between Los Ramos and El Palmar villages and controls the location of the Cresta del Gallo Range (Fig. 1C). The SW segment is more complex, represented by: i) the Carrascoy Fault, that uplifted the Carrascoy Range and ii) the Algezares-Casas Nuevas Fault, that runs for 23km from the villages of Algezares to Casas Nuevas (Fig. 1C). The Carrascoy Fault is a sub-vertical fault with a general left-lateral strike-slip kinematics and the Algezares-Casas Nuevas Fault develops as a low-angle oblique reverse fault (Martín- Banda et al., 2016). The surface geology and the information derived from geomorphological mapping, morphotectonic analysis and paleoseismological trenches performed along the Carrascoy Fault System, revealed the existence a system of thrusts related to the Algezares-Casas Nuevas Fault and transported to the north since 209.1±6.2ka (Martín- Banda et al., 2016;Martín-Banda, 2020;Martín-Banda et al., 2021). The Algezares-Casas Nuevas Fault constitutes the currently active front of the Carrascoy Range whose development progressively moves to the NW (Martín- Banda et al., 2016) forming a fault propagation anticline. The study profile, located at the southeastern border of the Guadalentín Depression (A-A' in Fig. 1D), samples the Algezares-Casas Nuevas Fault where it shows two strands, referred as F1 and F2.
The oldest rocks cropping out along the study profile ( Fig. 1D) correspond to the metamorphic basement of the Internal Zone of the Betic Cordillera, located in the highest part of the Carrascoy Range. Towards the north, bordering the main relief, upper Miocene marine sedimentary rocks (marl, biocalcarenite and limestone) crop out divided into two units of Tortonian and Messinian age. Along and parallel to the Algezares-Casas Nuevas Fault, a band of Pliocene-middle Pleistocene continental deposits (conglomerate, gravel, sand, marl, and limestone) is deformed in the hanging wall of this fault. These Pliocene-middle Pleistocene deposits are locally known as the Unidad Roja (Jerez et al., 2015) and are considered to be old alluvial fans coming from the erosion of the higher parts of the Carrascoy Range. The youngest materials are middle Pleistocene to Holocene alluvial deposits that fill the Guadalentín Depression. They were deposited during five generations of alluvial fan systems, intercalated with valley floor deposits, glacis and endorheic deposits (e.g. Martín- Banda et al., 2016), constituting an example of the transition from alluvial to fluvial sedimentary systems (Silva, 2014).

METHODOLOGY
The seismic data acquisition was designed to characterize the structure of the Algezares-Casas Nuevas Fault at depth. The profile studied has a NNW-SSE orientation, perpendicular to the fault (Fig. 1D). The seismic data recorded were processed to obtain the Vp and Vs models whose details are described below.

Seismic data acquisition
The seismic data acquisition consisted of a 240-channel system built up by connecting 10 GEODE recording units with 24 channels each, kindly provided by GIPP-GFZ (Potsdam, Germany). The seismic source used was a 200kg accelerated weight-drop provided by the Instituto Técnico Superior of the University of Lisbon (Lisbon, Portugal). The sample rate was 1ms and the total recording time was 4s. The data were recorded using conventional single vertical component exploration geophones with a natural frequency of 10Hz (Table 1).
The acquisition was carried out by rolling along the recording spread (cable). Once the source was located at the spread's center, the recording spread was moved just before the current source location and the process started again, following a leapfrog scheme. The source was fired every 6m moving along the spread until the source position reached, approximately, half of the spread. The geophone spacing was 2m, resulting in a maximum offset of nearly 500m. The total length of the profile is approximately 2496m and a total of 433 shots were recorded (Fig. 2).

P-wave Travel Time Tomography
P-wave travel time tomography is currently a wellestablished and broadly used inversion scheme to resolve Vp velocity structure. The reader can find comprehensive reviews of the methodology (Aki and Richards, 1980;Nolet, 1993;Thurber and Atre, 1993). Common approaches use the first arrival travel times and search for the most plausible velocity model that can reproduce the observables by minimizing the time difference between the estimated travel times. Theoretical travel times are thus calculated by using a ray tracing forward modelling scheme. In our study, the first arrivals were handpicked from the shot records ( Fig. 3A). A total of 84373 travel time picks from the 93528 available traces were picked, representing 90.2% of the total data available (Fig. 3).
To determine the subsurface velocity distribution from the first arrivals picked travel-times, we used the Delta-T-v inversion method (Rohdewald, 2011). This technique is based on the Common Mid-Point (CMP)-refraction concept (Gebrande and Miller, 1985), that considers the CMPtravel times can be seen as a function of the independent variables CMPx coordinates and the CMP-constant offsets (Fig. 3A). In that sense, two partial differentiations deliver the reciprocal CMP-apparent velocity. In principle, this technique is similar to a forward and reverse-shot analysis with the advantage that extrapolations to the shot points are not necessary, but local layer thicknesses H(X), are found at each CMP. The refractor is obtained as an envelope of circles with a radius around the CMPs at the surface.
The CMP scheme assumes a layered media characterized by a velocity step function (with positive or negative Vp increments) and each layer features a constant velocity gradient (Rohdewald, 2011). It starts with determining the velocity at the base of a layer from CMP travel times curves and then it inverts numerically the velocity at the top of the gradient layer. This approach known as "seismic stripping" can be understood as to physically lower the source and receiver for each ray to the top of the layer below. The algorithm is able to model the diffraction of seismic rays at the top of constrained low-velocity layers. It automatically identifies precise time delays on CMP curves transforming these delays into velocity-depth anomalies. A 1-D velocitydepth function is constrained beneath each CMP. All 1-D velocity-depth functions are integrated through a gridding scheme building up a final 1.5D velocity model. A simple and smooth 1D velocity model is needed to initialize the process. This is obtained by laterally extending a simple 1D layered model along the profile. Figure 3C shows the first-arrivals times in CMP scheme, in this way effects of dipping layers are averaged and minimized. The travel-times are smoothed by stacking CMP-sorted traveltime curves over a 40 adjacent CMP's. Subsequently, each curve is "Delta-t-v inverted". For the Vp final model (Fig. 4), FIGURE 2. Example of seismic shot records acquired along the profile. The shot gathers displayed are raw shots, the seismic traces have been normalized by the maximum amplitude for display purposes.
a Root Mean Square (RMS) error was of 9.09ms, resulting in a normalized RMS error of 3.1% (RMS error, divided by maximum pick time of all traces modelled). The relative misfit function was estimated to reach 41.35ms (squared error summed over all traces modelled, divided by the modelled trace count). The maximum absolute error is 120.43ms for trace 81 in shot 42. Thus, these quantitative indicators suggest that this Vp final model shows a relatively high degree of reliability. However, because of the method used assumes 1.5D geometry we cannot rule out the effects of the real 3D contributions to the final velocity anomaly model.

Surface-wave analysis
Surface waves have a lower velocity than body waves and usually present strong energy in reflection seismic records, commonly known as ground roll. While the high amplitude surface waves constitute a serious problem in conventional P-wave seismic reflection processing, where they need to be attenuated or even muted out, in the pursued surface wave analysis they constitute valuable data. The depth sampled by a particular frequency component of surface waves is proportional to its wavelength. This characteristic makes the surface wave velocity frequencydependent, i.e. dispersive. The dispersion diagram reveals the surface wave field in terms of frequency and phase velocity. This representation can be inverted to constrain the S-wave velocity (Vs) distribution as a function of depth, allowing to detect effectively Vs anomalies at shallow depths (between 1.5 and 100m) (Lai and Rix, 1999;Lai et al., 2002Lai et al., , 2005Rix and Lai, 2005;Xia et al., 2000). When this is applied to multichannel data, it is known as  Multichannel Analysis of Surface Waves (MASW). This methodology was introduced in the late 1990's and early 2000's Xia et al., 1999Xia et al., , 2000Xia et al., , 2002 and it is able to resolve 1D and 2D Vs models of the shallow subsurface. The conventional approach to address the shearwave analysis comprises two main steps (Foti et al., 2018). The first step is to determine the experimental dispersion curve; the second is to determine the parameters that define the layered model, i.e. the S-wave velocities and densities of each layer. To resolve the inversion problem, the model parameters are defined accordingly (Lai and Rix, 1999;Lai et al., 2002Lai et al., , 2005Mari, 1984;Xia et al., 2000).
The seismic data presented here features a high amplitude cone of dispersive surface waves, which is very prominent in the shot records (Fig. 5). The surface waves appear within a triangular area with the vertex of the triangle located at zero offset. The frequency content in this region is mostly distributed between 5-50Hz, with the higher amplitudes being within the 10-45Hz frequency band. The observer notes were used to incorporate the acquisition geometry that is the source and receiver UTM coordinates.
This information was used to obtain offset an azimuth for each trace in all shot records. Then, to obtain the dispersion characteristics of the surface waves, we used the frequency/ slowness (ray-parameter: f, p) derived from intercept time and ray parameter (τ-p slant-stacks in frequency domain) in different shot gathers (Fig. 5C, D). As a result, we obtained the dispersion diagrams (Fig. 5C, D). The f-Vph (frequencyphase velocity) picking was supervised by the operator, and only the fundamental model was considered.
The second step allows for more alternative procedures as it focuses on the inversion strategy. Although, in general, the inversion strategy is based on the minimization of the misfit between the observations and synthetic simulations (Menke, 1984;Snieder and Trampert, 1999;Tarantola and Vallete, 1982), the evaluation of this function can be achieved in different forms (Socco et al., 2010). The inversion problem is most commonly solved by linear approximation that uses a 1D forward modeling scheme and retrieve a 1D S-wave velocity depth function for each dispersion diagram (obtained from each shot record). The general assumption is that the subsurface is constituted Characterization of the shallow subsurface structure using P-wave tomography and MASW 8 by a number of layers with constant physical properties for each layer. The resulting velocity-depth function, and/or subsurface physical model can be considered to be represented by a series of "n" layers with shear wave velocities (Vs1, Vs2, …, Vsn), densities (ρ1, ρ2, …, ρn) and thicknesses (t1, t2, …, tn). Through inversion theory (Constable, 1987;Menke, 1984;Parker, 1994;Tarantola and Valette, 1982) the scheme aims to constrain a model (unknown) that is able to reproduce the observables as closely as possible. The inversion is nonlinear and thus an iterative approach is devised, guided by constrained optimization algorithms. The difficulties of such schemes are widely discussed in inversion contributions (Constable, 1987;Menke, 1984;Parker, 1994) and are mostly due to the non-linearities of the problem and/or the non-uniqueness character of its solution. In order to avoid such issues, we considered a strategy that aims to favor the simplest (i.e. the smoothest model) that is able to most closely reproduce the measured data. This approach is known as Occam's Inversion. The mathematical details on the differences between the approaches can be found in Haney and Qu (2010).
The Occam's Inversion approach was used here to address the site characterization of the fault in the study area. This algorithm has been first described and used extensively by the electromagnetic and magnetotelluric subsurface geophysical characterization community (Constable, 1987;Key, 2009;Werthmüller, 2017). In this approach the subsurface is over-parametrized by using a large number of relatively thin layers (Constable et al., 1987;de Groot-Hedlin and Constable, 1990). The strategy of the algorithm aims to find the smoothest velocitydepth model subject to the constraint of a specified misfit between measured and theoretically computed data, taking into account the uncertainties associated with the data measurements. It is an iterative inversion scheme that requires a starting model, which is consisted of the density, Vp, Vs, and layer thickness (Haney and Qu, 2010;Press et al., 1992;Xia et al., 1999). In our case, the starting model displayed a constant Vs of 500m/s and a density of 2.0g/ cm³. This constant velocity model features over 20 layers (all with identical physical properties), with a thickness of 0.5m. The input data consists of frequency and phase velocity picks and error estimates.
The output consists of velocity and standard deviation estimates as a function of depth. The standard deviation is understood as the estimated variability of the inverted velocity. The inverted velocity profiles reach depths down to 250m, and the resolution power decreases with depth, as indicated by the systematic increase of the standard deviation. Therefore, velocity estimates are considered resolved when their standard deviation is below 15%. On average, this assumption limited the resulting models down to a depth of 150m below the surface. The analysis and inversion were applied to individual shot gathers, and thus 433 analyses and their corresponding 1D velocity-depth functions were obtained. The analyzed shot records feature 6m spacing, and thus the 1-D velocity-depth functions are also assumed to be determined every 6m. Considering that the highest amplitude contributions to the signal are within the frequency band of 10-45Hz, the velocity estimates are assumed have averages over 25-35m length.
The inversion result should not depend on the a priori established information of the starting layered model (n, the number of layers; and their physical properties and geometry). The specific characteristics of the input model (physical properties and layer thickness) are not essential in order to minimize the misfit. Thus, to test this point (the dependence of the result with the starting model), three different starting models were used in which the layer thickness, the velocity, and the maximum depth of the model were changed (depending on the lowest frequency of the data). This was done in a way that allowed to statistically sample the model space and maximize the possibilities that the results do not correspond to a local minimum. In all inversions, a maximum smoothness and regularity in the solution are imposed. This scheme avoids resulting models that are too complex. The output consists of a smooth model and its standard deviation for each shot record. The starting models have included depth down to 400m and layer thicknesses that have ranged from 0.5 to tens of meters. Thus, the preferred starting model consisted of 600 layers, 0.5m thick layers which reached a depth down to 250m.
Once the corresponding 1D velocity-depth functions were obtained for each shot record, they were pasted together to build up a 2D velocity depth model (Fig. 6A). A spatial smoothing operator was applied to the 2D velocity model, in order to spatially balance small-scale sharp anomalies (Haney and Qu, 2010;Key, 2009;Werthmüller, 2017). This operation can be justified by considering that the inverted model is an average estimate over 25-35m wavelength combining the frequency of the seismic signal and the offset contributions. Furthermore, each velocity-depth function is the smoothest model that can reproduce the observations (dispersion diagrams) within a 10-12% standard deviation. A direct comparison between the composite velocity model and the smooth version is displayed in Figure 6B.

RESULTS
The P-wave velocity model The resulting Vp-depth model (Fig. 4) provides an image of the first 250m of the subsurface, allowing us to characterize the Unidad Roja (alluvial fans sediments), the weathered Messinian marlstones and the rest of Messinian and Tortonian rocks located underneath. Low Vp (<1500m/s) are modeled in the shallower part of the profile as a zone that increases slightly in thickness towards the north, into the Guadalentín Depression. Underneath this low velocity zone, intermediate Vp (1800-3500m/s) shows a variable lateral thickness implying changes in the Messinian rocks most probably due to the effect of the faults affecting the area. A Vp contrast of 3500-4500m/s is identified at 80-90m below sea level along the northern part of the profile and at 40-50m above sea level in the southern part. This Vp contrast has allowed us to interpret the position of the Messinian/Tortonian boundary. In general, Vp values of carbonate rocks (upper Miocene units) range from 3500-4500m/s, which is in the carbonate rock velocity interval in the range of around 3000-6000m/s (Anselmetti and Eberli, 1993;Eberli et al., 2003).

The S-wave velocity model
The resulted Vs model is shown in Figure 6B. This profile provides reliable information of the ~150m depth and reveals a low Vs (200m/s) layer subparallel to the surface, which can be correlated with the critical zone (Befus et al., 2011;Handoyo et al., 2022). This shallow layer is approximately 15m thick in the southern part of the profile and reaches a maximum thickness of 30-35m at a distance ca. 1400-1500m. Following the description for site characterization of the National Earthquake Hazard Reduction Program (NEHRP, Brown (1981), Table 2), the shallower levels of the profile (Vs ranging from 200 to 500m/s), would be consistent with moderately weathered to weathered sedimentary rocks domains (Types V and IV in Table 2), correlated with the upper Miocene to Holocene sediments cropping out along the profile (Fig.  1D). In addition, this layer thickens towards the north of the northern branch of the Algezares-Casas Nuevas Fault, where the topography is lower corresponding to the sedimentary infill of the Guadalentín Depression.
Under this layer, Vs are significantly higher with abrupt lateral changes. The Vs model derived in this analysis vary from 200m/s (near the surface) to about 1800m/s at depth and presents lateral changes with areas of relatively low-velocity values 200-750m/s, moderate-velocity 750-1300m/s, and high-velocity over 1500m/s. Towards the southern part of the profile, high Vs are located starting at approximately 50m below the surface with their value A B FIGURE 6. Vs model derived from the surface-wave analysis. A) Vs-depth model considering each profile as they were resolved by the inversion form in each of the shot records. B) Vs-depth model once it has been smoothed by considering 35m smoothing operator (see text for explanation). Thin black lines correspond to iso-velocity contours and the numbers indicate their value. The plot has a vertical to horizontal relation of 2:1. Black vertical lines in (A) and (B) indicate the position of fault strands F1 and F2 on the surface (Fig. 1D). o l o g i c a A c t a , 2 0 . 9 , 1 -1 9 ( 2 0 2 2  Characterization of the shallow subsurface structure using P-wave tomography and MASW 10 decreasing towards the southern branch of the Algezares-Casas Nuevas Fault at a distance ca. 900m (F2 in Fig. 6B). The area of Vs between 900-1200m/s, can be correlated to a slightly altered and fractured area, and/or to a minor degree of relatively unconsolidated rock (Type IV in Table  2). This could indicate an unconsolidated domain most probably a fractured zone associated to the Algezares-Casas Nuevas Fault. This Vs zone is laterally interrupted by a relatively broad slow domain (~800m/s) representing a relatively sharp decrease coinciding with the outcrops of the Guadalentín Depression and limiting with the northern branch of the Algezares-Casas Nuevas Fault (F1) at a distance ca. 1400m. Towards the northern end of the profile, at approximately distance ca. 2000m, velocities increase again, reaching the velocity >1000m/s, from a depth down to 50m. These velocities are above 1000m/s and most probably are indicative of harder rock or more consolidated domains in the subsurface (Adel et al., 2013;Kanli et al., 2006;Olona et al., 2010).

Contribution of Vp and Vs modelling to shallow subsurface characterization
The presence of fractures and/or faults greatly influence the physical properties of the shallow subsurface and, often, the effects of these changes in physical properties are observable by means of seismic methods. Fracturing involves a decrease in rigidity (bulk and shear modulus) of the faulted area and surrounding materials. In addition, the presence of groundwater (fluid content) can further reduce the value of Vp and Vs (Catchings et al., 2014(Catchings et al., , 2020. As an implication, the velocity values of Vp and Vs around the fault plane decrease with respect to the surrounding rock mass. Previous studies including empirical methods using laboratory data (Wang et al., 1978), seismic velocity from borehole geophysics (Boness and Zoback, 2004), controlled source seismic reflection (Alcalde et al., 2013;Catching et al., 2014Catching et al., , 2020Martí et al., 2002Martí et al., , 2006aMartí et al., , 2008, and earthquake data (Eberhart-Phillips and Michael, 1998;Thurber and Atre, 1993) reported a significant decrease in seismic velocity (Vp and Vs) within the fault zone, which in can reach up to 50%. In the subsurface, Vs is strongly influenced by the shear modulus and Vp is strongly influenced by the bulk modulus.
Fracturing and weathering reduce shear wave velocity and control surface wave dispersion Xia et al., 1999). Therefore, the velocity and dispersion characteristics of surface waves (MASW method) are used to determine the variation of shear wave velocities with depth (Foti et al., 2018;Park et al., 1999). The MASW technique assumes that the subsurface is a laterally homogeneous layered earth media, in which the average  Odum, et al., 2007;Sun et al, 2015]). The latter is included in the last two columns.  (Odum et al., 2007;Sun et al., 2015)). The latter is included in the last two columns velocity can be assigned to the middle of the receiver spread (Foti et al., 2018;Socco et al., 2010). The MASW method (acquisition, computation and analysis) allows us to investigate shallow subsurface structures in detail, particularly structural features characterized by lateral changes in the seismic velocities. The sensitivity of depth penetration is influenced by several parameters: i) the natural frequency of the observed seismic signal; ii) the array setup (geometry) used for collecting data; iii) the instrument frequency-bandwidth and iv) the velocity sediment/rock of the site (Foti et al., 2018;Park and Carnevale, 2010). Therefore, MASW can be used to characterize lateral shear wave velocity changes in the shallow subsurface.
Following the conventional MASW analysis, in this study we employed surface wave dispersion curves (a total of 416) derived from multiple shot gathers to infer 1D and 2D Vs-depth models. The signal frequencies ranged from 5-10Hz to 45-50Hz (Fig. 5C, D). These frequencies were shown to be sensitive to geological discontinuities at depth, with high frequencies placing constraints on the shallow layers and low frequencies providing information of deeper levels. In addition, the shallowest depth of investigation (Zmin) is resolvable between λ_min/3< Zmin < λ_min/2, where λmin~2dx (dx is the receiver spacing). Meanwhile, the maximum wavelength λmax corresponds to the array length L, where the value is about λ_max≈L (Park and Carnevale, 2010). In practice, the expected maximum depth of examination is about from λ_max/3< Zmax <λ_max/2 (Foti et al., 2018). According to the acquisition geometry of this profile (Table 1), the best resolution is obtained between ~2m and ~150m (Fig. 6).
The use of high-resolution Vs and Vp models, together with the surface geology information (trench and geological mapping) can help to increase confidence in the determination of the shallow subsurface structures and the characterization of the critical zone (Handoyo et al., 2022;Ivanov et al., 2006;Miller et al., 2000;Rempe and Dietrich, 2014). This work provides a successful example of combination of a P-wave tomography and MASW to unravel the shallow subsurface structure in an active tectonic setting. The Vp model provides a reliable image down to ca. 250m and allows to identify stratigraphic contacts and faults. The Vs model complements the Vp model by showing areas of damage due to intense tectonization, and subsequent more intense fluid circulation. The combined interpretation of the resulting Vs and Vp models allows us to characterize: i) the thickness of the critical zone, Unidad Roja and alluvial fans; ii) the depth of the contact between the Messinian and the Tortonian rocks; iii) the geometry of the Algezares-Casas Nuevas Fault; iv) blind faults and a damaged area associated to the Algezares-Casas Nuevas Faults and v) a new slip rate estimation for the Algezares-Casas Nuevas Fault. These outcomes are discussed in the following sections.

Thickness of the critical zone, Unidad Roja and alluvial fans
The high-resolution Vp model presented here has allowed us to accurately locate the contact of the shallower geological units: i) the Unidad Roja and the Quaternary alluvial fan deposits towards the north of F2 and ii) a meteorized layer of marlstone towards the south of F2 (Fig.  7A). The maximum thickness of the Unidad Roja and the Quaternary deposits is identified under the Guadalentín Depression, where they show a thickness of ca. 50m. Complementary, the Vs model provides information about the thickness of the critical zone and a damage area (Fig.  7B) (Handoyo et al., 2022).
Towards the south, F2 shows two strands that most probably merge at depth and form a monocline in the hanging wall, as observed in the surface. The Unidad Roja outcrops until ca. 650m along the study profile ( Fig. 1) and south of this position. The Messinian rocks show also a very low seismic velocity in the first tens of meters, which may correspond to a weathered layer of marlstone. There, the critical zone is formed by different layers of disaggregated materials and weathered bedrock, the so-called regolith, and jointed/fresh bedrock (Befus et al., 2011;Rempe and Dietrich, 2014) that in our profiles is as a layer subparallel to the surface in the Vs model.
Taking a closer look at this layer in the elevated area towards the south of F2, a small thickness difference at both sides of the point of highest elevation at a distance of ca. 600-1000m corresponds to the monocline on the hanging wall of F2. It is relevant to point out that this small thickness difference between the northern facing and south facing slopes is also detected in the un-smoothed velocity model (Fig. 6A). This thin layer of slow Vp and Vs to the south of F2 could be considered as a relatively thin blanket of heterogeneous deposits, loose and unconsolidated rocks (Type IV and V in Table 2) (Francisca and Bogado, 2019) and where fluids flow recharging and/or discharging the aquifers. This domain would be consistent with the regolith. Thicker regolith on north facing slopes compared with south facing slopes of the monocline has been reported from geophysical studies (Befus et al., 2011;Rempe and Dietrich, 2014;St Clair et al., 2015). In these areas, the interaction of tectonic stresses with topography influences the thickness of the regolith. In the study area, the thickness of this shallow layer is strongly influenced by the overall oblique reverse deformation that characterizes the surroundings of the Algezares-Casas Nuevas Fault.

Depth of the Messinian-Tortonian contact
The stratigraphic series in this profile is, from deepest to shallowest: i) Tortonian: limestone sandstones, marls and gypsum; ii) Messinian: conglomerates, sandstone and limestone, and marls in contact with Unidad Roja; iii) Pliomiddle Pleistocene-Unidad Roja: cemented conglomerates; iv) Quaternary: conglomerates less consolidated, the more recent generations are unconsolidated sediments (gravels, sands and silts).
The contact Tortonian-Messinian is marked by conglomerates and sandstones in the Messinian base and marls in the Tortonian top. We have interpreted that the contact between the Messinian and the Tortonian rocks is located along the Vp contrast of 3500m/s to 4500m/s. This model reveals a thickness for the Messinian rocks range from ~50m between F1 and F2, to 170m immediately south of F2 (Fig. 7A). The thickness of the Messinian inferred from geological mapping in the vicinity of the profile varies between 40m and 100m. This is coherent with the contrasts observed in the velocity models and the thickness obtained. This boundary is located at 80-90m below sea level in the southern part of the profile (number 1 in Fig. 7A) and at 172m depth to the northern end of the profile (number 2 in Fig. 7A). The lateral variation on the depth of this boundary may suggest the existence of blind faults that may disrupt it and would contribute to increase the damage reducing the Vs values under the Guadalentín Depression (Fig. 7B).

Geometry of the Algezares-Casas Nuevas Fault
The high-resolution seismic velocity models presented here allow us to identify the structure of the strands of the Algezares-Casas Nuevas Fault (Fig. 7). F1 and F2 correspond to relatively low-angle N-vergent faults which would merge at depth (as proposed by Martín- Banda et al., 2016), although the depth resolution of our models is not enough to image that geometry. In addition, we have interpreted two blind faults to the south of F2 and one blind fault between F1 and F2 based on the displacement of the Messinian-Tortonian contact at distances ca. 200m, 600m, and 950m at the depth between 50m to 0m (in Fig. 7A).
In Figure 6B, F1 is positioned at a distance of ca. 1400m and F2 at a distance of ca. 900m base on their outcrops. From Figure 7A and B, based on the Vp and Vs models, F1 is depicted following the outcrops data Characterization of the shallow subsurface structure using P-wave tomography and MASW 13 that dipping to the south until a depth of ca. 80-100m with the velocity Vp 1700-2000m/s and Vs 600-800m/s. Furthermore, in the deeper zone (Messinian-Tortonian contact), the F1 fault is not characterized by Vs model. The F1 fault is interpreted by velocity Vp 3000-3500m/s at the depth of ca. 0m to -20m. Meanwhile, near the surface (Unidad Roja-Messininian contact), the F2 fault is characterized by velocities of Vp 1500-2000m/s and Vs 600-800m/s at a depth of 150m to the surface. Then, at a depth of about 0-30m (Messinian-Tortonian contact), the fault F2 is represented by the velocity Vp 3000-3500m/s. In addition, the fault propagation anticline in the Unidad Roja associated to F1 could be extended deeper connecting with an area where a south vergent thrust is interpreted to uplift the Tortonian rocks between F1 and F2.
Indeed, F1 and F2 are relatively low-angle oblique reverse faults that displace the Unidad Roja and the middle-upper Pleistocene alluvial fans and, therefore, with activity during the Quaternary (Martín- Banda et al., 2016). Between F1 and F2, the rocks of the Unidad Roja form an anticline in the hanging wall of F1. In turn, the southern strand, F2, may be formed by several branches which may merge at ca. 30-40m above sea level, carrying in its hanging wall a monocline structure. This monocline could be linked to a blind fault on the Messinian-Tortonian contact at a distance of ca. 600m and at ca. 0m above sea level may have led to a change in the bedrock height levels, with implications for tectonic stress in the sub-vertical direction to the surface.

Blind faults and damage zone under the Guadalentín Depression
Fracture systems can play a key role in weathering processes, as they represent pathways for fluid flow, facilitating the interactions of corrosive fluids with minerals and allowing the drainage of neutral (chemically balanced) brines. There is an almost linear relationship between the density of cracks and the decrease in the elastic modulus and thus with the Vs. Theoretical studies by Budiansky and O'Connell (1980) and references therein, derive seismic velocity in terms of aspect ratios (thickness/length) considering the fractures as cracks geometrically described as oblate ellipsoids, and the elastic properties of the matrix and fluids. The larger the concentration of fractures, the higher the degree of weathering interaction and, therefore, the lower the Vs. The differences in the values of velocity anomalies between the interpreted fractures (F1 and F2), probably reflect the effect of the different degree of fracturing and, therefore, of weathering.
Under the Guadalentín Depression, the Vp model reveals up to two south-vergent blind thrusts that disrupt the Messinian/Tortonian contact. We interpret these structures as strands of the Alhama de Murcia Fault which crops out further north, in the northern border of the Guadalentín Depression. Indeed, several branches of the Alhama de Murcia Fault have been interpreted based on gravimetric data (Amores et al., 2002(Amores et al., , 2021, which characterized a highly tectonized area under the Guadalentín Depression. The gravimetric data (density models) obtained in this study was carried out with a trajectory stretching from Carrascoy Fault to AMF through the Guadalentín Depression in the middle. The interpretation of the geological structure from the density model shows that there are several faults under the Guadalentín Depression, which are probably the branches of the AMF. Base on the low-velocity under the Guadalentín Depression, we interpret this area as a damage zone (Fig.  7B). Therefore, under the Guadalentín Depression, a high density of blind faults would have intensely tectonized the area, allowing fluids to infiltrate lowering the Vs values (700-900m/s). The low velocities would then be related to fracturing and local weathering that would reduce the shear modulus value and associated to weathering. In the same area of the Vp profile, not a significant decrease of the Vp is identified (~3000m/s), suggesting that the damage area is caused by high fluid circulation that affected mainly the Vs. In addition to the outcropping strands of the Algezares-Casas Nuevas Fault, we interpret the presence of a blind south vergent thrust that plays most probably in F1 and that forms an antiform on its hanging wall at the distances ca.1650m and 2100m (Fig. 7A). Furthermore, two blind faults are mapped at 50m and 100m below sea level, slightly displacing at the Messinian-Tortonian contact/displacement. In the Vs model (Fig. 7B), the Algezares-Casas Nuevas Fault (area between F1 and F2) is interpreted as a contrast area where a reduction of the Vs may indicate a more tectonized or altered zone.

Slip rate estimations
Based on the above, the shallow subsurface structures in the southern Guadalentín Depression and northernmost edge of the Carrascoy Range have been unraveled. Here, we provide an image of the first 250m depth of the area allowing us to identify the structure at depth of outcropping faults, the locus of several blind thrusts, the thickness variations in the Upper Miocene units, as well as areas of intense fluid circulation.
In addition, the identification of the Messinian/ Tortonian and Unidad Roja/Messinian contacts in both blocks of the Algezares-Casas Nuevas Fault enables us to estimate an apparent vertical slip rate knowing that the onset of the formation of the fault was around 209.1±6.2ka (average value from Martín- Banda et al., 2016). Specifically, the Messinian/Tortonian contact apparent vertical displacement from one block to the other is quoted in 135±15m (points 1 and 1' respect to points 2 and 2' in Fig. 7A). This displacement is consistent and practically identical, 140±20m, taking as reference the Unidad Roja/ Messinian contact respect to the base of the outcrop of the Messinian on the upthrown block located southwards of the seismic profile. Considering these displacements accumulated since the onset of fault activity in 209.1ka, we calculate a vertical slip rate of 0.66±0.06m/ky considering the error propagation from the displacement measurements as well as for the time. This value is slightly higher than the 0.34±0.02m/ky (Ferrater et al., 2017) and 0.37±0.08m/ ky estimated by Martín- Banda et al. (2016), and lower than the 1.3±0.2m/ky, calculated based on GPS data in the Alhama de Murcia Fault (Echeverria et al., 2013). However, our estimated value is very similar to the 0.64±0.4m/ky established for the net slip rate along the NE Segment of the Carrascoy Fault for the last 220ky (similar period) (Martín- Banda et al., 2021). This difference in the slip rate calculated is not surprising, as previous values were calculated in a different part of the same fault, where the Algezares-Casas Nuevas Fault features one major strand, and according to the restitution of the top of the Unidad Roja based on data from trenches (Martín- Banda et al., 2016). Our slip rate value is estimated from the apparent displacement of the Messinian-Tortonian contact in an area where the Algezares-Casas Nuevas Fault depicts a complex geometry of fault strands that may have played a role displacing the Messinian-Tortonian contact.
On the basis of the above, we can conclude that Vp tomography provides a reliable image of the subsurface and that combined with MASW allow to image the geometry of the shallow structure (Anderson, 2015;Anderson et al., 2011Anderson et al., , 2013Parsekian et al., 2015;St Clair et al., 2015), where tectonized or altered areas represent preferential paths for fluids. Both methodologies combined are useful for the characterization of the critical zone that comprises the regolith and the jointed/fresh bedrock. To further study the stress interaction between strike-slip and oblique reverse tectonics, the kinematics of the blind faults that could represent a higher risk as previously though, and the lateral continuations of the active structures, a multidisciplinary coincident 3D geophysical data is required, including resistivity and Lidar, among others (Anderson et al., 2013;Leone et al., 2020). Further work will provide new Vp and Vs models along other faults of the EBSZ providing a larger picture of the shallow subsurface structure of this strike-slip system. The combination of methods presented in this article represent valuable high-resolution resources for the characterization of the shallow subsurface in tectonically active regions, improving seismic hazard models and reducing the potential risk caused by active structure zones in the other areas.

CONCLUSIONS
The P-wave velocity model and Multichannel Analysis of Surface Waves (MASW) have been applied to unravel the shallow subsurface structure of one of the most active fault systems in the Iberian Peninsula, the Eastern Betics Shear Zone (EBSZ). Specifically, we have sampled the Algezares-Casas Nuevas Fault, a relatively low-angle oblique reverse fault located in a left-lateral strike-slip fault system between the Carrascoy range and the southern border of the Guadalentín Depression. In this study, we provide new Pand S-wave velocity (Vp and Vs respectively)-depth models, reaching a maximum depth of 250m. P-wave velocity was determined from the Delta-t-v method and S-wave velocities are estimated from the surface waves recorded within the control source shot records of conventional seismic reflection data (MASW).
The interpretation of these models has allowed us to identify significant features of the subsurface: i) the critical zone and thickness of the Pliocene-Holocene sedimentary cover; ii) the depth of the contact between the Messinian and the Tortonian rocks; iii) the geometry of the strands of the Algezares-Casas Nuevas Fault and iv) the location of blind faults and a damaged area associated to Algezares-Casas Nuevas Faults and possibly to the Alhama de Murcia Fault. Additionally, it has allowed us to calculate a new vertical slip rate estimate for the Algezares-Casas Nuevas Fault. The critical zone along the study profile consists of a layer of low Vp and Vs parallel to the surface and slightly thicker to the north of the profile, the Guadalentín Depression. This layer is mainly formed by sediments of the Pliocene-Middle Pleistocene (Unidad Roja) and the Plio-Quaternary alluvial sediments in the Guadalentín Depression and towards the southern part, it corresponds to weathered Messinian marlstone. The depth of the Messinian/Tortonian boundary is identified along the Vp profile as a contrast of 3500 to 4500m/s. The apparent disruption of the interpreted surface allowed us to interpret at least two antithetic faults. The Algezares-Casas Nuevas Fault is interpreted to be formed mainly by two major N-verging strands (F1 and F2), although other related blind faults can also be interpreted. It is noteworthy the presence of a blind S-verging thrust between F1 and F2 that may form an antiform in its hanging wall. Under the Guadalentín Depression, two S-verging blind faults are interpreted disrupting the Messinian/ Tortonian boundary. Those structures could be related to strands of the Alhama de Murcia Fault, which crops out in the northern border of the Guadalentín Depression, about 5km from these S-verging structures. Specifically, the southernmost part of the Guadalentín Depression features very low Vs under the critical zone. We interpret this as a damage zone resulted from an intense fluid circulation in the aquifers. Finally, we propose a slip rate of 0.66±0.06m/ ky (in the last 209.1±6.2ky) for the Algezares-Casas Nuevas Fault based on the depth of both the Unidad Roja/Messinian and Messinian/Tortonian contacts south of F2 and to the north of this damage zone identified in the Vs model. This work represents a successfully combined MASW and P-wave tomography study to unravel the shallow subsurface structure of a tectonically active zone.