Optimal siting and sizing of renewable sources in distribution system planning based on life cycle cost and considering uncertainties

1.   Introduction

Numerical analysis has demonstrated considerable advantages over traditional methods in providing detailed insights and improved understanding of geotechnical stability, limit-state design, and soil-structure interaction problems. A critical factor for its success lies in accurately determining constitutive model parameters, which often relies on limited soil data and engineering judgment. Prediction competitions highlight the variability in results when the same problem is approached by different engineers. This raises the question of whether automating parameter determination could minimize this variability and improve the reliability of numerical analysis in geotechnical design ^[1]. Over time, soil constitutive models have advanced significantly. While more sophisticated models can represent soil behavior with greater accuracy than simpler ones, they also require a larger number of parameters. These parameters are typically derived from laboratory tests, such as triaxial and oedometer tests, which may not always be accessible for every project.

In situ tests offer an alternative approach for determining soil parameters. Compared to laboratory tests, they are generally faster, more cost-effective, and tend to cause less disturbance to the soil. However, soil parameters cannot be directly obtained from in situ test results. Instead, various empirical correlations have been developed to relate field measurements to soil parameters. Since multiple correlations often exist for the same parameter, the resulting values can vary widely, introducing uncertainty. This uncertainty arises from the fact that these correlations are typically developed for specific soil types or conditions, limiting their general applicability.

An ongoing research project is focused on developing an automated parameter determination (APD) framework to identify constitutive model parameters using in situ test data. This is particularly valuable in the early stages of a project when soil data is limited. At this stage, (relatively inexpensive) field tests are carried out prior to a full laboratory test campaign. The goal is not to replace laboratory tests with in situ methods but to complement them, as laboratory tests remain essential for refining soil and constitutive model parameters in the final design.

In APD, parameters are determined using a graph-based approach, which leverages certain principles of graph theory ^[2]. The primary goal of the project is to develop a parameter determination system that is both transparent and adaptable. Transparency is achieved by clearly demonstrating how parameter values are derived from the available data, while adaptability is ensured by enabling users to integrate their own experience and expertise into the system.

Section 2 presents the APD framework, followed by a case study in Section 3, where soil parameters for a sand site are assessed based on CPT data. In Section 4, these parameters are used to determine material sets for the Clay and Sand Model (CASM). The resulting constitutive model parameters are then evaluated through numerical simulation of a CPT, and the output is compared with in situ CPT data from the sand site. Finally, Section 5 discusses the conclusions and future research directions.

2.   Automated parameter determination framework

The parameter determination framework is developed in the programming language Python and is structured as a series of sequential modules that connect raw field measurements to constitutive model parameters. In the case of computing parameters from CPT, the framework comprises five modules, as illustrated in Figure 1. The first module handles the import of field measurements, which are then passed to the second module for CPT stratification. The stratified layers are sent to the third module, where the state of each layer is determined, including the overconsolidation ratio ($OCR$) and the coefficient of earth pressure at rest ($K_0$). In the fourth module, soil parameters are computed for each layer, while the fifth module evaluates the corresponding constitutive model parameters. The framework and the five modules are described in detail in ^[3].

Figure 1.  Overview of modules, used for automated parameter determination ^[3].

DownLoad: Full-Size Img PowerPoint

2.1.   Graphs as a parameter determination tool

Graph theory, a branch of discrete mathematics, explores relationships between objects within a network. A graph, in this context, serves as a mathematical representation of a network and consists of two key components: nodes, which represent entities, and edges, which define the connections between them. The APD framework utilizes a weighted directed graph, where connections between nodes have a defined direction, and edges can carry weights. A common example of such a graph is a one-way road network, where weights assigned to edges represent attributes like distance or travel time ^[2].

Figure 2 illustrates the fundamental concept of the graph-based approach used in APD, where source parameters (CPT raw data) are linked to destination parameters (soil or constitutive model parameters) through intermediate parameters. This connection is established using a predefined set of correlations. The system identifies all possible paths (sequences of correlations) that connect source parameters to destination parameters and calculates the values of both intermediate and destination parameters. Within the APD framework, the term "method" is used as a general descriptor, encompassing correlations, formulas, equations, and rules of thumb. The graphs are generated using the Python graph visualization library graphviz ^[4].

Figure 2.  Graph-based approach implemented in APD ^[5].

DownLoad: Full-Size Img PowerPoint

2.2.   Generating the graph

The graphs are constructed based on the parameters, methods, and their interconnections. This is accomplished by supplying the system with input files in comma-separated value (CSV) format. Two files are used: one containing the methods and the other containing the parameters. Each file includes specific properties that must be defined to establish links between methods and parameters, enabling the generation of the graph.

To demonstrate these properties, the correlation (Eq. 3.12) by ^[6], which estimates the relative density ($D_r$) from CPT data, is used as an example for creating CSV files for methods and parameters. The corresponding CSV files for this method, along with the generated graph, are shown in Figure 3.

Figure 3.  Format of methods and parameters CSV files ^[3].

DownLoad: Full-Size Img PowerPoint

Starting with the methods CSV file, several key properties must be defined: method_to, Formula, parameters_in, parameters_out, validity and Reference. The method_to field specifies the name of the method, which in this case is method_to_Dr. The Formula field contains the equation used in the method (Eq. 3.12); however, since "$\ln$" is not recognized in Python, "$\log$" (logarithm with base e) is used instead, as shown in Figure 3 (note that $\times100$ is removed for illustration purposes). The input parameters for the method ($q_c$ and $\sigma^{'}_{v}$) are listed under parameters_in, while the output parameter ($D_r$) is specified in parameters_out. The validity field defines the applicability of the method, as some methods are suitable for all soil types, while others are restricted to either coarse- or fine-grained soils. The validity classification follows Robertson's modified non-normalized soil behavior type (SBT) chart ^[7] and is outlined in Table 1. For example, if a method is applicable only to silts, the validity should be set to SBT(4). In the case of the example shown in Figure 3, since relative density applies to coarse-grained soils, the validity is assigned as SBT(5678) (refer to Table 1). Beyond SBT classification, additional criteria such as the normalized excess pore pressure ($B_q$) can also define validity. If a method is only applicable for $B_q$ values within the range [0.1, 1.0], this can be specified as Bq_min(0.1) Bq_max(1.0). If $B_q$ falls outside this range, the method will be excluded from consideration. Additionally, the method's author can optionally be cited in the Reference field.

Table 1.  SBT zones according to ^[7].

Zone	Soil behavior type (SBT)
1	Sensitive fine-grained
2	Clays: organic soil
3	Clays: clay to silty clay
4	Silt mixtures: clayey silt and silty clay
5	Sand mixtures: silty sand to sandy silt
6	Sands: clean sands to silty sands
7	Dense sand to gravelly sand
8	Stiff sand to clayey sand (overconsolidated)
9	Stiff fine-grained (overconsolidated)

---

 | Show Table

DownLoad: CSV

For the parameters CSV file, the following properties must be defined: Symbol, Description, Unit, Value, and Constraints. This file includes all parameters used in the methods CSV file, specifically those listed in the parameters_in and parameters_out fields. The Symbol field follows the same notation as in the methods CSV file to ensure consistency. The remaining fields are optional but provide additional clarity. The Description field allows users to specify the parameter's name, while the Unit field is recommended to prevent errors due to unit conversions. The value field is used to assign a specific value to a parameter when manually defined by the user (e.g., the unit weight of water). If the value is computed by the system, this field should remain empty. The Constraints field defines upper and lower limits for a parameter. If the system calculates a value exceeding these limits, it is discarded. For instance, in Figure 3, the relative density is constrained within a lower bound of 0% and an upper bound of 100%.

There are two types of nodes in the generated graph: parameter nodes and method nodes. The parameter nodes are defined in the parameters CSV file, while the method nodes are specified in the methods CSV file. The connections between these nodes—i.e., the edges—are established based on the parameters_in and parameters_out fields in the methods CSV file. Specifically, parameters_in defines the incoming edges from the required parameters to the corresponding method node, while parameters_out defines the outgoing edge from the method node to the resulting parameter node.

2.3.   Correlations database

The correlations database serves as an external input to the system, with users responsible for selecting the methods used to compute parameters. However, APD includes a standardized, validated database of methods and parameters, which is continuously updated and refined ^[3,8,9,10]. Even when utilizing the provided database, users should apply their expertise and judgment when interpreting the system's output. Ultimately, the system aims to generate "reasonable" parameter values while allowing for user oversight and refinement.

Currently, APD comprises three main workflows, allowing parameter assessment based on CPT, dilatometer test (DMT), and shear wave velocity measurements. As a result, three validated databases of methods and parameters are provided, corresponding to these workflows. These databases collectively include over 200 methods. It is important to note that both soil parameters and constitutive model parameters are included within the same database.

3.   Application of APD on the Øysand site

This section describes the use of the APD system to calculate soil parameters for a test site, employing the CPT workflow. The test site is described in Section 3.1, while some soil parameters are presented in Section 3.2. These soil parameters are then compared to "reference" values obtained from laboratory tests, as illustrated in Section 3.3. Additional soil parameters are discussed in Section 3.4. These additional soil parameters are not compared to reference values; however, they are needed for the numerical simulation of the CPT as described in Section 4.

3.1.   Test site

Data for this analysis was gathered from one of the Norwegian GeoTest Sites (NGTS) ^[11] using the web-based application Datamap.

3.1.1.   Datamap

The web application Datamap has been developed to collect, classify, and manage geotechnical data. It serves as a platform for making geotechnical information accessible while enabling researchers to create and share their projects. The application can be accessed at www.geocalcs.com/datamap ^[12].

3.1.2.   Norwegian GeoTest Sites (NGTS)

Between 2016 and 2019, the Norwegian Geotechnical Institute (NGI), the Norwegian University of Science and Technology (NTNU), SINTEF Building and Infrastructure, the University Centre in Svalbard (UNIS), and the Norwegian Public Roads Administration (NPRA) collaborated to establish five GeoTest Sites (NGTS) across Norway. These sites represent various soil types, including clay, silt, quick clay, sand, and permafrost ^[11].

3.1.3.   Øysand test site

The Øysand test site is located approximately 15 km south west of Trondheim, Norway. Site characterization at the Øysand research facility started in 2016 as part of the NGTS project. Since 2016, a wide range of in situ tests, geophysical techniques, sampling techniques, and laboratory tests have been used to assess the geological history and geotechnical properties of the sand deposits ^[13].

Figure 4 shows the test layout as well as the cone tip resistance ($q_c$) for two CPTs, namely OYSC09 and OYSC20. OYSC09 is located in the area of the initial sampling campaigns where several other CPTs were conducted. The stratification in this area is divided into two main units. The first unit (upper 6 to 10 m) consists of coarse to gravelly sand. The lower unit consists mainly of fine silty sand with some layers of clay and silt ^[13]. The two units are highlighted in the OYSC09 results as shown in Figure 4b. The soils at the Øysand site are normally to slightly over-consolidated, and the laboratory triaxial compression tests show a strong dilative response ^[13].

Figure 4.  (a) Layout of the test site, (b) results of OYSC09 and (c) results of OYSC20.

DownLoad: Full-Size Img PowerPoint

Sampling the gravelly sandy soil using piston technology was extremely difficult at the Øysand site ^[14]. Piston sampling techniques disturb saturated, uncemented coarse-grained soils. Therefore, ground freezing for sampling was considered as an alternative option to obtain high-quality sand samples ^[14]. As the ground freezing sampling campaign required access to a road, a different site location was selected (compared to the initial sampling campaign). The location of the ground freezing is shown in Figure 4a. The conditions under which the ground freezing occurred were different from those encountered during the initial sampling campaign. Gravel was present in the upper 7 m, resulting in a need for pre-drilling ^[14]. Therefore, the results of OYSC20 shown in Figure 4c are only available from 7 m to 20 m.

The stratigraphy in the ground freezing zone consists of an upper layer extending from 0 to 6 m, composed of agricultural silty sand and poorly graded sand and gravel. Beneath this, a thick layer of silty sand and sandy silt is present from 6 to 10.5 m. The main sand layer is found at depths of 12 to 18.2 m ^[14]. This sand layer was the focus of the ground freezing campaign and is the primary subject of this study. This layer is highlighted in the OYSC20 results as shown in Figure 4c (sand layer). The groundwater level is 2 m below the ground surface ^[13].

3.2.   Soil parameters

In this section, the methods selected for determining the unit weight ($\gamma_t$), void ratio ($e_0$), relative density ($D_r$), peak friction angle ($\varphi_p^\prime$), and shear wave velocity ($V_s$) by CPT are presented.

As previously outlined ^[3], it is necessary to initially estimate the unit weight in module 1. A selection of methods for determining the unit weight are outlined in Table 2. The initial unit weight can be defined using any method or value. In this study, Eq. 3.4 (in Table 2) was used to calculate the initial unit weight.

Table 2.  Methods for the unit weight ($\gamma_t$).

Method		Author
$\gamma_w [0.27(\log{R_f}) + 0.36(\log{q_t/p_a}) + 1.236]$	(3.1)	[15]
$19.5 - 2.87 [\frac{\log(\frac{9000}{q_t})}{\log(\frac{20}{R_f})}]$	(3.2)	[16]
$26 - \frac{14}{1+[0.5\log f_s+1]^2}$	(3.3)	[17]
$0.254 \cdot \log (\frac{q_t-u_2}{p_a})+1.54$	(3.4)	[18]

---

 | Show Table

DownLoad: CSV

When the soil is fully saturated, the saturated unit weight of soil can be determined as follows:

Equation 3.5 can be rewritten to determine the void ratio from the saturated unit weight:

In this study, the specific gravity ($G_s$) is set to 2.65, and the unit weight of water ($\gamma_w$) is set to 10 kN/m³. The saturated unit weight ($\gamma_{sat}$) can be determined using the methods presented in Table 2. Equation 3.4 was used to compute $e_0$.

The methods selected for determining $D_r$, $\varphi^{\prime}_p$, and $V_s$ are highlighted in Tables 3, 4, and 5, respectively.

Table 3.  Methods for the relative density ($D_r$).

Method		Author
$100\times\sqrt{\frac{1}{305 OCR^{0.2}} Q_{tn}}$	(3.7)$^{a}$	[19]
$100\times(0.268\times\ln{\frac{q_t/p_a}{(\sigma_v^\prime/p_a)^{0.5}}}-0.675)$	(3.8)	[20]$^{b}$
[0.3cm] $\frac{0.87(q_t/p_a)}{(\sigma_v^\prime/p_a)^{0.5}}$	(3.9)	[17]
$-98+66\times\log{\frac{q_c}{\sigma_v^\prime}}$	(3.10)	[21]$^{c}$
$100\times\frac{1}{2.41}\ln{\frac{Q_{tn}}{15.7}}$	(3.11)	[22]$^{d}$
$100\times\frac{1}{2.91}\ln{\frac{q_c}{61\sigma_v^{0.71}}}$	(3.12)	[6]
$100\times\frac{1}{3.10}\ln{\frac{q_c/p_a}{17.68(\sigma_v^\prime/p_a)^{0.5}}}$	(3.13)	[20]
$^a$ $Q_{tn}$ is the normalized cone resistance corrected for stress level and the normalized SBT index ($I_{cn}$). OCR is discussed in section 3.4. $^b$ as cited in [18]. $^c$ as cited in [23]. $^d$ as cited in [24].

---

 | Show Table

DownLoad: CSV

Table 4.  Methods for the peak friction angle ($\varphi^{\prime}_p$).

Method		Author
$17.6+11\log(Q_{tn})$	(3.14)	[19]
$\tan^{-1}[\frac{1}{2.68}\log\frac{q_c}{\sigma_v^\prime}+0.29]$	(3.15)	[25]
$25(Q_{tn})^{0.10}$	(3.16)	[26]$^{a}$
$53-6.9I_{cn}$	(3.17)	[18]
$^a$ as cited in [27].

---

 | Show Table

DownLoad: CSV

Table 5.  Methods for the shear wave velocity ($V_s$).

Method		Author
$51.6 \ln f_s +18.5$	(3.18)	[28]
$277(\frac{q_c}{1000})^{0.13}(\frac{\sigma_v^\prime}{1000})^{0.27}$	(3.19)	[29]
$13.18q_c^{0.192}(\sigma_v^\prime)^{0.179}$	(3.20)	[30]
$12.02q_c^{0.319}f_s^{-0.0466}$	(3.21)	[30]
$(10.1\log q_c-11.4)^{1.67}(f_s/q_c \times100)^{0.3}$	(3.22)	[30]
$(10^{(0.55I_{cn}+1.68)}\times\frac{q_t-\sigma_v}{p_a})^{0.5}$	(3.23)	[24]
$2.27q_t^{0.142}I_{cn}^{0.989}z^{0.033}$	(3.24)	[31]
$10^{(0.8I_{cn}-1.17)}Q_{tn}/(\frac{p_a}{\sigma_v^\prime})^{0.25}$	(3.25)	[32]
$0.0831q_{c1N}(\frac{\sigma_v^\prime}{p_a})^{0.25}e^{1.786I_{cn}}$	(3.26)	[33]
$18.4q_c^{0.144}f_s^{0.0832}z^{0.278}$	(3.27)	[34]
$135(q_c/1000)^{0.23}(\frac{\sigma_v^\prime}{p_a})^{0.135}$	(3.28)	[35]
$^a$ $z$ is the depth below the ground surface in m. $^b$ $q_{c1N}$ is the normalized cone tip resistance (for $I_{cn} < 2.6, q_{c1N}=(\frac{q_c}{p_a})(\frac{p_a}{\sigma_v^\prime})^{0.5}$).

---

 | Show Table

DownLoad: CSV

3.3.   Results

Soil parameters are determined based on OYSC20 (see Figure 4c). The detailed results of OYSC20 in terms of $q_c$, sleeve friction $f_s$, and pore pressure measurements $u_2$ are shown in Figure 5(a-c). The different soil units are also highlighted in the figure. The focus is on the third unit, namely the sand layer.

Figure 5.  Results of OYSC20, (a-d): profiles of $q_c$, $f_s$, $u_2$, and $I_{cn}$.

DownLoad: Full-Size Img PowerPoint

Section 3.2 presented the methods for determining soil parameters for the CPT test illustrated in Figure 5. In this study, the CPT results were averaged at 4 cm intervals (manual layering, see also ^[3,8,9,10]), with these values then being utilized as the input for the analysis. The rationale behind this methodology (i.e., manual layering and averaging) is to generate a sufficient number of layers that encompasses the CPT sounding, facilitating comparison with the reference values of various parameters. The APD predominantly calculates parameters based on layers; consequently, selecting layers with reduced thickness would result in a greater number of layers and, consequently, a larger number of graphs. Other intervals, rather than 4 cm, could be used; the choice of 4 cm in this study was made to demonstrate the system's ability to handle very fine stratification.

Given that the layers are determined manually, the user is required to provide the SBT for each layer. This SBT will act as a validity criterion for the methods CSV file. In this study, as the focus is on the sand layer, all the layers had SBT(6) corresponding to sand, according to ^[7] (refer to Table 1). The averaging process resulted in 337 layers for the entire CPT sounding, out of which 148 are within the sand layer of interest. The APD system finally generated 337 graphs corresponding to the 337 layers.

The comparison between the output of APD (using the methods presented in Section 3.2) and the reference values interpreted at the Øysand test site is presented in Figure 6. The results are presented for the sand layer (between 12 and 18 m). The blue shaded areas represent the range of values obtained from the CPT workflow, while the blue lines with circle markers correspond to the average value of the methods. Thus, the circle markers indicate the mid-depth of the 148 thin (4 cm thick) layers.

Figure 6.  Comparison between APD and interpreted values at the Øysand test site (the blue shaded areas represent the range of values obtained from the CPT workflow): (a) unit weight, (b) void ratio, (c) relative density, (d) peak friction angle, and (e) shear wave velocity.

DownLoad: Full-Size Img PowerPoint

The total unit weight was determined indirectly from water content measurements of the samples assuming full saturation ^[13]. However, the water content measurements obtained from disturbed samples were found to be lower than expected, presumably due to the disturbance of the samples ^[14]. Water migration can be expected to occur during sampling, transportation, and storage ^[14]. Consequently, the in situ unit weight is expected to be lower than the reported ones. The performance of the four methods illustrated in Table 2 is shown in Figure 6a. Equation 3.4, which acts as the initial unit weight, represents the upper bound of the blue area. The overall tendency of the methods underestimates the reference values; however, as previously stated, these reference values are estimated based on disturbed samples, and it is expected that the in situ values will be lower. The influence of the unit weight on the results obtained was also investigated for another test site in ^[3]. The initial unit weight influences several components of the APD from the outset of the analysis. Therefore, it is advisable to establish a representative unit weight using cost-effective tests that can be performed during the early phases of a project. In cases where such tests are not available, the selected unit weight should be cross-checked against values typically observed at the specific site, based on prior experience. To address this, a feature is currently being developed that allows users to rerun the analysis with an updated initial unit weight determined in module 4 of the various APD workflows.

Five isotropically consolidated triaxial tests—three undrained (IUC) and two drained (IDC)—were conducted on undisturbed (frozen) samples at depths of 13.3, 13.8, 14.1, 14.3, and 14.6 m ^[14]. Consequently, five values are reported as reference values for $e_0$, $D_r$, and $\varphi_p^\prime$ as highlighted in Figure 6(b-d). The determination of the void ratios of frozen undisturbed samples is discussed in ^[14]. Figure 6b shows the values obtained through Eq. 3.6 using the four unit weight methods outlined in Table 2. The void ratio derived from using the initial unit weight (Eq. 3.4) is employed in the subsequent numerical analysis (refer to Section 4). This unit weight represents the lower bound of the range shown in Figure 6b.

The reference relative density values are found to be significantly influenced by the method employed to determine $e_{min}$ and $e_{max}$ in the laboratory ^[14]. The performance of the seven methods presented in Table 3 is shown in Figure 6c. It is evident from the figure that the upper bound of the range successfully captures the first three reference values, while the lower bound of the range captures the lower two reference values. The average of all methods underestimates the first three reference values, while it overestimates the lower two reference values.

The peak friction angles were determined from three undrained isotropically consolidated (IUC) and two drained isotropically consolidated (IDC) triaxial tests as highlighted in Figure 6d. The performance of the four methods presented in Table 4 is also highlighted in the figure. With the exception of the fourth reference value (at depth of 14.3 m), the overall tendency of the methods is to predict the peak friction angle in a reasonable manner.

The reference values of the shear wave velocity were determined from the seismic SCPT20 (OYSC20 in the project). The eleven methods presented in Table 5 predict the in situ $V_s$ reasonably well.

3.4.   Other parameters

This section presents the derivation of other soil parameters that are required for the numerical simulation of the CPT presented in Section 4. A comprehensive summary of these parameters is provided in Table 6.

Table 6.  Soil parameters used in the numerical simulation of the CPT.

Parameter	Symbol	Method		Author
Slope of the compression line in the ($\ln p', v$) plane	$\lambda$	$\frac{F_r}{10\times2.303}$	(3.29)	[36]$^{a}$
State parameter	$\psi$	$0.56-0.33\log Q_{tncs}$	(3.30)	[37]
Permeability	$k$	$10^{(0.952-3.04I_{cn})}$	(3.31)	[38]
Critical state friction angle	$\varphi_{cv}^\prime$	$\varphi_p^\prime + 48\psi$	(3.32)	[39]
Slope of the critical state line in $p^{\prime}$ - q space	$M$	$\frac{6\sin\varphi_{cv}^\prime}{3-\sin\varphi_{cv}^\prime}$	(3.33)
Overconsolidation ratio	$OCR$	$\frac{0.33(q_t-\sigma_v)^{m^\prime}}{\sigma^{'}_{v}}$	(3.34)	[40]
Coefficient of earth pressure	$K_0$	$(1-\sin \varphi_{cv}^\prime) OCR^{\sin \varphi_{cv}^\prime}$	(3.35)	[41]$^{f}$
$^a$ as cited in [39], $F_r$ is the normalized friction ratio ($F_r=f_s/(q_t-\sigma_v)$). $^b$ $Q_{tncs}$ is the clean sand equivalent normalized cone resistance as described in [37]. $^c$ for $1.060 I_{cn}\leq3.27$. $^d$ by rearranging ($\varphi_p^\prime= \varphi_{cv}^\prime-48\psi$). $^e$ $m^\prime$ is the yield stress exponent that increases with fine content and decreases with mean grain size ($m^\prime = 1 - \frac{0.28}{1+(\frac{I_{cn}}{2.65})^{25}}$) [42]. $^f$ as cited in [19].

---

 | Show Table

DownLoad: CSV

The results for the parameters presented in Table 6 are shown in Figure 7. $\lambda$ is determined using Eq. 3.29, and the output of APD is highlighted in Figure 7a. The average value of $\lambda$ for the sand layer is 0.01967. The negative values of the state parameter presented in Figure 7b indicate a dense state for the sand layer (average value of $\psi$ for the sand layer is –0.07251). The permeability ($k$) of the sand layer is demonstrated in Figure 7c, with an average value of $4\times10^{-5}$ m/s. The critical state friction angle ($\varphi_{cv}^\prime$) is computed using Eq. 3.32. The peak friction angle used in Eq. 3.32 is the average of the four methods shown in Table 4 (blue line in Figure 6d). The average critical state friction angle for the sand layer is $35.18$°. $M$ (illustrated in Figure 7e) is derived from the critical state friction angle presented in Figure 7d. The mean value of $M$ is 1.426. The $OCR$ is computed through Eq. 3.34, and the average value of the sand layer is 1.6. $K_0$ is determined from Eq. 3.35, and the mean value of the sand layer is 0.555. The reliability of the approach for different soil types is ensured through the use of a correlations database that includes methods specifically developed for both granular and fine-grained soils. In this study, only correlations applicable to granular soils were used, as the site consists predominantly of sand. More generally, the APD system automatically selects suitable methods based on the soil behavior type (SBT) identified at each depth, ensuring the appropriateness of the CASM parameters for varying soil types.

Figure 7.  Results of other parameters: (a) $\lambda$, (b) $\psi$, (c) $k$, (d) $\varphi_{cv}^\prime$, (e) $M$, (f) $OCR$, (g) $K_0$.

DownLoad: Full-Size Img PowerPoint

4.   Numerical simulation of CPT

In this study, parameters derived from the APD framework based on in situ CPT data will be used to numerically simulate CPT. The goal is to assess the parameters presented in Figures 6 and 7 by comparing the simulation results with in situ measurements. The simulation will be conducted using the CASM constitutive model, which is discussed in Section 4.1. The CPT will be simulated using the geotechnical particle finite element method (G-PFEM), with the methodology detailed in 4.2.

4.1.   Clay and Sand Model (CASM)

CASM is an isotropic elastoplastic constitutive model developed by ^[43], featuring a single yield function and plastic potential to describe the behavior of both fine-grained and coarse-grained soils. It was formulated to unify critical state theory by incorporating the state parameter concept. CASM was later extended to finite deformations by ^[44], integrating the hyperelastic framework proposed by ^[45]. A broad range of geomaterial behaviors can be captured by adjusting the stress-state coefficient or shape factor ($n$) and the spacing ratio ($r$) in the yield function:

where $q$ represents the deviatoric stress, $p'$ is the mean effective stress, $M_\Theta$ is the slope of the critical state line, which varies with the Lode angle ($\Theta$), and $p'_0$ denotes the preconsolidation pressure.

Experimental observation shows that $n$ typically ranges from 1 to 5, while $r$ for sands varies between 1 and infinity $(\infty)$. The spacing ratio in sands is generally much higher compared to clays, which typically ranges from 1.5 to 3.0 ^[43].

CASM employs a non-associated flow rule, where the plastic potential was originally derived from Rowe's stress dilatancy relation ^[46]. However, an alternative plastic potential function (Eq. 4.2) proposed by ^[47] is adopted in this study to address the limitations of the original equation for stress paths with low stress ratios ^[43]. The parameter controlling the shape of the plastic potential function ($m$) is estimated from $K_0$, following the strategy proposed by ^[47].

4.1.1.   Calibration of r and n

The calibration of $r$ and $n$ was carried out by simulating an IDC test and comparing it to one of the reported tests in ^[14]. The CASM material set used for the calibration process is presented in Table 7. These values represent the mean values described in Section 3.4. The value of $e_0$ corresponds to the IDC test at depth of 14.6 m, as reported in ^[14]. The swelling parameter ($\kappa$) is set to 0.005, a typical value for sands ^[43]. Alternatively, it could be evaluated from the shear wave velocity measurements through the small-strain shear modulus ($G_0$) and the bulk modulus ($K$). Poisson's ratio $\nu'$ is set to 0.3. The simulation of the IDC test was conducted in the finite element code PLAXIS ^[48] using CASM as a user-defined model (USDM) ^[49]. The calibration results are presented in Figure 8, employing n and r values of 5 and 1000, respectively.

Table 7.  CASM material set used for calibration of r and n.

$\lambda$	$\kappa$	$\nu'$	$\varphi_{cv}^\prime$	$e_0$	$\psi$	$m$
0.01967	0.005	0.3	35.18	0.55	-0.07251	2
-	-	-	°	-	-	-

---

 | Show Table

DownLoad: CSV

Figure 8.  Calibration of n and r (a) $\varepsilon_1-q$, (b) $\varepsilon_1-\varepsilon_v$, isotropic compression $p' = 150$ kPa.

DownLoad: Full-Size Img PowerPoint

As previously discussed, CASM is an elastic-plastic hardening/softening model that defines a yield surface expanding in response to $p_0'$. At any given moment, it inherently divides the stress space into a high-stiffness elastic region and an elastic-plastic region where stiffness experiences a sudden reduction. Consequently, this class of models inevitably introduces abrupt stiffness transitions in dense specimens, as observed in the stress-strain curve in Figure 8a. Such behavior deviates from the typical stress-strain response of soils. However, this does not imply that elastic-plastic models are unsuitable; rather, it highlights that certain aspects of soil behavior may not be fully captured by these models ^[50].

4.2.   Geotechnical particle finite element method (G-PFEM)

The CPT was simulated using the geotechnical particle finite element method (G-PFEM). G-PFEM is an implementation of the particle finite element method (PFEM) designed to solve large deformation penetration problems in geotechnics, initially proposed by ^[51]. The framework was developed by the International Center for Numerical Methods in Engineering (CIMNE) and the Polytechnic University of Catalonia (UPC). G-PFEM has been integrated into the Kratos Multi-Physics framework ^[52], an object-oriented platform for multi-disciplinary numerical analysis.

G-PFEM utilizes an updated Lagrangian formulation, which is particularly effective for the numerical analysis of coupled problems involving large deformations. To address numerical challenges associated with large deformation problems, continuous remeshing is employed. While this approach increases computational costs, it is mitigated by the use of low-order elements in combination with a mixed, stabilized formulation. This strategy not only reduces computational time but also alleviates the locking phenomenon typically encountered with low equal-order approximations for displacement and water pressure fields. In ^[53], various mixed formulations were implemented and assessed, introducing an additional degree of freedom. Specifically, the determinant of the deformation gradient ($J$) or the mean pressure ($p$) is incorporated alongside the displacement ($u$) and water pressure fields ($p_w$). It was found that the mixed ($u$-$J$-$p_w$) formulation performed well, and as a result, this approach is adopted in the present study to solve the coupled hydro-mechanical problem in cone penetration scenarios.

This tool has been previously utilized for CPT simulations across various soil types, including clay ^[54], organic clay ^[55], intermediate soil ^[56], and liquefiable soil (loose sand, mine tailings, hydraulic fills, and quick clay) ^[57,58,59]. In this study, it will be applied to simulate penetration in medium-dense sand exhibiting dilative behavior.

4.3.   Numerical model

A fully coupled hydro-mechanical model is employed to simulate cone penetration under axisymmetric conditions. The model boundaries are defined 300 cm away from the CPTu tip to minimize potential boundary effects. The bottom boundary is constrained in both lateral and vertical directions, while the side boundary is laterally fixed. For hydraulic conditions, free drainage is permitted along the upper and lower boundaries, whereas the lateral boundary is impervious.

The cone is initially embedded 20 cm into the ground and then pushed at a constant velocity of 2 cm/s. A time step of 0.0022 seconds is used, with re-meshing performed at each step. In this study, the minimum mesh size is set to 0.4 cm in the interior and 0.8 cm along the boundaries. The robustness of the adopted numerical framework, particularly regarding convergence behavior and mesh sensitivity, has been thoroughly demonstrated by ^[60], while the stabilization, re-meshing, and mapping procedures are described in ^[61].

In this simulation, the soil is assumed to be weightless, resulting in uniform in situ effective stresses ($\sigma_v^\prime$, $\sigma_h^\prime$) and pore water pressure ($u_0$). To maintain equilibrium, a total stress boundary condition is applied at the top surface of the model (see Figure 9a). Under the assumption of a homogeneous stress state, it is expected that the penetration responses ($q_t$, $f_s$, and $u_2$) will reach a stationary equilibrium value (convergent) after sufficient penetration depth.

Figure 9.  Axisymmetric model for cone penetration test (a) geometry and boundary condition, (b) initial mesh before penetration, and (c) updated mesh at the end of penetration.

DownLoad: Full-Size Img PowerPoint

4.4.   Simulation input

When conducting CPT in sandy soils, the measured values of $q_t$ and $f_s$ tend to exhibit greater variability compared to penetration tests performed in fine-grained soil. This results in a range of $q_t$ and $f_s$ values within the same sand layer (see Figure 5), and is also reflected in the fluctuations of soil parameters derived from APD correlations, shown in Figures 6 and 7.

In the current state of practice, numerical CPT simulations cannot account for varying soil layers within a single model. To address this limitation, numerical simulations will be performed at multiple depths, incorporating distinct input parameters determined by the APD framework for each respective depth. A summary of the input parameters derived from the APD framework for each simulated depth is presented in Table 8.

Table 8.  Input parameters for different depths for the numerical simulations.

Depth (m)	$\lambda$ (-)	$\varphi_{cv}^\prime$ (°)	$k$ (m/s)	$e_0$ (-)	$\psi$ (-)	$\sigma_v^\prime$ (kPa)	$\sigma_h^\prime$ (kPa)
12.36	0.02605	34.75	1.40E-05	0.61	-0.07895	121	71
12.52	0.01120	36.28	6.01E-05	0.60	-0.06211	122	72
12.76	0.00796	36.36	7.73E-05	0.60	-0.05908	125	72
13.36	0.01609	35.70	6.07E-05	0.56	-0.08581	131	82
13.52	0.01467	35.51	3.98E-05	0.59	-0.07312	133	77
13.76	0.01525	35.39	8.14E-05	0.55	-0.09949	135	87
14.04	0.02808	34.77	9.09E-06	0.62	-0.06645	138	75
14.36	0.01068	35.44	1.31E-04	0.54	-0.10056	142	90
14.52	0.02793	34.21	4.11E-05	0.52	-0.12845	143	98
14.76	0.05238	33.57	1.98E-06	0.64	-0.07015	146	76
15.04	0.02513	35.01	8.17E-06	0.62	-0.05220	149	75
15.36	0.01927	35.31	1.13E-05	0.63	-0.04592	152	75
15.52	0.00927	36.02	2.28E-05	0.64	-0.03077	153	72
15.76	0.01368	35.62	4.61E-05	0.57	-0.07131	156	86
16.04	0.03434	34.54	2.44E-06	0.66	-0.03543	159	73
16.36	0.02718	34.35	2.29E-05	0.55	-0.10173	162	96
16.52	0.02030	35.00	1.70E-05	0.59	-0.06677	164	85
16.76	0.01846	35.14	1.68E-05	0.60	-0.05993	166	84
17.04	0.01592	35.29	4.24E-05	0.55	-0.08040	169	94
17.36	0.02837	34.74	7.77E-06	0.60	-0.06059	173	86
17.52	0.01361	35.64	5.07E-05	0.55	-0.07264	174	94
17.76	0.01903	34.99	2.14E-05	0.57	-0.07123	177	92
12.68	0.01082	36.03	4.48E-05	0.61	-0.05705	124	70
13.08	0.01178	36.17	5.89E-05	0.59	-0.06451	128	75
13.6	0.01477	35.51	2.22E-05	0.62	-0.05488	134	71
14.28	0.01205	35.60	9.36E-05	0.55	-0.09061	141	87
14.64	0.03793	33.83	1.39E-05	0.55	-0.11503	145	92
14.88	0.02949	34.58	1.02E-05	0.60	-0.07609	147	81
15.24	0.01536	35.62	1.27E-05	0.64	-0.03493	151	72
15.48	0.01246	35.86	1.40E-05	0.65	-0.02665	153	70
15.88	0.02536	34.73	1.36E-05	0.59	-0.07419	157	85
16.16	0.02773	34.58	1.22E-05	0.58	-0.07838	160	87
16.68	0.01667	35.21	2.30E-05	0.59	-0.06445	165	86
17.08	0.02268	34.63	2.50E-05	0.55	-0.09138	170	96
17.44	0.02158	34.94	1.54E-05	0.58	-0.06631	173	89
17.84	0.01675	35.40	1.18E-05	0.61	-0.04051	178	82

---

 | Show Table

DownLoad: CSV

During cone penetration, the contact friction angle at the soil-cone interface depends on the cone roughness ($R_a$) and soil particle size. Currently, CPT friction sleeves are manufactured following the standards set by ^[62], with a roughness of $R_a = 0.4 \pm 0.25$ $\mu$m, measured in the longitudinal direction. In this study, the weighted average particle size ($D_{av}$), estimated from the grain size distribution results presented by ^[14], is approximately 1000 $\mu$m. Based on these characteristics, the interface efficiency ($\delta_{cv} / \varphi_{cv}$) is determined to be 0.3 ^[63], resulting in an interface friction angle ($\delta_{cv}$) of 10$^{°}$.

4.5.   Simulation results

The results of the CPT simulations for each individual depth are presented in Figure 10. For a total penetration depth of 0.5 m, the $q_t$ values reached a stationary state at approximately 0.25 m (see Figure 10a), and the actual simulation results are based on the average values obtained beyond this depth. However, a different trend is observed for $f_s$, as most simulations do not reach a stationary state even after 0.5 m of penetration (see Figure 10b). These findings highlight the need for further investigation into the $f_s$ results, potentially by extending the penetration depth beyond 0.5 m with additional stabilization methods to reduce accumulated error or refining the existing soil-cone interaction formulation. This aspect remains an ongoing area of research.

Figure 10.  GPFEM simulation output result (a) $q_t$ and (b) $f_s$.

DownLoad: Full-Size Img PowerPoint

The simulation results over depth (green) are presented in Figure 11. These results are compared with in situ measurements at point OYSC20 (black) and the APD framework input, which represents a thin 4 cm layer (red). While the simulated tip resistance is slightly overestimated, with an average deviation of 3.5 MPa, the overall pattern and magnitude are well captured. Additionally, the simulated pore water pressure during the cone penetration generally aligns with the theoretical hydrostatic line and remains within the range of in situ measurements.

Figure 11.  Comparison between OYSC20 and the results obtained from G-PFEM (a) $q_t$, (b) $u_2$.

DownLoad: Full-Size Img PowerPoint

Most of the simulations indicate fully drained conditions. However, two simulation points (z = 16.04 m and 14.76 m) exhibit slightly undrained behavior during penetration, with $u_2$ values slightly exceeding the theoretical hydrostatic line and lower $q_t$ and $f_s$ values. This discrepancy arises because, in the simulation, the cone penetrates a uniform soil layer with constant parameters and permeability throughout the process. As a result, pore water pressure accumulates around the cone and dissipates slowly due to the lower permeability.

As penetration progresses and reaches a stationary state, the significant stress zone beneath the cone tip extends approximately five times the cone diameter (±18 cm), as illustrated in Figure 12. This suggests that the soil within this zone continues to influence cone resistance measurements. However, accurately capturing this behavior is challenging due to the limitations of the current simulation technique, which cannot incorporate soil layers with varying properties within a single simulation.

Figure 12.  Influence zone during penetration (a) mean effective stress (p'), and (b) deviatoric stress (q).

DownLoad: Full-Size Img PowerPoint

4.6.   Comments related to Section 4

In Section 4, an attempt was made to evaluate the constitutive model parameters obtained from APD. The study's objective was to employ a CPT from the test site in APD to determine constitutive model parameters, and subsequently utilize these material sets to simulate CPT measurements at several depths. The simulated CPT was then compared to the CPT at the site, enabling the evaluation of the material sets obtained from APD.

The present study encountered several challenges. First, the constitutive model was limited to CASM, as it is the only model available within the framework of G-PFEM. Second, determining CASM parameters for sands using CPT was not straightforward, as there are not a lot of linkages between the parameters and CPT results, which is a requirement for computing parameters in APD. Additionally, the simulation results for $f_s$ exhibited uncertainties concerning the convergence of the results.

In view of the challenges described above, it is evident that further investigations and studies are required into these aspects. The ongoing research includes the following: the employment of other constitutive models (for example, an extension to CASM that is suitable for dense sands); the investigation of other methods for determining CASM parameters for sands; the exploration of the possibilities of determining the parameters using other workflows of APD (DMT and shear wave velocity measurements); and the investigation of the simulation results for $f_s$ for sands. Despite the challenges, it is believed that this approach could enhance the validation of determining the constitutive model parameters from CPT results.

One example of the ongoing investigations is the use of APD's shear wave velocity ($V_s$) workflow to determine CASM parameters. As previously discussed (see Section 4.1.1), $\kappa$ can be derived from $V_s$ measurements based on the ratio between the shear modulus ($G$) and the small-strain shear modulus ($G_0$). This approach yields different $\kappa$ values instead of assuming a typical 0.005 for sands. From APD's $V_s$ workflow output, two $\kappa$ values—0.0065 and 0.008—were selected for further comparison using a reduced number of simulated depths (13 instead of 36). The results of this analysis are presented in Figure 13. Figure 13(a-b) show the triaxial test results (n = 5, r = 1000) for varying $\kappa$, while Figure 13c illustrates the obtained $q_t$ profile, demonstrating that increasing $\kappa$ leads to a further reduction in $q_t$ values, thereby improving the match with in situ measurements. This study also underscores the benefits of expanding the APD system to incorporate additional in situ tests.

Figure 13.  (a) $\varepsilon_1-q$, (b) $\varepsilon_1-\varepsilon_v$, isotropic compression $p' = 150$ kPa, (c) $q_t$ profile.

DownLoad: Full-Size Img PowerPoint

5.   Conclusions and next steps

This paper presents an automated parameter determination (APD) framework that employs a graph-based methodology to derive parameters from in situ tests. This approach is particularly beneficial in the early stages of a project when soil data availability is limited. At this stage, cost-effective field tests like CPT and DMT are typically conducted before comprehensive laboratory investigations. By integrating APD into the preliminary design phase, users can obtain more detailed insights efficiently. The objective is not to replace laboratory testing but to complement it, as laboratory analyses remain crucial for refining soil properties and constitutive model parameters in the final design. APD is characterized by transparency and adaptability, enabling users to contribute their expertise and expand the system's database by integrating additional methods and parameter correlations.

Section 3 illustrated the process of determining soil parameters using the CPT workflow. The results obtained from various methods were compared with reference values derived from laboratory tests conducted at the Øysand test site. This comparison serves to validate individual methods and refine the compiled database of correlations. The results were presented in terms of lower and upper bounds, as well as the average. However, using the average as a representative value for different parameters is problematic, as it may include contributions from less accurate methods. The main challenge lies in the large number of methods available within APD and the significant variability in the estimated values for the same parameter. To address this, a statistical module is currently being developed to support the selection of a more reliable representative value. To further support the selection of reliable profiles for soil characterization, the statistical module will also enable users to simulate element-level tests using different material sets derived from the available statistical measures, thereby allowing an assessment of how parameter selection influences mechanical response.

Section 4 explored the evaluation of constitutive model (CASM) parameters obtained in Section 3 by numerically simulating a CPT and comparing the results with in situ CPT data, which served as input for APD. Several challenges were encountered during this study, including the limited number of methods available for deriving CASM parameters for sands from CPT and uncertainties in the simulation results for $f_s$. Despite these challenges, the initial results for $q_t$ at different depths, as shown in Figure 11a, demonstrate that this approach can be used to assess the derived material parameters. It is anticipated that addressing the aspects discussed in Section 4.6 will enhance the accuracy of the simulations and further validate the obtained material sets.

Author contributions

Conceptualization, I.M., A.E.W., H.F.S. and F.T.; methodology, I.M., A.E.W., H.F.S. and F.T.; software, I.M. and A.E.W.; validation, I.M., A.E.W., H.F.S. and F.T.; formal analysis, I.M., A.E.W., H.F.S. and F.T.; investigation, I.M. and A.E.W.; resources, I.M., A.E.W., and F.T.; data curation, I.M.; writing---original draft preparation, I.M. and A.E.W.; writing---review and editing, I.M., A.E.W., H.F.S. and F.T.; visualization, I.M. and A.E.W.; supervision, H.F.S. and F.T.; project administration, F.T.

Use of AI tools declaration

The authors declare they have not used Artificial Intelligence (AI) tools in the creation of this article.

Acknowledgments

The authors extend their gratitude to the members of the APD group, Arny Lengkeek, and Ronald Brinkgreve, for their valuable contributions. They also sincerely appreciate the support and collaboration of the research team at CIMNE. Special thanks go to Dr. Laurin Hauser for providing access to G-PFEM and for the insightful discussions.

Conflict of interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.