Despite inequalities in spatial resolution between stochastic geological models and flow simulator models, geostatistical algorithms are used for the characterisation of groundwater systems. From available data to grid-block hydraulic parameters, workflows basically utilise the development of a detailed geostatistical model (morphology and properties) followed by upscaling. This work aims to design and test a two-step methodology encompassing the generation of a high-resolution 3D stochastic geological model and simplification into a low-resolution groundwater layer-type model. First, a high-resolution 3D stochastic model of rock types or hydrofacies (sets of rock types with similar hydraulic characteristics) is generated using an enhanced version of the sequential indicator simulation (SIS) with corrections for local probabilities and for two- and three-point template statistics. In a second step, the high-resolution geological model provided by SIS is optimally simplified into a small set of layers according to a supervised simulated annealing (SA) optimisation procedure and at the end equivalent hydraulic properties are upscaled. Two outcomes are provided by this methodology: (1) a regular 2D mesh of the top and bottom limits of each hydrogeological unit or layer from a conceptual model and (2), for each layer, a 2D grid-block of equivalent hydraulic parameters prepared to be inputted into an aquifer simulator. This methodology was tested for the upper aquifer area of SPEL (Sociedade Portuguesa de Explosivos), an explosives deactivation plant in Seixal municipality, Portugal.