Modeling the Relationship Between SPT-N Value and Compression Index (Cc) Using Copula Theory for Thi-Qar Clay Soil

Understanding geotechnical properties such as the coefficient of compressibility (C_c) is essential to foundation design. The calculation of the consolidation settlement depends on the quality of information available from the site and laboratory tests. The standard penetration test (SPT) is one of the most widely used field tests in Iraq and around the world to determine soil resistance. This is because it is simpler and less expensive compared to laboratory tests like the uniaxial compression test. The SPT test was developed in the United States in 1927 [1]. Because of its antiquity and widespread use, studies have presented numerous experimental equations based on its results to evaluate soil properties, such as the correlation of standard penetration number (SPT-N) with relative density (Dr) of sand, undrained cohesion (cu) with SPT-N for clay soil [2], the consistency index (CI), and SPT-N [3]. In many studies, correlations between SPT-N values and sand friction angle have been proposed; a study [4] provided a critical review of these correlations, and a new correlation was proposed based on data collected from the literature. Correlation between SPT-N and the modulus of elasticity (E) has been performed for granular soils [5]. In addition, several studies have examined the correlations between shear wave velocity and SPT-N [6], reviewing the most important of these relationships.

Similarly, the C_c is a crucial metric for assessing the settlement of foundations in clay soils subjected to stresses. Precise and expensive tests in the lab are usually needed to calculate it. Many researchers have taken numerous approaches worldwide to create empirical correlations for compressibility with easily obtainable index properties like void ratio, liquid limit (LL), plasticity index (PI), and natural moisture content [7]. A helpful study for Nasiriyah soil also correlates the compressibility index with Aterberg limits (LL, plastic limit (PL), and PI) [8].

When resources are limited for many projects, it is important to have equations that can estimate soil properties without complex laboratory tests. Empirical relationships about the compressibility index of soil or SPT parameters are thus essential. The compressibility factor has rarely been estimated directly from field test data such as the SPT test. However, some attempts have been made to estimate C_c from SPT data [9]. The researcher used linear regression analysis to determine the relationship between the SPT and the clay soil C_c. The correlation between the soil compression index and SPT-N value of stiff clay soil was linear, as the compression index decreases with increasing SPT-N value, and an unreliable correlation for soft clay. A contrary conclusion was observed for samples included that were collected from different coastal areas, where higher SPT values were associated with higher C_c values [10].

In the domain of probabilistic modeling of geotechnical variables, there exists a necessity to denote the interdependence between parameters, particularly when these variables are not entirely autonomous. Some traditional methods have been applied to accomplish this purpose, such as the Rosenblatt transform [11] and the Nataf transform [12]. However, both approaches face significant restrictions that constrain their proficiency and adaptability. The Rosenblatt transform, alternatively, demands exact foreknowledge of the joint probability density function (PDF), a task that is often extremely complex or impractical, particularly when addressing incomplete data [11]. The Nataf method assumes only a linear relationship between variables and can only be applied to distributions belonging to the normal distribution family, like the Gaussian distribution [13]. This renders it ineffective when the relationship between variables is nonlinear or asymmetrical, as is regularly the case in natural soil properties.

In contrast, the copula theory is a more advanced and flexible option, as it allows for the detachment of marginal distributions from the dependency structure between variables, permitting the precise representation of intricate statistical relationships. The copula does not necessitate a linear relationship and allows the joint cumulative distribution function (CDF) to be constructed independently of marginal distributions [14]. For this reason, the copula has become a widely used instrument in the fields of risk examination, reliability analysis, and geotechnical engineering, particularly when addressing incomplete data. Numerous studies have used copula theory to generate the joint probability distribution of interrelated geotechnical parameters with limited probability distribution information. Uzielli and Mayne [15] used Copula theory to analyze and represent the dependence of load-displacement model parameters to generate samples that account for parametric uncertainties in model inputs, which were then used in a Monte Carlo Simulation (MCS) to assess the impact of these uncertainties on settlement estimation. Wu and Xin [16] used a bivariate copula to simulate the cement fly ash gravelly (CFG), in pile load-settlement behavior regression curve parameter dependence. When probability information is missing, the copula function creates a joint distribution of rockfill material nonlinear strength coefficients [17]. The vine copula was employed in a previous study [18] to model the dependence structure between consolidation parameters such as initial void ratio, compressive index, and saturated density of soft clayey soil, which contributed to providing a more realistic and reliable analysis of consolidation settlement. The copula theory was used in a previous study [19] to model the dependency between uncertain geotechnical design parameters in reliability-based shallow foundation design.

Although the SPT is a common in situ test employed worldwide, and although the C_c is a primary parameter needed for settlement prediction, the literature reveals several critical limitations that motivate the present study. First, SPT–C_c correlations are very scarce and tend to vary with soil type, with some studies noting a decrease in C_c with increasing SPT-N values, and vice versa in other geological settings. Second, past efforts mainly employed simple linear regression, which is an inappropriate tool for modeling linear or asymmetric dependence, which is often the dominant feature of natural soil behaviour, limiting the generalizability of such correlations. Third, to the best of the authors’ knowledge, no prior research has applied copula theory to model the SPT–C_c relationship, whether in Iraq or internationally for regionally representative clay soils, leaving a methodological gap in probabilistic geotechnical modeling. Fourth, in Thi-Qar Governorate, the need for reliable C_c estimation is especially pressing due to the scarcity of undisturbed sampling and the high cost or limited availability of laboratory consolidation tests. Therefore, a probabilistic framework capable of leveraging inexpensive field data such as SPT-N₆₀ to characterize soil compressibility is of practical importance for both routine geotechnical design and reliability-based settlement assessment.

This research aims to develop a copula-based statistical model to describe the relationship between SPT-N60 values and the C_c for clay soils in Thi-Qar Governorate, by integrating real field data with copula dependence structures. The research consists of three parts: the first is collecting the necessary data, then the Copula approach is used to get the dependency, and then simulating the SPT and C_c. The research discussed the marginal distributions of C_c and SPT and the effect of value on the Gaussian copula. Important comparisons with other models are also stated in this research. This research contributes to developing a framework to construct the dependency of C_c and SPT, then includes it in a simulation of the probabilistic state of the two parameters. It is very helpful in examining the reliability analysis and design of shallow foundations and pile foundations.

3.1 Marginal distribution of the N₆₀ and C_c

Initially, it is essential to identify the appropriate marginal distributions. Table 3 includes the candidate marginal distributions. The lowest values of BIC and AIC indicate the best probability distribution that simulates the data [40]. As shown in Table 4, based on the lowest AIC and BIC values, the marginal distributions are estimated to be a Weibull distribution for N₆₀ and a TruncGumbel distribution for C_c. The value of AIC for all the probability distributions is approximately similar to each other. It ranges from 100.9926 to about 103.2645. The narrow range of AIC and BIC may refer to the fact that all the probability distributions fit the data well. However, the lowest value, even with slightly different forms, serves as a basis for comparison between them. The BIC for the second parameter (C_c) was very low values, which refer to a very good fit of the data.

Table 3. Probability density functions and domains of distribution parameters for candidate marginal distributions [14]

Table 4. Values of AIC and BIC for different marginal distributions of data

Figure 3. Histograms and PDFs of the best-fitting marginal distributions of the measured dataset

The PDFs of the four candidate marginal distributions are plotted together with the histograms of the measured data in Figure 3. It may visually evaluate the fitting of the probability distribution to the data. This strengthens the capability of the probability distribution to simulate the data.

3.2 General correlation between SPT and C_c

Before modeling the correlation between SPT-N₆₀ and C_c, the measured data are transformed from the original space to a uniform space to facilitate understanding of the relationship and to provide an idea of the appropriate copula. Figure 4 contains the scatter plot in the uniform space obtained after applying Eq. (24), where the effect of the bounded distribution is isolated. The scatter plot in Figure 4 shows a negative correlation between SPT-N₆₀ and C_c. The negative correlation is also verified by Pearson's rho and Kendall's $\tau$ coefficients, which are calculated using Eqs. (1) and (7).

The results showed a negative correlation between SPT-N₆₀ and C_c Pearson’s rho coefficients, and Kendall’s $\tau$ are obtained as −0.5591 and −0.4280, respectively. These findings are consistent with the results reported by Alam et al. [9]. This study showed that an increase in SPT-N values is accompanied by a decrease in C_c values.

The reason for this difference is that Pearson's rho measures linear dependence influenced by the actual values of the data; thus, it is more sensitive to contrarian behavior and outlier values. By contrast, Kendall's $\tau$ captures dependence based on ranks, providing an indication of the general monotonic trend, as opposed to strict linearity.

Figure 4. Scatter plots of the measured dataset

Due to the observed negative correlation, the Gaussian, Plackett, Frank, and No.16 copulas are suitable for describing the dependence between SPT₆₀ and C_c. Thereafter, the copula parameters θ associated with the four prospective copulas, utilizing Eqs. (4) or (15) for the Gaussian copula and Eq. (19) for the Archimedean Copulas are listed in Table 5. After obtaining the copula parameters, the AIC and BIC values are computed for each copula by using Eqs. (22) and (23). Table 5 summarizes the calculated results, where the Gaussian copula is the most suitable because it has the smallest AIC and BIC values.

Table 5. Correlation coefficients of data, copula parameters, AIC, and BIC values of fitted copulas

Given that the available dataset for the variable pair N₆₀–C_c consists of only 17 samples, relying solely on the original sample to estimate probabilistic model parameters or to select the appropriate copula may be unreliable. Small samples are highly sensitive to data fluctuations, as highlighted by Tang et al. [41]. Their study demonstrated that sample size is the most influential factor in identifying the correct copula or marginal distribution, more so than the COV or the distributional shape. Small datasets tend to produce large variability in AIC values, increasing the likelihood of selecting an incorrect model, particularly when the sample size is less than 30. The study also showed that a stronger true dependence reduces the amount of data required for reliable model identification.

To address uncertainty due to a small sample size, Tang et al. [41] recommended using the bootstrap resampling method. The bootstrap method is based on the idea of random sampling with replacement from the original dataset. For the N new samples we created from sample n. This may cause some of the values to appear in more than one row, but also some not to appear on any row. Then we compute the respective statistics for each of the resampled datasets, e.g., mean, standard deviation, copula dependence parameter, and AIC values. By repeating this process many times, every metric gets a sampling distribution that quantifies the uncertainty associated with it and allows for more robust evaluation of probabilistic models. Bootstrap sampling and coupling fit evaluation procedures are a methodology that can be generalized to various probabilistic relationships between geotechnical variables when data is limited [13].

Based on the recommendations of Sepulveda-Garcia and Alvarez [13], a total of N = 10,000 bootstrap samples were adopted to evaluate the uncertainty in the dependence parameter of each copula. For every resampled dataset, the dependence parameter is computed using the method of moments based on Kendall’s $\tau$ for each copula type. After completing all iterations, a set of 10,000 estimated values of the dependence parameter is obtained, forming a sampling distribution that characterizes the uncertainty associated with each copula.

Table 6 shows bootstrap results from 10,000 resampled datasets, indicating unambiguous best-fit results for the Gaussian copula relative to other models. It had the lowest mean AIC value (−6.07) and was selected as the best-fitting model in about 75% of all iterations. This further implies its objectivity and reliability to map the dependency between N₆₀ and C_c, regardless of the limited amount of data points used in the proposed model (17 samples). On the other hand, the Frank and No.16 copulas provided fair performance (selection rates of 11.5% and 12%, respectively), indicating that they are capturing some of the dependence structure, but with less stability than the Gaussian copula. Among them, the Plackett copula demonstrated the poorest performance since it was only selected in 1.57% of cases, having the highest values of AIC, which mirrors its reduced capacity to fit the data.

Table 6. Bootstrap AIC statistics and selection frequencies for candidate copulas

Figure 5. Bootstrap PDFs of AIC scores for fitted copulas

Figure 5 shows the AIC distributions for four candidate copulas based on the bootstrap analysis. The Gaussian copula has the lowest AIC values, for most of its distribution tilted to the negative zone, thus confirming its clear temporal superiority as a goodness-of-fit over other models. In comparison, the Frank and No.16 copulas have distributions that are more Gaussian-like; however, they lean toward higher AIC values, indicating reasonable but less stable performance. In contrast, the AIC values for the Plackett copula are consistently lower and more tightly concentrated, indicating a more limited capacity of the Plackett copula to accommodate the dependence between N₆₀ and C_c.

3.3 Generation of data based on a copula

Figure 6 compares the scatter plots of 1000 simulated samples from 17 measured samples using a Gaussian copula with the appropriate marginal distributions. It is worth noting that these simulations were performed using the algorithms presented by Li et al. [23].

Figure 6. Scatter diagrams showing both measured and simulated data using a fitted copula with regression analysis

The scatterplot shows a clear negative, monotonic relationship between N₆₀ and C_c. Linear regression on the measured data yields R2 ≈ 0.31, whereas the copula-based simulated samples reproduce a very similar slope and intercept C_c = −0.007 N₆₀ + 0.288 with a slightly higher explanatory power R2 ≈ 0.36. Since R2 only evaluates how well a straight line fits the data, it can under- or over-state dependence when the true relationship is only partially linear—An effect also reflected by the gap between Pearson’s ρ and Kendall’s $\tau$. Most of the data points are clustered in the mid-range of N₆₀ (8–15) and C_c (0.20–0.30); this clustering increases the scatter around the regression line and explains the relatively low coefficient of determination (R2). Nevertheless, the Gaussian copula was able to successfully reproduce this mid-range concentration, confirming its suitability for capturing the dependence structure between the measured data.

3.4 Comparison of copula shapes with the optimal Gaussian copula

The probability density curves of the proposed copulas were plotted and compared to the optimal copula (Gaussian) with the same margins (Weibull distribution for C_c and Truncgumbel for N₆₀) with $\tau$ = −0.4280, as shown in Figure 7. The results indicate that the Plackett and the Frank fit the Gaussian very closely in the middle of the data, while copula No.16 shows a clear departure from the other copulas, with noticeable variation at the upper and lower ends of the relationship, reflecting its nonlinear nature and asymmetric end dependence when the negative correlation is not strong. This suggests that cupula 16 highlights nonlinear dependence effects that the Gaussian cupula cannot accurately represent.

Figure 7. Iso-density contours for copula models (Gaussian, Plackett, Frank, No.16)

3.5 Effect of t value on the Gaussin copula

Figure 8 demonstrates the variation in the dependence structure of the Gaussian copula for different values of Kendall’s $\tau$, ranging from strong positive to strong negative correlation. Each surface represents the bivariate PDF of the copula in the uniform domain (u₁, u₂) given the observed data.

The upper row shows when $\tau$ is positive ($\tau$ = 0.90, 0.428, and 0.016) and the corresponding Gaussian ρ = 0.988, 0.623 and 0.025, respectively. When $\tau$ increases, the joint PDF is more peaked along the main diagonal (u₁ ≈ u₂), with the implication of strong positive dependency: Large (or small) values of the one variable tend to occur at the same time that large (or small) values of the other variable occur. As $\tau$ tends to zero, the surface becomes flatter, which suggests a near independence between the two uniform variables; conversely, in the bottom row, where the correlation is negative.

3.6 Nataf versus Gaussian and No.16 copulas

For the purpose of comparison between the Nataf distribution and both the Gaussian copula and Copula No.16, 1,000 samples were simulated based on the dependence structure inferred from the measured data and the best-fitting marginal distributions, as shown in Figure 9. Modeling results indicate that the Nataf transformation exhibits behavior very similar to that of the Gaussian copula when the dependence between the variables is linear or quasi-linear.

When comparing the Nataf transform with copula No.16, which exhibits asymmetry when the negative correlation is weak, the difference becomes clear, while Nataf maintains a symmetric elliptic dependence, copula No.16 exhibits an asymmetric tail behavior, and highlights the nonlinearity between the two variables. Therefore, the Nataf transform is equivalent to the Gaussian copula when the relationship is approximately linear, but it loses accuracy when modeling nonlinear or tail-dependent variables.

The comparison shows that the Nataf transformation cannot describe other nonlinear or tail-dependent behavior between soil parameters since it inherently imposes a symmetric, Gaussian-type dependence structure. On the other hand, copula models like those based on Gaussian and No.16 have a parametric shape of dependence, but this can be learned from data. For the current N₆₀-C_c dataset, it is nearly monotonic and symmetric, so Nataf and the Gaussian copula behave similarly. However, in datasets where soils exhibit stronger nonlinearities or asymmetric tail dependence, the Nataf model may misrepresent the joint variability, whereas copula-based approaches provide a more reliable characterization of the underlying dependence.

Figure 8. 3D surface plot of bivariate density for the Gaussian copula with different values of $\tau$

Figure 9. Simulation of measured data from Nataf distribution, Gaussian, and No.16 copulas

Table 7. Comparison of simulation outcomes related to different dependence models

Both the Pearson correlation coefficient (ρ) and the Kendall correlation coefficient ($\tau$) were calculated for the simulated correlated samples in the uniform (U) space to assess the dependency structure in the copula space. These coefficients were then recalculated in the physical space of variables (X) defined by marginal distributions. The results are summarized in Table 7, which also includes the correlation coefficients measured from the data (N₆₀, C_c).

It is noted that the Kendall coefficient ($\tau$) values between X₁ and X₂ remained almost constant compared to those calculated between U_I and U₂ for all the reliability models studied. This is attributed to the stability of the Kendall coefficient against the monotonic transformations used to move from the uniform distribution space U to the physical space X. In contrast, the ρ shows significant variation between the two spaces due to the nonlinear nature of the inverse cumulative distribution transformations used to generate the values of X. In comparison, the Nataf model produced correlation values close to the Gaussian correlation, demonstrating similar linear dependence.

In general, these results validate that $\tau$ is invariant under monotonic transformation, but ρ is influenced by the shape of the marginal distributions. The Nataf model also exhibits similar behavior to the Gaussian copula in representing linear dependence, but differs from copula (16) when describing nonlinear relationships, pointing to the significance of model choice in this case concerning some types of geotechnical data in the non-Gaussian case.