Abstract An innovative methodology for diagenesis characterization and quantification is presented. It includes different geostatistical modeling workflows applied to a partially dolomitized carbonate platform. The case study consists of a Lower Cretaceous (upper Aptian) shallow-water carbonate platform from the Basque–Cantabrian basin (northern Spain), in which a widespread burial dolomitization occurs. Previous studies at basin scale suggested that the flow of dolomitizing fluids through the carbonate succession was channeled by regional faults and that subsequently the dolomite distribution was partially controlled by depositional facies and their modifications after early meteoric diagenesis. Here, at reservoir scale, several carbonate facies were differentiated and grouped in five depositional environments. Two depositional sequences corresponding to transgressive–regressive cycles and three stages of the platform evolution were distinguished. The statistical data treatment indicated that the dolomitization is mainly concentrated in the regressive part of the first sequence, corresponding to the second stage of the platform evolution. The most dolomitized environments are the inner platforms and the shoal. Facies from these shallower/proximal depositional environments were more exposed to early meteoric diagenesis, possibly controlling later dolomitization. The total macroscopic porosity is directly proportional to the degree of dolomitization: pores are most abundant in fully dolomitized portions of the succession, particularly in the rudist-bearing and grain-dominated facies. Abundant aragonitic shells (rudists, corals), easily leached or recrystallized during early meteoric diagenesis, could justify the higher moldic porosity in these facies. For geostatistical modeling purposes, several statistical rules were elaborated in order to associate to each depositional environment, in each of the three platform stages, different proportions of dolomitization and related pore abundance. A direct simulation of the distribution of depositional environments, degree of dolomitization, and pore abundance was achieved using a bi-plurigaussian simulation (PGS) algorithm. A nested-PGS algorithm was used to simulate the same parameters independently: dolomite and pore abundance were distributed within each depositional environment, based on the statistical rules previously defined. These simulations allowed three-dimensional (3D) visualization of the original depositional facies and textures affecting the distribution of dolomitization and pore abundance. Modeling using both bi-PGS and nested simulations accounted for the 3D dolomite body extension: the dolomitized succession is thicker in the north and thins toward the south, in agreement with evidence from mapping of the dolomite geobodies.