In Iceland, due to the prevalent surface rock condition, whether it is older bedrock or more recent lava-rock, and the relatively thin and thus easily removable topsoil layer, site effects have largely been assumed to be small and insignificant in different engineering applications. The most commonly used empirical earthquake ground motion models (GMMs) and consequently all past seismic hazard studies in Iceland largely ignore site effects. However, recent comprehensive site effect studies using advanced Bayesian statistical methods (Rahpeyma et al., 2016, 2019, 2022) have revealed that (i) the site effects can vary significantly, even over relatively short distances;(ii) some rock sites exhibit considerable and frequency-dependent site-effects; and (iii) site effects on variable geology are both much stronger and more variable than for rock. In general, 𝑉S30 (ie, the time-averaged shear-wave velocity in the uppermost 30 m) is known as the most common indicator to quantify site effects and used for site classifications. However, systematic estimation of 𝑉S30, mapping the velocity profile with depth, or that of other proxies for the purpose of quantifying site effects has not been carried out in Iceland. To overcome this limitation, we implement a new Bayesian Hierarchical Model (BHM), that enables us to partition the ground motion residuals into various terms associated with the source, propagation path, and station terms (Rahpeyma et al., 2018). We use 83 strong motion data from 6 strike-slip Icelandic earthquakes recorded by the Iceland Strong motion Network (ISMN) and the first small-aperture urban array (ICEARRAY I) strong-motion stations in the South Iceland Seismic Zone (SISZ). The data were recorded by 34 strong motion stations (ie, 25 ISMN and 9 ICEARRAY I stations). The new Bayesian GMM for peak ground acceleration (PGA) and spectral responses of various simple oscillators (PSA at different oscillator periods of 𝑇= 0.01-3.0 sec) can be developed as: log𝑌𝑒𝑠= log𝜇𝑒𝑠+ 𝛿𝐵𝑒+ 𝛿𝑆2𝑆𝑠+ 𝛿𝑊𝑆𝑒𝑠+ 𝛿𝑅𝑒𝑠 for 𝑒= 1,…, 𝑁 & 𝑠= 1,…, 𝑄 (1) where log𝜇𝑒𝑠 is a predictive model (ie, GMM) that provides median ground motion in terms of independent variables. We use the predictive functional form of the latest GMM developed for SISZ proposed by Kowsari et al.(2020) as follows: log𝜇𝑒𝑠= 𝛽1+ 𝛽2𝑀𝑒+ 𝛽3 log10√ 𝑅𝑒𝑠 2+ 𝑍 (𝑀𝑒) 2 𝑍 (𝑀𝑒)= 𝛽4+ 𝛽5 (𝑀𝑒− 𝛽6) 2𝐻 (𝑀𝑒− 𝛽6)