1 1 Introduction Spatial uncertainty is endemic in geospatial data due to the imperfect means of recording, processing, and representing spatial information (Zhang and Goodchild, 2002; Shi et al., 2016). As geospatial data often serve as inputs in models with spatially distributed parameters, such as physically-based flow simulation models, the propagation of spatial data uncertainty to uncertainty in model predictions is a critical requirement in GIScience and related fields (Heuvelink, 1998; Caers, 2011). Although analytical and/or quasi-analytical uncertainty propagation methods have been developed in the literature; see, for example, Şalap-Ayça et al. (2018), Monte Carlo simulation is rather routinely used for uncertainty propagation purposes, as it does not call for, often limiting, assumptions regarding the form of the spatial model itself. In a nutshell, Monte Carlo simulation consists of generating alternative samples (realizations) from the input parameters, evaluating the model response for each of these realizations, and constructing the corresponding distribution of model predictions. The spatial distribution of input parameters is often modeled within a geostatistical framework, and spatial Monte Carlo simulation is performed within the context of geostatistical simulation (Goovaerts, 1997). Any realistic uncertainty analysis, however, calls for the availability of a representative distribution of model outputs, and can become expensive in terms of both time and computer resources in the case of complex models (Helton and Davis, 2002; Caers, 2011). This problem is far more pronounced in earth and environmental sciences applications, where, in hydrogeology for example, three dimensional grids of hydraulic conductivity values are used along with other parameters to simulate flow and transport in porous media (Gutjahr and Bras, 1993; Chilès and Delfiner, 2012). The computational cost associated with classical Monte Carlo spatial uncertainty propagation calls for the development of more efficient geostatistical simulation methods. This paper proposes key modifications to classical geostatistical simulation to render it more efficient in terms of generating more representative attribute realizations that better span the range of possible realizations corresponding to a geostatistical specification. Computing model predictions using a set of fewer, yet representative, input parameter realizations, is illustrated to reproduce model output sampling variability corresponding to a much larger input parameter set, thus reducing significantly the computational cost associated with Monte Carlo based spatial uncertainty propagation. 2 Efficient geostatistical simulation 2.1 Geostatistical simulation In a geostatistical context, the spatial distribution of values of attributes serving as inputs for spatial models is typically conceptualized via a random field; that is, a set of spatially correlated random variables, *()   +, one per location (Goovaerts, 1997), where () denotes a random variable (RV) defined at a location with coordinate vector . Geostatistical simulation aims at generating multiple (a set of ) simulated attribute values at a set of  locations *      +, typically coinciding with the nodes of a regular grid discretizing the study area; i.e., joint realizations from the M respective RVs *(  )      +. Those realizations are often constrained by (or reproduce) N attribute values Efficient geostatistical simulation for spatial uncertainty propagation Stelios Liodakis University of the Aegean University Hill Mytilene, Greece stelioslio@geo.aegean.gr Phaedon Kyriakidis Cyprus University of Technology 2-8 Saripolou Str., 3036 Lemesos, Cyprus phaedon.kyriakidis@cut.ac.cy Petros Gaganis University of the Aegean University Hill, 81100 Mytilene, Greece gaganis@aegean.gr Abstract Spatial uncertainty and error propagation have received significant attention in GIScience over the last two decades. Uncertainty propagation via Monte Carlo simulation using simple random sampling constitutes the method of choice in GIScience and environmental modeling in general. In this spatial context, geostatistical simulation is often employed for generating simulated realizations of attribute values in space, which are then used as model inputs for uncertainty propagation purposes. In the case of complex models with spatially distributed parameters, however, Monte Carlo simulation could become computationally expensive due to the need for repeated model evaluations; e.g., models involving detailed analytical solutions over large (often three dimensional) discretization grids. A novel simulation method is proposed for generating representative spatial attribute realizations from a geostatistical model, in that these realizations span in a better way (being more dissimilar than what is dictated by change) the range of possible values corresponding to that geostatistical model. It is demonstrated via a synthetic case study, that the simulations produced by the proposed method exhibit much smaller sampling variability, hence better reproduce the statistics of the geostatistical model. In terms of uncertainty propagation, such a closer reproduction of target statistics is shown to yield a better reproduction of model output statistics with respect to simple random sampling using the same number of realizations. The need to process fewer (yet more representative) input realizations implies that the proposed method could contribute to the even wider the application of Monte Carlo-based uncertainty propagation in practice. Keywords: Monte Carlo simulation, random field models, conditional simulation, Latin hypercube sampling, environmental modeling.