Co-simulation of continuous and categorical variables: application in the Shuiyindong gold deposit modeling

Introduction: Accurate characterization of concealed orebodies in complex geological environments remains a principal challenge in deep mineral exploration. In particular, the strong spatial interdependence between continuous (e.g., grade) and categorical (e.g., lithology) variables greatly increases the complexity of three-dimensional modeling and uncertainty assessment.

Methods: In this study, we implement a unified co-simulation workflow integrating Plurigaussian simulation and the Turning Bands algorithm to co-simulate ore grade and rock type data from drillholes. The approach explicitly captures both auto-correlations and cross-correlations between variables, enabling realistic spatial reproduction of orebody geometry and internal heterogeneity.

Results: The resulting three-dimensional stochastic models effectively represent the spatial anisotropy and variability of mineralization, accurately delineating orebody continuity and its lithological controls.

Discussion: These results not only enhance the understanding of ore distribution at the Shuiyindong deposit but also identify zones with high spatial variability that warrant further deep exploration. Although sparse drilling data introduce local uncertainty, the workflow demonstrates a practical framework for three-dimensional geological modeling and uncertainty management in similarly heterogeneous ore deposits.

1 Introduction

In recent years, advances in three dimensional (3D) geological modeling have strengthened Earth science research and provided practical support for mineral exploration (He et al., 2015). Beyond visualization, modern 3D models are increasingly expected to quantify spatial uncertainty and to integrate heterogeneous information, thereby optimizing target selection and risk management (Li et al., 2018; Chen et al., 2019; Ciazela, 2025; Guo et al., 2025). This transition has been significantly underpinned by methodological innovations, most notably implicit modeling and probabilistic simulation. Implicit modeling facilitates rapid scenario testing by representing geological units as continuous scalar fields (Li et al., 2024; Adoko and Madani, 2024). While probabilistic approaches utilizing ensembles help evaluate decision robustness by quantifying grade uncertainties (Bolarinwa et al., 2024).

A core task in mineral deposit evaluation is the construction of 3D models of grade and rock types, which underpin resource estimation and geological characterization, respectively (Emery, 2007; Pyrcz et al., 2015). These variables are often interdependent: grade varies systematically with lithology and lithological boundaries. Neglecting their spatial interdependence may bias estimates of ore tonnage and average grade; furthermore, it may misrepresent transitions between geological units. (Emery and Silva, 2009). Consequently, a robust framework must simultaneously capture the continuity of individual variables and the dependence among them, while honoring the conditioning data.

Geostatistical frameworks offer various strategies to address the co-simulation of continuous and categorical variables. Traditional cascade workflows, which simulate geology prior to grade, remain common in practice. However, this approach often imposes artificial “hard” boundaries, effectively severing the spatial cross-correlation between variables (Deutsch and Journel, 1998; Emery and Silva, 2009). While alternatives such as Sequential Indicator Simulation or Multiple Point Geostatistics improve categorical modeling, they present their own challenges. These algorithms are often sensitive to sparse conditioning data or training image representativeness, and they frequently struggle to explicitly reproduce joint spatial dependence in complex geological settings (Goovaerts, 1997; Hu and Chugunova, 2008).

Gaussian random field frameworks offer a rigorous mathematical foundation for the co-simulation of mixed variables (Maleki and Emery, 2015; Silva and Deutsch, 2019). Plurigaussian simulation has proven effective for modeling categorical variables. This technique transforms multiple Gaussian random fields into lithological domains via truncation rules, effectively reproducing both simple and complex geological contact relationships (Matheron et al., 1987; Galli et al., 1994; Beucher and Renard, 2016; Veliz et al., 2023). For continuous variables, the Turning bands algorithm represents an efficient method for non-conditional co-simulation, which can be subsequently conditioned using co-kriging, thereby simultaneously honoring available data and capturing the cross-correlations between variables (Emery and Lantuéjoul, 2006).

However, the application of these frameworks to complex Carlin-type environments remains limited. Unlike the stratiform copper deposits analyzed by Maleki and Emery (2015), the Shuiyindong deposit presents challenges due to its discontinuous mineralization and heterogeneous alteration. Specifically, the deposit exhibits composite mineralization controls: while orebodies are predominantly stratabound within specific strata, they are locally complicated by mineralization associated with faults (Zhang et al., 2010; Li et al., 2025). This approach establishes a dependence between grade and lithological units along with non-stationary spatial continuity. Furthermore, the geological history involves multiple episodes of hydrothermal alteration, creating complex enrichment patterns across multiple scales (Su et al., 2008; Tan et al., 2019). This intricate geological setting poses challenges of accurately simulating the spatial distribution of orebodies at Shuiyindong and places strict demands on conventional modeling methodologies.

This study addresses these challenges by applying an integrated Plurigaussian and Turning Bands co-simulation framework to the Shuiyindong gold deposit. We extend the scope of previous research by not only reconstructing the realistic heterogeneity and the interdependence between grade and rock types within the deposit, but also incorporating a systematic uncertainty analysis to delineate high probability targets for deep exploration. Validation against independent modeling results confirms the framework’s proficiency in capturing specific anisotropy and soft-boundary transitions, offering technical insights for the exploration of Carlin-type gold deposits.

The paper has been organized as follows: Section 2 outlines the methodology that has been used to undertake the study. Section 3 provides a description of the Shuiyindong gold deposit and the data processing techniques. Section 4 outlines the results obtained, with the associated discussion outlined in Section 5, and finally the conclusions of the study have been provided in Section 6.

2 Methodology

The methodology employed in this study adopts an established geostatistical framework to address the research problem: the co-simulation of continuous and categorical variables within the Shuiyindong gold deposit. As an application-focused paper, this section details the workflow and justification for the selected standard methods, rather than their theoretical derivation.

2.1 Workflow

The proposed co-simulation workflow (Figure 1) builds on established mixed variable geostatistical frameworks for jointly modeling a continuous grade variable and rock types (Emery and Silva, 2009), and follows recent practical implementations (Valakas and Modis, 2023).

1. Data Preparation and Transformation: The continuous grade data were transformed using a normal score transformation to a standard normal distribution. The categorical rock type data were converted into indicator variables.

2. Spatial Structure Analysis: Experimental direct variograms for the transformed grade and each rock type indicator, as well as cross-variograms between grade and indicators, were computed to quantify spatial auto-correlation and cross-correlation.

3. Co-regionalization Modeling: A linear model of co-regionalization (LMC) was fitted jointly to all direct and cross-variograms to ensure a valid multivariate covariance structure (Wackernagel, 2003). The fitted LMC provides the correlation model used for co-simulation.

4. Geostatistical Simulation: The Plurigaussian simulation method was selected to model the complex geometries of the rock types. The Turning bands algorithm was used for the conditional simulation of the grade.

5. Post-processing and Validation: The simulated Gaussian values were back-transformed to their original scales, generating multiple equiprobable realizations of the deposit. A comparison between these simulation results and observed geological profiles verifies the plausibility of the model.

6. Uncertainty Analysis: The resulting ensemble of realizations was used to perform uncertainty analysis and delineate areas with high mineralization potential.

Figure 1

Figure 1. Flow diagram of co-simulation of continuous and categorical variables.

2.2 Data transformation and spatial modeling

A prerequisite for many geostatistical algorithms is that the data follow a multivariate Gaussian distribution. The raw grade data from the deposit exhibited significant positive skewness, which violates this assumption. Therefore, a normal score transformation was applied to the grade data to obtain a standard Gaussian random field.

The rock types, being categorical variables, were modeled using indicator codes $I_k\left(x\right)$ for each category $K$ (Equation 1). The spatial relationships between the rock types were observed to be non-sequential and exhibited complex contact relationships. Categorical variable modeling is achieved through Plurigaussian simulation, which requires converting rock types into indicator data. The conversion is performed as follows:

$I_k\left(x\right)=\left\{\begin{array}{l}1\hspace{0.33em}if\hspace{0.33em}x\hspace{0.33em}belongs\hspace{0.33em}to\hspace{0.33em}the\hspace{0.33em}kth\hspace{0.33em}rock\hspace{0.33em}type\\ 0\hspace{0.33em}otherwise\end{array}\right.(1)$

To capture the spatial continuity of all variables and their interdependence, a Linear Model of Co-regionalization was constructed (Wackernagel, 2003). This matrix model ensures that the variance-covariance structure defined by all direct and cross-variograms is mathematically valid and provides the necessary input for co-simulation.

2.3 Co-simulation for grade and rock type

Given the complex geological contacts, the Plurigaussian simulation was chosen for modeling the rock types (Armstrong et al., 2011; Dowd et al., 2003). This method utilizes a set of underlying Gaussian random fields and a truncation rule, derived from the observed geological relationships and global proportions, to define the spatial distribution of the categorical variables. For the grade variable, the Turning bands simulation was selected (Paravarzar et al., 2015). This method is computationally efficient and robust for simulating large 3D stationary Gaussian fields (Maleki and Emery, 2015).

Crucially, these algorithms were implemented within a co-simulation framework (Emery, 2008). The Turning bands Simulation algorithm was used to generate all Gaussian random fields (both $Y_0$ for grade and the underlying fields $\left(Y_1,Y_2,...,Y_n\right)$ for Plurigaussian simulation) simultaneously, conditional to the co-regionalization model. This approach ensures that the simulated grade values honor the simulated lithological boundaries, and that both variables respect the spatial auto-correlations and cross-correlations defined in the co-regionalization model. Following co-simulation, the Gaussian fields were back-transformed to produce the final realizations of grade and rock types. This ensemble of realizations provides a robust basis for quantifying resource uncertainty and identifying exploration targets.

3 Study area and dataset

3.1 Geological characteristics of the Shuiyindong gold deposit

The Shuiyindong gold deposit is located in southwestern Guizhou Province, within the Youjiang Basin at the southwestern margin of the South China Block (Figure 2a). The basin has experienced multiple tectonic events, which are considered fundamental to the regional metallogenic evolution (Wang and Groves, 2018). Two major metallogenic episodes at approximately 210 Ma and 140 Ma have been proposed for Carlin-type gold mineralization in the basin, corresponding to the late Triassic Indochina orogeny and the early Cretaceous subduction of the paleo-Pacific plate, respectively (Gao et al., 2022; Jin et al., 2021). Carlin-type gold deposits in this region are predominantly hosted in Permian–Triassic sedimentary rocks and are commonly classified into strata-bound and fault-bound orebodies (Hu et al., 2024). Fault-bound orebodies typically occur in calcareous siltstones and silty mudstones and controlled by high-angle reverse faults. In contrast, strata-bound mineralization mainly occurs in limestones or along unconformities, controlled by anticlines or domes.

Figure 2

Figure 2. (a) Regional geological map of the Shuiyindong gold deposit (modified from Hu et al., 2024); (b) Profile of exploration line 7 in Shuiyindong Gold deposit (modified from Liu et al., 2017)

In the Shuiyindong gold deposit, the exposed strata are Permian to Triassic and include from older to younger, the Middle Permian Maokou Formation (P₂m), the Longtan Formation (P₃l), the Changxing Formation (P₃c), the Dalong Formation (P₃d), and the Triassic Yelang Formation (T₁y) (Figure 2b). Mineralization is closely related to carbonate-bearing stratigraphic units and their boundaries (Liu et al., 2017a). In particular, an important ore-hosting architecture is developed near the karst unconformity between the Maokou Formation (P₂m) and the Longtan Formation (P₃l), where tectonic deformation and hydrothermal alteration formed a structural alteration body (Zhang et al., 2010). Within the study area, the dominant orebody is a silicified breccia–claystone type hosted in the P₃l–P₂m formation, and its geometry is broadly consistent with the unconformity surface (Li S. et al., 2025). In addition, strata-bound orebodies occur within the Longtan Formation, especially in bioclastic limestone interbeds, where the orebody geometry is generally consistent with the host lithology (Li et al., 2023).

The primary ore-controlling structures in Shuiyindong are folds and faults. The dominant fold is the E–W trending Huijiabao Anticline, which exerts a strong control on orebody distribution (Kang et al., 2024). Orebodies are mainly located around the anticline core and on the limbs within an approximately 800 m corridor from the core toward the flank, and their attitudes are broadly concordant with bedding (Li et al., 2025). Along strike, orebodies can exhibit undulating trends with an overall eastward plunge and may occur as vertically stacked and overlapping bodies in 3D space. Ore-bearing faults trend approximately E–W, are distributed along and generally parallel to the anticline axis, and locally provide additional structural preparation for mineralization (Liu et al., 2017a; Liu et al., 2017b).

3.2 Dataset

This study collected data from 97 boreholes in the Shuiyindong gold deposit. The data comprise downhole Au assays and logged rock types. These data provide the conditioning constraints for the subsequent co-simulation. We filtered the dataset to match the target mineralized interval. Available logs interpret samples above elevation 1,400 m as non-mineralized, so we excluded the 1,400–1,500 m interval. We retained samples within elevations of 1,080–1,400 m to capture the main stratabound mineralization. The final filtered dataset contains 13,288 samples.

We restricted the simulation domain to the sedimentary sequences and intentionally excluded the structural alteration zones to account for the dual style mineralization characteristic of the Shuiyindong deposit. As documented by Su et al. (2009), this deposit comprises two distinct mineralization regimes: a deep-seated structural system within alteration zones and a stratabound disseminated system hosted in the Upper Permian strata. These systems represent contrasting fluid dynamic and geochemical precipitation environments, resulting in disparate grade distributions. Mineralization occurring in altered zones typically exhibits continuity of high grade and discrete geometric boundaries, whereas the style bound to strata is characterized by a diffuse distribution. The geostatistical workflow applied in this study, which incorporates variogram analysis and the LMC, rests upon the assumption of spatial stationarity. Aggregating these two populations into a single domain would introduce significant non-stationarity and statistical bias (Madani and Emery, 2017; Mantoglou and Wilson, 1982). Furthermore, the sparse sampling density within the alteration zones limits the ability to constrain such features. The exclusion of these zones therefore ensures the statistical stability and geological realism of the simulations within the primary stratabound host domain.

3.3 Data preprocessing

3.3.1 Processing of rock type data

The rock type data for the Shuiyindong gold deposit were classified into five categories for analysis: (1) non-mineralized rocks (siltstone and limestone), (2) bioclastic rocks, (3) calcareous sandstone, (4) breccia, and (5) claystone. These categories were denoted as lithologies 1–5 for ease of reference (Table 1).

Table 1

Table 1. Basic statistics of grade data contained in different rock types.

3.3.2 Gaussian transformation of grade data

Co-simulation methodology relies on a Gaussian random field assumption. Given that grade distributions often vary significantly between lithological units, grade data within each stationary geological domain (rock type) were transformed separately. This transformation was executed using the nscore function in Geostatistical Software Library (Deutsch and Journel, 1998) to conform each domain’s data to a standard Gaussian distribution, ensuring consistency with the assumptions for spatial variability analysis.

3.4 Determination of rock boundaries

Accurately characterizing rock boundaries as either “hard” (marked by abrupt grade transitions across contacts) or “soft” (characterized by gradual grade transitions) is critical for effective geological modeling. In the case of hard boundaries, grades can be modeled independently within each lithology, assuming no spatial cross-correlation between different rock types. Conversely, for soft boundaries, it is necessary to account for spatial correlations of grades across rock type contacts. Properly classifying these boundaries thus directly influences the modeling approach: hard boundaries imply independent grade distributions among lithologies, while soft boundaries suggest significant spatial dependencies that must be incorporated into the model (Emery and Silva, 2009).

To characterize the behavior of grades across lithological interfaces, two distinct analytical techniques were employed. First, contact analysis was performed, following the methodology of Rossi and Deutsch (2013), to examine the trend of average grade as they approached and crossed the interfaces. Second, this analysis was supplemented by lag scatter plots to assess the spatial correlation of grades between adjacent rock types (Emery and Silva, 2009; Glacken and Snowden, 2001). This technique uses sample pairs selected from different lithologies on opposite sides of a contact, plotted against each other at varying lag distances. A boundary is inferred to be ‘soft’ when the scatter plot exhibits a high correlation coefficient and its center of gravity aligns closely with the $y=x$. This alignment indicates a significant spatial dependency and a smooth grade transition across the contact (Maleki and Emery, 2020). Conversely, a ‘hard’ boundary is indicated by a low correlation coefficient, where the plot is disorganized or non-linear, suggesting the grades are spatially independent across the interface.

The results from both analyses were consistent, identifying the boundaries between rock type 1 and rock type 5, as well as between rock type 3 and rock type 5, as soft (Figure 3). This finding confirms a significant spatial dependency of grades across these interfaces, thereby underscoring the necessity of co-simulation to honor this cross-boundary correlation.

Figure 3

Figure 3. (a,b) Scatterplot of lag of boundary properties; (c,d) Mean grade graphs analysis between rock types.

3.5 Rock type modeling

First, vertical proportionality curves were used to capture the vertical variations in rock proportions within the study area. These curves represent the average areal proportion of each rock type calculated within horizontal strata relative to a common reference datum. They were graphically represented to illustrate the changing proportions of distinct rock types with elevation (Figure 4a).

Figure 4

Figure 4. (a) Vertical scale curve of rock types. (b) Truncation rules in study area.

Second, we determined the lithologic contact rules of the study area based on the drilling data. Each subdomain $\left(D_{1}to{D}_n\right)$ is assumed to be rectangular, with contact relationships graphically represented. Borehole lithological data revealed the following contact relationships: (1) Rock 1: Contacts 2, 3, 4, and 5; (2) Rock 2: Contacts 1, 3, 4, and 5; (3) Rock 3: Contacts 1, 2, 4, and 5; (4) Rock 4: Contacts 1, 2, and 3. Given this complexity, a 3D truncation rule was adopted, represented as [5, 1, 1, 5, 1, 1, 1, 5, 2, 2, 5, 3, 4]. The first Gaussian random field was truncated by two thresholds, while one threshold was used for the second and third fields (Figure 4b).

Third, determine cutoff threshold based on relative proportions of rock types. The proportions were calculated as follows: Rock 1: 55.8%; Rock 2: 7.5%; Rock 3: 4.6%; Rock 4: 1.9%; and Rock 5: 30.2%. Figure 4 depicts the relationship between the rock types of indicator data and truncation thresholds. The truncation thresholds for the study area were calculated as −0.5417, 0.7091, 0.1517, and 0.8426 (Equation 2).

$I_1\left(x\right)=\left\{\begin{array}{l}1\hspace{0.33em}{Y}_1\left(x\right)>t_1\&Y_3\left(x\right)\le t_3\\ 0otherwise\end{array}\right.,$

$I_2=\left\{\begin{array}{l}1\hspace{0.33em}{Y}_1\left(x\right)>t_1\&Y_2\left(x\right)\le t_2\&Y_3\left(x\right)>t_3\\ 0otherwise\end{array}\right.,I_3=\left\{\begin{array}{l}1\hspace{0.33em}{Y}_1\left(x\right)\le {t_1}^{\prime }\&Y_2\left(x\right)>t_2\&Y_3\left(x\right)>t_3\\ 0otherwise\end{array}\right.,(2)$

$I_4=\left\{\begin{array}{l}1\hspace{0.33em}{Y}_1\left(x\right)>{t_1}^{\prime }\&Y_2\left(x\right)>t_2\&Y_3\left(x\right)>t_3\\ 0otherwise\end{array}\right.,$

$I_5=\left\{\begin{array}{c}1\hspace{0.33em}{Y}_1\left(x\right)\le t_1\\ 0\hspace{0.33em}otherwise\end{array}\right..$

For the rock type data, cross-variograms were fitted iteratively to describe the relationships between the lithologic indicator variables and the Gaussian random fields $\left(Y_0,Y_1,Y_2,Y_3\right)$. The first and second spherical structures of $Y_0$ were assigned to $Y_1$, the third spherical structure to $Y_2$, and the fourth and fifth exponential structures to $Y_3$. Nested variogram structures were determined, and goodness of fit was evaluated using the weighted sum of squares (WSS) by trial and error. A smaller WSS indicates a better fit, and a final value of 0.00343 indicates satisfactory results calculated from (Equation 3).

${\gamma }_{11}=0.06gau\left(30,30,8\right)m+0.04gau\left(300,300,300\right)m+0.9sph\left(800,800,245\right)m,$

${\gamma }_{22}=0.2sph\left(50,50,50\right)m+0.8sph\left(100,100,100\right)m,(3)$

${\gamma }_{33}=0.15\mathit{exp}\left(200,200,200\right)m+0.85\mathit{exp}\left(500,500,500\right)m.$

3.6 Grade modeling

The study area was approximately rectangular, measuring about 1,200 m in the east–west direction, 700 m in width, and 300 m in vertical height. Drill holes were spaced approximately 50–60 m apart, and the ore body stroke east–west with a dip either to the south or north at an angle of 5°–10°. The variogram was calculated based on these parameters, and the parameter settings are listed in Table 2. The horizontal and vertical gold grade variograms were computed using GSLIB (Figure 5), and the fitted variogram model consisted of five nested structures, effectively capturing the geometry and spatial anisotropy of the orebody (Equation 4).

$\begin{array}{l}gamm{a}_{00}=0.3gau\left(30,30,8\right)m+0.26sph\left(800,800,245\right)m\\ +0.1sph\left(100,100,100\right)m\\ +0.3\exp \left(200,200,200\right)m+0.04\exp \left(500,500,200\right)m\end{array}(4)$

Table 2

Table 2. Variation function setting parameters.

Figure 5

Figure 5. (a) Fit of the horizontal grade variogram to the nesting structure. (b) Fit of the vertical grade variogram to the nesting structure.

After defining the variogram structures, the spatial cross-correlations between the Gaussian random fields were modeled using a co-regionalization matrix. In this matrix, the diagonal elements represented the direct variogram coefficients of each Gaussian random field, whereas the non-diagonal elements captured the cross-variogram coefficients. The co-regionalization matrix was constrained to be positive and semi-definite, requiring all eigenvalues to be non-negative to ensure mathematical validity. The five off-diagonal parameters were determined through trial and error fitting of the rock types and grade cross-variogram (Figure 6), and the results are presented in Table 3, with specific relationships defined as follows: a and b for $Y_0$ and $Y_1$; c for $Y_0$ and $Y_2$; and d and e for $Y_0$ and $Y_3$.

Figure 6

Figure 6. (a–e) Cross-variogram models fitted between ore grade and rock types 1–5.

Table 3

Table 3. Co-regionalization model parameters.

4 Results

4.1 Co-simulation results

Twelve simulation results were produced. For clarity and focus, Figure 7 presents the realization that most closely resembles the exploration profiles, alongside the mean model for grade and the modal model for rock types. The simulations revealed gradual or relatively smooth transitions in gold grade between rock types. Results from multiple stochastic simulations indicate significant spatial variability in gold grade within the deposit. While uncertainty is inherent in this method, each realization exhibits a distinct distribution of grade values.

Figure 7

Figure 7. (a) One modeling result of the grade. (b) One modeling result of the rock type. (c) Mean model of grade. (d) Modal model of rock types.

4.2 Interpretation of the results

Variograms of the simulated gold grade were compared with those derived from borehole data. Figures 8, 9 present the variogram of the drill holes and simulated grade data, respectively. The variograms in the horizontal and vertical directions indicate the preservation of the spatial properties of the original data in the simulation. Table 4 compares the rock type proportions of the original data with the ensemble mean proportions derived from all simulated realizations. The statistics demonstrate that the simulation reproduces the lithological distribution characteristics of the original dataset.

Figure 8

Figure 8. (a) Horizontal variogram of drill hole grade. (b) Vertical variogram of drill hole grade.

Figure 9

Figure 9. (a) Horizontal variogram of the simulated grade results. (b) Vertical variogram of the simulated grade results.

Table 4

Table 4. Comparison of simulated and original rock type proportions.

Furthermore, the simulation results were overlaid with the profile map of the mine area to verify the correlation between the high-grade regions identified in the simulation and the ore bodies depicted in the profile map, thus validating the accuracy of the simulation outcomes. The simulation results were compared with the exploration line profiles to analyze the causes of any inconsistencies in high-grade areas. Two exploration lines, 7 and 8, were selected for the overlay comparison analysis.

The simulation results within the orebody area of the profile also reflect higher grade values, as illustrated in the overlay map of exploration line 7 (Figure 10a). However, additional high-grade zones were observed in the P₃l³ and P₃l² strata between drill holes ZK708 and ZK701, which is inconsistent with the ore body profile. This discrepancy can be attributed to the larger spacing between ZK708 and ZK701 as compared with the other drill holes, resulting in weaker constraints in the simulation. In drill hole ZK709, the stratigraphic boundary between P₃l¹ and P₃l² yielded high-grade simulation results. This was because of the limited sample data in drill hole ZK709, which reduced the ability to accurately represent the grade distribution. In the overlap map of exploration line 8 (Figure 10b), the simulation results within the ore body area of the profile also exhibited higher grade values. Some simulation outcomes reflected high-grade zones in the stratigraphy of drill holes ZK813, ZK817, ZK825, and P₃l³, which were attributed to the limited sampling data available for this region.

Figure 10

Figure 10. (a) Comparison of the ore body and co-simulation result for exploration line 7. (b) Comparison of the ore body and co-simulation result for exploration line 8.

In summary, the simulation results reasonably reproduced the spatial variability of gold grade and rock type proportions in the study area, showing a general correspondence between high grade zones and the delineated ore body. Moreover, the simulation estimated rock type and grade at locations with missing drill hole data, providing insight into the grade distribution and mineralization in unsampled areas.

4.3 Uncertainty of results

To assess the reliability of the simulation results, we employed cardinality and Shannon entropy as practical uncertainty indicators (Rényi, 1961). The computations were efficiently carried out using the LoopUI 1.0 Python library (Pirot et al., 2022). Both metrics produced highly consistent spatial distributions of uncertainty (Figures 11, 12). The results demonstrated generally low uncertainty for the simulated rock types, with elevated levels confined to edge areas characterized by sparse data and geological complexity, thereby underscoring the overall stability of the simulation. A key insight from the entropy analysis (Figure 12b) was the emergence of specific uncertainty hotspots. These are interpreted as a consequence of the high variance in grade values and the dispersed nature of the ore bodies at the Shuiyindong gold deposit, challenges that persist despite the presence of complete drill hole data. Ultimately, this application of uncertainty quantification serves to validate the model’s consistency and provides critical insights for targeting areas where future data collection or model refinement could be most beneficial.

Figure 11

Figure 11. (a) Results of rock type model uncertainty quantification using cardinality indicators. (b) Results of rock type model uncertainty quantification using Shannon entropy. Map represents the top surface of the voxet, WE, NS represent east-west and north-south profile.

Figure 12

Figure 12. (a) Results of grade uncertainty quantification using cardinality indicators. (b) Results of grade uncertainty quantification using Shannon entropy. Map represents the top surface of the voxet, WE, NS represent east-west profile and north-south profile.

4.4 Implications for further mineral exploration

This study analyzed drilling data and exploration profiles to constrain mineralization primarily within the stratigraphic interval from P₃l¹ to P₃l². To assess the validity of the co-simulation results, we built a grade-only base model that simulates gold grade independently without lithological co-simulation, using the same conditioning data and variogram settings. Figure 13 compares the experimental results under identical grade thresholds. The co-simulation model highlights two deep regions with sparse drilling data. However, the correlation between variables indicates a high probability of mineralization in these areas. Two high probability mineralization zones were identified through uncertainty quantification and comparison with a stratigraphically unconstrained gold element interpolation model (Figure 13). Exploration engineering has confirmed the presence of multiple gold orebodies at depths of 800–1,400 m, with proven reserves amounting to 28 t (Li et al., 2025). These findings are expected to guide mineral exploration within the Shuiyindong gold deposit.

Figure 13

Figure 13. (a) Predicted prospecting area for mineralization of co-simulated grade model. (b) Au grade interpolation model for comparison.

5 Discussion

Treating grade and rock type independently is widely used in practice. This approach may reproduce individual spatial continuity, but it often neglects cross-correlation. Consequently, unrealistic grade patterns arise near lithological boundaries. To address this, we compared our proposed co-simulation with independent simulations. The co-simulation reduces boundary artifacts because the cross-variograms constrain grade variability to follow the same structural and stratigraphic directions that control rock type continuity in the dataset.

High grade anomalies in sparsely drilled volumes do not arise from spacing alone, because the model produces them only where the co-simulated rock type realizations and the fitted cross-correlation jointly favor grade values. In other words, the anomalies reflect an extrapolation of coupled geological controls, and future drilling can test them by targeting intersections of favorable host horizons with mapped structural corridors at depth. The co-simulation identifies deep favorable mineralized zones within the P₃l¹ stratigraphic horizon. This spatial constraint relies on the fitted cross-correlation between grade and lithological categories.

Several limitations remain. First, the LMC assumes stationarity within modeling domains, and strong heterogeneity may require locally varying coregionalization or other non-stationary multivariate formulations. Second, the Turning bands simulation is efficient for large 3D grids, but it can introduce artifacts when anisotropy is poorly constrained, so we recommend sensitivity tests on anisotropy and nested structures.

Our current framework couples grade with rock type, but it does not explicitly encode hydrothermal alteration, even though alteration assemblages at Shuiyindong are closely linked to mineralization in many studies. In future work, we intend to incorporate alteration related indicators or continuous variables into the multivariate framework. Furthermore, we will evaluate Intrinsic random functions or Multiple point geostatistics to better characterize non-stationary patterns and complex connectivity in structurally controlled ore systems (Madani and Emery, 2017). This direction aligns with recent integrative exploration workflows that combine geological knowledge, quantitative modeling, and uncertainty-aware targeting in deposit studies (Ciążela et al., 2024; Ciążela, 2025).

6 Conclusion

This study demonstrated the effective application of a joint simulation framework, integrating Plurigaussian and Turning bands techniques, to co-simulate grade and rock type in the complex Shuiyindong gold deposit. A core-regionalization model was successfully implemented to capture the spatial cross-correlations between these variables, an approach that proved crucial for minimizing error propagation and enhancing simulation accuracy. The results show a strong spatial agreement between the simulated models and the known ore body, while systematic uncertainty quantification highlighted areas of high reliability as well as zones requiring further verification. This study integrated co-simulation results with uncertainty characterization to delineate prospective mineralization zones, providing technical support for the deep exploration of the Shuiyindong gold deposit. Future research could explore incorporating soft constraints or alternative algorithms to model alteration zones, thereby enhancing the applicability of the joint simulation framework.

Data availability statement

The datasets presented in this article are not readily available because No right to share the data. Requests to access the datasets should be directed to CWd2FuZ2NoYkBjdWcuZWR1LmNu.

Author contributions

JZ: Conceptualization, Methodology, Writing – original draft. YH: Conceptualization, Methodology, Validation, Writing – review and editing. CW: Conceptualization, Funding acquisition, Project administration, Supervision, Writing – review and editing. JC: Methodology, Project administration, Supervision, Writing – review and editing. JL: Resources, Writing – review and editing.

Funding

The author(s) declared that financial support was received for this work and/or its publication. This study was supported by the National Key R&D Program of China (grant numbers 2022YFF0801202, 2022YFF0801200, and 2024ZD1001205) and the Key Program of the Geological Bureau of Hunan Province (grant number HNGSTP202401).

Conflict of interest

Author YH was employed by Wuhan Zhi Bo Chuang Xiang Science and Technology Co., Ltd.

The remaining author(s) declared that this work was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Generative AI statement

The author(s) declared that generative AI was not used in the creation of this manuscript.

Any alternative text (alt text) provided alongside figures in this article has been generated by Frontiers with the support of artificial intelligence and reasonable efforts have been made to ensure accuracy, including review by the authors wherever possible. If you identify any issues, please contact us.

Publisher’s note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

References

Adoko, C. G., and Madani, N. (2024). Stochastic modeling of geological domains using a truncated Gaussian collocated co-simulation approach. Stoch. Environ. Res. Risk Assess. 38 (5), 2081–2094. doi:10.1007/s00477-024-02670-x

CrossRef Full Text | Google Scholar

Armstrong, M., Galli, A., Beucher, H., Loc’h, G., Renard, D., Doligez, B., et al. (2011). Plurigaussian simulations in geosciences. Springer Science and Business Media.

Google Scholar

Beucher, H., and Renard, D. (2016). Truncated Gaussian and derived methods. Comptes Rendus. Géoscience 348 (7), 510–519. doi:10.1016/j.crte.2015.10.004

CrossRef Full Text | Google Scholar

Bolarinwa, S., Femi, B. O., Obinna, J. O., and Henpy, O. I. (2024). Techniques for improved reservoir characterization using advanced geological modeling in the oil and gas industry. Int. J. Appl. 6 (9), 2060–2088. doi:10.51594/ijarss.v6i9.1542

CrossRef Full Text | Google Scholar

Chen, Q., Liu, G., Ma, X., Zhang, J., and Zhang, X. (2019). Conditional multiple-point geostatistical simulation for unevenly distributed sample data. Stoch. Environ. Res. Risk Assess. 33, 973–987. doi:10.1007/s00477-019-01671-5

CrossRef Full Text | Google Scholar

Ciążela, J. (2025). Editorial for the special issue “Recent Advances in Exploration Geophysics”. Appl. Sci. 15 (3), 1251. doi:10.3390/app15031251

CrossRef Full Text | Google Scholar

Ciazela, J., Junge, M., Helmy, H. M., Chen, L., and Nadeau, O. (2024). Ore formation and critical metal deposits: geological contribution to the clean energy transition. Front. Earth Sci. 12, 1371997. doi:10.3389/feart.2024.1371997

CrossRef Full Text | Google Scholar

Deutsch, C. V., and Journel, A. G. (1998). GSLIB: geostatistical software library and user’s guide. Oxford University Press.

Google Scholar

Dowd, P. A., Pardo-Igúzquiza, E., and Xu, C. (2003). Plurigau: a computer program for simulating spatial facies using the truncated Plurigaussian method. Comput. Geosciences 29 (2), 123–141. doi:10.1016/s0098-3004(02)00070-5

CrossRef Full Text | Google Scholar

Emery, X. (2007). Simulation of geological domains using the Plurigaussian model: new developments and computer programs. Comput. Geosciences 33 (9), 1189–1201. doi:10.1016/j.cageo.2007.01.006

CrossRef Full Text | Google Scholar

Emery, X. (2008). A turning bands program for conditional co-simulation of cross-correlated Gaussian random fields. Comput. Geosciences 34, 1850–1862. doi:10.1016/j.cageo.2007.10.007

CrossRef Full Text | Google Scholar

Emery, X., and Lantuéjoul, C. (2006). TBSIM: a computer program for conditional simulation of three-dimensional Gaussian random fields via the turning bands method. Comput. Geosciences 32 (10), 1615–1628. doi:10.1016/j.cageo.2006.03.001

CrossRef Full Text | Google Scholar

Emery, X., and Silva, D. A. (2009). Conditional co-simulation of continuous and categorical variables for geostatistical applications. Comput. Geosciences 35 (6), 1234–1246. doi:10.1016/j.cageo.2008.07.005

CrossRef Full Text | Google Scholar

Galli, A., Beucher, H., Le Loc’h, G., Doligez, B., and Group, H. (1994). “The pros and cons of the truncated Gaussian method,” in Geostatistical simulations: proceedings of the geostatistical simulation workshop, fontainebleau, France, 27–28 may 1993 (Dordrecht: Springer Netherlands), 217–233.

Google Scholar

Gao, W., Hu, R., Mei, L., Bi, X., Fu, S., Huang, M., et al. (2022). Monitoring the evolution of sulfur isotope and metal concentrations across gold-bearing pyrite of Carlin-type gold deposits in the Youjiang Basin, SW China. Ore Geol. Rev. 147, 104990. doi:10.1016/j.oregeorev.2022.104990

CrossRef Full Text | Google Scholar

Glacken, I. M., and Snowden, D. V. (2001). “Mineral resource estimation,” in Mineral resource and ore reserve estimation—the AusIMM guide to good practice. Editor A. C. Edwards (Melbourne: The Australasian Institute of Mining and Metallurgy), 189–198.

Google Scholar

Goovaerts, P. (1997). Geostatistics for natural resources evaluation. New York: Oxford University Press. doi:10.1093/oso/9780195115383.001.0001

CrossRef Full Text | Google Scholar

Guo, J., Zheng, Y., Liu, Z., Wang, X., Zhang, J., and Zhang, X. (2025). Pattern-based multiple-point geostatistics for 3D automatic geological modeling of borehole data. Nat. Resour. Res. 34 (1), 149–169. doi:10.1007/s11053-024-10405-6

CrossRef Full Text | Google Scholar

He, W., Mo, X., He, Z., White, N., Chen, J., Yang, K., et al. (2015). The geology and mineralogy of the Beiya skarn gold deposit in Yunnan, Southwest China. Econ. Geol. 110 (6), 1625–1641. doi:10.2113/econgeo.110.6.1625

CrossRef Full Text | Google Scholar

Hu, L. Y., and Chugunova, T. (2008). Multiple-point geostatistics for modeling subsurface heterogeneity: a comprehensive review. Water Resources Res. 44 (11). doi:10.1029/2008wr006993

CrossRef Full Text | Google Scholar

Hu, N., Hu, R., Chen, H., Fu, S., Yang, J., Xu, L., et al. (2024). Lithium and oxygen isotopic constraints on the source and evolution of ore-forming fluids: a case study from the Shuiyindong Carlin-type gold deposit, SW China. Miner. Deposit. 59 (2), 313–328. doi:10.1007/s00126-023-01211-w

CrossRef Full Text | Google Scholar

Jin, X. Y., Zhao, J. X., Feng, Y. X., Hofstra, A. H., Deng, X. D., Zhao, X. F., et al. (2021). Calcite U-Pb dating unravels the age and hydrothermal history of the giant Shuiyindong Carlin-type gold deposit in the golden triangle. South China. Econ. Geol. 116 (6), 1253–1265. doi:10.5382/econgeo.4870

CrossRef Full Text | Google Scholar

Kang, H., Liu, Y., Hu, K., and Han, S. (2024). Geochemical evidence of ore-forming processes in the shuiyindong gold deposit of Southwest Guizhou Province, China. ACS Omega 9 (38), 39365–39386. doi:10.1021/acsomega.4c02165

PubMed Abstract | CrossRef Full Text | Google Scholar

Li, N., Song, X., Xiao, K., Li, S., Li, C., and Wang, K. (2018). Part II: a demonstration of integrating multiple-scale 3D modelling into GIS-based prospectivity analysis: a case study of the Huayuan-Malichang district, China. Ore Geol. Rev. 95, 292–305. doi:10.1016/j.oregeorev.2018.02.034

CrossRef Full Text | Google Scholar

Li, R., Tan, Q., Wang, X. L., Sun, X., Yang, T., Xia, Y., et al. (2023). A metasedimentary origin for gold deposits in the Dian-Qian-Gui “Golden Triangle” of Southwest China. Ore Geol. Rev. 159, 105560. doi:10.1016/j.oregeorev.2023.105560

CrossRef Full Text | Google Scholar

Li, M., Chen, C., Liang, H., Han, S., Ren, Q., and Li, H. (2024). Refined implicit characterization of engineering geology with uncertainties: a divide-and-conquer tactic-based approach. Bull. Eng. Geol. Environ. 83 (7), 282. doi:10.1007/s10064-024-03765-z

CrossRef Full Text | Google Scholar

Li, G., Xue, C., and Pi, Q. (2025). Occurrence Characteristics and Enrichment Mechanism of Au in Shuiyindong Gold Deposit in Southwestern Guizhou Province, Southwest China. Geol. J., gj.5244. doi:10.1002/gj.5244

CrossRef Full Text | Google Scholar

Li, S., Yang, C., Xie, Z., Tan, Q., Wang, Z., Chen, F., et al. (2025). Ore-prospecting breakthrough and research advances on the fully concealed super-large Shuiyindong Carlin-type gold deposit in Guizhou Province, China. Ore Geol. Rev. 183, 106674. doi:10.1016/j.oregeorev.2025.106674

CrossRef Full Text | Google Scholar

Liu, J., Xia, T., Fu, Y., Wang, Z., Wu, W., Tan, Q., et al. (2017). Study on the SBT in Southwest Guizhou. Wuhan: China University of Geosciences Press.

Google Scholar

Liu, J., Yang, C., and Wang, Z. (2017). Geological research of Shuiyindong gold deposit in Zhenfeng county, Guizhou. Geol. Surv. China 4 (3), 32–41.

Google Scholar

Madani, N., and Emery, X. (2017). Plurigaussian modeling of geological domains based on the truncation of non-stationary Gaussian random fields. Stoch. Environ. Res. Risk Assessment 31 (4), 893–913. doi:10.1007/s00477-016-1365-9

CrossRef Full Text | Google Scholar

Maleki, M., and Emery, X. (2015). Joint simulation of grade and rock type in a stratabound copper deposit. Math. Geosci. 47, 471–495. doi:10.1007/s11004-014-9556-8

CrossRef Full Text | Google Scholar

Maleki, M., and Emery, X. (2020). Geostatistics in the presence of geological boundaries: exploratory tools for contact analysis. Ore Geol. Rev. 120, 103397. doi:10.1016/j.oregeorev.2020.103397

CrossRef Full Text | Google Scholar

Mantoglou, A., and Wilson, J. L. (1982). The turning bands method for simulation of random fields using line generation by a spectral method. Water Resour. Res. 18 (5), 1379–1394. doi:10.1029/wr018i005p01379

CrossRef Full Text | Google Scholar

Matheron, G., Beucher, H., Fouquet, C., Galli, A., Guerillot, D., and Ravenne, C. (1987). Conditional simulation of the geometry of fluvio-deltaic reservoirs. Dallas: Society of Petroleum Engineers, 27–30.

Google Scholar

Paravarzar, S., Emery, X., and Madani, N. (2015). Comparing sequential Gaussian and turning bands algorithms for cosimulating grades in multi-element deposits. Comptes Rendus. Géoscience 347 (2), 84–93. doi:10.1016/j.crte.2015.05.008

CrossRef Full Text | Google Scholar

Pirot, G., Joshi, R., Giraud, J., Lindsay, M. D., and Jessell, M. W. (2022). loopUI-0.1: indicators to support needs and practices in 3D geological modelling uncertainty quantification. Geosci. Model Dev. 15 (12), 4689–4708. doi:10.5194/gmd-15-4689-2022

CrossRef Full Text | Google Scholar

Pyrcz, M. J., Sech, R. P., Covault, J. A., Willis, B. J., Sylvester, Z., Sun, T., et al. (2015). Stratigraphic rule-based reservoir modeling. Bull. Can. Petroleum Geol. 63 (4), 287–303. doi:10.2113/gscpgbull.63.4.287

CrossRef Full Text | Google Scholar

Rényi, A. (1961). On measures of entropy and information. Proc. Fourth Berkeley Symposium Mathematical Statistics Probabilit 4, 547–562. Available online at: https://projecteuclid.org/euclid.bsmsp/1200512181.

Google Scholar

Rossi, M. E., and Deutsch, C. V. (2013). Mineral resource estimation. Springer Science and Business Media.

Google Scholar

Silva, D., and Deutsch, C. (2019). Multivariate categorical modeling with hierarchical truncated Pluri-Gaussian simulation. Math. Geosci. 51 (5), 527–552. doi:10.1007/s11004-018-09782-5

CrossRef Full Text | Google Scholar

Su, W., Xia, B., Zhang, H., Zhang, X., and Hu, R. (2008). Visible gold in arsenian pyrite at the Shuiyindong Carlin-type gold deposit, Guizhou, China: implications for the environment and processes of ore formation. Ore Geol. Rev. 33 (3-4), 667–679. doi:10.1016/j.oregeorev.2007.10.002

CrossRef Full Text | Google Scholar

Su, W., Hu, R., Xia, B., Xia, Y., and Liu, Y. (2009). Calcite Sm-Nd isochron age of the Shuiyindong Carlin-type gold deposit, Guizhou, China. Chem. Geol. 258 (3-4), 269–274. doi:10.1016/j.chemgeo.2008.10.030

CrossRef Full Text | Google Scholar

Tan, Q., Xia, Y., Xie, Z., Wang, Z., Wei, D., Zhao, Y., et al. (2019). Two hydrothermal events at the Shuiyindong Carlin-Type Gold deposit in Southwestern China: insight from Sm–Nd dating of fluorite and calcite. Minerals 9 (4), 230. doi:10.3390/min9040230

CrossRef Full Text | Google Scholar

Valakas, G., and Modis, K. (2023). GeoSim: an R-package for plurigaussian simulation and Co-simulation between categorical and continuous variables. Appl. Comput. Geosciences 19, 100130. doi:10.1016/j.acags.2023.100130

CrossRef Full Text | Google Scholar

Veliz, V., Maleki, M., Madani, N., Soltani-Mohammadi, S., Mery, N., and Emery, X. (2023). Plurigaussian modeling of non-stationary geological domains to assess geological uncertainty in a porphyry copper deposit. Ore Geol. Rev. 162, 105707. doi:10.1016/j.oregeorev.2023.105707

CrossRef Full Text | Google Scholar

Wackernagel, H. (2003). Multivariate geostatistics. Berlin, Heidelberg: Springer Berlin Heidelberg.

Google Scholar

Wang, Q., and Groves, D. (2018). Carlin-style gold deposits, Youjiang Basin, China: tectono-thermal and structural analogues of the Carlin-type gold deposits, Nevada, USA. Miner. Deposita 53 (7), 909–918. doi:10.1007/s00126-018-0837-x

CrossRef Full Text | Google Scholar

Zhang, Y., Xia, Y., Su, W., Tao, Y., Zhang, X., Liu, J., et al. (2010). Metallogenic model and prognosis of the Shuiyindong super-large strata-bound Carlin-type gold deposit, Southwestern Guizhou province, China. Chin. J. Geochem. 29, 157–166. doi:10.1007/s11631-010-0157-z

CrossRef Full Text | Google Scholar