Improved implementation of travel time randomization for incorporating Vs uncertainty in seismic ground response

1Department of Civil, Architectural, and Environmental Engineering, University of Texas at Austin, USA 2Department of Statistics and Data Science, University of Texas at Austin, USA 3Department of Civil and Environmental Engineering, Utah State University 4Department of Structural, Geotechnical and Building Engineering, Politecnico di Torino, Italy 5The Charles E. Via, Jr. Department of Civil & Environmental Engineering, Virginia Tech, USA *Corresponding author: mhallal@utexas.edu


Introduction
Predicting the influence of subsurface geomaterials on earthquake ground shaking, known as site effects, is a crucial aspect that controls the seismic hazard for many engineering projects. In current practice, it has become the standard to conduct seismic hazard analyses for critical facilities and major infrastructure within a probabilistic framework, and to model site effects using predominantly one-dimensional (1D) ground response analysis (GRA). Owing to the flexibility of probabilistic seismic hazard analysis (PSHA), different sources of uncertainty can be integrated into the hazard calculations using procedures that continue to evolve. However, recent studies have shown that the uncertainty in 1D GRA and its propagation into PSHA can significantly drive the hazard estimates. Thus, proper quantification and modeling of this uncertainty is of paramount importance (Passeri et al., 2019;Pehlivan et al., 2016;Rathje et al., 2010).
Uncertainties in 1D GRA namely stem from: (1) selection of the input rock motions, (2) characterization of the shear wave velocity (V s) profile, (3) characterization of the nonlinear dynamic material properties, (4) specification of the analysis method, and (5) suitability of the underlying 1D site assumptions (Passeri et al., 2019;Rathje et al., 2010;Thompson et al., 2012). In this article, we focus on the V s profile, which is a fundamental component of any site-specific GRA. The V s profile can be characterized using invasive/borehole methods (e.g., downhole, crosshole, and P-S suspension logging) or noninvasive/surface wave methods (e.g., active-and passive-surface wave measurements). However, regardless of the dynamic site characterization method, uncertainty in the derived V s profile is inevitable, and the purpose of this article is to highlight challenges that are encountered in current and newly proposed procedures that attempt to model this uncertainty.
Invasive and noninvasive V s measurements are influenced by two types of uncertainty: (1) aleatory variability, and (2) epistemic uncertainty. Aleatory variability is related to the inherent spatial variation of subsurface properties. This variability corresponds to natural changes in the layering and stiffness of subsurface materials, and it can be quantified by testing at multiple locations across the site. Epistemic uncertainty, on the other hand, stems from uncertainties in the measured data and those introduced during analysis. This uncertainty is attributed to our inability to perfectly measure the V s profile and it can be quantified by using multiple characterization techniques and/or different data processing and analysis methods. Although this distinction of uncertainties is conceptually evident, their separation in practice is nontrivial.
For many years, both variability and uncertainty have been modeled as a single random variable, in which a reference base-case V s profile is randomized using a statistical model, such as that developed by Toro (1995). More recently, the state-of-practice in PSHA has moved toward separating aleatory variability from epistemic uncertainty. It is becoming increasingly common to model the latter using alternative base-case profiles, typically a best-estimate, upper-case, and lower-case V s profiles. Even when using alternative base-case profiles to account for epistemic uncertainty in V s, aleatory variability within PSHA is additionally accounted for by randomizing each of these alternative V s profiles (EPRI, 2012).
A growing body of literature has investigated the efficacy of these commonly used approaches in capturing V s aleatory variability and epistemic uncertainty. For example, several researchers have compared site response predictions from V s randomization with ground motion observations at vertical seismometer arrays (Field and Jacob, 1993;Kaklamanos et al., 2020;Passeri et al., 2020;Teague et al., 2018). Besides these studies, others have performed theoretical assessments (without comparison to observations) to examine the influence of V s randomization (Pehlivan et al., 2016;Rathje et al., 2010;Barani et al., 2013;Bazzurro and Cornell, 2004;Griffiths et al., 2016a;Tran et al., 2020) and alternative V s profiles (Passeri et al., 2019;Griffiths et al., 2016a) on site response predictions. An unsettling conclusion from these studies is that, although it is common for engineers to use randomized and/or alternative V s profiles to account for V s uncertainty in GRAs, these procedures, as currently implemented, often generate unrealistic V s profiles that do not represent well the actual site conditions. More specifically, both approaches commonly generate V s profiles that fail to match other site-specific experimental data (e.g., layering from boring logs, fundamental site frequency, surface wave dispersion data) and result in uncontrolled and potentially unconservative median site response predictions, which can be detrimental for seismic hazard estimates. These challenges have motivated the development of alternative approaches to effectively model general V s uncertainty. One such development is that proposed by Passeri et al. (2020) based on randomizing the cumulative shear wave travel time (tts).
Despite several potential improvements gained through tts randomization, more validations are required to ensure its successful implementation in practice.  assessed multiple approaches to account for V s variability in 1D GRAs, including tts randomization. The authors concluded that although tts randomization likely offers an improvement over traditional V s randomization, some outstanding challenges remain that require further attention before adopting this approach in practice. The aim of this article is to scrutinize these challenges and to provide guidelines which will work toward better modeling of V s variability in 1D GRA. The findings of this study will enhance our understanding of the tts randomization framework, and the recommendations we provide will allow researchers and engineers to adopt this framework more confidently.
The article begins with a presentation of the statistical framework of V s randomization and a discussion of its challenges. Then, the adjustments introduced by Passeri et al. (2020) to tackle these challenges using the tts randomization framework are described. Despite the potential benefits of the tts randomization framework, we argue, based on real experimental data and Monte Carlo simulations, that some of the introduced adjustments may pose additional challenges, which might in turn lead to bias and unrealistic results if not implemented appropriately. For each challenge, we examine revised implementation approaches and demonstrate that they provide a clear improvement on the current framework. The recommendations proposed in this article do not require any additional model parameters beyond those specified by Passeri et al. (2020), but rather focus on methodologies for improved implementation. Finally, we close with a discussion on future directions and potential improvements in need of further investigation.
2 Vs randomization and its challenges 2.1 Statistical framework of Vs randomization V s randomization is currently performed within a Monte Carlo framework, predominantly using the Toro (1995) model, which is one of the earliest published approaches and remains one of the few systematic statistical models available for this purpose. This framework operates on three components: (1) a layering model to randomize layer thicknesses, (2) a bedrock model to vary bedrock depth, and (3) a velocity model to simulate V s variability. For these various components, Toro (1995) developed input parameters based on a total of 557 V s profiles, most of which are from various California sites, and these input parameters are applicable to specific site classes (across an entire seismic site class) or applicable to a generic site (across all seismic site classes).
The layering model attempts to reflect the trend observed by Toro (1995) in layer thicknesses, which is thin layers near the surface transitioning to thick layers at depth. This is achieved by modeling the occurrence of a layer interface as a non-homogeneous Poisson process with a depth-dependent occurrence/interface rate parameter (λ), as shown in Fig. 1a. Toro (1995) defined λ using a modified power-law function: where c 1 , c 2 , and c 3 are fitting parameters and z is depth in meters. Toro (1995) recommended generic values of 10.86, 0.89, and 1.98 for c 1 , c 2 , and c 3 , respectively, which result in a decreasing value of λ with depth ( Fig. 1a).
Because the occurrence of a layer interface is modeled as a non-homogeneous Poisson process, the distance between two interfaces, which is layer thickness, follows an exponential distribution with the same parameter λ. The theoretical mean layer thickness is the reciprocal of the occurrence rate parameter (i.e., 1/λ), as shown in Fig. 1b. A sample exponential probability density function for thickness at a particular depth (around 18 m) with λ = 0.1 is shown in Fig. 2. By plotting similar probability density functions for layer thickness with depth, a twodimensional thickness probability density function can be visualized, as shown in Fig. 1b. Additional information contained in Fig. 1 will be discussed later in the paper.
Figure 1: The layering model developed by Toro (1995) to randomize V s profile layer thicknesses. Shown are: (a) a depth-dependent rate parameter (λ) using the generic parameters recommended by Toro (1995), and (b) a two-dimensional probability density function (pdf) of the layer thicknesses, which is defined by the same depthdependent λ. Also shown are a set of 275 simulated layers from 50 V s profiles. Each circle symbol corresponds to a simulated layer plotted at its mid depth. Bedrock depth in the Toro (1995) model is simulated separately from the soil layer interfaces. However, unlike the layering model, little guidance is provided. Toro (1995) utilized a uniform distribution but did not specify how the parameters were set. In research (Passeri et al., 2020;Teague et al., 2018;Griffiths et al., 2016a,b) and practice, uniform, normal, and lognormal distributions have been used to model the depth to bedrock, and their parameters are either assumed or calculated when multiple V s measurements are available.
One-dimensional V s profiles in the Toro (1995) model are simulated from a non-stationary lognormal random field, which is defined by a depth-dependent V s lognormal standard deviation (σ lnV s ) and an inter-layer correlation model (ρ). Toro (1995) observed that if the V s in a layer is above average, then the V s in the layer directly below it has a similar tendency. This tendency is modeled using a positive ρ, which is defined as a function of depth-dependent and thickness-dependent correlations (ρ d and ρ t , respectively): where ρ 200 , d 0 , b, ρ 0 , and ∆ are parameters for which Toro (1995) provided site class generic values, d is the average of midpoint depths of successive layers, and t is the difference between these midpoint depths. Generally, this model gives higher ρ for deeper and/or thinner layers.
Finally, a first-order autoregressive function is used to define V s in each layer. This model predicts the V s of a layer as a function of the V s of the layer above it and ρ. The V s of layer i (V s i ) is simulated from a lognormal distribution with a median equal to the base-case V s at its mid-depth ( V s i ) and a specified σ lnV s , as given by: where Z i and Z i−1 are standard scores (number of standard deviations from the mean, i.e., ln( V s i )) for layers i and i − 1, respectively, and ε i is a standard normal random variable for layer i. Toro (1995) provided site class generic values for σ lnV s that range between 0.27-0.37. Stewart et al. (2014) considered a subset of the V s profiles from Toro (1995) and recommended "site-specific" σ lnV s values of 0.15 for depths ≤ 50 m and 0.22 for depths > 50 m, which are preferable whenever a site-specific V s profile is available and when spatial variations are not expected to be extreme.

Challenges of Vs randomization
Most researchers agree that V s profiles developed from blind application of V s randomization using generic parameters are often unrealistic, even when using the lower, "site-specific" σ lnV s from Stewart et al. (2014). The unsatisfactory results from blind application of V s randomization can be attributed to a combination of factors. First, the assumed non-homogeneous Poisson model proposed by Toro (1995) results in layer thicknesses that are highly variable, even if the mean values are reasonable. For example, and as discussed in the next paragraph, Fig. 1b shows simulated layers that display extreme variability. Second, the median V s of each randomized layer is set based only on the base-case V s at its mid-depth (Equation 5), which is a discrete point. While the mid-depth V s might be appropriate in some cases, it provides very limited information that ignores potential V s variations at shallower and/or deeper depths along the simulated layer thickness. Third, separate randomization of layer thicknesses and V s may result in double counting of uncertainty (Passeri et al., 2020). These factors combined may yield shallow, thick layers with unrealistic velocities that represent deeper/stiffer material.
For example, Fig. 1b shows a set of layer thicknesses obtained from 50 randomized V s profiles using the generic λ provided by Toro (1995) (Fig. 1a). Note that each circle symbol corresponds to a simulated layer plotted at its middepth. Because of the assumed statistical layering model, randomized V s profiles might have a different number of layers. In Fig. 1b, most of the randomized V s profiles needed 5-6 layers to reach a depth of 50 m, resulting in a total of 275 simulated layers for all 50 randomized V s profiles. As evident, the simulated layer thicknesses show extreme variability. For example, at depths between 20-40 m, the simulated layer thicknesses range between 1-40 m, even though the theoretical mean values are between 10-15 m. While this variability is not surprising given the nature of the assumed exponential distribution for layer thicknesses, it might not be realistic for many sites. In addition, the simulated top layer in some of the randomized profiles has a thickness as high as 40 m. These thick near-surface layers will consequently be assigned V s values solely based on their mid-depth, which will be as deep as 20 m in this case, and likely inconsistent with the near-surface V s.
Another challenge with V s randomization is the need to specify several input parameters which are site-specific, but not typically available. Assumptions are typically required to set parameters needed to define the rate parameter controlling the layering randomization (Equation 1), to define the statistical distribution of bedrock depth and its parameters, to calculate inter-layer correlations (Equations 2-4), and to set the σ lnV s (Equation 5). Because a wealth of site characterization information is not typically available, these parameters are often guessed or assumed based on generic relations derived for other sites, which introduces additional uncertainty into the site response predictions.
To investigate the consequences of these challenges, site response predictions with V s randomization have been compared with observations at vertical seismometer arrays by, for example, Field and Jacob (1993) at the Turkey Flat array in California, Teague et al. (2018) at the Garner Valley array in California, Kaklamanos et al. (2020) at 10 arrays from the Kiban-Kyoshin network (KiK-net) in Japan, Passeri et al. (2020) at the Mirandola array in Italy, and  at the Treasure Island and Delaney Park arrays in the United States. The collective results of these studies indicate that many of the randomized V s profiles represent unrealistic site conditions, which in turn produces site response predictions that are different from observations. In addition, Griffiths et al. (2016a,b), Teague et al. (2018), and Teague and Cox (2016) noted that many randomized V s profiles are inconsistent with the experimental site signature, which is defined by experimental dispersion data and estimates of the fundamental site frequency (f 0 ) from horizontal-to-vertical (H/V) spectral ratios of noise measurements (f 0,H/V ).
This extreme V s variability generally decreases the amplitude of the predicted site response (Rathje et al., 2010), causing a broadening and underprediction of site amplification in comparison to observations at downhole arrays, which is unconservative. Several procedures have been recommended to mitigate some of the challenges of V s randomization, such as screening randomized V s profiles based on the experimental site signature Teague et al., 2018;Griffiths et al., 2016b) or randomizing tts instead of V s (Passeri et al., 2020). We focus on tts randomization and its challenges next. Passeri et al. (2020) introduced tts randomization, which included several adjustments meant to resolve some of the challenges of traditional V s randomization. First, they proposed a method to mitigate unrealistic layer thicknesses by defining acceptable layer thickness bounds. Second, tts includes cumulative information about the V s profile, which avoids potential problems arising from randomizing about a discrete point in the V s profile (i.e., mid-depth of each layer). Third, tts randomization deals with a single random variable, which eliminates double counting of uncertainty that results from separate layering and V s randomization. However, despite potential improvements gained by applying these adjustments, tts randomization may pose additional challenges that require due care to avoid bias and unrealistic conditions in the simulated V s profiles. Below, we summarize the statistical framework of tts randomization and then discuss challenges and recommendations for improved implementation in the following section.

Statistical framework of tts randomization
Randomization can be performed on tts, which deals with a single random variable that is a function of both V s and layer thicknesses: where n is the number of layers, h i is the thickness of layer i, and V s i is its shear wave velocity. The proposed framework by Passeri et al. (2020) adopted most of the formulations and equations from Toro (1995), but applied them to tts. The proposed framework operates on three components: (1) a column model to randomize soil layer thicknesses and tts, (2) a bedrock model to vary bedrock depth and its V s, and (3) merging of the column and bedrock models. Passeri et al. (2020) discuss different approaches for handling V s profiles developed from invasive and noninvasive methods. When multiple V s profiles from inversion of surface wave data are available, site-specific input parameters can be estimated, but we only discuss the generic application herein, which requires a single base-case V s profile and aligns more closely with the Toro (1995) approach.
Although randomization is performed on tts, it is essential to represent the randomized tts profiles as V s profiles, which are the input to GRAs. For this purpose, the same layering model from Toro (1995) is adopted for tts randomization, but it is recommended that the parameters in Equation 1 (i.e., c 1 , c 2 , and c 3 ) be specified such that they result in a theoretical mean layer thickness with depth that is consistent with the measured thickness profile. In addition, Passeri et al. (2020) recommended bounding the simulated layer thicknesses by defining an acceptable range determined by upper and lower factors (e.g., ±50%) applied to λ, as shown in Fig. 3. Imposing these acceptable bounds is intended to eliminate unrealistic layer thicknesses. The same 275 simulations of layer thickness from Fig. 1b are also shown in Fig. 3b, but now color coded based on their location within or outside the acceptable range bounds. The results are shown prior to imposing the bounds to highlight unacceptable simulated layers. Passeri et al. (2020) imposed the acceptable bounds by truncating the simulated layer thicknesses at the bounds. Results and challenges from imposing the bounds in this way are discussed later in this article.
After randomizing layer thicknesses, and instead of modeling the variability in V s, the model operates on tts using the inter-layer correlation model (Equations 2-4) and first-order autoregressive function (Equations 5 and 6) developed by Toro (1995). Thus, the tts of adjacent layers are assumed to be positively correlated and the tts of each layer is assumed to be lognormally distributed with a median and a specified tts lognormal standard deviation (σ lntts ). Passeri et al. (2020)  Randomized tts profiles are then converted to V s profiles, which are referred to as column models. While this information is sufficient to define the column models, one other important detail is worth discussing. Passeri et al. (2020) did not explicitly detail how the median tts should be obtained. For example, V s randomization models commonly set the median V s for each randomized layer equal to the base-case V s at the mid-depth of that layer. However, if tts randomization is similarly performed at the mid-depth, any V s (or tts) information in the second half of each randomized layer will not be accounted for. Therefore, for tts randomization, the median tts used for the randomization of a specific layer should be set equal to the base-case tts at the bottom-depth of that layer. This is simply because tts includes cumulative information about the base-case V s profile, which guarantees capturing the full properties.
Bedrock properties are randomized separately using a joint distribution, where bedrock depth and V s are assumed to be lognormally distributed and positively correlated (i.e., higher bedrock V s tends to be associated with deeper bedrock), with the latter having a recommended generic correlation coefficient of 0.508. Instead of arbitrarily imposing a bedrock depth on a soil column model, an equal number of bedrock depths and soil columns are generated, which are then merged by pairing each bedrock depth (from a lognormal distribution) with the soil column that has the closest bottom depth (from the non-homogeneous Poisson process). This merging is Figure 3: The layering model adopted by Passeri et al. (2020) for tts randomization, which is the same as that developed by Toro (1995), but with the addition of enforcing acceptable bounds, as shown for: (a) interface rate, and (b) layer thickness. These bounds are simply defined as factors applied on λ (e.g., ±50% λ) to bound the layering randomization and avoid generating layers with unrealistic thicknesses. The same 275 layer simulation results from Fig. 1b are shown here, but color coded based on their location within or outside the acceptable range bounds. We note that results are shown prior to imposing the bounds.
intended to avoid unrealistic profiles, which could arise from enforcing a bedrock depth that is inconsistent with the layer interfaces.

Challenges and recommendations for implementing tts randomization
The work by Passeri et al. (2020) is indeed a step in the right direction of addressing some of the limitations of traditional V s randomization. However, in scrutinizing some of the adjustments introduced by Passeri et al.
(2020), we show that they give rise to three additional challenges that require further improvements. Particularly, these challenges pertain to: (1) bounding simulated layer thicknesses, (2) merging column and bedrock models, and (3) randomizing tts. In the next subsections, we use Monte Carlo simulations and/or real data to demonstrate these challenges and recommend improved implementations that are grounded in fundamental principles of statistics and probability.

Bounding simulated layer thicknesses
Despite the potential for unrealistic variability in layer thicknesses resulting from the statistical model developed by Toro (1995) (i.e., non-homogeneous Poisson model), as demonstrated in Fig. 1, this model continues to be commonly used in practice and research, and has also been directly incorporated into at least one site response software package, Strata (Kottke and Rathje, 2009). In addition, Passeri et al. (2020) adopted the same model but with modifications, as previously discussed. Thus, rather than developing a different statistical model for randomizing layer thicknesses, we strictly adhere to the work of Toro (1995) and Passeri et al. (2020), but with improved implementation. Passeri et al. (2020) proposed truncating simulated layers at acceptable bounds set as a symmetric factor on λ (e.g., ±50% λ). We show that imposing and defining the acceptable bounds in this way can lead to bias by clustering layer thicknesses at the upper and lower acceptable boundaries. Therefore, to mitigate this bias, and rather than incorporating additional parameters, we simply provide recommendations to impose and define the acceptable bounds in a different, unbiased manner.

Imposing acceptable layer thickness bounds
The resulting layer thicknesses from applying the method adopted by Passeri et al. (2020), which we refer to as Method 1, are shown in Fig. 4a. Note that these data points are based on the 275 simulated layers shown in Fig. 1 and Fig. 3, although the number of layers will differ after imposing the bounds. In addition, a histogram of the final layer thicknesses after imposing the bounds is shown in Fig. 4c. Because layer thicknesses vary with depth, the histogram values are normalized by the theoretical mean thickness for better illustration. It is evident that forcing unacceptable layers to be equal to the bounds causes many of the simulated layer thicknesses to be clustered at the bounds, and more so at the lower thickness bound (Fig. 4c). This is an artifact of the assumed distributions and how the bounds are imposed. Recall from Fig. 3b that many of the simulated layer thicknesses fell outside the acceptable bounds, particularly below the lower thickness bound, which is due to the high probability density function of the exponential distribution at low values (Fig. 2). Consequently, Method 1 adopted by Passeri et al. (2020) could be improved to produce a better distribution of the final layer thicknesses.
To overcome this challenge, we propose a simple fix by rejecting simulated layers that are outside the acceptable bounds and repeating their simulation until they fall within the bounds. Results from applying this method, which we refer to as Method 2, are shown in Fig. 4b and d. In contrast to the results discussed earlier, the final layer thicknesses in this case are not clustered at any of the bounds, but rather have a good coverage range that is more consistent with the assumed exponential distribution. Consequently, care should be taken in how the acceptable layer thickness bounds are imposed on the randomizations to avoid any bias in the final layer thicknesses. It is highly recommended that unacceptable simulated layer thicknesses be rejected and re-simulated, as proposed in this article, rather than forced to the boundaries. We note that, on average, unacceptable layers in Fig. 4b had to be re-simulated twice until an acceptable thickness was obtained, but this number will vary depending on the width of the bounds.
Despite the better results achieved by Method 2, it is evident that the final layer thicknesses are biased to thicker layers on aggregate. As shown in Fig. 4c, the mean of the simulated layer thicknesses is higher than the theoretical mean when imposing the bounds using Method 1. This bias remains even when Method 2 is applied (Fig. 4d).
The overall bias to thicker layers is governed by the boundary values used by Passeri et al. (2020) (i.e., ±50% λ), which lead to unsymmetrical thickness bounds, as discussed next. (1) truncating unacceptable layers at the boundaries (i.e., forcing all simulations that are outside the bounds equal to the values of the bounds), which is adopted by Passeri et al. (2020), and (2) re-simulating unacceptable layers, which is the method proposed in this article for improved implementation. For each method, we show: (a, b) the final layer thicknesses with depth, and (c, d) a histogram of the final normalized layer thicknesses. Because the thickness increases with depth, the histogram values are normalized by the theoretical mean thickness for better illustration.

Defining acceptable layer thickness bounds
The assumed statistical distribution of layer thickness (exponential) is unsymmetrical and long tailed (Fig. 2), which makes defining appropriate bounds nontrivial. Passeri et al. (2020) used symmetric bounds on interface rate, as defined using factors on λ (e.g., ±50% λ), which we will refer to as rate-symmetric bounds. To demonstrate the effect of these bounds on the acceptable range, Fig. 5 shows the statistical distribution of layer thickness at a particular depth with λ = 0.1 (refer to Fig. 1 and Fig. 2). The acceptable thickness range based on ±50% λ at this depth, as highlighted in the green region, is 6.7-20 m, which is unsymmetrical relative to the theoretical mean thickness of 10 m. This asymmetry causes the mean thickness of the acceptable range (green dotted line in Fig. 5) to be higher than the theoretical mean, introducing an overall bias to thicker layers. Thus, if engineers develop a site-specific λ based on available data, and then bound the simulations using rate-symmetric bounds, the final layer thicknesses, on aggregate, will be greater than, and inconsistent with the site-specific λ.
To ameliorate this bias, the acceptable bounds should be specified such that the mean thickness of the acceptable range is consistent with the theoretical mean, which we refer to as mean-consistent bounds. Since the ratesymmetric bounds are biased to thicker layers, then the acceptable range should either be extended further toward thinner layers, or reduced to limit thicker layers. Fig. 5 shows, in red, the acceptable range based on limiting thicker layers (i.e., using rate bounds of -30% to +50% λ instead of ±50% λ). By appropriately limiting thicker layers from the rate-symmetric bounds, the mean thickness of the acceptable range becomes consistent with the theoretical mean (red dotted line in Fig. 5).  Exponential ( = 0.1) Figure 5: Alternative approaches to define acceptable bounds for the layering model and their corresponding thickness ranges at a particular depth with λ = 0.1 (refer to Fig. 2). The wider acceptable range shown in green is based on symmetric interface rate bounds, as defined using symmetric factors on λ (e.g., ±50% λ), which were used by Passeri et al. (2020); we refer to these as rate-symmetric bounds. The narrower range shown in red, which we propose as an improved approach for defining acceptable layer thickness bounds, has a mean thickness that is consistent with the theoretical mean; we therefore refer to these bounds as mean-consistent bounds. Fig. 6 shows the resulting layer thickness simulations from using the rate-symmetric ±50% λ bounds from Passeri et al. (2020) (same results as Fig. 4b and d) and the mean-consistent -30% to +50% λ bounds proposed herein. Both bounds were imposed on the simulations by re-simulating unacceptable layers (i.e., Method 2 in Fig. 4). The layer thicknesses based on the mean-consistent bounds ( Fig. 6b and d) have a simulated mean that is in perfect agreement with the theoretical mean, unlike the results based on the rate-symmetric bounds, which are biased to thicker layers ( Fig. 6a and c). This provides further evidence supporting that rate-symmetric bounds result in an overall bias to thicker layers, which can be corrected by using the proposed mean-consistent bounds.
The -30% to +50% λ is just one set of mean-consistent bounds, and Fig. 7 shows other combinations of upper and lower mean-consistent bounds. Also shown are the rate-symmetric bounds used by Passeri et al. (2020), with the circle symbols being the specific cases demonstrated in Fig. 5 and Fig. 6. Fig. 7 thus provides engineers with an easy way to define the acceptable bounds recommended in this article, by simply choosing any point along the relationship for mean-consistent bounds. The -30% to +50% λ bounds (red circle symbol in Fig. 7) are believed to be appropriate for sites with moderate spatial variability in layer thicknesses, and wider bounds could be used for sites with high spatial variability in layer thicknesses. However, we do not recommend using bounds beyond -50% to +150% λ, because otherwise, the benefit of applying the bounds is negligible. More research on this topic is required to establish typical bound values for different site conditions.

Merging column and bedrock models
After generating soil layer thicknesses that are within the acceptable bounds, which are referred to as column models, an equal number of bedrock depths are independently simulated using an assumed lognormal distribution. Then, each bedrock depth is merged with a column model that has the closest bottom depth. This merging is intended to avoid enforcing a bedrock depth that is inconsistent with the layer interfaces. However, an exact match between the simulated bedrock depths (from a lognormal distribution) and layer interfaces in the column model (from a Poisson model) is virtually impossible. Bedrock depth is enforced on the column model after merging, and thus, the thickness of the deepest column layer will be altered. Depending on the magnitude of this alteration, the thickness of the deepest column layer might become outside the acceptable bounds, which defies the goal of bounding layer thicknesses.
To demonstrate this undesirable effect, Fig. 8a, c, and e present simulation results after merging 50 column and 50 bedrock models. To investigate the effect of the width of the acceptable soil layer thickness range on merging, results are shown for three different sets of mean-consistent bounds, with Fig. 8a having tight bounds (less thickness variability) and Fig. 8e having wide bounds (more thickness variability). Each circle symbol corresponds to a layer in a column model, with each column model having multiple layers. All simulated layers fell entirely within the acceptable bounds prior to merging with bedrock models. Bedrock models were simulated using a lognormal distribution with a median of 45 m and a depth to rock lognormal standard deviation of 0.05.   Figure 8: Resulting layer thicknesses after merging column and bedrock models using two different methods: (a, c, e) Method 1 simulating column models only once, which is implemented by Passeri et al. (2020), and (b, d, f) Method 2 re-simulating column models with unacceptable deep layers, which is proposed in this article. Results are shown for three sets of mean-consistent bounds to investigate the effect of the width of the acceptable layer thickness range on merging. Note that although shallower layers are unaffected by merging and remain within the acceptable bounds, all layers for models with unacceptable deep layers are highlighted in magenta.

Rate-Symmetric Bounds ( 50% )
It is evident from Fig. 8a, c, and e that merging column and bedrock models without checking the thickness of the deepest layer after merging, which is the method implemented by Passeri et al. (2020) and which we refer to as Method 1, could potentially cause some of the deepest layers to fall outside of the acceptable thickness bounds. Although shallower layers are unaffected by merging and remain within the acceptable bounds, all layers for models with unacceptable deep layers are highlighted in magenta (Fig. 8a, c and e). The results show that merging is particularly more critical for tighter bounds (e.g., -30% to +50% λ in Fig. 8a), simply because even small changes in the thickness of the deepest layer after merging could cause it to fall outside the bounds. In Fig. 8a, 12 out of 50 models had unacceptable deep layers after merging. However, with wider bounds (e.g., -50% to +150% λ in Fig. 8e), there is a greater margin within which the thickness of the deepest layer above bedrock could be altered and remain acceptable. In this case, only a single model had an unacceptable deep layer after merging.
In an effort to ensure full consistency between the randomized layer thicknesses, their acceptable bounds, and the randomized bedrock depths, we simply propose re-simulating a new set of column models and then merging them only with the bedrock models that resulted in unacceptable deep column layers, which we refer to as Method 2. Results from applying this method are shown in Fig. 8b, d, and f, in which the re-simulated models are highlighted in dark blue. These results assert that Method 2, proposed herein, yields layers that are entirely within the acceptable bounds, even after merging. Therefore, it is recommended that column model with unacceptable bottom layer thicknesses after merging be completely rejected (i.e., all layers), entirely re-simulated, and then re-merged with the remaining bedrock models. We note that it might be necessary to re-simulate unacceptable column models more than once until the resulting model, after merging, is entirely within the acceptable bounds. This can be easily automated.

Randomizing tts
The final step of tts randomization is simulating tts for each layer in the soil column from an assumed lognormal distribution, which is defined by a median tts, σ lntts , and ρ. As discussed earlier, the median tts of each layer should be set equal to the bottom-depth base-case tts. The simple rationale behind this is that tts is a cumulative parameter, and thus the bottom-depth value captures the base-case properties to the depth of interest. In this subsection, we focus on challenges and recommendations for specifying appropriate σ lntts and ρ values.
The extent of material property variability incorporated into GRAs for tts randomization is controlled by σ lntts . This is analogous to σ lnV s for V s randomization, which has been shown to have a significant effect on the predicted seismic site response (Rathje et al., 2010;Boushehri et al., 2020). Passeri et al. (2020) suggested that σ lntts is depth independent and recommended values between 0.01-0.02 based on V s profiles derived from the inversion of active-source surface wave data. Thus, this recommended σ lntts range reflects uncertainties associated with the non-uniqueness of surface wave inversion. However,  hypothesized, based on comparisons of predicted site response from tts randomization with recorded ground motions at several borehole array sites, that using the recommended σ lntts of 0.02 underestimates spatial variability. In addition, some researchers have hesitated using values as low as 0.01-0.02 in practice. In this article, we use numerous borehole and surface wave V s measurements from a real dataset along with Monte Carlo simulations to further investigate appropriate σ lntts values.  presented preliminary results on the influence of σ lntts from Monte Carlo simulations. While our results indicate the need for higher σ lntts , we show that this might yield unrealistic V s profiles, with excessively high or even negative V s values, which must be mitigated by increasing ρ.

Investigations of σ lntts Based on Real Data
We first compare the recommended σ lntts range of 0.01-0.02 with σ lntts associated with in-situ measured V s profiles from a site located in the East Gulf Coastal Plain of western Georgia, USA (Damm et al., 2021), which we refer to as the Greenfield site. A comprehensive deep V s profile site characterization project was conducted at the Greenfield site using several invasive and noninvasive tests, resulting in a total of 27 V s profiles. The invasive methods included 10 seismic cone penetration tests (mostly in the top 30 m) and 8 PS suspension logging tests, whereas the noninvasive tests resulted in 5 V s profiles from spectral analysis of surface waves (SASW) surveys and 4 median V s profiles from combined multichannel analysis of surface waves (MASW) and microtremor array measurements (MAM). All tests were performed within an area of approximately 1 km 2 .  (Fig. 9b), which is consistent with low-to-moderate variability (Stewart et al., 2014). The low spatial variability has also been confirmed through consistency between H/V curves from ambient noise recordings at 38 distinct locations across the site. While not shown here for brevity, the H/V curves are nearly identical, providing further evidence of the site's low spatial variability. However, despite the apparent low spatial variability at the Greenfield site, the σ lntts of the 27 V s measurements is in the range of 0.05-0.10 across the entire depth of 150 m, which is 5 times as much large as the range recommended by Passeri et al. (2020). Again, we would like to emphasize that the 0.01-0.02 range recommended by Passeri et al. (2020) mostly reflects uncertainties associated with the non-uniqueness of surface wave inversion. However, in the context of PSHA and ground response analyses, randomization approaches are used to account for V s aleatory variability. While we only present measurements from a single dataset, albeit a high quality one, the considerably higher σ lntts at a site where there is evidence of low spatial variability provides strong experimental support to the conclusions of , that to account for spatial variability, a σ lntts of 0.01-0.02 is perhaps too low, at even relatively homogeneous sites.

Investigations of σ lntts based on simulations
Based on this field evidence, we next assess the effects of using higher σ lntts on tts randomization results. We investigate the influence of three different σ lntts values of 0.02, 0.05, and 0.10 on the simulated variability in V s and on subsequent site response predictions. The assessments are performed at the Treasure Island Downhole Array (TIDA) so that the GRA predictions from tts randomization can be compared to actual ground motion observations. The TIDA site is located in a seismically active region of northern California, approximately 5 km northeast of San Francisco, 12 km west of the Hayward fault, and 20 km east of the San Andreas fault. The site has been characterized through several testing programs, and existing data show that the near subsurface is horizontally stratified, with relatively low spatial variability (Cheng et al., 2021;Papadopulos and Eliahu, 2009;Tao and Rathje, 2019). The site has a V S30 of 165 m/s, which classifies as Site Class E, soft soil (ASCE, 2022).
An invasively measured V s profile from PS suspension logging at TIDA (Graizer et al., 2004) was used to develop a base-case tts profile for randomization. Then, for each investigated σ lntts (i.e., 0.02, 0.05, and 0.10), a set of 250 V s profiles were simulated using the recommended tts randomization procedures discussed above for bounding layers and merging column and bedrock models (i.e., Methods 2). The simulated V s profiles from tts randomization are shown in Fig. 10. To better interpret the significance of different σ lntts values, the simulated V s profiles from tts randomization were used to compute σ lnV s values. As shown in Fig. 10, the σ lnV s of the simulated V s profiles are compared to a typical range of 0.10-0.30, which was chosen based on common values used in practice and in research (Stewart et al., 2014;Rodriguez-Marek et al., 2014;Garofalo et al., 2016a).
It is evident from Fig. 10a, c, and e that the median V s profiles are similar to the base-cases, indicating that the V s profiles from tts randomization capture the assumed baseline models, on average, regardless of σ lntts . However, increasing σ lntts causes an increase in the variability of the simulated V s profiles and their resulting σ lnV s . This observation is not surprising because σ lntts controls the magnitude of the modeled variability. However, Fig. 10b shows that using σ lntts of 0.02 results in σ lnV s generally below 0.10, which is less than the lower bound of the typical range used in practice and research. These simulation results lend further support to the hypothesis that σ lntts larger than that recommended by Passeri et al. (2020) might be required to more appropriately model subsurface variability. For example, using a σ lntts of 0.05 at the TIDA site yields σ lnV s values below 0.20, but mostly within the typical range (Fig. 10d). Thus, values closer to 0.05 might be more appropriate for sites with low-to-moderate variability, like the TIDA site. Although a σ lntts of 0.10 (Fig. 10f) appears to yield higher σ lnV s values that are still consistent with high variability sites, this value of σ lntts results in a few unrealistic V s models with very high V s in some cases (Fig. 10e); this aspect is further addressed later in the article.  To assess the influence of σ lntts on the predicted site response, linear-viscoelastic, theoretical transfer functions (TTFs) were calculated for each simulated V s profile between the surface and deepest triaxial ground motion accelerometer at the TIDA site, which is located at 122 m. TTFs were calculated using the closed-form solution for damped, layered soil over elastic bedrock. The individual TTFs were then averaged to obtain a representative mean TTF. The required input soil properties include V s, mass density, and small-strain damping ratio (D min ). The damping profiles were assigned using the laboratory-based D min (Darendeli, 2001) using estimates of in-situ confining stresses obtained from typical soil unit weights and assumed plasticity indices of 25 for clay layers and 0 for sand layers, similar to Tao and Rathje (2019).
Next, to compare the site response predictions from different σ lntts levels with observations, a suite of 36 recorded low-amplitude ground motions (peak ground accelerations < 0.01 g) was used to compute a median empirical transfer function (ETF) between the deepest accelerometer and the ground surface. For each horizontal component of all selected ground motions, an ETF was computed as the ratio of the Fourier amplitude spectra (FAS) of the horizontal surface accelerations divided by the horizontal rock accelerations. Then, the lognormal me-dian and the lognormal standard deviation (σ lnET F ) of all individual ETFs from the two horizontal components of recorded earthquake events were computed. By including only low-amplitude motions, we ensure that challenges associated with modeling nonlinear soil behavior do not interfere with our attempts to understand the effects of spatial variability. It should be noted that since the ETFs represent the surface-to-rock ratio in which the rock time histories include both upgoing and downgoing waves, the TTFs are modeled to include those effects using what is referred to as the within boundary condition at depth. For more details, please refer to ; Hallal and Cox (2021). Fig. 10g compares the mean TTFs computed from tts randomization using the three different levels of σ lntts with the median ETF. Increasing levels of variability cause a decrease in the predicted mean TTF amplitude, which is in line with observations by Rathje et al. (2010) for traditional V s randomization. While the mean TTF from using σ lntts = 0.02 is the only one that clearly captures all amplification modes observed in the ETF, there is persistent amplification overprediction at frequencies below 5 Hz (Fig. 10g). Using a σ lntts of 0.05 gives a better match to the ETF at TIDA over the entire frequency range, whereas a higher value of 0.10 significantly over-flattens amplification at higher frequencies. As was noted above, the TIDA site is characterized by low spatial variability. In fact, due to the nature of its hydraulic fill construction, it is known as a prototypical 1D site. Even at this low spatial variability site, a σ lntts = 0.02 value appears to be too low to accurately capture the amplitude range of the fundamental mode, and a higher value of approximately 0.05 appears to better represent spatial variability within the tts randomization framework. While a σ lntts of 0.10 might not be appropriate in this case, as revealed by the underprediction of the observed site response at all amplification modes (Fig. 10g), this higher value might be appropriate for sites with higher spatial variability. However, Fig. 10e shows that this also generates some unrealistic V s profiles, which are examined hereafter.
A concerning observation when implementing tts randomization with high σ lntts is generating unrealistic V s models. This is evident in several of the V s profiles in Fig. 10e that have soil layers in the top 60 m with V s > 800 m/s. Additionally, some V s profiles at depths between 60-80 m have layers with V s > 2000 m/s, despite the fact that the base-case profile has V s < 400 m/s over this same depth range. One of the simulated V s profiles even has a layer with negative V s between approximately 65 and 85 m.
To help understand the reason underlying such highly variable and unrealistic V s profiles, Fig. 11 is used to more closely examine one of the most erratic tts-randomized V s profiles from Fig. 10e. In addition to the V s profile, Fig. 11 shows the base-case and randomized tts profile and the inter-layer correlation coefficients. Recall that tts is randomized at each layer interface from a lognormal distribution, such as those shown in gray in Fig. 11b, with the assumption that successive randomized tts values are positively correlated. Each assumed tts distribution has a median equal to the base-case tts (black circle symbols in Fig. 11b) and a specified σ lntts (0.10 in this case). A simulated V s profile from tts randomization is obtained from the slopes of the lines connecting successive randomized tts values (e.g., the red square symbols in Fig. 11b). Because tts is cumulative (Equation 7), it should strictly increase with depth. However, this condition is not enforced in the tts randomization framework. For example, consider the region highlighted in orange in Fig. 11b, within which the successive tts distributions at approximately 65 m and 88 m overlap. Overlapping distributions at successive depths allow for simulated tts values at these depths to possibly be equal to one another, or even decrease with depth; the latter would result in a negative V s. For the example V s profile under consideration, the randomized tts values at approximately 65 and 88 m are quite similar, resulting in a nearly vertical tts slope (Fig. 11b), and consequently an unrealistically high V s, in excess of 3000 m/s (Fig. 11a).
Problems like this are a consequence of two main factors: (1) wider tts distributions at depth (i.e., higher absolute variance for the same σ lntts due to larger base-case tts at greater depths, as evident in Fig. 11b), and (2) generally closer tts distributions at depth (i.e, lower inter-layer travel times due to typically higher velocities at depth). Both of these factors generally occur simultaneously at greater depths, exacerbating the possibility of unrealistic results. While all log tts distributions in Fig. 11 have the same σ lntts , the variance of the tts distributions (not in log space, but similar to those plotted in Fig. 11b) changes. Particularly, it is a fact that the variance of lognormal distributions is proportional to the mean value. Consequently, because the mean tts increases with depth, the assumed lognormal tts distributions become wider, even when using the same σ lntts . This is clearly demonstrated in Fig. 11b, which shows that the tts distribution below 80 m has a higher variance than those at shallower depths, although all have the same σ lntts . If the tts distributions have very different base-case values (i.e., are far apart from each other in time space), then using a σ lntts = 0.10 will not necessarily result in an overlap of successive distributions, as evident at shallow depths in Fig. 11b. However, the results are exacerbated by the fact that the subsurface gets stiffer with depth, meaning that tts increases at a much slower rate. Consequently, successive tts distributions overlap because they are wider and their locations are closer in time, both of which generally occur simultaneously at greater depth. This overlap poses a challenge to tts randomization. The unfavorable results specifically arise when the simulated tts at successive depths have very different Z-scores (i.e., number of standard deviations from the base-case). For example, the randomized tts at approximately 65 m in Fig. 11b is comparable to its base-case (Z-score ≈ 0), but at approximately 88 m it is considerably lower (Zscore ≈ -2). This causes the randomized tts slope to be near vertical, resulting in an unrealistic V s layer. Recall that the tts of successive layers are assumed to be positively correlated, meaning that the simulated Z-scores of successive layers have similar tendencies. The value of ρ (Equation 2) controls the magnitude of this tendency and thus, it is likely that with higher ρ, successive Z-scores would be more similar, and as a result, unrealistic V s values would be mitigated. As shown in Fig. 11c, the parameters recommended by Passeri et al. (2020) for the correlation model result in ρ below 0.6 for the sample simulated profile. However, as we argue next, because tts is a cumulative parameter, values at successive depths should be more highly correlated. Below, we revisit the Greenfield site and then use Monte Carlo simulations to examine appropriate ρ values.

Investigations for ρ based on real data
Fig. 12 compares the ρ of the tts profiles from the V s measurements at the Greenfield site to those calculated using the model adopted by Passeri et al. (2020) (Equations 2-4). The functional form of this model is based on Toro (1995), but we use the parameters recommended by Passeri et al. (2020) for tts randomization. As shown in Equations 2-4, ρ is assumed to be depth-and thickness-dependent. Consequently, to investigate the influence of inter-layer separation distance (i.e., t in Equation 4) on the Greenfield site data, ρ was computed for thicknesses of 2, 10, and 20 m, as shown in Fig. 12. We note that ρ relates the tts Z-scores at successive depths (similar to Equation 6). Thus, to back-calculate the empirical ρ values from the field V s profiles, tts Z-scores were first computed for each tts profile. Then, ρ values were computed between all Z-scores using separation distances of 2, 10, and 20 m (Fig. 12b). The 2 m results are better suited for the near surface, where layers are generally thinner, whereas those from either 10 m or 20 m are better suited for greater depths.
The high tts inter-layer correlation is strikingly evident in Fig. 12b, which shows that a tts profile that is above the median (i.e., positive Z-score), remains above the median with almost the same Z-score, and vice versa. Quantitatively, the ρ for 2 m separation distances based on the Greenfield site data is above 0.85 at all depths, whereas that based on Passeri et al. (2020) remains below 0.8 over nearly the entire depth (Fig. 12c). For higher separations distances, ρ is low in the near surface for the Greenfield site data, which is not surprising because thicker layers generally do not occur at shallow depths. However, as depth increases, the Greenfield site data ρ values for both 10 and 20 m rapidly increase beyond 0.8. All three separation distances show similar trends of high ρ at depths greater than 40 m that approach 1.0, which is much greater than the ρ values from Passeri et al. (2020). Fig. 12 strengthen our argument that, since tts is cumulative, ρ is higher than results from the generic model. Therefore, in what follows, we will investigate the effect of using higher ρ on tts randomization results.

Investigations for ρ based on simulations
Fig . 13 shows the V s profiles simulated from tts randomization at the TIDA site with σ lntts = 0.10 and using the generic ρ model and higher ρ obtained by adding fixed values of 0.2 and 0.4 to the generic ρ. We note that we believe σ lntts = 0.10 might be too high for the TIDA site, but is used herein to examine whether increasing ρ mitigates unrealistic V s profiles from tts randomization. The ρ values from the generic model ranged between 0.45-0.65, depending on the layer thicknesses and depth, and those increased by 0.2 thus ranged between 0.65-0.85, whereas those increased by 0.4 ranged between 0.85-1.00. Fig. 13 demonstrates that increasing ρ mitigates unrealistic V s profiles, and consequently, decreases σ lnV s , even though all results are based on the same σ lntts . The decrease in σ lnV s with increasing ρ is a manifestation of not allowing profiles with unrealistic V s to be generated from overlapping tts distributions.
The most intriguing observation from Fig. 13 is that although ρ influences the simulated V s profiles and their σ lnV s , the results suggest that it has a minor effect on the mean TTFs (Fig. 13g). Except for slight differences in the mean TTFs near the fundamental mode, the results are comparable for all ρ and still provide poor site response estimates at higher frequencies. One possible explanation is that all simulations were based on the same σ lntts , which is the key parameter controlling the extent of variability incorporated into GRA. Nonetheless, even if the mean TTF appears to be unaffected by ρ for this particular case, it is important that this parameter be appropriately specified to avoid generating unrealistic V s profiles. The appropriate value of ρ is site-specific, but it appears that values in excess of 0.7-0.8 are more appropriate when implementing tts randomization, particularly for deeper layers. In order to maintain depth-dependent values, engineers can either increase ρ from the generic model by approximately 0.2-0.3 or adjust the model parameters in Equations 2-4 to guarantee a higher ρ. Sitespecific tts ρ values can also be developed if multiple V s measurements are available (refer to Fig. 12).
While using higher ρ helps mitigate unrealistic V s profiles with excessively high or negative V s values, we strongly recommend that when performing any statistical randomization, a robust screening criteria be invoked. For example, Griffiths et al. (2016a,b), Teague et al. (2018), and Teague and Cox (2016) proposed and examined a strategy to exclude randomized V s profiles that fail to reasonably fit the experimental "site signature" (i.e., broadband experimental surface wave dispersion data and fundamental site frequency estimates from H/V curves).  implemented similar site signature screening for V s randomization results at two vertical seismometer arrays and demonstrated the significant improvements that can be achieved to site response predictions. Screening unrealistic profiles could also be potentially performed using the fundamental site frequency estimated from H/V curves and V S30 values, which are more readily available than broadband experimental surface wave dispersion data. However, this requires further examination to establish best screening practices.

Discussion and future developments
The assessments presented above demonstrate some of the pressing challenges associated with tts randomization, for which we provided guidelines and recommendations for improved implementation. A summary of the recommendations and improvements introduced in this article is provided in Table 1. We would like to emphasize that these recommendations do not require any additional model parameters or constraints. For the sake of simplicity, we opted to strictly adhere to the work of Passeri et al. (2020), but with improved implementation. While the framework presented herein forms the basis for an improved methodology to incorporate variability into site response, this research raises additional questions in need of further investigation.  (2020), helps mitigate unrealistic layers, this requires two additional input parameters: upper and lower acceptable bounds. In addition, the assessments undertaken in this article have revealed that the procedures through which these bounds are imposed and defined can introduce unintended bias into the distribution and mean of the final simulated layer thicknesses. While our recommendations for re-simulating unacceptable layers and using mean-consistent bounds overcome this bias, this is still based on the adjustments of Passeri et al. (2020). The prospect of being able to randomize layer thickness more simply and realistically using an alternative statistical model that does not need to be bounded in a specific way should serve as an incentive to future research.
Second, while our assessments at the Greenfield site and the TIDA site provide evidence of the need for using higher σ lntts than 0.02 when accounting for spatial variability, even for sites with low variability, both evaluated sites have V S30 values consistent with soil conditions. Future work should focus on developing appropriate σ lntts recommendations for different site conditions and for sites with differing degrees of variability. On a wider level, research is also needed to determine whether σ lntts is depth independent, as suggested by Passeri et al. (2020). Results from the Greenfield site (Fig. 9d) suggest that σ lntts is higher at shallow depths, which is consistent with the higher near-surface V s variability reported in other studies (Garofalo et al., 2016a,b).
Third, we have shown that using higher ρ for tts is critical to avoid unrealistic V s profiles, particularly at depth and at stiffer sites, where tts distributions might overlap. Because the data we presented herein and all models adopted in the literature argue that ρ is likely depth-dependent, we sought to achieve a higher depth-dependent ρ by simply increasing the values computed from the generic model from Passeri et al. (2020) by a fixed value between 0.2-0.3. However, developing a new set of correlation model parameters is an important aspect to pursue. While recommendations for higher σ lntts and a higher ρ values are provided, it is clear that engineers must put more thought into the site-specific determination of σ lntts and ρ when implementing tts randomization.
Fourth, based on preliminary assessments the authors have performed, randomizing the time-averaged V s instead of tts is another potential improvement. Passeri et al. (2020) argued that randomization can be performed on tts, or equivalently, on the time-averaged V s. However, the authors discussed their framework only in terms of tts. Nonetheless, the time-averaged V s can help mitigate unrealistic V s profiles resulting from overlapping tts distributions without using higher ρ. This is because overlapping time-averaged V s distributions are not necessarily unrealistic, and the slope of the time-averaged V s between two depths can be vertical or backward. These would correspond to a V s equal to the time-averaged V s or lower, respectively, rather than excessively high or negative V s when performing tts randomization. More research on this topic is needed.

Conclusions
It is standard practice to model the variability in V s profiles in 1D GRA for critical facilities and major infrastructure using a stochastic framework, such as that developed by Toro (1995). However, several studies have expressed concerns about the blind application of this approach. Randomization can be alternatively performed on tts, which deals with a single random variable with depth, as proposed by Passeri et al. (2020). Nonetheless, despite potential improvements from implementing tts randomization, adoption of this new framework in practice is limited due to the fact that it has not been thoroughly validated yet.
In this article, we scrutinized the outstanding challenges of tts randomization and devised strategies for improved implementation that are grounded in fundamental principles of statistics and probability. These challenges and recommendations particularly pertain to: (1) methods used for bounding simulated layers, (2) approaches for merging column and bedrock models, and (3) appropriate values of σ lntts and ρ for randomizing tts. The positive impacts of this work are better understanding of tts randomization and more confidence in adopting this approach. The guidelines for improved implementation of tts randomization, combined with the original framework of Passeri et al. (2020) and the future directions that were outlined in the Discussion and Future Developments section, will work towards better modeling and propagation of V s variability in 1D GRA.