Assessing the complexity of subsurface facies heterogeneity and patterns, and evaluating their impact on fluid flow dynamics, are two of the primary objectives of reservoir modeling. Geostatistical tools have traditionally been applied to model petroleum reservoirs, significantly contributing to decision-making by supporting well placement strategies aimed at optimizing hydrocarbon recovery under uncertainty. The construction of high-resolution models that integrate well and seismic data, the generation of uncertainty maps, and the creation of heterogeneity models for input into flow simulations are among the most common geostatistical practices in reservoir characterization. In parallel, physics-based approaches have also been employed to evaluate facies heterogeneity. One such approach is Chevron’s novel modeling framework, CompStrat, which simulates high-resolution stratigraphic realizations of deltaic deposits based on physical parameters governing water flow and sediment transport. These digital analogs represent grain size distributions and are characterized by high fidelity and detail, with exhaustive coverage consisting of millions to billions of data points. The models are structured on a uniform grid in the X and Y directions and a non-uniform grid in the Z direction, preserving even the thinnest stratigraphic layers. The application of Geostatistics to CompStrat-driven models is a relatively unexplored area. Bridging Geostatistics with physics-based modeling presents an opportunity to enhance current methodologies and broaden the toolkit available for reservoir modeling. This Thesis explores this integration and investigates its potential benefits. From this integration, new programs and methodologies were developed. The first main contribution of this Thesis is the development of a new GSLIB (Geostatistical Software Library)-like program, sgsim_covtab, for performing sequential Gaussian simulation without the need for variograms. The only input required by sgsim_covtab is a data-driven, positive definite covariance lookup table, enabling the generation of both conditional and unconditional realizations. The exhaustive nature of the grain size CompStrat models made this development feasible by allowing the computation of sufficiently large covariance maps and volumes. Within the sgsim_covtab framework, two additional programs -covamap and varmap_upd- were also developed for computing covariance and variogram maps and volumes. These serve as more efficient alternatives to the conventional varmap GSLIB program. Additionally, a novel methodology is proposed for correcting data-driven covariance maps and volumes to ensure positive definiteness. The efficiency of the new program is evaluated through comparisons with the conventional approach using sgsim and variograms, based on unconditional and conditional simulations involving 2D and 3D reference models extracted from CompStrat realizations. The second major contribution of this Thesis is the development of a novel workflow for constructing n-variate distributions using geologically reasonable synthetic bivariate distributions and grainsize data derived from CompStrat models. A new algorithm is introduced for assigning weights to bivariate distributions to match the marginal distributions of shared variables across different pairs, while preserving the original bivariate structures. These constructed multivariate distributions enable the simulation of additional variables based on the CompStrat grain size models using another novel methodology proposed in this Thesis. Detailed step-by-step examples are provided for constructing both trivariate and quadrivariate distributions, followed by generating realizations of all included variables. Extensive validation is conducted to confirm the accuracy and reliability of all the new tools and methodologies developed in this Thesis.