Modelling principal stress orientations in the Arabian Plate using plate velocities Open Access

The understating of the present-day in situ stress field in the Earth's crust and its variability is critical for earth science disciplines that require an understanding of geodynamic processes (e.g. Townend and Zoback 2004; Heidbach and Ben-Avraham 2007; Gaucher et al. 2015). In addition, it is important for a range of practical engineering applications, including the production of hydrocarbons and geothermal energy, mine safety, seismic hazard assessment, and underground storage of CO₂ (e.g. Cornet et al. 1997; Fuchs and Müller 2001; Heidbach et al. 2019; Zoback 2010; Moeck et al. 2011; Altmann et al. 2014; Ziegler et al. 2016). The importance of mapping the maximum principal horizontal stress (S_Hmax) has been demonstrated by the World Stress Map (WSM) Project (Heidbach et al. 2018), which shows that lithospheric in situ stress is controlled by the forces exerted at tectonic plate boundaries, drag forces from slow mantle convection acting at the base of plates and gravity-induced deformation. Following on from this work, we focus this study on the Arabian Plate and similarly model the azimuth of S_Hmax. This is particularly pertinent given the region's numerous prolific hydrocarbon resources, associated subsurface developments and engineering projects. Presently, assessing stress magnitudes within the plate remains challenging due to the inaccessibility of most of the existing downhole data and a relative seismic quiescence in the plate's interior (discussed in the section ‘The present-day stress pattern in the Arabian Plate: compilation of published data to augment the WSM database'.).

The Arabian Peninsula is part of a small tectonic plate that is characterized by active deformations along its boundaries with: (1) extension and plate divergence in the Red Sea and Gulf of Aden oceans; (2) strike–slip transform along the Dead Sea and Owen Faults; and (3) compression due to plate collision along the Zagros fold belt range (Fig. 1). These tectonic processes acting at the plate boundaries are accompanied by appreciable seismic activity. Furthermore, recent studies have highlighted the occurrence of small earthquakes within the plate (Al-Shijbi et al. 2019). However, the interior of the Arabian Plate is largely aseismic, suggesting rigid plate motion at lithospheric stress levels which are not sufficient to activate pre-existing faults in the interior of the plate (Viltres et al. 2022). Thus, analysis and inversion of seismic focal mechanisms for present-day in situ stress is not possible. The eastern region of the plate is characterized by substantial horizontal compression, as evidenced by some S_Hmaxin situ stress indicators estimated in boreholes (Ahmed et al. 2007, 2014; Ameen et al. 2010; Ameen 2014), which impact wellbore integrity, hydrocarbon production and well-stimulation processes. Unfortunately, the WSM database does not contain any present-day stress indicators (Rajabi et al. 2014), even though thousands of oil and gas wells have been drilled in eastern Saudi Arabia, Iraq, Kuwait, SW Iran, Abu Dhabi and Oman.

We use numerical simulations to deepen our understanding of the present-day in situ stress state in the Arabian Plate. Specifically, we predict S_Hmax azimuths for the whole plate using numerical finite element modelling tools together with a plate velocity model. A simplified, yet accurate, mesh adequately discretizes the brittle lithosphere. The model contains the peninsula's tectonically active boundaries, gravity-induced stresses, other prominent structural elements, and enough lithological layers to capture mechanical stratigraphies such as the sedimentary cover, the crystalline basement and the upper mantle. We use a 3D unstructured mesh (with finer resolution at critical locations) and solve the linear momentum equation. Boundary conditions are imposed based on the GPS velocity model proposed by Viltres et al. (2022). To model the initial stress state, we follow a two-step approach by first applying gravity followed by lateral displacement, as explained in detail in Hergert et al. (2015).

For simplicity, we provide simulation results at 4 km depth, which enable us to capture the sedimentary layers in the eastern part of the plate where most of the petroleum activities occur. Thus, we provide a possibility to verify our model against any proprietary data which may later become available in the public domain. However, comparing the S_Hmax azimuths recorded in the WSM database (Heidbach et al. 2018), and also listed in a few publications (as discussed in the section ‘Numerical model’), our simulation results agree well with the available data.

The Arabian Plate originated around 30 Ma ago during the late Oligocene by rifting the NE African Plate to form the Red Sea and the Gulf of Aden oceans (e.g. Bosworth et al. 2005). The NE drift of the plate is accommodated in the north and east by the Zagros and Makran fold-and-thrust belts (Fig. 1). Laterally, the plate is flanked by two major transform fault zones: the left-lateral Dead Sea Fault to the west and the right-lateral Owen Fault to the east (De Lamotte and Leturmy 2014). GPS measurements show that the Arabian Plate is slowly moving north relative to Africa, with spreading rates of 16–22 mm yr⁻¹ in the Gulf of Aden and decreasing northward from 18 to 8 mm yr⁻¹ in the Red Sea (Viltres et al. 2022)

The Arabian Plate is a structurally and compositionally complex plate, measuring c. 2600 km north–south and c. 3000 km east–west, and tilted toward the east. The Arabian Shield on the western side is composed of Neoproterozoic basement rocks (1000–540 Ma old) formed by amalgamated magmatic arcs overlain by sedimentary and volcanic basins and intruded by granitoids. Locally, Archean and Paleoproterozoic rocks are structurally intercalated in the Neoproterozoic in the southern and eastern parts of the shield (Johnson and Woldehaimanot 2003). The Arabian Platform on the eastern side predominantly accumulated siliciclastic sediments during the Paleozoic and then transitioned to largely marine carbonates during the Mesozoic and Paleogene. The maximum thickness of this sedimentary basin reaches 10 km (e.g. Murris 1980; Geert et al. 2001; Konert et al. 2001; Ziegler 2001). The onset of the collision between the Arabian and Eurasian Plates at the end of the Oligocene allowed the return of clastic sedimentation at the front of the future Zagros Mountain Range (e.g. Mouthereau 2011), and also in the young rift basins of the Red Sea and the Gulf of Aden (e.g. Leroy et al. 2012).

Owing to the fact that the WSM database contains only very sparse information on in situ stress from the Arabian Plate, we conducted a public-domain literature review to scout for and extract such information. We found 23 publications from which we generated a regional database consisting of 143 in situ stress data points extracted from a variety of publications addressing topics and themes that require information on the present-day in situ stress state. These stress indicators include wellbore breakouts, drilling-induced fractures and recent geologic structures. They reflect the present-day in situ stress orientation. Table 1 provides an overview of the references from which the stress information was extracted, categorized by country, the number of data points and the corresponding WSM quality.

Useful stress indicators include drilling-induced wellbore breakouts and tensile fractures, together with geologic structures (igneous intrusions) as described in detail by Zoback (2010). Geological indicators include structures from the Quaternary age. Tensile fractures, igneous dykes and volcanic alignments grew parallel to the S_Hmax direction at the time when these features formed (Sperner et al. 2003).

The extracted stress data in Table 1 come from circular histograms (rose diagrams), tabular data and azimuths plotted on maps of the respective area. Herein, we list and explain the procedure used to select the data (Fig. 2):

1. Acoustic image log data from wells in North Kuwait fields (Abdali, Bahrah, Dhabi, Mutriba, Raudhatain, Ratqa, Sabriyah and Umm Niqa) show a general trend of N40°E–N50°E (Cretaceous formation) and N20°E–N90°E (Jurassic formation). The S_Hmax azimuths variability from Cretaceous to Jurassic formation is explained by the presence of the Gotnia Salt that, due to its rheologic nature, decouples in situ stress below from that above (Al-Ammar et al. 2012; Perumalla et al. 2014).

2. The stress state in Qatar's North Field located in the Arabian Gulf is constrained using data from three gas wells. The S_Hmax azimuth trends NE–SW (Jonaghani et al. 2019).

3. Borehole image analysis from 16 wells of the Ahmadi carbonate reservoir, located in the Bahrain Field, indicates that the S_Hmax has a NW–SE orientation (Ali et al. 2009).

4. The present-day stress condition of Oman has been extensively studied in the north-central region. The overall trend of the S_Hmax is NE–SW, with slight variations in regions with active faults. In these regions, the in situ stress state is reoriented from NW–SE to east–west. In the Khazzan gas field, the regional S_Hmax trend is NE–SW (Filbrandt et al. 2006).

5. The geomechanical model in the Raoq Field in North Yemen suggests NE–SW S_Hmax azimuths. Furthermore, studies from the Say'um Masila Basin region indicate N55°E S_Hmax azimuths (Barton et al. 2010; Murray and Montgomery 2014).

6. Stress inversion from 277 earthquake focal mechanisms in NW Iran–SE Turkey results in S_Hmax trending N158°E (Mokhoori et al. 2021).

Note that we assessed the extracted in situ stress data's quality according to the rigorous ranking system of the WSP, excluding data ranked as qualities D and E. Consequently, we based our numerical validation on only 143 data points corresponding to qualities A–C.

Figure 3 shows the study area of this work and the plate boundaries. We simplify the fault systems in the Red Sea and Gulf of Aden areas as the spacing between them is below the model resolution (Fig. 3b). The Red Sea and Gulf of Aden areas are mainly related to model discretization that allows accurate and affordable computational simulation runs. Although these fault systems might significantly contribute to the deformation process, their mechanical responses have a local impact only. Consequently, in this work, which considers present-day in situ stress patterns on a plate scale, we do not include fault systems. Furthermore, we model the Arabian Plate as a dry, elastic rock solid. Therefore, we do not account for pore pressure changes (e.g. resulting from petroleum activities in oil or gas fields) to model the overall mechanical behaviour of the Arabian Plate. We acknowledge that there may be such localized effects within certain reservoirs, but for the purposes of this paper, we chose to focus on the plate-scale deformation. Finally, the area encompassed by the boundaries in Figure 3b is extruded vertically to create the 3D model volume.

To reduce the Arabian Plate geological complexity, we build a 3D watertight subsurface model that includes the oceanic domains (i.e. the Red Sea, the Gulf of Aden and the Arabian Sea oceans). It is composed of three main compartments: the sedimentary cover, the crystalline basement and the upper mantle. In Figure 3a, we use the freely available GEBCO grid, which provides the topography and the bathymetry surfaces of the model (Kapoor 1981). Moreover, we use updated versions of the top crystalline basement and crust–mantle (Moho) discontinuity depth maps. To build crystalline basement and the Moho discontinuity depth maps (Fig 3c, d), we compiled isopach and isobath maps, geological sections, seismic refraction profiles, seismic stations and proprietary depth-migrated seismic reflection profiles available from the literature (Kaban et al. 2016, 2019). These data were geo-referenced, digitized and interpolated using QGIS 3.22 software and Move 2020.1 software: Arabian Platform, Geert et al. (2001), Kaviani et al. (2020); Gulf of Oman, Ali et al. (2020), Ninkabou et al. (2021), Tarapoanca et al. (2010), Weidle et al. (2022); Arabian Sea: Exxon (1985); Gulf of Aden, Autin et al. (2008), Cochran (1981), d'Acremont et al. (2005), Gillard et al. (2021), Laske et al. (2013), Leroy et al. (2012), Nonn et al. (2017, 2019), Watremez et al. (2011); Red Sea, Al-Damegh et al. (2005), Blanchette et al. (2018), Davison et al. (1998), Doornenbal et al. (1991), Gaulier et al. (1988), Hansen et al. (2007), Makris et al. (1991). In line with Goteti et al. (2021), we adopted a model configuration with the bottom layer situated at a depth of 200 km, which is positioned beneath the isostatic compensation level for the Arabian Plate. This choice ensures that mechanical loading resulting from topographic and bathymetric changes becomes negligible at the fixed bottom layer.

In order to capture the Arabian Plate subsurface geometry with sufficient resolution in our numerical model, we use a finite element approach to solve the force balance partial differential equation (see equation A8). The finite element method has been extensively used to model plate-scale deformation, providing accurate numerical results for complex geometries (Buchmann 2008; van-der-Zee et al. 2011; Reiter and Heidbach 2014; Ahlers et al. 2021). Thus, we discretize the subsurface model, as illustrated in Figure 4a, using tetrahedral finite elements. Tetrahedral finite element meshes allow modelling of complex geometries while having affordable computational runs. Our finite element mesh is composed of 3.1 million C3D4 Abaqus general-purpose tetrahedral elements. The discretization resolution plays a major role in the accuracy of the finite element solution. When using linear functions for regions with steep elevation gradients, sufficiently small discretization is required to capture these variations (e.g. the transition from the Arabian Shield region with exposed basement to the Arabian Platform sedimentary cover). We refine our finite element mesh using a central node remeshing strategy to capture the thickness variations in the sedimentary basin. Figure 4b portrays our model discretization. This mesh geometry and configuration do not possess several distorted elements, which would result in unrealistic numerical solutions. Figure 4c–e shows the basin thickness variations from three distinct perspectives alongside the basement structure (Fig. 4f).

The mesh quality is assessed using established metrics, including aspect ratio, shape factor and tri-faced corner angles, to ensure the reliability and accuracy of the finite element simulations. The colour-coded representation in Figure 5 highlights the mesh quality across the model. Unshaded areas indicate elements that conform to the defined quality criteria, while areas with deviations from the criteria are depicted in yellow. Notably, only few elements exhibit these deviations, predominantly localized near the corners of the model. These minor deviations from the quality metrics do not compromise the overall integrity of the finite element results. It is important to emphasize that most of the mesh elements meet the quality standards, underscoring the robustness of the mesh for simulating deformation processes. The observed deviations in mesh quality can be attributed to the inherent challenges posed by the transition zones, especially structural heterogeneities. However, we emphasize that the impact of these deviations on the accuracy of the simulations is minimal, given their localized nature.

Figure 6 outlines the key stages of a comprehensive stress modelling workflow. Beginning with the Arabian Plate subsurface model and material properties, it then defines initial stress state and boundary conditions. These inputs feed into a geomechanical model, followed by its discretization into finite elements. The process involves solving force balance equations, culminating in the post-processing of stress results for the stress modelling analysis.

We use the ABAQUS™/Standard implementation to solve the linear static and elastic analysis (Abaqus 2011). Material properties for the linear elastic compartments are listed in Table 2. Without loss of generality, we consider isotropic materials throughout the Arabian Plate. Therefore, the material properties do not change spatially. Furthermore, we adopt a stress notation where compressional stresses are positive and negative stresses indicate tension. To calculate principal stress orientations, the linear static analysis following equation (A8) must be solved using either displacement or force boundary conditions.

An essential aspect of modelling in situ stress azimuths in rock volumes is to account for as many sources of deformation as possible. In this study, we consider the following:

1. Lateral boundary forces between different plates. In the case of Arabian Plate, we find six main plate interactions: Arabian–Anatolian, Arabian–Sinai, Arabian–Nubian, Arabian–Somalian, Arabian–Indian and Arabian–Eurasian (see Fig. 1 to locate these plates). We excluded the Arabian–Danakil interaction as their boundary is not yet well constrained.

2. Gravitational potential energy differences caused by lateral density variations, e.g. transition from the Arabian Shield to the Arabian Platform sedimentary cover.

3. Equilibrium between model geometry and state of initial geostatic stress.

Figure 7a shows the angular rotations $ω˙$ together with their locations. Based on the various velocities around the plate, we divide the plate boundary into six segments, each with a constant Euler vector and location. Figure 7b shows the GPS velocity model proposed by Viltres et al. (2022). Using the linear relation between velocity and time, we calculate the respective displacement boundary conditions. We choose a time span of 100 kyr to estimate cumulative deformation. Assuming a linear relation for velocity changes over 100 kyr is feasible since the plate has not moved substantially over this time period (ArRajehi et al. 2010; Viltres et al. 2022). For the uppermost boundary model, we do not constrain the displacement field since the Earth's surface serves as a free boundary. For the bottom layer, roller boundary conditions in the vertical direction are applied. Vertical roller boundary conditions constrain the model from vertical movement while allowing horizontal motion.

Accounting for the initial stress state in the rock volume is challenging. We assume the rock volume is in quasi-static equilibrium with in situ stresses, gravitational loading, tectonics and stress-assisted thermal deformation, to name a few. In our framework, we estimate the initial stress state using the model geometry, gravity and in situ stresses. Thus, before imposing displacement boundary conditions, we initialize the stress state in the model using the geostatic analysis procedure in ABAQUSTM/Standard. This step involves the roller boundary conditions and gravity to seek equilibrium with the model geometry and initial in situ stresses, preventing unrealistic elastic compaction caused by the gravitational loading. As a result, we can model the kinematic behaviour of the rock volume during a time span while accounting for a pre-stressed solid. Our numerical simulations use a vertical stress gradient of 24.5 MPa km⁻¹ and ratio of horizontal to vertical effective stresses of 0.8 (Goteti et al. 2021). We post-process our simulation to showcase the maximum horizontal stress azimuths. Furthermore, we assume that the vertical stress, S_V, is a principal stress direction.

Figure 8 presents the initial S_Hmax azimuths, located at 4 km depth below sea-level, before imposing boundary conditions. The non-uniform pattern arises from the combined effects of the applied initial stress state and gravitational forces. However, to accurately replicate the current S_Hmax azimuths within the Arabian Plate, additional adjustments involving pushing and pulling at the six boundaries of the domain are necessary to correctly modulate the displacement field. This geostatic procedure is typically employed as a first initial step in geomechanics analysis, during which gravity loads are applied to arrive at an initial stress state. In this step, the loads and initial stresses must equilibrate, resulting in zero deformations, which is the case with our model. Figure 9 shows the initial model displacements after the geostatic step.

In Figure 10, S_Hmax azimuths result from the combination of the initial stress state (initial stress state, gravitational loading and model geometry) with the stress induced by applying displacement boundary conditions (derived from Fig. 7). Figure 11 shows the displacement field magnitude that leads to the observed S_Hmax azimuths. We conducted a sensitivity analysis to pinpoint the model parameters that significantly influence S_Hmax azimuths in the model. Our findings indicate that the most critical parameter is the displacement field applied across the six different boundaries. Figure 12 presents the influence of three displacement fields on S_Hmax azimuths, employing black arrows to represent the displacement field orientations. In Figure 12a, the displacement field orientation is predominantly north–south in the Red Sea, Gulf of Aden, and Owen Fault areas. Conversely, in the Dead Sea Fault region, the orientation mainly follows a NNW–SSW pattern. These orientations result in an east–west orientation of the S_Hmax azimuths across most of the plate. In the Iraq region, there is a notable shift in the S_Hmax azimuths, transitioning from a NE–SW to an NNE–SSW orientation. Close to the boundary of the Zagros fold belt, the S_Hmax azimuths take on a NW–SE orientation. In Figure 12b, the displacement field orientations are primarily north–south in the Red Sea, NNW–SSE in the Gulf of Aden and east–west in the Owen Fault region. Conversely, within the Dead Sea Fault region, the orientation shifts from an east–west to a WNW–ESE pattern, resulting in an overall east–west alignment of the S_Hmax azimuths across most of the tectonic plate. In the Iraq region, there is a notable shift in the S_Hmax azimuths, transitioning from NE–SW to NNE–SSW. Near the boundary of the Zagros fold belt, the S_Hmax azimuths take on a NW–SE orientation. Additionally, a stress rotation is observed in the transition from the Red Sea to the Gulf of Aden, shifting from north–south to NW–SE, and subsequently returning to north–south in the transition from the Gulf of Aden to the Owen Fault region. It is important to note that the earlier displacement field does not accurately represent plate motion, leading to the observed deviations in stress data. In contrast, as depicted in Figures 10 and 11, the displacement field calculated from GPS velocities yields a stress azimuth pattern that closely aligns with the observed data, highlighting a robust agreement between simulation results and real-world observations.

Notably, different applied displacement fields can yield the same S_Hmax azimuth pattern. Therefore, to constrain the model and its displacement boundary conditions further, it is necessary to calibrate stress magnitudes. However, given the lack of published stress magnitude data (and those we found are associated with large uncertainties) across the Arabian region (in particular in Saudi Arabia) such a calibration process is currently not feasible. Previous work on plate-scale deformation modelling constructs the solid geometry, such that the horizontal boundary's faces are perpendicular to the azimuths of the maximum principal compressive stresses (Altmann et al. 2014; Reiter and Heidbach 2014; Ziegler et al. 2016; Ziegler and Heidbach 2021). Consequently, applying perpendicular displacement boundary conditions captures the stress azimuths naturally. In our approach, we model the Arabian Plate using the real plate boundary shape and related crustal deformation (using the most current GPS data) to capture the impact of the plate geometry in the stress field orientations.

GPS velocities used in this work suggest counterclockwise rotation for the Arabian Plate (see e.g. fig. 7b in; Viltres et al. 2022). When translating from linear velocities to displacements, the deformation leads to slightly deviated north–south azimuths of the maximum principal horizontal compressive stresses for most of the Arabian Plate. As previously noted, the resulting pattern considers the initial stress state, which is influenced by gravitational loading as well. We notice that by increasing the number of points to evaluate the linear velocities along the plate boundary, namely the sampling, the simulation results replicate the observed data more accurately, suggesting the dependency on both the imposed amount of deformation and their spatial distribution along the plate boundary. To emphasize again, the sampling also depends upon the resolution of the finite element mesh, since these points must coincide with nodal points. In this study, we use all nodal points for imposing boundary conditions as inferred from their corresponding angular rotations and Euler poles. The data points for stress azimuths used in this work are within the range of 0–12 km. We chose to represent the simulated stress azimuths at 4 km depth since the model within the 0–12 km range does not undergo significant changes. Furthermore, this depth falls within the range accessed by industry for energy extraction processes.

The simulation results show NE–SW S_Hmax azimuths in northeastern Saudi Arabia and Kuwait. NE S_Hmax azimuths dominate the Arabian Gulf region. The collected data, as depicted in Figure 2b, reveal S_Hmax azimuths that closely align with the calculated S_Hmax azimuths in Figure 10. It is worth noting that minor deviations between the calculated S_Hmax azimuths and the observed data fall within the range of uncertainties associated with breakouts and drilling-induced fractures. Farther north, from the Dead Sea transform area to the Arabian­–Anatolian boundary in Iraq's western area, S_Hmax azimuths rotate counterclockwise from NNW–SSE to NW–SE. We observed the same pattern between the calculated S_Hmax azimuths in Figure 10 and the collected S_Hmax azimuth data in Figure 2b. This behaviour highlights the influence of the pinned boundary along the collision zone and the displacement field variations along the Arabian–Anatolian boundary. S_Hmax azimuths in the Dead Sea transform region trend NW–SE and NNW–SSE.

Again, a pinned boundary, but also gravitational loading, in the Zagros Mountain Range, results in a S_Hmax azimuths shift. In the southwestern region of the plate, azimuths change due to the high-relief topography in conjunction with large displacement imposed in the eastern part of the Gulf of Aden relative to the southwestern part of the Red Sea boundary (intersection between two different plate boundaries, Arabia–Africa and Arabia–Somalia). Figure 11 illustrates that the applied displacement field exhibits the greatest magnitude within this region. Therefore, the azimuths rotate counterclockwise from ENE–WSW to north–south.

Our simulation results allow the understanding and study of the current S_Hmax azimuths in the Arabian Plate. To a large extent, the model results suggest NE–SW S_Hmax azimuths in northeastern Saudi Arabia and Kuwait, while the Dead Sea transform areas show NW–SE to NNW–SSE azimuths, and the rest of the plate is characterized by predominant north–south S_Hmax azimuths.

The modelling of the current in situ stress azimuths in the Arabian Plate presented in this work uses a simplified geometry of the brittle lithosphere in the sense that the present model does not account for fault systems, sub-basin boundaries, lateral facies changes, and changes in material properties and rheologies. Such discontinuities contribute significantly to the local perturbations of the in situ stress field. Nevertheless, the overall stress field is mostly controlled by gravitational potential energy arising from density variations together with tectonic processes described in the ‘Numerical model' section.

We apply the displacement field perpendicular to the boundary faces, meaning that the imposed displacement field is not a function of depth. The latter is relevant to capture deformation processes at subduction zones and seafloor spreading. Moreover, in our simulations, we do not include the slow creeping motion of the continental lithosphere base caused by convection currents carrying heat from the interior to the Earth's surface.

Our modelling highlights how the stress orientation pattern in the Arabian Plate results from the superposition of gravitational loading and tectonism. We acknowledge that various factors, including mechanical discontinuities (e.g. faulting at various scales), basin morphology, human activities (such as production or/and injection) and salt domains, can markedly affect stress conditions locally. Consequently, strategic and accurate stress measurements are required to calibrate deformation models at a plate-scale level, albeit they are currently unavailable for the Arabian Plate. Despite this limitation, our numerical analysis demonstrates that, for most of the plate, stress patterns are primarily governed by tectonic interactions along plate boundaries and gravitational loading.

Using a finite element model of the Arabian Plate, Goteti et al. (2021) also derived an in situ stress map of the Arabian Plate. Their model results agree in some measures with the available stress data for the Arabian region. As we find in our study, the model geometry plays a significant role when applying boundary conditions inferred from GPS velocities. In this work, we provide a comprehensive geometry of the Arabian Plate, which favours capturing the stress behaviour for essential regions such as Western Iraq, Oman, Kuwait and the Dead Sea transform. We use the most recent GPS velocity model available for the Arabian region, which allows us to constrain the model more accurately (Viltres et al. 2022). The validation process involved in this work compares simulated stress azimuths with stress azimuth information obtained from various sources in the literature, including drilling-induced fractures, volcanic vents and breakout data, as highlighted in the section ‘The present-day stress pattern in the Arabian Plate: compilation of published data to augment the WSM database’. Although not ubiquitous across the modelled region, the data nonetheless validate our modelling framework and make the resulting stress azimuths more robust. However, it is essential to acknowledge the need for additional measurements – not only present-day in situ stress orientation but also magnitudes. Despite the abundance of drilling projects conducted throughout the Arabian region, we found that there was a marked lack of publicly available stress magnitude data. This represents a significant gap in our current understanding. Should such data become available, they not only have the potential to refine our understanding of the current in situ stress and associated models (regionally and locally), but also reduce uncertainties and manage risks in order to optimize resource exploration, exploitation and other subsurface activities.

We present a static present-day in situ stress – S_Hmax azimuth – model of the Arabian Plate to elucidate ongoing tectonic processes and plate boundary interactions. The modelling framework is based on a finite element method in conjunction with a plate velocity model. The finite element method allows us to build a comprehensive 3D representation of the Arabian Plate, focusing on the Arabian Shield, sedimentary basin and upper mantle geometries. Further, we use a recently published plate velocity model to evaluate the current deformation process at plate boundaries. The combination of the finite element model and the plate velocity model facilitates the study of the deformation field in the Arabian Plate at plate and regional scales. Finally, we conduct an extensive literature review on available in situ stress data, providing new data points that complement those already in the WSM Project. The main results are summarized as follows:

1. The plate velocities used in this work constrain the amount of deformation at plate boundaries when imposing accumulative displacement boundary conditions.

2. The Arabian Plate shows appreciable stress azimuths variability caused by topography changes, plate boundary intersections and plate boundary interactions.

3. Across the transition from the Arabian Shield to the basin, the azimuths rotate from north–south to NNE–SSW and north–south to NW–SSE in several regions, particularly those close to the model boundaries.

4. The effect of the Zagros collision zone, modelled as pinned boundary conditions, has a significant impact on the stress field since the azimuths rotate toward the western part of the plate.

5. We found 143 in situ stress azimuth data points extracted from wellbore breakouts, drilling-induced fractures and volcanic vents from a variety of publications. These new data points serve as validation points in this work as well as in future numerical studies.

6. It is crucial for future research to prioritize acquiring and releasing already existing stress azimuth and magnitude data to calibrate and further refine our understanding of the Arabian Plate's present day in situ stress field.

The authors wish to thank King Abdullah University of Science and Technology (KAUST) for supporting this work and are grateful to Prof. Oliver Heidbach and Dr Moritz Ziegler for their fruitful comments and insights. The authors thank Sander De Jong from Baker Hughes for providing feedback on implementing the finite element mesh. The authors also thank Baker Hughes for providing licenses to the subsurface modelling software JewelSuite™.

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

SPC: conceptualization (lead), data curation (lead), formal analysis (lead), investigation (lead), methodology (lead), software (lead), validation (lead), visualization (lead), writing – original draft (lead); AD: data curation (supporting), investigation (supporting), visualization (supporting); GB: investigation (supporting), visualization (supporting), writing – review & editing (supporting); AMA: supervision (supporting), writing – review & editing (supporting); TF: conceptualization (supporting), formal analysis (supporting), investigation (supporting), supervision (lead), writing – review & editing (lead).

This project was funded by Ali Al-Naimi Petroleum Engineering Research Center of King Abdullah University of Science and Technology.

Data sharing is not applicable to this article as no datasets were generated or analysed during the current study.