Earthquake recordings on two small-aperture (covering km2 each) strong-motion arrays in Iceland (ICEARRAY I and II) exhibit considerable variations in the spatial distribution of ground-motion amplitudes. To better understand this spatial variability, we use a Bayesian Hierarchical Model (BHM) that incorporates ground motions models (GMMs) for peak ground accelerations (PGA) developed from ground motion databases recorded by the two arrays, respectively. The posterior distributions of the model parameters are then determined using Markov Chain Monte Carlo simulations in the context of Bayesian statistical methods. The BHM allows the partitioning of a GMM into event, station, and event-station terms, which in turn allows the relative contributions of source, path, and site effects to be quantified. The results indicate that site effects can dominate the spatial distribution of ground-motion parameters (e.g., PGA) observed across both ICEARRAY I and II. Although the site conditions across ICEARRAY I have been classified as uniform (i.e., “rock” with a relatively flat topography), station terms contribute to the total variability in the amplitudes of predicted ground motions across the array. In contrast to ICEARRAY I, the variation of the geologic profiles and topography is much greater across ICEARRAY II. As a result, the inter-station variability is shown to contribute up to of the total variability in the amplitudes of predicted ground motions across the array, with the contributions being less constrained for ICEARRAY II than ICEARRAY I due to the relative sizes of the recorded ground motion databases. These results facilitate our understanding of the key factors that affect the variation of seismic ground motions across a relatively small area. Such a detailed microzonation is of great importance for earthquake hazard assessment on a local scale and has practical implications for engineering decision making.