Integration of Soft Data Into Geostatistical Simulation of Categorical Variables

Uncertain or indirect “soft” data, such as geologic interpretation, driller’s logs, geophysical logs or imaging, offer potential constraints or “soft conditioning” to stochastic models of discrete categorical subsurface variables in hydrogeology such as hydrofacies. Previous bivariate geostatistical simulation algorithms have not fully addressed the impact of data uncertainty in formulation of the (co) kriging equations and the objective function in simulated annealing (or quenching). This paper introduces the geostatistical simulation code tsim-s, which accounts for categorical data uncertainty through a data “hardness” parameter. In generating geostatistical realizations with tsim-s, the uncertainty inherent to soft conditioning is factored into both 1) the data declustering and spatial correlation functions in cokriging and 2) the acceptance probability for change of category in simulated quenching. The degree or sensitivity to which soft data conditions a realization as a function of hardness can be quantified by mapping category probabilities derived from multiple realizations. In addition to point or borehole data, arrays of data (e.g., as derived from a depth-dependency function, probability map, or “prior realization”) can be used as soft conditioning. The tsim-s algorithm provides a theoretically sound and general framework for integrating datasets of variable location, resolution, and uncertainty into geostatistical simulation of categorical variables. A practical example shows how tsim-s is capable of generating a large-scale three-dimensional simulation including curvilinear features.


INTRODUCTION
In many hydrogeological modeling applications, much of the available characterization data for categorical variables, such as lithology, texture, or hydrofacies, are uncertain or less than 100% accurate. In some geostatistical applications, the categorical data such as soil texture are treated as 100% accurate or "hard" data despite the fact these data are uncertain for various reasons (Carle, 1996;Burow et al., 1997;Carle et al., 1998;. In other applications, uncertain or indirect "soft" data such as geophysical imaging are available but found difficult to apply as "soft conditioning" to geostatistical simulations of categorical variables (Falivene et al., 2007;Koch et al., 2014). Categorization of hydrogeological variables often has uncertainty attributable to sample quality, geologic interpretation, or indirectness of measurement (e.g., geophysical logs, cone penetrometer data). For example, so-called "driller's logs" or lithologic descriptions by well drillers based on interpretations of drilling conditions, cuttings, and limited core samples are, understandably, uncertain (Oatfield and Czarnecki, 1989;Smith, 2002;Dumedah and Schuurman, 2008;Arihood, 2009;Tsai and Elshall, 2013). Yet driller's logs may provide the most detailed characterization information available for many local site to basin-scale hydrogeological modeling applications (Ezzedine et al., 1999;Carle et al., 2006;Fleckenstein et al., 2006;Elshall et al., 2013). Alternatively, hydrogeological categories can be inferred from geophysical measurements or imaging, but the resulting inferences are inherently uncertain. Paradis et al. (2015) found that even after application of sophisticated machine learning techniques to cone penetrometer test and soil moisture and resistivity probe data, classification errors persist in hydrofacies identification. The resolution of geophysical images varies as a result of petrophysical relationships and, therefore, provides indirect contraints on lithology, texture, or hydrofacies, as further discussed in application of three-dimensional resistivity maps derived from airborne electromagnetic surveys to characterization of spatial occurrences of sand-clay textures (Koch et al., 2014;Hoyer et al., 2015) and electrical resistance tomographs to determine spatial distributions of alluvial or fluvial hydrofacies Carle and Ramirez, 1999;Hermans and Irving, 2017). It is a well-accepted fact in the Earth sciences that uncertainty is common to subsurface data.
A general framework for integrating soft data into categorical geostatistical simulation should be useful to stochastic subsurface characterization and inversion of discrete heterogeneity within hydrogeologic systems. One, two, or three-dimensional (3-D) spatial information derived from geophysical imaging or logging (e.g., seismic, electrical) or hydrogeologic interpretation (e.g., cross-sections or calibrated flow models) could be treated as a type of soft data available for conditioning geostatistical simulation of categorical variables. In stochastic inversion, a general framework for integrating prior information could be used, for example, to selectively manipulate discrete heterogeneity structure, such as the interconnectivity of permeable units or the continuity of impermeable units (Carle and Ramirez, 1999;Aines et al., 2002;Wainwright et al., 2014). A "Markov-Bayes" approach has been proposed to transform the soft data to prior probability distributions (Zhu, 1991;Deutsch and Journel, 1998), but this approach has been deemed intractable because of nonlinear relationships and large volumetric scale of the soft data (Deutsch and Wen, 2000).
In this paper, a simple theoretical framework is developed for incorporating uncertain, indirect, or soft categorical data into categorical geostatistical simulation. Geostatistical realizations will honor or be conditional to both hard and soft data. The theory considers that soft data should not be treated the same as hard data in formulating (co)kriging equations and objective functions in simulated annealing (or quenching) as implemented in the original categorical stochastic simulation codes using bivariate spatial statistics such as isim3d (Gomez-Hernandez and Srivastava, 1990), tsim (Carle 1996;Carle et al., 1998), sisim and anneal (Deutsch and Journel, 1998), and iksim (Ying, 2000).
A new version of tsim, called tsim-s, has been coded to enable incorporation of soft categorical data or prior information. The tsim-s algorithm was originally conceived to enable tsim to produce and perturb stochastic realizations for Monte Carlo Markov chain inversion Carle, 2003;Glaser et al., 2004). The new capabilities in tsim-s have more general applicability and flexibility to handle conditioning data of variable quality or uncertainty. tsim-s will be made available by request as the open-source fortran code tsim has been distributed in the past. The development of the equations necessary for implementation of the tsim-s algorithm are included in this paper to fully document the methods and to facilitate coding of the tsim algorithms in higher-level languages such as R (Sartore, 2013;Sartore et al., 2016). As will be seen in the equations, the computational overhead for tsim-s is not signifcantly different from tsim because the only modifications are to the entries in the cokriging matrices and the parameters of the quenching objective function.
This paper provides the theory behind the tsim-s algorithms and results from example applications. The paper first reviews transition probability-based indicator geostatistical theory implemented in the tsim algorithms. Next, the paper derives the equations used for implementing the new soft data capabilities to account for uncertainty in categorical variables using the "hardness" concept previously introduced for continuous variables (Deutsch and Wen, 2000). Cokriging equations and simulated quenching objective functions are re-ormulated to account for soft data using the hardness concept. Example applications are given for hard and soft conditioning at boreholes and from prior conditioning by an array of data such as another realization or a depth-dependent uncertainty function. Results of a 3-D simulation by Carle et al. (2006) are included to demonstrate that tsim-s can be used to produce large-scale stochastic realizations useful for investigation of flow and transport processes in hydrogeologic systems. These examples are provided to show the flexibility of the algorithm and by no means cover all potential applications, a topic that is beyond the scope of this paper. Further discussion is added to clarify the capabilities and limitations of the tsim and tsim-s algorithms in relationship to variogram-based and multi-point geostatistical methods and the current understanding of how the simulation algorithms can and should be implemented in a geological context.

Transition Probability-Based Indicator Geostatistics
Transition probability-based indicator geostatistics is a categorical geostatistical approach where the transition probability bivariate statistic is used to analyze spatial variability and formulate cokriging equations (Carle and Fogg, 1996). The transition probability approach enables consideration of spatial cross correlations (e.g., how different facies tend to locate in space relative to each other) and facilitates a Markov chain modeling framework that can be linked to geologic interpretation (Krumbein and Dacey, 1969;Miall, 1973;Carle and Fogg, 1997). Indicator variogram-based geostatistical approaches do not fully consider spatial cross-correlations and rely on data-intensive empirical curve fitting for model development (Deutsch and Journel, 1992;Deutsch and Journel, 1998;Goovaerts, 1996;Goovaerts, 1997).
In a categorical geostatistical approach, an indicator variable is defined with respect to mutually exclusive or discrete categorical variables (e.g., lithofacies, hydrofacies) by where x is location, and K is the number of categories. The probability that category k occurs at x is equivalent to the expected value of the indicator variable: In transition probability-based indicator geostatistics (Carle and Fogg, 1996), the transition probability bivariate spatial statistic is used to quantitatively describe spatial variability of the discrete categorical variables, which we will generally refer to as "facies." Assuming second-order stationarity, the transition probability t jk (h) is defined as a conditional probability that depends on a lag separation vector h by t jk (h) Pr k occurs at x + h j occurs at x .
Applying Bayes theorem and Eqs. 1 and 3 is formulated with respect to indicator variables by t jk (h) Pr j occurs at x and k occurs at x + h Pr j occurs at x Pr I j (x) 1 and I k (x + h) 1 The transition probability entries, t jk (h), form the transition probability matrix, T(h), as Other bivariate statistics, such as the indicator (cross-) variogram or covariance can be used to implement indicator geostatistical techniques. However, the transition probability has several advantages: • The transition probability is defined as a conditional probability, which facilitates the connection of statistical measures to geologic interpretation of facies architecture (Miall, 1973;Carle et al., 1998). • The geologically observable and interpretable parameters of proportions, mean length, and juxtapositional tendencies can be used to develop Markov chain models (Carle and Fogg, 1996;Carle and Fogg, 1997). • Non-symmetric juxtapositional tendencies can be considered (Carle and Fogg, 1996). • Three-dimensional (3-D) transition probability models of spatial variability are readily developed from 1-D Markov chains along principal stratigraphic directions (Carle and Fogg, 1997). • Continuous-lag Markov chains have been found suitable for 3-D modeling of vertical and lateral spatial transitioning among geo-or hydro-facies (Carle 1996;Carle and Fogg, 1996;Carle and Fogg, 1997;Carle et al., 1998;Fogg et al., 1998;Zhang and Fogg, 2003;Proce et al., 2004;Ye and Khaleel, 2008;Engdahl et al., 2010a;Bianchi et al., 2011;Pozdniakov et al., 2012;Purkis et al., 2012;Bakshevskaia and Pozdniakov, 2016;Krage et al., 2016;Sartore et al., 2016;Zhu et al., 2016a;Meirovitz et al., 2017;Guo et al., 2019b).

Simulation with Hard Data Only
The tsim computer code is used to generate geostatistical "realizations" of categorical variables such as lithology, soil texture, or hydrofacies. The realizations generated by tsim consist of a rectangular block of regularly-spaced grid cells. The conditional simulation algorithm consists of two steps: (1) cokriging-based "sequential indicator simulation" (SIS), and (2) simulated quenching.

Sequential Indicator Simulation Step
The algorithm and code used by the SIS step in tsim is modified from the sisim code (Deutsch and Journel, 1992;Deutsch and Journel, 1998). The SIS step in tsim traces a random path through every grid cell in the realization. At each grid cell tsim uses cokriging (instead of repeated kriging steps as in sisim) to estimate conditional probabilities that a facies occurs at a grid cell given surrounding conditioning data, which are typically facies occurrences located at nearby grid cells. Initially, hard data are the only conditioning information. In the process of completing the simulation, nearby simulated values serve as hard conditioning data for the future cokriging estimates along the random path. Based on the cokriging estimates of the conditional probability that a facies occurs at a particular grid cell given facies occurrences at other nearby cells, a uniformlydistributed random number is used to select the category that occurs at a grid cell in the realization. This process continues one grid cell at a time until all cells have been reached by the random path.
The indicator cokriging estimate at a location, x 0 , is formulated as a weighted sum by (6) where i j (x α ) are indicator data values, N is the number of data, K is the number of categories, and w jk,α represent weighting coefficients. The weighting coefficients w jk,α are computed by the transition probability-based cokriging system of equations (Carle, 1996;Carle and Fogg, 1996) where N is the number of data and The indicator cokriging estimate is only an approximation of the conditional probability on the left side of Eq. 6. The SIS step provides the "initial configuration" for the next step in tsim, simulated quenching (Carle, 1997;Carle et al., 1998).

Simulated Quenching Step
The SIS step alone does not ensure that the realization will honor the model of spatial variability including all spatial auto-correlations [t jk (h) for j k] and cross-correlations The simulated quenching step is used to improve the match between modeled and simulated spatial variability by attempting to minimize an objective function, O, defined by where h l denote l 1, . . . , M specified lag vectors and "SIM" and "MOD" distinguish simulated and modeled transition probabilities, respectively (Aarts and Korst, 1989;Deutsch and Journel, 1992;Deutsch and Journel, 1998;Deutsch and Cockerham, 1994;Carle, 1997;Deutsch and Journel, 1998). Simulated quenching is implemented by cycling through every grid cell of the realization several times along a random path and querying whether a change in facies will decrease O; if so, the category is changed. Conditioning of hard data is maintained during quenching by not allowing changes of categories at conditioning locations. The quenching step continues until a specified number of iterations through the every grid cell is reached or O is reduced below a specified minimum threshold Carle, 1999. Simulated quenching is the "zero temperature" form of simulated annealing, where an "annealing schedule" determines a probability of acceptance for changes that increase O to avoid high-valued local minima in the solution space for O (Deutsch and Journel, 1992;Carle, 1997). The main advantages of using simulated quenching over annealing are 1) the difficulty of designing and implementing an annealing schedule is avoided and 2) quenching is much faster. In tsim, the cokriging-based SIS step avoids high-valued local minima of O by providing a spatially-correlated initial configuration prior to quenching. The quenching step simply modifies existing spatial structures in the initial configuration to be consistent with the transition probability model.

Soft Data Conditioning in Categorical Geostatistical Simulation
Two concepts are presented here to enable location-specific soft data conditioning for mutually exclusive categories: • "prior probability," which assigns probabilities between zero and one to each category (Deutsch and Journel, 1998;Ying, 2000), and • "hardness," which assigns one uncertainty measure to account for overall uncertainty of the data at a given location (Deutsch and Wen, 2000).
A common practical situation is that the data are categorized (e.g., as textures in driller's logs or by interpreted hydrofacies) but known to be uncertain (e.g., because the driller's logs or geological interpretations are not 100% accurate). In this situation, the latter approach is more straightforward to apply, although both approaches can be applied simultaneously.

Prior Probabilities
Uncertainty in the indicator data or the "softness" of categorical variables can be accounted for in tsim by assigning prior probabilities between zero and unity to the indicator values, which  implemented in application use of driller's logs. This prior probability approach has also been implemented in integration of geophysical imaging to condition stochastic simulations of hydrofacies architecture Carle and Ramirez, 1999;Hermans and Irving, 2017). However, the prior probability approach can be problematic in many practical situations for several reasons: • The most readily available categorical data (e.g., from driller's logs or geological interpretations) are discretely categorized and not presented in the form of prior probabilities. • The practitioner may find difficulty, tedium, and confusion in assigning multiple prior probabilities between zero and unity to the different facies categories over hundreds, thousands, or more data points. • Uncertainty in indicator data is not accounted for in the cokriging Eq. 7 and objective function Eq. 9 used the SIS and simulated quenching steps, respectively, of tsim.
The latter reason is a persistent theoretical shortcoming of past and current bivariate geostatistical methods, wherein the (co) kriging or simulated annealing equations are weighting the soft data in the same manner as the hard data. On this topic, Deutsch and Wen (2000) stated that "a significant problem with krigingbased approaches (to stochastic simulation) is that there is no convenient way to handle the fact that the soft data have locally variable precision (or accuracy)" and, as a result, proposed a simulated-annealing approach. However, we believe there is a convenient way to handle soft data in both cokriging and simulated quenching of categorical variables.

Categorical Data with Hardness
We introduce a simple alternative approach to integration of uncertain or soft categorical data into stochastic simulation by extending the concept of data hardness to categorical variables, as previously proposed for continuous variables (Deutsch and Wen, 2000). This approach requires assignment of a hardness value ranging between zero and unity to the set of facies (or indicator) probabilities given for each data location. On the extremes, a hardness value of 1.0 represents hard data, and 0.0 represents data that provide no additional information. To incorporate hardness into transition probability-based indicator geostatistics, a soft datum is assumed to consist of a weighted sum of both certain and uncertain information, with weights that sum to unity. The certain portion is a set of indicator values represented in binary form (e.g., the presence or absence of a certain lithology) or as prior probabilities, and the uncertain portion is, in effect, a state of no useful information.
Assuming stationarity, the condition of complete uncertainty for the expected value of the indicator variable, I k (x), for a facies at location x is the marginal probability or proportion, p k , such that The soft indicator value, denoted byĩ k (x), consists of a weighted sum of a hard indicator value i k (x) and the marginal probability, p k , according toĩ where the weights, α(x) and β(x), indicate hardness and softness, respectively, at location x. The hardness and softness weights are complementary to each other as Values of hardness or softness weights are assumed to depend only on location and not on individual categories. A soft indicator variableĨ k (x) is defined with respect to a hard indicator variable In practice, a single hardness value is assigned to the set of indicator values (i.e., facies probabilities). Compared to the original tsim code, hardness values are the only additional conditioning data information needed to implement the soft data approach described herein for tsim-s.

Transition Probabilities and Cokriging with Soft Data
The transition probability values used in the cokriging Eq. 7 were originally formulated under the assumption of hard data (Carle and Fogg, 1996). To incorporate soft data, transition probability values in Eq. 7 will need to be modified to reflect the uncertainty of the data used to formulate the cokriging estimate. For example, if a datum has zero hardness [α(x) 0 or β(x) 1] and, therefore, provides no additional information, that datum should not impact the cokriging estimate. In particular, the left hand side matrix of Eq. 7, which accounts for the "declustering" of the data, and the right hand side of Eq. 7, which accounts for the statistical closeness (spatial correlation) of the data with respect to the estimation location (Isaaks and Srivastava, 1989;Journel, 1992, Deutsch andJournel, 1998), should be modified to account for data uncertainty.
To account for soft data, the transition probability entries in Eq. 7 must be modified to account for decreased spatial correlation of soft data relative to hard data. This decrease in spatial correlation is derived below using the transition probability as the bivariate spatial statistic. Substituting Eq. 2 into Eq. 4, the transition probability is defined with respect to hard indicator variables by Substituting soft indicator values as defined by Eq. 12 into Eq. 13, a "soft transition probability"t jk (h) is formulated bỹ Expanding (as shown step-by-step in the Appendix), assuming stationarity, applying Eqs. 10 and 11, and combining terms, Eq. 14 reduces to Assuming that hardness values are independent and again applying Eq. 10, Eq. 15 reduces tõ Applying Eq. 13 and simplifying the right hand side, Eq. 16 reduces tõ According to Eq. 17, the soft transition probability is a weighted sum of the transition probability t jk (h) and the marginal probability p k . The weight for t jk (h) is the product of the hardness values α(x) and α(x + h) at the two datum locations, and the weight for p k is the complement to the weight for t jk (h).
To consider soft data, the transition probability-based indicator cokriging equations are simply modified by substitutingT(h) for T(h) in Eq. 7 as follows If α(x m ) 1 and α(x n ) 1, the soft transition probability matrix T(x m − x n ) is identical to the hard transition probability matrix T(x m − x n ) as defined by Eq. 5. If either α(x m ) 0 or α(x n ) 0, where the entries in each column k are the proportions p k . Assuming stationarity and ergodicity, for lags h ϕ in any direction ϕ. Eq. 18 indicates that for large lags (beyond the range of spatial correlation), the transition probabilities converge on the marginal probabilities. Thus, the matrix values forT(x m − x n ) when α(x m ) 0 and α(x m ) 0 are identical to the case for hard data where the lag (x m − x n ) is large enough such the two data are not spatially correlated. Importantly, the right hand side entries,T(x 0 − x n ) for n 1, . . . , N, should be formulated assuming α(x 0 ) 1 because: • Cokriging is used to estimate the probability that a category exists at location x 0 assuming that I k (x 0 ) for k 1, . . . , K are hard indicator variables, as defined in Eq. 1.
• The estimation location, x 0 , is known.
However, if the location of x 0 is not certain, the hardness value α(x 0 ) could be used to account for estimation location uncertainty.

Use of Acceptance Probability
In tsim, the simulated quenching step enforces hard conditioning by not allowing any changes in categories at grid cell locations with hard data. At grid cell locations with no data, categories are changed along the random path wherever O can be reduced. One can view the quenching algorithm in terms of a bi-modal acceptance probability for changing the category on the basis of reducing O: 0.0 if the category is determined by data and 1.0 if the category is not determined by data.
In tsim-s, categories are also queried for change at each grid cell, but with a lesser probability of acceptance at grid cell locations containing soft data as compared to cells with no data. The probability of accepting a change of categories that reduces O is set at β(x); corresponding to Eq. 11, the probability of rejecting the change is α(x). Thus, if hardness equals unity [α(x) 1], categories are not allowed to be changed at location x. For soft data with hardness less than unity [α(x) < 1], changes in categories that reduce O are accepted with a probability of β(x) along the random path implemented by the quenching algorithm until O is sufficiently reduced in Eq. 9. This approach enables the simulated quenching step to impart the most change in the realizations at the locations where data are least certain, and the least change where data are most certain. This algorithm is, in effect, a location-dependent simulated annealing schedule where the acceptance probability is proportional to the softness of the data.

Use of the Joint Probability
In the example discussed later in Section 3.1, one category (gravel) has a very low proportion of 0.006. In application of tsim, lowproportion facies could be problematic in the simulated quenching step because of the small amount of sample statistics for matching the measured and modeled transition probabilities. Alternatively, the objective function can be reformulated with respect to the joint probability, p jk (h) p jk (h) Pr k occurs at x + h and j occurs at x In tsim-s an option is available for using the joint probability defined in Eq. 19 to formulate the simulated quenching objective function O as Frontiers in Earth Science | www.frontiersin.org November 2020 | Volume 8 | Article 565707 for the lags l 1, . . . , M included in the quenching process. The joint probability formulation of O using Eq. 20 de-emphasizes matching of measured and modeled transition probabilities for low-proportion categories compared to Eq. 9.

Initial Characterization of Savannah River Site
The shallow subsurface beneath the Savannah River Site (SRS) in South Carolina consists of Tertiary siliciclastic sediments deposited in shoreline and nearshore depositional environments (Aadland et al., 1995;Falls et al., 1997). Characterization of heterogeneity at SRS is of interest to improve understanding of vadose zone, groundwater flow, and contaminant migration processes (Miller et al., 2000). A primary concern is characterization of vertical and lateral extent of clay lenses within the sand-dominated flow and transport regime. A two-dimensional (2-D) analysis (for lateral x and vertical z directions) is performed for the following SRS example, with the goal of generating 2-D realizations of lithofacies heterogeneity conditioned by both hard and soft data. Figure 1 shows lithologic and geophysical log data from SRS treated as hard and soft data, respectively, to condition realizations generated by tsim-s. Four texturally-based lithofacies are distinguished, with proportions in parenthesis: gravel (0.006), sand (0.735), clayey sand (0.156), and clay (0.103). The hard data (labeled "Hard") are derived from core descriptions in the borehole data, and the soft data (labeled "Soft-1" and"Soft-2") are inferred from resistivity logs from two boreholes. In some sedimentary environments, resistivity log data are not necessarily highly correlated with texturallyderived facies (Burow et al., 1997).
The lithologic data for two boreholes in the vicinity of the crosssection were used to calculate matrices of transition probability measurements with dependence on vertical lag shown by circles in Figure 2. The transition probability measurements associated with the gravel category are somewhat erratic because of the very low proportion of gravel. As a practical matter, this sort of measurement variability caused by data sparseness should not be unduly fitted in the transition probability modeling process. A Markov chain model, shown by the solid lines in Figure 2, was fitted the calculated vertical transition probabilities through use of a matrix exponential function where the rate coefficients in R z are given in units of m −1 . As shown by Agterberg (1974) and Carle and Fogg (1997), the entries in T(h z ) for Eq. 21 are computed by an eigenvalue decomposition, so that each entry, t jk (h z ), is a weighted sum of the k (column) category proportion and three exponential functions: For the lateral x direction, Markov chain transition probability models were developed for the SRS example based on prior geological estimates of lithofacies mean lengths and juxtapositional tendencies indicated by geologic cross-sections. Such geological information can be converted into lateral transition rates (e.g., R x ) to enable development of 2-and 3-D Markov chain models and demonstrated by Carle and Fogg (1997), Carle et al. (1998), and.  Figure 1 at z 98-102.5 m, more clay and clayey sand is indicated on the right. The realizations in Figure 3 honor these data, but to a lesser extent than between z 87-90 m where lateral correlation of the soft data is stronger. This comparison is a simple example of how tsim-s can be used to further constrain the realizations with soft data as compared to Frontiers in Earth Science | www.frontiersin.org November 2020 | Volume 8 | Article 565707 10 conditioning with only hard data, as could otherwise be implemented with tsim.

Effect of Soft Data
Another way to visualize the impact of the soft data is to average the indicator values over many realizations to produce a "probability map" for each lithofacies. Figure 4 shows probability maps for each lithofacies derived from 100 realizations with and without soft data conditioning. In case (A) with hard data only, the lithofacies probabilities approach marginal probabilities (proportions) toward the left portion of the realizations. In case (B) with hard and soft data, the soft data further constrain the probability structure toward the left side of the realizations. However, less-refined "gray areas" remain between the borehole data because of the limited lateral correlation of the lithofacies units. Overall, the soft conditioning at 0.80 hardness imparts a strong influence on the realizations.
Through the hardness parameter, the degree of influence by the soft data can be controlled as needed in application of tsim-s. Figure 5 shows probability maps for the sand, clayey sand, and clay lithofacies where hardness values for the soft data are reduced to 0.5 and 0.2. Reducing the hardness produces less contrast or Frontiers in Earth Science | www.frontiersin.org November 2020 | Volume 8 | Article 565707 more gray areas in the probability map, which reflects the increased uncertainty of the soft data.

Use of Realizations as Soft Conditioning
There is increasing interest in use of stochastic inversion approaches to modify heterogeneity structures within geostatistical realizations to be more consistent with geophysical or hydraulic testing data Harp et al., 2008;Wainwright et al., 2014;Berg and Illman, 2015;Wang et al., 2017). To implement these approaches, there can be a need to modify an initial heterogeneity structure in a controlled or incremental manner. Another application of the soft data capability in tsim-s is to use all or part of a realization as soft conditioning to exert control on the modification of heterogeneity structures from one realization to the next. One realization (or any available field of categorical values) can be used as prior information for conditioning of a new realization, which makes several new capabilities available in tsim-s: • The degree of correlation between a series of realizations or can be controlled. • The degree of variation from one realization to the next can be controlled at different locations within each realization. • By exerting control on the difference between one realization and the next, Monte Carlo Markov chain algorithms can be implemented as a Bayesian inverse approach to optimization of local heterogeneity structure Wainwright et al., 2014;Wang et al., 2017).
In practice, we refer to a realization used for soft conditioning as the "prior realization" and a subsequent realization that is produced as the "posterior realization." Use of prior realizations as soft conditioning in the SIS step is implemented by adding only one additional soft datum to each cokriging estimate Eq. 6; that additional soft datum consists of the indicator values from the grid cell of the prior realization corresponding to the cokriging estimation location for the posterior realization. More soft data can be used, but one prior realization soft datum at the cokriging estimation location itself provides sufficient conditioning for generating correlated realizations without adding much more computational burden. The degree of correlation between the prior and posterior realizations is controlled by setting the hardness values, which may vary with location. Figure 6 shows a sequence of six realizations where each posterior realization is soft-conditioned to the prior realization for cases of (A) hardness 0.5 and (B) hardness 0.9. Because the degree of hardness controls the degree of similarity (or rate of change) between one realization and the next, the realizations in case (A) are less similar (or more different) from one realization to the next. Comparing the left side (away from the hard conditioning) of realizations 1 and 6, similarities have largely disappeared for case (A) but remain for case (B). Thus, the higher rate of change (lesser hardness) shortens the memory of heterogeneity patterns within a sequence of realizations softconditioned by prior realizations. Figure 7 shows an example using the same realization as soft data (upper left) with hardness varying with depth from 0.8 at top to 0.0 at bottom. The probability maps for gravel, sand, clayey sand, and clay facies, derived from 100 realizations generated by tsim-s, show increasing uncertainty with depth. In particular, the location of individual clayey sand and clay lenses in the soft data become less distinct with depth in the probability maps. Likewise, the sand probabilities become less distinct with depth and approach the proportion of 0.735. The influence of a single gravel lens in the soft data is evident in the gravel probability map. This example illustrates how 2or 3-D geophysical images or geological interpretations of categorical data might be used to soft-condition geostatistical realizations with consideration of variable spatial resolution such as decreasing resolution with increasing depth.

Computational Aspects
Like tsim, tsim-s is readily extended to large-scale threedimensional (3-D) applications. From a computational standpoint, the main difference between running tsim and tsim-s will be in memory use. For very large grid cell counts (tens of millions or more), the memory use in tsim is mostly taken up by integer arrays storing the grid cell category and anisotropy direction information (if used). In the late 1990s when memory was quite restricted compared to the year 2020, a 45-million cell 3-D realization was produced on a computer with 48 megabytes (not gigabytes!) of memory by modifying tsim to use a 1-byte integer format for the grid cell array and implementing an analytical function to define variable stratigraphic directions (Carle, 1996;Carle et al., 1998;Carle, 1999;Tompson et al., 1999). The open-source nature of the T-ProGS package of codes enables the user to make similar modifications to conserve memory. Depending on the application, tsim-s will require a factor of as much as ten or more times the memory use as tsim, but computational time tsim-s will not be significantly higher because the corresponding arrays for the cokriging equations and quenching objective function are of the same dimensions under similar model parameter settings. Considering that computer memory and computational speed are orders of magnitude higher today and into the future as compared to the late 1990s, we do not anticipate large-scale simulation of categorical heterogeneity with tsim-s to be highly constrained by the current or future computational technology. Frontiers in Earth Science | www.frontiersin.org November 2020 | Volume 8 | Article 565707 16

Llagas Basin Example
As discussed in the introduction, the use of driller's logs in categorical geostatistical simulation is of interest in hydrogeology, yet presents some uncertainty in how to treat driller's logs as conditioning data. As a demonstration of the use of driller's logs and tsim-s to large-scale simulation, Figure 8A shows an interpreted version of a driller's log dataset, and Figure 8B shows a 3-D simulation of hydrostratigraphic architecture for the Llagas groundwater subbasin south of San Jose, California Carle et al., 2006).
In this application, the driller's log data were categorized into aquifer, interbedded, and aquitard hydrofacies, where the interbedded category represents relatively thin interlayers of aquifer and aquitard materials. The addition of the interbedded hydrofacies addresses unresolvable fine-scale heterogeneity of coarse-and fine-grained textural classifications and adds flexibility to development of the Markov chain model. If only two hydrofacies, aquifer and aquitard, had been distinguished, the 3-D Markov chain model would consist of only four parameters -proportion and mean length in the three principal depositional directions for one of the hydrofacies; the remaining transition probabilities are all determined by probability law (Carle and Fogg, 1996). A twocategory characterization of heterogeneity presents no distinct advantage in the Markov chain spatial variability modeling framework; the spatial variability model is equivalent to a twocategory indicator variogram-based approach modeled by an exponential variogram. With three categories, the number of parameters for the 3-D Markov chain is raised to 14, allowing for development of more complexity in the heterogeneity structure including asymmetry such as fining-upward and outward tendencies (Miall, 1973;Carle and Fogg, 1997;Fogg et al., 1998).
This simulation was generated by tsim-s with soft conditioning data divided into three sets of driller's logs to which hardness levels were set to 0.3, 0.7, and 1.0 based on data quality. Two separate tsim-s simulations of 162,500,000 and 117,000,000 cells were generated for the upper and lower portions of the final realization of hydrofacies architecture to address differences in the spatial structure of deeper and shallower alluvial hydrofacies evident in the driller's logs. Carle et al. (2006) provides further detail on the process of selection of the hardness parameter in the practical hydrogeological situation of using driller's logs for conditioning data. This example confirms that large-scale 3-D stochastic simulation using tsim-s is feasible. The stochastic analysis was further applied to investigate permeability heterogeneity effects on nitrate transport from agricultural sources toward municipal wells .

DISCUSSION
The discussion focuses on capabilities and limitations of t-sim and tsim-s with attention to the current literature on comparison and evaluation of geostatistical methods for subsurface characterization of categorical variables.

Curvilinear Features
Both tsim and tsim-s have the capability to produce curvilinear features by specifying azimuthal and dip angles local to each grid cell in an a priori manner (Carle, 1999; that can be deterministic Carle et al., 2006) or stochastic . These angles may be derived or inferred from prior geological knowledge or geologically reasonable interpretation of the depositional or stratigraphic architecture, surface mapping of the deep soil horizons, interpretation of seismic or surface geophysical data, or stochastic modeling (e.g., by modeling variation in the major axis of deposition due to meandering by a gaussian random field). Local anisotropy directions are implemented in tsim and tsim-s in a manner similar to "local anistropy kriging" (te Stroet and Snepvangers, 2005). This is not a coordinate transformation approach, as applied to the variogram-based isim3D (Gomez-Hernandez and Srivastava, 1990).
The Llagas basin example application of tsim-s uses the local anisotropy direction option to impart variable angles of dip and principle direction of deposition into the geostatistical realizations. Figure 9 shows vertical and horizontal slices through the 3-D Llagas basin realization to better reveal the nature of the curvilinear features. These include variable dip angles evident in longitudinal and transverse-plane cross-sections (A) and (B) and variable direction in the major axis of deposition in the horizontal-plane cross-section (C).
The simulated aquifer/interbedded facies architecture for the Llagas basin example does not show pronounced sinuousity, which is consistent with how channel belt deposits were conceptualized in this alluvial depositional setting by the California Department of Water Resources (1981). This is in contrast to the "true" or "training image" concept of continuous "channels" worming their way through homogeneous lowpemeability media, a common argument posed for replacing bivariate statistical methods with multi-point statistical (mps) methods (Strebelle, 2000;Caers, 2001;Krishnan and Journel, 2003;Feyen and Caers, 2006;Ronayne et al., 2008;Li et al., 2014a;Li et al., 2015;Mariethoz and Caers, 2015;Zovi et al., 2017;Ramgraber et al., 2020). Such surficially-based conceptual models gloss over fundamental geologic concepts showing how depositional processes produce amalgamations of channel and adjacent sediments that are broader, less sinuous, and more variable in lateral extent as compared to an active fluvial channel viewed on the Earth's surface (Galloway and Hobday, 1996;Miall, 2013). It is a well-known fact in interpretation of borehole data that sedimentary features on the surface are not necessarily preserved in the subsurface (Smith, 2002).
We further discuss methods comparison below because tsim-s carries forward several capabilities of tsim that appear not to be recognized in the methods comparison literature. Flawed methods comparison causes a trickle-down effect of selective references that propogate misleading appraisals of the capabilities of the available methods for stochastic simulation of categorical variables.
In a methods comparsion of sisim, tsim, and mps, dell'Arciprete et al. (2012) did not consider variable anistropy directions in either conceptualization or parameterization of their analysis of spatial variablity and applications of tsim despite obvious dipping structures in their data. They chose not to apply Markov chain modeling concepts in the structural framework of a stratigraphic coordinate system relevant to the hydrofacies of a sedimentary depositional system Tompson et al., 1999;Carle et al., 2006). The example of a tsim realization shown in Figure 4 Carle and Fogg (1997). To "reveal" limitations of Markov chain models, Langousis et al. (2018) execute a "simple test-case" of a 2-D dipping layer with a coordinate system anchored in the vertical and horizontal directions of their statistical analysis. Neither dell' Arciprete et al. (2012) nor Langousis et al. (2018) appeared to grasp the concept that bivariate geostatistical analysis and simulation should be applied in a geologically-based coordinate system, as demonstrated for at least 30 years (e.g., Gomez-Hernandez and Srivastava, 1990). The geological realism of the mps application by He et al. (2017) stands wide open to geological criticism too, with a mps-generated "real world buried valley system" showing unrealistic topography and isolated occurrences of the valley fill buried beneath pre-valley-fill strata. Geostatistical analysis should recognize and make adjustments to account for geological slopes, directions, and depositional ordering and not be strictly anchored in a rectilinear cooordinate system pertinent only to data locations.
Another pitfall of methods comparison is use of oversimplified example applications. In Kessler et al. (2013), tsim was compared to mps in application to stochastic simulation of sand lenses within clayey till. This application of tsim adhered to an oversimplified two-category Markov chain to model "a complex type of heterogeneity" exhibiting a bi-modal distribution in the size of sand lenses in 2-D training and reference images. A simple way to address this facies size distribution complexity would be to define two types of sand lenses (perhaps "small" and "large") and model a three-category system. Instead of adhering to a textural basis for defining categories (e.g., "sand" and "clay"), a more geological approach is to interpret hydrofacies appropriate to the geometric framework of the depositional system .
The impact of methods comparison studies can get more convoluted by selective referencing of the literature, particularly on the topic of curvilinear features. For example, in promoting mps, Barfod et al. (2018) dissmiss applicability of tsim by stating that ". . .T-ProGS also has difficulties in reconstructing curvilinear geological features" without any reference to previous work in which the variable anisotropy direction capability of tsim was implemented. Barfod et al. (2018) refer to the oversimplified two-category sand-clay Markov chain analysis of Kessler et al. (2013) as "revealing a sub-optimal pattern reproduction, in comparison to other simulation tools such as multiple-point statistics. . . " Linde et al. (2015) in reviewing "variogram based models" state that "transition probability techniques such as T-ProGS . . . cannot properly produce curvilinear features. . ." This is after referencing only dell'Arciprete et al. (2012), Falivene et al. (2007), Lee et al. (2007), andReffsgaard et al., 2014) as T-ProGS applications, all four of which did not employ variable anisotropy direction capability in tsim. The misconception that all bivariate statistical methods for stochastic simulation of categorical variables lack the ability to produce curvilinear features appears to derive from selective referencing and comparison of mps to only variogram-based methods (Strebelle, 2002;Caers, 2001;Krishnan and Journel, 2003;Feyen and Caers, 2006;Linde et al., 2015;Barfod et al., 2018) irrespective of previous hydrogeologic studies Tompson et al., 1999;Carle et al., 2006;Green et al., 2010;Engdahl et al., 2012) and the T-ProGS user manual (Carle, 1999;.
If tsim or tsim-s is producing spatial structure that falls well short of the intended model for represention of the spatial heterogeneity, the conceptualization, implementation, or utilization of the capabilities of tsim or tsim-s may be at fault. Methods comparison studies in geostatistics should fully investigate the capabilities of each method including both the geological and statistical conceptual underpinnings before making sweeping judgments. The introduction of this paper provides references to varied applications of T-ProGS that should be useful for methods comparison and capabilities assessment.

Method Limitations
All models have limitations. The T-ProGS package was originally conceived in 1996 (Carle, 1996;Carle, 1997;Carle and Fogg, 1996;Carle and Fogg, 1997) to improve or add to the capabilities of the then state-of-the art variogram-based geostatistical methods of the Geostatistical Software Library (Deutsch and Journel, 1992) and subsequently analyze pumping test data in a highly heterogeneous alluvial aquifer system (Carle, 1996;Lee et al., 2007). It is easy to argue that natural systems are geometrically inconsistent with certain geostatistical representations or exhibit far more complexity than any bivariate geostatistical model could ever characterize. Different practitioners will be more comfortable with different methods or levels of complexity in the statistical or bio/chemo/hydrogeological components of their models. There will always be an open question as to what levels of complexity are necessary to sufficiently analyze the subsurface processes of interest.
The stationarity assumption in regard to the spatial continuity model and category proportions can be limiting. Early applications of tsim recognized this issue and incorporated "nonstationary" qualities into the geostatistics through data conditioning, geologically-based zones based on stratigraphic analysis or sedimentary environments, and spatially variable angles of anistropy (Carle, 1996;Carle, 2000;Carle et al., 1998;Tompson et al., 1999;. There continue to be applications of tsim that recognize and address nonstationarity Traum et al., 2014;Weissmann et al., 2015;Meirovitz et al., 2017;Zhu et al., 2016a;Zhu et al., 2016b;Zhang et al., 2018;Liao et al., 2020;Maples et al., 2020). As discussed previously for application of tsim-s, Carle et al. (2006) employed two-zones of differing spatial statistics and apply variable anisotropy directions in their application of tsim-s. To address geological realism, adding a modicum of geological insight can be more effective than adding more statistical complexity.
The model of spatial variability is another limitation to any geostatistical approach. A methods limitations analysis by Langousis et al. (2018) criticizes the interpretive framework of transition probability-based Markov chain model development as "based on unverified/untested simplifying assumptions" and "adhoc manipulations." Langousis et al. (2018) further contend that.
. . .stochastic modeling of actual geologies using the [T-ProGS] approach of Carle and Fogg (1997), is characterized by simplifying assumptions and theoretical limitations, with the simulated random fields exhibiting statistical structures that strongly depend on the problem under consideration and the modeling assumptions made, leading to increased epistemic uncertainties in the obtained results.
We offer some perspectives on assessing the limitations of T-ProGS that apply to tsim-s as well: • Use of geological concepts in categorical geostatistical simulation may involve subjectivity, which can be viewed by some as either a strength or a limitation to reducing uncertainty. • The implication that injection of subjective geologic interpretation increases epistemic uncertainties in stochastic modeling would appear to expose a lack of understanding of subsurface geology and it's role in modeling the subsurface. • As referenced in the introduction of this paper, many applications in hydrogeology and related fields have found the Markov chain modeling framework to be useful to characterization of bivariate spatial statistical cross-relationships (i.e., juxtapositional tendencies), which can be related in an interpretive manner to geological concepts such as Walther's Law in the stratigraphic context of sedimentary depositional environments (Leeder, 1982;Doveton, 1994). • As evident in the T-ProGS manual (Carle, 1999;, the Markov chain is not actually required to run tsim (or tsim-s). So one could investigate epistemic aspects of other transition probability or indicator covariance models using tsim, tsim-s, or the various variogram-based simulation methods, if so desired. However, such methods comparison exercises, which can certainly be expanded to all stochastic models, will not prove that a spatial Markov chain is not useful to modeling subsurface hydrofacies heterogeneity when applied in the appropriate geologic context.
In every subsurface heterogeneity modeling project we have encountered, the heterogeneity contains both deterministic and stochastic aspects that should be treated differently. For example, conventional geologic stratigraphic analysis is quite effective for identifying and mapping the major formations, depositional systems, and the bounding unconformities and structural discontinuities (e.g., faults). Accordingly, such features can typically be treated deterministically and then used as the basis for appropriate geologic zoning of the system into quasistationary subdomains. The stochastic aspect of subsurface characterization typically lies within those geologically defined subdomains. If one, however, lumps together both the deterministic and stochastic aspects of the heterogeneity and calls on the stochastic geostatistical algorithm to sort out those spatial patterns, one only invites naive mischaracterization of the heterogeneity that produces unnecessary uncertainty and unrealistic results. An effective way to reduce or moderate uncertainty is to recognize and separate out deterministic and stochastic parts of the problem.

Stochastic Methods Evaluation
Stochastic methods evaluation in hydrogeology should strive to determine appropriate levels of complexity necessary for gaining insight to 3-D subsurface flow and transport processes at scales relevant to measurement diagnostics used in decision-making. An example of a process-oriented methods comparison, Damico et al. (2018) compared dynamics and trapping metrics for 3-D carbon-dioxide plume simulations using subsurface representations of heterogeneity derived from tsim and a more rigorous model for representation of complex features in fluvial architecture. They concluded: . . .in the context of representing plume dynamics and residual trapping within fluvial deposits, and within the scope of the parameters used here, the simpler geostatistical model of braided fluvial deposits appears to give an adequate representation of the smaller scale heterogeneity. The depositional-and geometric-based benchmark models represented more features of the fluvial architecture, including variability in the dip of cross sets, variability in the geometry and orientation of unit bars, and the occurrence of channel fills. Depending on context, representing those features may be quite important to understanding some multiphase flow processes in aquifers and reservoirs. However, the simpler geostatistical model (tsim) is able to capture the important aspects of fluvial architecture within the context of understanding the general effect of smaller scale heterogeneity on residual trapping of CO 2 in geosequestration reservoirs, within the scope of the parameters used here.
The worthiness of a stochastic methods to subsurface applications does not necessarily depend on statistical rigor or geological detail, it also depends on the perception of relative value, including both benefits and costs, in the application (Ginn, 2004). The relative value of a model is not necessarily in its complexity, given that calibration of more complex models may erode rather than enhance predictive ability (Doherty and Christensen, 2011). This paper is offering an approach to address uncertainty in data conditioning of categorical geostatistical models. It would be quite straightforward to add more statistical complexity to transition probability or hardness concepts. However, in our nearly 25-year experience with using transition probabilitybased geostatistics, we find simpler and more interpretetable geologically-based tools to be quite useful in the study of the effects of various scales and types of heterogeneity on subsurface flow and transport processes.

CONCLUSIONS
Many Earth science applications would benefit from increased ability to incorporate "soft" (uncertain or indirect) data to further constrain subsurface models of heterogeneity. In categorical geostatistical simulation applications, often abundant soft data on lithology or hydrofacies (e.g., geophysical logs and imaging, geological interpretations, driller's logs, etc.) offer opportunity for imposing increased or relaxed model constraint.
A soft data capability has been incorporated into the categorical geostatistical simulation code tsim-s. Soft data for categorical variables are input either as indicator values or prior probabilities, and a "hardness" value accounts for uncertainty in the data. This approach is particularly conducive to soft data that is already categorical, such as texture inferred from driller's logs, hydrofacies interpretations, or electrofacies based on resistivity cutoffs. In generating realizations with tsim-s, the impact of uncertainty in the soft data is factored into formulation of both the cokriging and simulated quenching geostatistical simulation steps. The extent to which the realizations honor the soft data is balanced by the values of hardness, the model of spatial variability, and the values of other nearby hard and soft data.
The degree to which soft data reduces variability in simulation outputs can be quantified by mapping facies probabilities derived by averaging indicator values from many realizations. Example applications in this paper using different values and spatial distributions of hardness illustrate how the impact of data uncertainty can be controlled in the stochastic realizations. Such control will be useful for assimilating different data sets of variable resolution and accuracy. The soft conditioning can be arrays of data, including "prior realizations," to incrementally adjust or evolve the spatial heterogeneity structure of the realizations. The ability to manipulate localized heterogeneity structure or rate of change in a sequence of realizations should be useful for flow and transport model calibration, inverse approaches, or sensitivity analysis. The tsim-s algorithm is ammenable to large-scale 3-D simulation including curvilinear features.
Overall, the tsim-s code more rigorously integrates data uncertainty and prior information into the categorical stochastic simulation algorithm as compared to previous indicator-based geostatistical simulation codes including its direct predecessor tsim. However, users and evaluators of bivariate geostatistical models should become familiarized with capabilities, limitations, and varied uses of T-ProGS or other geostatistical software packages before applying or evaluating tsim or tsim-s.

DATA AVAILABILITY STATEMENT
The data analyzed in this study is subject to the following licenses/ restrictions: must be obtained by permission. Requests to access these datasets should be directed to carle1@llnl.gov.

AUTHOR CONTRIBUTIONS
SC developed code, executed applications, and wrote first draft of article. GF aided in theoretical development, outside application, tie-in to other research, and assisted SC in revision of the manuscript.