New Numerical and Theoretical Model to Characterize the Upper Crustal Structure of the Moroccan Atlas from Wide-Angle Seismic Reflection Data

The upper crust beneath the Moroccan Atlas has been modeled using the travel time of wide-angle arrivals from the SIMA (Seismic Imaging of the Moroccan Atlas) dataset. The detailed knowledge of the internal structure of this orogen allows understanding its uncommon characteristics, featuring high topography, moderate tectonic shortening and moderate crustal thickening. The ~700 km long SIMA wide-angle reflection seismic profile has provided a high resolution geophysical data set to obtain a shallow P-wave velocity model along the transect. The seismic data processing has enabled us to accurately pick the first arrivals of the seismic records. Subsequently, the development of a numerical code to mathematically model the hodochrones defined by the first arrivals, has provided us the P-wave velocity structure of the crust down to 10 km. The resulting model shows a detailed image of the Atlas upper crust and reveals several relevant features that help to understand the structure of the orogen and its composition.


Introduction
The great development of the seismic wave methods was led by the oil and gas exploration industry, but these methodologies were quickly adopted in other exploration applications such as mineral, water, geothermal energy or waste disposal sites.These methods provide detailed information about the internal structure, the physical properties or the composition of the subsurface and also give insights about the geological evolution and the geodynamics of the target area.This field of study covers a very broad spectrum of ground motions, from earthquakes to very weak seismic pulses.
Two main groups of seismic methods can be differentiated according to the nature of the seismic source: a) passive seismic, that uses natural sources, mainly earthquakes, and also in recent studies the seismic noise has revealed as a powerful characterization tool; b) controlled source seismic that uses artificial seismic sources and is focusing on reflection and refraction seismic studies at different depths, with different resolutions and using distinctive acquisition geometries.
Since the 80's, seismic methods have been one of the most important tools to study the main orogens around the world, providing unique information about their geodynamics and the deep crustal structure.Nowadays, one of the most challenging orogens to study is the Moroccan Atlas.Located at more than 600 km of the Africa-Europe plate boundary, the High and Middle Atlas intracontinental mountain chains feature a really high topography (up to 4000 m).However, shortening, which is basically related to the convergence between Africa and Europe, is moderate, generally less than 20%, [1]- [6].
Multidisciplinary geophysical studies suggest that the maximum crustal thickness in the High Atlas is ~40 km, and accordingly, a general state of Airy isostatic under-compensation is inferred.A Bouguer anomaly of about -120 mGal [7]- [9], implies that the observed crustal thickness is not sufficient to compensate the topography.
All these features, together with the presence of an intriguing alkaline volcanism and a high heat flow [5], make the understanding of the geological evolution of the Atlas Mountains a matter of great interest in the Earth Science research community [1]- [6].
Modern geological, geodynamical and geophysical studies are proposing as explanation to this apparent contradiction, the existence of dynamic topography and the existence of a thin and hot lithosphere beneath the High Atlas, thus contributing to the relief by thermal doming [1], [2].In other words, the mantle is actively contributing to the topography.
In 2010, a North-South high-resolution wide-angle seismic reflection transect was carried out across the three major geological zones in Morocco, the High Atlas, the Middle Atlas, and the western edge of the Rif Mountains.The purpose of this transect was to image the crustal structure and the Moho geometry beneath these mountain chains obtaining a 2D, P-wave seismic velocity model.
The analysis of the data reveals the presence of several conspicuous crustal phases, including those from the Moho.The current manuscript contribution is only

Fadila Ouraini et al.
focused in the study of the first arrivals from the upper part of the crust.The approach used consists on the analysis of the hodochrones, which are the representation of the relationship between travel times and the offset observed in the seismic phases.Additional knowledge of the structure of this area will help us to understand the variations in elastic wave speed and/or reflector depth and geometry.To achieve our goal, we have developed a numerical code, capable to calculate the properties of each layer, using the intercept method, and focusing on the shallowest part (upper 10 km) of the subsurface layering.The resulting 2D model obtained reveals information about the sedimentary (Mesozoic) layering and exhibits a noticeable variation in the velocity gradient along the cross section.

Geological setting
From a geological point of view, the Atlas Mountains are an inverted Rif (reactivated by compressional tectonics) developed over a continental basement.The plateaus bordering or included within the belts are also bound to the Atlas system, e.g., the Middle Atlas "Causse" to the north, the "High Plateaus" (Oran Meseta) to the east.These young mountains were uplifted during the Cenozoic as a result of the Alpine orogeny, such as the Rif Mountains to the north.However, the High and Middle Atlas thrust-fold belts were created by the tectonic inversion of pre-existing extensional basins of Triassic-Jurassic age [4], [5], [21]- [27].Compared to the Alpine-type, collisional Rif belt, the Atlas system is an intracontinental, autochthonous system, developed over a continental crust, which was only slightly thinned during its pre-orogeny evolution.
The Atlas mountains upper crust consists of a few triassic rocks, mostly Jurassic, with the deposition of marine carbonates and shales capped by continental red beds and the inverted Triassic and Jurassic basins of the High Atlas are thrust outward over Cretaceous postrift rocks which overlie basement or very reduced Triassic-Jurassic remains [5].
The geodynamic evolution of the Atlas system comprises two major periods.The pre-orogeny period is characterized by rifting, which affected the Variscan crust, and then by the filling of Mesozoic basins.The orogeny period is characterized by the basin inversion, the shortening of the basement and cover units, and the formation of syn-orogenic basins.
As a compressional belt, the High Atlas exhibits a relatively small crustal root that contributes to the mountain topography, as revealed by gravity [7], lowresolution [8], [11] and a high resolution wide-angle reflection seismic experiment [28].However, the total elevation of the Atlas system (where the mean altitude of the mountain belts exceeds 2000 m over large areas and some of the undeformed foreland basins stand above 1200 m) is only partially explained by shortening and crustal thickening, being the system in a state of isotactic under-compensation at the crustal scale [5], [7], [15], [29].

Data acquisition and processing
A controlled source seismic acquisition experiment consists on an array of sensors (geophones) usually deployed at the surface that record the seismic energy generated by an artificial source (for instance, explosives, Vibroseis trucks, accelerated weight drop or sledge hammer).Every geophone records, during a listening time window that depends on the target depth of the study, the reflected and refracted waves.These are generated by the changes of impedance in the subsurface.The movement of the ground is sampled, digitized and stored on magnetic media.The recorded seismic traces are gathered by shot point and sorted by offset (distance between receiver and shot point) to make up the data set in a seismic survey.
Seismic waves are a perturbation of the environment that propagates in space and time.The propagation depends on the elastic properties and the density of the medium and it can be described by Hooke's law, which relates the stress to deformation, and the second Newton's law that relates the force to the acceleration.These seismic waves recorded on the seismic traces includes wanted seismic energy, like body waves (P-and S-wave) and unwanted energy such as surface waves, Airy waves or guided waves.The study of the desired seismic energy, reflected or refracted, determine the kind of seismic experiment carried out, affecting the geometry of the acquisition, the seismic instrumentation and the later processing of the data set.

Database acquisition
The North-South SIMA (Seismic Imaging of the Moroccan Atlas) wide-angle seismic transect was acquired in 2010 across three major geological zones in Morocco, the recently uplifted High Atlas, the Middle Atlas, and the western edge of the Rif Mountains.This experiment was designed to determine the internal structure of the crust, the Moho geometry and the seismic velocity structure of this complex orogenic system (see Fig. 1.).
To achieve these objectives, 934 Reftek 125a (TEXANS) digital seismic recording instruments from the IRIS (Incorporated Research Institutions for Seismology)-PASSCAL Instrument Pool were deployed along a transect of almost 700 km from the Sahara desert to Tangier.The average distance between the stations was 400 m between Merzouga to El-Hajeb (central Morocco) increasing to 1000 m in the Rif domain.A total of six shot points were distributed along the transect (see TABLE I ) with a distance interval from 40 to 90 km.Every shot point was charged with a 1000 kg of explosives in variable depth shot holes (30-60 m).The shot points were carefully chosen based on the field investigations.The seismic data was recorded in SEGY format (conventional seismic format for reflection and refraction surveys).SEGY data begins with 3200-bytes ascii header providing a human-readable description of the seismic data followed by a 400-bytes Binary File Header that contains binary values that affect the whole SEG Y file.Certain values in this header are crucial for the processing of the data in the file, particularly the sampling interval, trace length and format code.After that a 240-bytes trace header containing the attributes of the seismic trace and after this header, the amplitudes values for every sample are included for every trace.One of the main values is the coordinates of every receiver (seismic trace) and the shot point locations, which were obtained by handheld GPS.
An internationally diverse field crew of 75 faculty and students from more than a dozen institutions in Africa, Europe, and North America conducted the 2-week long experiment in April and May 2010.The project, Seismic Investigations of the Moroccan Atlas (SIMA), is affiliated with the PICASSO program in Spain and Morocco.

Method -Travel time modeling
The seismic energy propagates from the seismic source (shot), in a relatively homogeneous layered media.The seismic waves behave as light does when travelling through materials of varying indices of refraction.In fact they are governed by Huygens and Fermat's principles, and will follow a minimum-time path.
We remember the equation for a plane wave: where k is a vector that points in the direction of propagation and thus, by definition, is a ray and A is the amplitude.For homogeneous material, k does not change.For a variable seismic velocity media, which is the more realistic situation, we must solve the equation: (2) , where E is the Young modulus and ρ the density, and represent the velocity of P-waves.We solve this partial differential equation by assuming a functional form: where W(x).ω/c0, which replaces k.x, is a function of position, and c0 is a reference velocity.Substitution yields:

TABLE I LISTING OF THE SHOT TIMES, POINT COORDINATES AND EMPLACEMENT OF EXPLOSIVES
After computations of spacial derivatives we can rearrange the obtained equation as: For high frequencies (small wavelengths), the term in the right is small, approximately 0, and equation 5 becomes: This equation is called the Eikonal equation, which is a partial differential equation that relates rays to the seismic velocity distribution.We consider the 3-dimensional wave surface shown in Fig. 2. The ray W(x), is characterized by traveling an arc length, s, in a time t.The direction cosines associated with the ray are given by dx1/ds, dx2/ds, dx3/ds, and must satisfy: Now we consider the physical connection between s and W(x): W(x) α s, which is just the statement that the gradient of a function (surface) is oriented normal to that function (surface).Thus we can see that dxi/ds must be proportional to W(x)/xi.This implies that we can rewrite 7 as: where a is the constant of proportionality.After comparing this equation with the eikonal equation we can notice the similarity if: a = c(x)/c0.The value a -1 = n = c0/c(x) is commonly called the index of refraction.
From the combination of the previous equations, we can verify how the normal equation changes along the path of the ray, in addition to taking the derivative of the normal equations with respect to ds: The generalized form of this equation is called the raypath equation: If we follow a ray through a material that has a change in velocity in only one direction, (depth in our case), then  = ( 3 ), and thus  = ( 3 ).Thus / 1 = / 2 = 0. Then: The ratio  1 / 2 confines the raypath to a normal plane to the  1  2 plane.In other words, the projection of the ray into the  1  2 plane is a straight line (see Fig 3).
For convenience, we can choose this plane to coincide with the  1  3 plane, reducing the previous equations to: At a given point the direction cosine of the ray is given by: Thus: The constant p is called the ray parameter, or horizontal slowness.The angle i is the angle of incidence and c is the velocity.We can see clearly that the resulting equation is the Snell's law, which can also be derived from Fermat's principle.
The horizontal distance x at which the wave reaches a depth z can be expressed in terms of the angle i as follows: The time t required for the wave to travel along the raypath from the surface to the point (x,z) within the subsurface medium is given by : As mentioned before the seismic rays outgoing from a source act as light waves, so on an interface, if we consider an incident wave with an angle θ1 respect to the normal, there will be a reflected wave with an equal angle and refracted wave, satisfying: 2 is the angle of refraction.We suppose that the velocity increases with depth.
The maximum value of  2 is 90°, in this case the wave travel horizontally in the second medium, with a velocity V2, and the incidence angle will be called the critical angle θc, which is sin   =  1  2 .
When the wave travels parallel to the interface, it radiates energy into the upper medium at an angle   .It is also called a head wave (Fig. 4) [31], [32].
In the following, the interpretation of time curves will be formulated, for an earth model composed by three layers supposed parallels and a condition governing the velocities:  1 <  2 <  3 (see Fig. 5).From the refraction theory [33], [34], the velocities correspond to the inverse of slopes of each segment in the time curve plot, representing the recorded arrival time on geophones Gi, of an explosion source at S. We will consider the refraction recorded on geophone G7 to calculate the thickness of each layer hi and the displacement pi, to provide to the points at each depth pi (see Fig. 5).In (Fig. 5) Vn are the velocities of each layer n, inm are the angles of critical refraction, ti the intercept times, tci and xci the times and positions of the crossover points, where the refracted waves reach the sensor/geophone first time.
 The first segment with a slope  The critically refracted wave on the first interface has as time equation: , and after computations we get: with  1 the intercept time.
 The refracted wave on the second interface separating mediums with the velocities  2 and  3 has as time curve equation: From where we can deduct after computation that: This leads to deduce the general form to calculate the thickness of each layer as: The displacement pk can be calculated with the formula: )) 1 = (25) From where we can calculate for example  2 : Even though our calculations, assume that the layers are parallel between two successive shots, we know that this criteria is not met because of the important dimension of the profile and the inhomogeneity of the structure along the transect.Accordingly, for a direct and its corresponding inverse shot, we are calculating the velocities and thicknesses independently instead of calculating the average velocity and thickness between the two shots.Consequently, to the North and the South of every shot, the calculations are done separately.
In order to model the picked first arrivals and following the formulation presented in this section, a Matlab code was developed with the aim to identify the characteristics of every lithological unit that features differentiated physical properties.
The first step consisted on the identification of the different slopes characterizing the first arrivals by selecting a number of points that define independent events representing single seismic velocities relative to the subsurface (see Fig. 5).We selected the group of points satisfying the dipping variation in the first arrivals in sections represented with a velocity reduction.Those selected points (travel-times), were converted back to real time (reduction velocity cannot be used) and became the input data of the code.
The code, uses a least-square method, and defines the intercept time and the different slopes observed by fitting the selected points to a line (see Fig. 5).The inverse of the calculated slopes provide the different seismic velocities observed in the subsurface and the intercept time, i.e. the intersection of the fitting line with the time axis.By means of the equations formulated (Eq.24 and 25), the code provides the position relative to the normal intersecting at each shot positions (displacement pk) and the thickness (hk) of this differentiated lithological unit; those locations together with the appropriate velocities are defined as velocity nodes.Using all this information, a velocity grid is constructed showing separated regions.To construct the velocity model we need to have an uniform distribution of velocity nodes.
The subdivision into different regions is a result of the model parameterization where a linear interpolation is applied between velocity nodes in the limits of each layer, and laterally between adjacent velocity values.

Data processing
The raw seismic data set need some previous processing to enable a precise picking of the first arrivals.First of all, the basic information characterizing every seismic trace needs to be introduced in the SEGY headers.This part is mandatory because information about the geometry of the acquisition experiment is needed in order to correctly sort the data by shot number and offset.The resulting shot gathers were edited to remove the noisy traces.To improve the signal-to-noise ratio and the continuity of the first arrivals to large offsets, seismic signal treatment has to be carried out.In the SIMA data set, spherical divergence and offset amplitude recovery were applied to account for the loss of energy depending on the distance.The analysis of the frequency content of the data by means of frequency filters is the key factor to improve the S/N ratio.This processing step tries to remove the unwanted energy keeping the frequency range of the seismic energy.Several frequency filters were tested to determine which frequency range was richer in seismic signal energy (Fig. 6).Similar data has been processed in the same way [35], [36].The band-pass filter applied to the shot gather f = 0.25-2-6-8 Hz basically retains all the original frequencies present in the data.The resulting processed shot gathers clearly imaged the first arrivals and even other seismic phases corresponding to intra-crustal and Moho reflections and refractions, allowing the picking up to large offsets.For imaging purposes, amplitude balance functions were applied and other gain corrections were also tested to assure a more homogeneous distribution of the energy of adjacent seismic traces to make easier the first arrival picking.Fig. 6, shows shot S1 plotted with reduced travel time (velocity reduction = 8 km/s) in order to show the wave propagating at 8 km/s as a horizontal hodochrone separating the phases characterized by lower and higher velocities.
The picking of first arrivals and the imaging of the shot gathers (Fig. 7 and Fig. 8) were carried out using academic software: Seismic Unix, developed at Colorado School of Mines, Generic Mapping Tools and other Linux open source New numerical and theoretical model 293 utilities.Observe that the distances up to which we identify first arrivals is almost 150 km.In Fig. 9, we have an example of the determination of the fitting lines that the code is using to calculate the slopes and the intercept times.From these results, the velocities, the thicknesses and the displacements (Eq.24 and 25) are obtained to build, by interpolation of the different calculated velocity nodes, a final homogeneous sampled grid (Fig. 10).

Results
Wide-angle seismic reflection/refraction data experiments are designed to image deep into the subsurface, typically crustal scale to upper mantle.To achieve these targets the acquired offsets need to be 4 or 5 times these depths to assure the recording of the diving waves that travels through these deep discontinuities [37].Nevertheless, in this contribution, the study of the first arrivals will focus on the upper part of the crust (up to 10 km).The study of the whole crust and the Moho discontinuity has been published elsewhere [28].
In the previous section, the description of the travel time modeling was formulated.In order to better determine the main points defining the different slopes in the first arrivals, the shot gathers with their corresponding picks were plotted in two reduced velocities, 8 km/s and 6 km/s (Fig. 8).The use of both velocities facilitates the identification of the changes in slopes, therefore the selection of input data to the code.
The resulting velocity model (Fig. 11) is presenting a detailed image of the upper crust up to 10 km corresponding to the region comprised between the first and the last shot (see Fig. 1).The main reason for that is because to the North of the shot S1 only two velocity nodes were obtained, located at 0.3 km and 8 km depth respectively.These two nodes are not enough to model the hodochrones, because the first one is very close to the surface and the second one is almost in the bottom of the model, which means that they are not uniformly located.In addition    New numerical and theoretical model 295 surface just in the High Atlas, slightly decreasing to the North.It is also noticeable the presence of another low velocity zone around S2 with an important velocity gradient to depth.

Discussion
The hodochrone mathematical modeling method developed in this paper is using first arrivals, to perform an accurate description of the slopes behavior, identifying their values and the intercept times to create a smooth velocity model.The SIMA data has provided us with good quality picks but with a limited amount of them.Nevertheless, they have enabled us to outline the geometry and the location of the main changes observed in the velocity model.The image obtained for the upper 10 km of the model is essential in the future inversion approaches that might be applied to the whole dataset, including the inversion of different phases focused on the study of the middle and lower crust and the Moho discontinuity.
From a geological point of view, the final velocity model is showing two main differentiated units with depth corresponding to the sedimentary cover and the basement, characterized by their seismic velocities.The maximum thickness of the sediments is clearly located beneath the Middle Atlas and the High Atlas featuring velocities around 4.5 km/s and 5.2km/s respectively at the surface, increasing to values of 5.3-5.4km/s at the bottom of the layer reaching depths of 3 to 4 km.The highest velocities (around 5.2 km/s) are located beneath the High Atlas (at distances of 150-190 km) and match the geology mapped at surface and recent studies carried out in the area [39], [4], [5].The High Atlas is mostly composed of Jurassic rocks, including carbonates, gabbros and Triassic evaporites and basalts that usually feature high P-wave velocities under these conditions of depth, pressure and temperature [40].The total thickness of the sedimentary layer seems to be thinner beneath shots S2 and S3 (210-260 km distance) featuring lower velocities at surface (around 4 km/s) regarding the velocities observed in High Atlas.The velocities decrease laterally until they get their minimum underneath S2.Again, surface geology is explaining this fact: Tertiary and Cretaceous rocks outcrop in the area between S2 and S3 characterized by lower velocities in laboratory measurements and other seismic experiments in orogenic belts [40].
The other main feature observed at surface is located south of High Atlas, beneath shot S5.This area features really low velocity values at the surface.The velocity gradient observed there is a really high, featuring velocity at surface of 2.9 km/s or less arriving to velocities of 5.4 km/s at 2.5-3.0 km depth.The surface geology shows that this area is composed by Cretaceous red beds that include porous clastic rocks and carbonates.These kind of sediments are characterized by little compaction featuring low P-wave velocities in the range that are observed in this area.The strong gradient could be explained by the presence of volcanic rocks at depth and/or by the proximity of the basement.The Mesozoic rifting episode positioned the basement at a lower level in the central part of the High Fadila Ouraini et al. and Middle Atlas.Maybe later inversion was not important enough to bring it to higher depths than to the North and South so velocities may be higher at depth near shot points S1, S2, S5 and S6.The area around S4 has maximum sediment thickness (Jurassic mostly).
The anomalies highlighted in the P-wave velocity model beneath the shots S2 and S5 are characterized by a small thickness, the low velocity zones beneath S2 is reaching 2.5 km thick and 1.5 km under S5.

Conclusion
A wide-angle transect was acquired, crossing the major geological structures in Morocco.The shot records reveal the presence of crustal arrivals, reflected and refracted waves from the Moho, that have allowed imaging the crust using only first arrivals, and focusing on the shallower 10 km.The method based on the development of a numerical code is capable to evaluate the properties of the subsurface underneath the Middle and High Atlas.The velocity structure shows two apparent low-velocity zones reaching 2.5 km and 1.5 km respectively under shot points S2 and S5, that are followed by an important velocity gradient beneath those shots.This feature, correlates well with the geological information about the filling of the Mesozoic basins.In the central part of the model, coinciding with the High Atlas, a relatively high surface velocity coincides with denser Jurassic layers.In that area, the High Atlas region, a moderate gradient was observed.
The comparison of the results with the geological sections confirms the major part of the model.
Our results are providing important details about the weathered layer along the profile, necessary in the oil prospection, being the first stage in a geophysical study in this purpose, and also as input for a future inversion tomography.

Fig. 1 .
Fig. 1.Tectonic map of the studied area, design of the wide-angle seismic refraction experiment and the main geological units [30].Stripe pattern indicates the Rifean plate-boundary orogenic belt, the grayed areas represent the Atlas chains and the extent of the Cenozoic basins is indicated by the dotted pattern.The location of the wide-angle refraction transect is indicated in a white line.Red stars indicate the shots positions, and the black squares the cities crossed by the profile.

Fig. 2 .
Fig. 2. Raypath for a medium in which the velocity is independent of the x2 and x1 directions Fig. 3. Three-dimensional wavefront with a normal ray with length

Fig. 7 .Fig. 8 .
Fig. 7. Observed seismic record sections plotted as a function of offsets, with two kind of reduction velocities, 6 and 8 km/s.The six shots are plotted, with the same processing applied.

Fig. 10 :
Fig. 10: The upper part of the velocity model created; the black line is: the topography; red crosses are: the shot points; green points are: velocity nodes (displacement pk, thickness hk); the values are the calculated velocities.

Fig. 9 :Fig. 11 .
Fig.9: Identification of the points corresponding to different slopes for a near offsets portion of the shot S1, red points are the picks of first arrivals in reduced time by 8 km/s, and the green points with a reduction velocity of 6 km/s, and the blue lines are the best fitting lines.

Fig. 12 :
Fig. 12: Comparison between the obtained model and the crustal-scale section of the High and Middle Atlas; (a) Hodochrones modeling velocity model; (b) Interpretation of geological and geophysical data according to Arboleya et al. [4] and the rectangle represent the area modeled above; ANMA: North Middle Atlas fault; ASMA: South Middle Atlas fault.