Beyond Surface Modulus: Shape-Aware Crack-Index Assessment of Mass Concrete Foundations Under Early-Age Thermal Stress

2.1. Transient Heat Conduction with Internal Heat Source

The three-dimensional transient temperature field T(x, y, z, t) in hydrating concrete is governed by the Fourier heat equation with an internal heat-generation term (Carslaw & Jaeger, 1959):

• ρ - the density (kg m⁻³),

• c - the specific heat (J kg⁻¹ K⁻¹),

• k - the thermal conductivity (W m⁻¹ K⁻¹),

• q - the volumetric heat-generation rate.

The volumetric heat q is obtained from an adiabatic temperature-rise curve of exponential form, widely adopted in MIDAS Civil and in the Korean Concrete Institute specifications (Korea Concrete Institute, 2009):

• K - the ultimate adiabatic temperature rise (°C),

• α - the reaction-rate coefficient (day⁻¹),

• t - the time (day).

In this study, K = 51.0 °C and α = 1.516 day⁻¹ are adopted, consistent with a mix containing 350 kg m⁻³ of ordinary Portland cement (OPC). For the cement-content sensitivity study (Section 3.7), K is scaled proportionally to the binder content following K = K_ref(C/C_ref), with α held constant. This proportional scaling is consistent with the maturity-based formulation of Ballim and Graham (2004), who derive total heat release from adiabatic-calorimeter binder fingerprints and report peak heat rates of approximately 3 W kg⁻¹ at t₂₀ ≈ 12 h for typical Portland cements. The adopted α gives a heat-release time scale of 1/α = 0.66 day, so its maximum heat-production window falls within the same first-day hydration interval rather than implying a later or slower heat source.

2.2. Boundary Conditions

Convective heat exchange with the surrounding environment is imposed on each exposed face through Newton's law of cooling:

Two distinct convection regimes are prescribed in the present model. The first (Type-1) corresponds to free surfaces that exchange heat directly with the ambient air, while the second (Type-2) corresponds to surfaces in contact with formwork or curing mats. The reference convection coefficients are sourced from the Korean Concrete Institute specifications (Korea Concrete Institute, 2009): h_c = 12 kcal m⁻² h⁻¹ °C⁻¹ for steel formwork in contact with water and h_c = 4.5 kcal m⁻² h⁻¹ °C⁻¹ for curing mats. The latter is equivalent to h_c ≈ 5 W m⁻² °C⁻¹, in close agreement with the timber-formwork value reported by Ballim and Graham (Ballim & Graham, 2004) under light-wind conditions. For surfaces exposed to direct air, Ali and Urgessa (Ali & Urgessa, 2012) proposed a wind-dependent expression h_c = 20 + 14w (kJ m⁻² h⁻¹ °C⁻¹) for wind speeds w ≤ 5 m s⁻¹, which is consistent with the Korean reference range used here.

2.3. Mechanical Properties and Thermal Stress

To capture the strong age dependence of early-age concrete mechanical properties, the compressive-strength gain follows the ACI 2092R formula (ACI 209.2R-08, 2008):

The elastic modulus and the splitting tensile strength are derived from the ACI 318 (ACI 318-19, 2022) correlations:

Figure 1:

Figure 1 illustrates the two principal time-dependent relationships underlying the hydration-heat and mechanical-property sub-models. Subplot (a) shows the exponential adiabatic temperature-rise curve T_ad(t) = 51(1 − e^−1.516t), which reaches 90 % of the ultimate heat release K = 51 °C within approximately 1.5 days. Subplot (b) shows the normalised evolution of compressive strength f_c(t)/f_c28, elastic modulus E_c(t)/E_c28, and splitting tensile strength f_sp(t)/f_sp28 per Equations (4) – (6); all three properties develop rapidly in the first three days, confirming that the thermal-cracking window (typically 12–120 h) coincides with the steepest gradient of mechanical-property gain.

2.4. Thermal Crack Index and Acceptance Criteria

The KCI Standard Specifications (Korea Concrete Institute, 2009) define the thermal crack index as:

• f_sp(t) – the age-dependent splitting tensile strength from Eq. (6) (or, when available, from a best-fit regression of laboratory data),

• σ_th(t) ≡ σ_t,max(t) – the maximum principal tensile stress at age t extracted directly from the 3-D thermal–stress analysis.

The Korean code adopts the three-band acceptance criterion summarised in Table 1.

The same three-band acceptance scheme has been adopted by the Japan Concrete Institute (Japan Concrete Institute, 2016; Sato et al., 2013), whose Guidelines further translate the crack index into a quantitative thermal-cracking probability through the Weibull-type expression P(I_cr) = [1 − exp{−(I_cr/0.92)^−4.29}] × 100, calibrated on 65 mass-concrete structures and 728 individual members. The 5% probability threshold corresponds approximately to I_cr ≈ 1.85, broadly consistent with the I_cr ≥ 1.5 “prevent cracks” band employed in the Korean code (Korea Concrete Institute, 2009), and is supported empirically by Ashizawa et al. (Ashizawa et al., 2019), whose 130-section field survey of box-culvert side walls reported a mean I_cr of 1.07 in cracked sections versus 1.10 in uncracked sections, with overall coefficients of variation of 13.4% for tensile stress and 7.4% for tensile strength. Based on project requirements and engineering judgement for WTG foundations, the following combined acceptance limits are adopted herein:

This acceptance table establishes the interpretation framework used throughout the Results and Discussion sections. In particular, the distinction between the prevent-cracks, limit-cracks, and allow-cracks bands provides the basis for evaluating whether the simulated responses remain acceptable under variations in geometry, placement temperature, ambient conditions, convection, and binder content. The subsequent discussion of I_cr,min therefore refers to Table 1 as the common acceptance benchmark.

2.5. Surface Modulus Definition

The surface modulus M is defined as the ratio of the area of surfaces subjected to environmental cooling A_cool (m²) to the total concrete volume V (m³) (Flaga & Klemczak, 2016; Klemczak & Żmij, 2021):

The cooling surface area is the area of faces that can exchange heat with the surrounding environment under the adopted boundary conditions. Faces cast against soil, against previously hardened concrete, or against heavily insulated boundaries are excluded from A_cool. For a cube of side length L resting on soil with five exposed faces, M = 5L² / L³ = 5/L. The same convention is applied consistently in Table 3: the bottom face is excluded for the slab and pile-cap cases (G5, G6 and G8), while the column case G7 includes the four vertical faces and one exposed end face, giving A_cool = 52 m². This convention keeps A_cool tied to the thermal boundary condition rather than to the geometric face count alone. The metric allows cubes, slabs, pile caps and columns to be compared on a common scale and is consistent with the CIRIA C766 mass-concrete classification M < 2 m⁻¹ (Flaga & Klemczak, 2016). Physically, M represents the heat-dissipation area available per unit concrete volume; lower M indicates a thermally massive member with greater core heat retention.

The enhanced massivity index of Klemczak and Smolana (Klemczak & Smolana, 2024), m_cor = m_s/(k_f·k_b·k_ΔT), augments Eq. (9) with binder-, content- and temperature-correction factors k_f, k_b and k_ΔT. The present study deliberately retains the simpler geometric form of M because T₀, T_amb and cement content are swept explicitly in separate parametric series, allowing their influence on I_cr,min to be observed rather than embedded in correction coefficients. Consequently, the proposed M_crit should be read as a condition-specific threshold for the stated mix, T₀, T_amb and convection assumptions, not as a universal constant. In applications where binder type, binder content or initial temperature depart substantially from the baseline, the m_cor framework provides a useful complementary screening variable.

2.6. Post-processing Pipeline and Reference Setup

All MIDAS Civil analyses reported in this paper use SI units throughout: lengths in metres (m), forces in kilonewtons (kN), temperatures in degrees Celsius (°C), thermal conductivity in W m⁻¹ K⁻¹, and convection coefficients in kcal m⁻² h⁻¹ °C⁻¹ as required by the Korean Concrete Institute input format (1 kcal m⁻² h⁻¹ °C⁻¹ ≈ 1.163 W m⁻² K⁻¹). The Python post-processor converts all MIDAS Civil outputs to SI before computing derived quantities. The post-processing routine, implemented in Python, parses each MIDAS Civil output dataset and reports four quantities at every solver step: (a) the temperatures recorded at the seven monitoring nodes, (b) the maximum principal tensile stress within the concrete domain, (c) the instantaneous core-to-surface temperature differential ΔT_cs(t_i) and its temporal maximum ΔT_cs,max, and (d) the crack index I_cr(t) evaluated using the time-dependent material models in Eqs. (4)–(7). The notation ΔT_cs,max is adopted to distinguish the governing core-to-surface temperature differential from a generic thermal gradient.

A 3 × 3 × 3 m concrete block resting on a 9 × 9 × 1.5 m sub-soil layer is modelled as the reference mock-up, representative of WTG foundations constructed in Vietnam (V. Q. Nguyen et al., 2022). Exploiting symmetry, only one quarter (1.5 × 1.5 × 3 m concrete and 4.5 × 4.5 × 1.5 m soil) is meshed into 2376 eight-node solid elements at a 250 mm cubic reference element size. Table 2 reports the baseline material and thermal properties.

Mechanical restraint and soil-foundation interaction are represented through an explicit sub-soil layer beneath the concrete block so that early-age stress develops from the coupled stiffness of the concrete and its support rather than from a free-standing body. Symmetry planes are restrained in the normal direction, the base of the soil layer is restrained vertically, and the lateral soil boundaries are constrained only to the extent required to eliminate rigid-body motion. Perfect thermal contact and bonded mechanical interaction are assumed at the concrete-soil interface; interface slip, uplift, creep relaxation, and construction-joint debonding are not modelled explicitly. These assumptions are appropriate for comparative screening, but they may influence σ_t,max and, consequently, I_cr in design cases involving weak subgrade restraint or intentionally introduced sliding layers.

The baseline properties in Table 2 combine a moderate-strength concrete, explicit sub-soil support, and reference thermal boundary conditions representative of practical foundation casting. Because all later sensitivity studies modify only one parameter group at a time relative to this baseline, Table 2 is essential for maintaining internal consistency when comparing Tables 4–10 and for ensuring that observed changes in T_max,core, ΔT_cs,max, and I_cr,min can be attributed to the intended variable rather than to simultaneous shifts in the reference setup.

Seven monitoring nodes are extracted from each MIDAS Civil run: N295 (core centre), N300 (side centre), N335 (corner centre), N540 (top centre), N545 (side-top), N580 (corner-top) and N50 (base centre). At every time step t_i, the instantaneous maximum gradient is computed as:

For the baseline case, the post-processor returns T_max,core = 65.16 °C at t = 62 h (node N295) and ΔT_cs,max = 26.17 °C at t = 86 h, governed by the core-to-corner-top node pair (N295–N580). The seven-node check confirms that the corner-top node, rather than the geometric side face, governs the critical surface temperature, consistent with the bell-shaped temperature distribution reported by Klemczak and Żmij (Klemczak & Żmij, 2021) for restrained foundation slabs.

2.7. Surface-modulus Parametric Matrix

Eight representative geometries spanning M = 1.00–2.50 m⁻¹ are defined in Table 3. Each geometry uses the same baseline material properties, boundary conditions, and mesh density (250 mm reference element size), so differences in T_max,core, ΔT_cs,max, and I_cr,min can be attributed to the combined effects of geometric scale, aspect ratio, and cooling-surface layout. The updated matrix provides overlapping-M comparisons rather than exact same-M pairs: G4 and G6 both have M = 1.00 m⁻¹, G3 and G5 are close at M = 1.25 and 1.33 m⁻¹, G8 and G2 lie at M = 1.50 and 1.67 m⁻¹, and G7 and G1 occupy the high-M range at M = 2.17 and 2.50 m⁻¹. These comparisons allow a direct test of whether the surface modulus alone captures geometry effects or whether an additional shape- or aspect-ratio-related descriptor is required.

The G-series matrix in Table 3 is structured to test whether geometries with similar or identical surface modulus values also produce similar thermal responses. The matrix is therefore not intended as a simple size sequence, but as a deliberate cross-geometry comparison set in which cubes, slabs, columns, and pile caps partially overlap in M. This design is important because it allows the later analysis in Section 3.8 to separate the effect of massivity from the effect of shape and heat-flow path length.

3.1. Baseline Thermal Response

The baseline configuration (T₀ = 15 °C, T_amb = 27.1 °C, h_c = 4.5 kcal m⁻² h⁻¹ °C⁻¹, Δt = 2 h, element size 250 mm) predicts T_max,core = 65.16 °C at t ≈ 62 h and ΔT_cs,max = 26.17 °C at t ≈ 86 h. Figure 2 shows the temperature–time histories at the seven monitoring nodes, and Figure 3 illustrates the temperature contour of the quarter model at the instant of peak core temperature.

The interior node N295 attains its maximum at approximately 62 h and then cools gradually, whereas the top-face and corner-top nodes heat up and cool more rapidly because of convective heat loss to the ambient air. The separation between the warmest curve (N295) and the coolest curve (N580) represents the instantaneous core-to-surface temperature differential defined in Eq. (10), and its maximum at approximately 86 h corresponds to ΔT_cs,max. Figure 3 confirms the expected bell-shaped temperature distribution: the highest temperature is concentrated in the core of the block, while the steepest gradient develops towards the exposed top corner. This response pattern is broadly consistent with the experimental observations of Le et al. (H. H. Le et al., 2020) on a 2.5 m cubic mock-up and with the field-instrumented diaphragm-wall study of Grabowski and Mitew (Grabowski & Mitew-Czajewska, 2021).

Figure 2:

Figure 3:

3.2. Effect of Time Step Δt

To assess temporal convergence, the analysis was repeated with five-time steps from 0.5 h to 8 h (cases A1–A5). Table 4 reports T_max,core, ΔT_cs,max, and I_{cr, min.}.

Results with Δt ≤ 2 h are effectively converged: cases A1 and A2 differ by only 0.09 °C in T_max,core. Coarser time steps overpredict the peak response; the 8 h case (A5) overestimates T_max,core by almost 9 °C and reduces I_cr,min from 1.793 to 1.284, thereby shifting the assessment from the “prevent cracks” band to the “limit cracks” band. A time step of 2 h is therefore recommended for this class of problem and is consistent with the implicit-explicit stability condition reported by Ballim and Graham (Ballim & Graham, 2004), who require Δt ≤ ρc_pΔ²/[4k(1 + Bi)] at surface nodes.

3.3. Effect of Element Size

Spatial-discretisation convergence was investigated using six characteristic element sizes ranging from 500 mm to 150 mm (cases B1–B6; Table 5). The peak core temperature T_max,core is nearly insensitive to mesh size, with only 0.5 °C separating the 500 mm and 150 mm meshes. By contrast, ΔT_cs,max and I_cr,min remain mesh-sensitive because the governing corner-top gradient is highly localised. Mesh refinement increases ΔT_cs,max from 17.50 °C (B1) to 29.38 °C (B6), while I_cr,min decreases from 1.793 at B4 to 1.396 at B6. The 250 mm mesh should therefore be regarded as a practical reference mesh for the comparative parametric campaign rather than as a fully converged crack-index solution. The B6 result indicates that final design checks for cases close to the I_cr = 1.2 or 1.5 acceptance limits should be repeated with an element size of approximately 150 mm, or with an explicit mesh-convergence correction applied to ΔT_cs,max and I_cr,min. Comparable mesh-density requirements were also reported by Ahorro and Estores (Ahorro & Estores, 2024) for probabilistic thermal-cracking simulations, who used 100 mm hexahedral elements for 0.5 m cube analyses in DIANA 10.3.

3.4. Effect of Convection Coefficient h_c

The convection coefficient on the exposed top face was varied from 1.5 to 12 kcal m⁻² h⁻¹ °C⁻¹ (cases C1–C4; Table 6). The results show an inverse relation between surface-cooling intensity and thermal-crack safety. Increasing h_c reduces T_max,core only marginally, by 0.96 °C between C1 and C4, because peak core temperature is governed primarily by internal heat generation. In contrast, ΔT_cs,max increases from 16.55 °C (C1) to 31.03 °C (C4), while I_cr,min decreases from 2.471 to 1.537, shifting the response from the “prevent cracks” band to the “limit cracks” band. This indicates that passive thermal insulation, achieved for example through insulated formwork or curing mats, is a more reliable crack-control measure than aggressive surface convection. Su et al. (Su et al., 2018) similarly recommended limiting the core-to-surface temperature differential to approximately 25 °C through a combination of insulation, controlled cooling, and low-heat materials.

3.5. Effect of Placement Temperature T₀

The placement (fresh) concrete temperature was swept from 15 °C up to 35 °C in five increments (cases D1–D5), with T_amb = 27.1 °C fixed (Table 7).

Increasing the placement temperature by 20 °C produces an almost one-to-one increase of 19.12 °C in T_max,core, because the adiabatic hydration contribution remains approximately unchanged while the initial temperature shifts the entire thermal history upward by T₀. For T₀ ≥ 25 °C, the 70 °C DEF threshold is exceeded and I_cr,min falls below 1.5. At T₀ = 35 °C, I_cr,min = 1.302, only marginally satisfying the project requirement I_cr ≥ 1.2. Pre-cooling of aggregates or the use of ice-water mixing should therefore be regarded as necessary measures for hot-weather mass-concrete pours in tropical climates. The sensitivity obtained here is also consistent with Le, Vu and Ho (V. M. Le et al., 2025), whose 57-case FE regression identified placement temperature as the dominant first-order term for T_max and the second most influential term for ΔT_{cs, max.}. Tasri (Tasri, 2025) independently reached a comparable placement-temperature limit from a V/A-based analysis.

3.6. Effect of Ambient Temperature T_amb

Ambient temperature was varied from 20 °C to 35 °C while keeping T₀ = 15 °C constant (cases E1–E5; Table 8). The influence of T_amb on T_max,core is negligible, with a total variation of only 0.77 °C between E1 and E5. Its effect on ΔT_cs,max, however, is substantial: the governing temperature differential decreases from 31.02 °C at 20 °C ambient temperature to 20.84 °C at 35 °C. The corresponding I_cr,min increases from 1.614 to 2.033. Cold-weather pours therefore remain a non-trivial cracking risk even when T_max remains acceptable, because the larger ambient-to-core differential induces higher tensile stress near the surface. This outcome is consistent with Abdel et al. (Abdel-Raheem et al., 2018), who recommended additional surface protection under cooler ambient conditions to prevent thermal-shock cracking caused by rapid surface cooling.

3.7. Effect of Cement Content

To probe the role of internal heat-generation intensity, the binder dosage was swept across four levels covering 300 to 450 kg m⁻³ (cases F1–F4, Table 9). The adiabatic temperature-rise parameter K was adjusted proportionally according to K = K_ref(C/C_ref) with α = α_ref and C_ref = 350 kg m⁻³.

The F-series demonstrates that binder content is among the strongest physical drivers of thermal-cracking risk in mass concrete. Over the parametric range investigated, a near-linear increase in T_max,core is observed, with a sensitivity of approximately 7.1 °C per 50 kg m⁻³ increment. At C = 350 kg m⁻³, T_max,core = 65.16 °C is comfortably below the 70 °C DEF threshold; at C = 400 kg m⁻³ the peak reaches 72.25 °C, exceeding the threshold even at T₀ = 15 °C; at C = 450 kg m⁻³ it reaches 79.35 °C. The minimum crack index follows a corresponding decrease, from I_cr,min = 2.16 at C = 300 kg m⁻³ down to 1.33 at C = 450 kg m⁻³. Linear interpolation places the I_cr,min = 1.5 limit at C ≈ 405 kg m⁻³ and the I_cr,min = 1.2 limit at C ≈ 470 kg m⁻³.

This finding is highly consistent with the case-study evidence reported for high-performance mass concrete columns and slabs. Kim and Lee (J.-H. Kim & Lee, 2004) documented two 2.4 m diameter Petronas Twin Tower column mock-ups: a binary OPC + silica-fume blend at 535 kg m⁻³ produced T_max = 91.6 °C and a core–corner gradient ΔT = 57.5 °C that resulted in visible thermal cracks, whereas a ternary OPC + 20 % PFA + silica fume blend at the same total binder content produced T_max = 87.0 °C, ΔT = 52.9 °C and no cracking. Woo et al. (Woo et al., 2010) developed a low-heat ternary cement (OPC + GGBS + FA) for a 4 × 4 × 4 m mock-up that achieved T_max = 60.0–64.0 °C and ΔT = 13.7–15.0 °C — both below the JCI 20 °C limit. For thermally critical mass-concrete pours, the present F-series therefore provides a quantitative upper bound: cement content should be capped at C ≤ 400 kg m⁻³ unless paired with low-heat binders or active cooling.

3.8. Effect of Surface Modulus M

All eight geometries in Table 3 were analysed using the baseline material properties and boundary conditions. Table 10 reports the key thermal outputs, and Figure 4 plots M against T_max,core, ΔT_cs,max, and I_cr,min. The 3 × 3 × 3 m reference cube (G2) reuses the baseline values from Section 3.1.; the remaining cases (G1, G3–G8) are additional MIDAS Civil runs generated through the automated framework described in Section 2.6.

Figure 4:

Dashed lines represent simple one-variable M trends fitted to the updated full G-series dataset in Table 10. The low R² values indicate that M alone is not a reliable shape-independent predictor across mixed geometries.

The G-series results in Table 10 link the cooling-surface definition in Table 3 to the stress-based crack index in Eq. (7), and the P_cr,JCI column converts each I_cr,min value into the JCI thermal-cracking probability using P(I_cr) = 1 − exp[−(I_cr/0.92)^−4.29] × 100. M-only regressions were examined first because the surface modulus is widely used as a massivity descriptor; however, their R² values are low for the mixed G-series dataset (0.227 for T_max,core, 0.080 for ΔT_cs,max, and 0.003 for I_cr,min). This indicates that higher-order polynomials in M would merely fit scatter without capturing the missing physical controls. A more defensible set of response-specific regressions is therefore adopted.

RMSE, MAE, and leave-one-out cross-validation errors are reported because the G-series contains only N = 8 geometries; accordingly, the models should be interpreted as empirical screening relationships rather than universal design equations. In Equations (12) – (13), d_min denotes the smallest member dimension and represents the shortest characteristic heat-flow path. The quadratic term in Eq. (12) captures the expected saturation of T_max,core with increasing thickness, whereas Eq. (13) combines heat-flow length (d_min) with cooling-area density (M). Equation (14) adopts an inverse form in ΔT_cs,max because I_cr is defined as a strength-to-stress ratio and the thermally induced tensile stress is approximately proportional to the governing temperature differential. This form is also physically more defensible than a purely empirical second-order polynomial fitted to I_cr,min versus M and provides the strongest interpretable fit for the present G-series (R² = 0.951).

The model-performance comparison in Table 11 confirms that the best-performing models are also the most physically interpretable. In particular, the strong performance of the d_min-based and ΔT_cs,max-based relations, together with the weak LOOCV performance of the M-only fits, reinforces the central argument of the paper: the surface modulus is useful for preliminary screening, but design-relevant crack assessment must remain response-specific and shape-aware. The consistency between statistical performance and mechanistic interpretation strengthens the credibility of Eqs. (12) – (14) as screening tools.

The non-cubic shapes deviate markedly from the M-only cube trend even when their M values overlap the cubic range. At the same M = 1.00 m⁻¹, the very large cube G4 yields I_cr,min = 1.095, whereas the thick slab G6 yields I_cr,min = 2.553; the corresponding ΔT_cs,max values are 33.67 °C and 23.00 °C, respectively. A similar divergence is observed between the large cube G3 (M = 1.25 m⁻¹, I_cr,min = 1.295) and the thin slab G5 (M = 1.33 m⁻¹, I_cr,min = 4.523). The stronger I_cr,min–ΔT_cs,max relation in Eq. (14) therefore indicates that the crack index is governed more directly by the realised thermal differential than by M alone.

These updated regressions refine the earlier A/V-based interpretation of Le, Vu and Ho (V. M. Le et al., 2025). Their conclusion that T_max and ΔT_cs,max are governed primarily by A/V is consistent with the cubic subset; however, the full G-series shows that combining cubes, slabs, columns, and pile caps require at least one additional geometric scale descriptor such as d_min. When the response variable changes to I_cr,min, the controlling correlation becomes gradient-based rather than M-based. This observation is consistent with the JCI/KCI crack-index framework and with enhanced massivity approaches such as that of Klemczak and Smolana (Klemczak & Smolana, 2024), in which the simple surface modulus is used as a pre-design screening parameter rather than as a complete predictor of cracking risk.

The inverse model I_cr,min = −1.253 + 77.323/ΔT_cs,max gives R² = 0.951 and is physically consistent with the crack-index definition as a strength-to-thermal-stress ratio.

This relationship further clarifies why the crack index is more consistently interpreted through the realised thermal differential than through the surface modulus alone (Figure 5). The data follow a clear monotonic inverse trend: geometries with larger ΔT_cs,max consistently exhibit lower I_cr,min, indicating lower thermal-cracking safety. From a practical standpoint, once the governing core-to-surface temperature differential has been obtained from three-dimensional analysis, Equation (14) offers a more direct and physically interpretable screening estimate of crack risk than an M-only rule for mixed geometries.

Figure 5:

4.1. Overall Sensitivity Ranking

The sensitivity ranking in Table 12 summarises the influence of all investigated parameters on T_max,core and I_{cr, min.}. The metrics ΔT_max,core and ΔI_cr,min denote the maximum absolute variation in peak core temperature and minimum crack index over each parameter range. The G-series is reported separately for the cubic-only M trend and for the full geometry/surface-modulus portfolio so that the screening regressions are not overgeneralised.

Among the physically controllable inputs, binder content and placement temperature emerge as the dominant drivers of T_max,core, whereas the convection coefficient and geometry-related effects dominate I_cr,min. The large variation, ΔI_cr,min = 3.428, in the full G-series is driven mainly by the non-cubic slab response, particularly G5. The revised regressions clarify this behaviour: T_max,core is controlled primarily by d_min, ΔT_cs,max by the combined M–d_min descriptor, and I_cr,min by the realised ΔT_cs,max rather than by M alone.

This sensitivity ranking aligns broadly with Le, Vu and Ho (V. M. Le et al., 2025), who reported F-values for T_max in the order T_p > r_t > T_c > W_c,u in their two-layer DOE, and with Zhang et al. (Zhang et al., 2024), who showed that moulding temperature, wind speed and insulation thickness control thermal stress in super-long mass slabs. The present ranking extends those findings by adding a 3-D geometry sweep at fixed material parameters and by separating the cubic M trend from the full cross-geometry response of I_{cr, min.}.

The ranking in Table 12 also helps reconcile the different roles of M, d_min, and ΔT_cs,max across the manuscript. No single scalar descriptor captures every response equally well: T_max,core is governed mainly by heat-generation intensity and characteristic thickness, whereas I_cr,min is linked more directly to the realised temperature differential and the associated restraint-sensitive stress state. This interpretation is consistent with the regression hierarchy in Table 11 and supports the use of M as a screening parameter rather than as a stand-alone crack-control criterion.

Accordingly, Tables 11 and 12 should be read as evidence for the relative usefulness of different screening descriptors, not as a substitute for explicit three-dimensional thermal-stress analysis when acceptance decisions depend on project-specific geometry, restraint, and material behaviour. Physically, the contrasting roles of M, d_min, and ΔT_cs,max can be understood from the governing heat-transfer and stress mechanisms. The peak core temperature T_max,core is set by the balance between internal heat generation (proportional to binder content and K) and heat dissipation through the minimum dimension d_min; once the characteristic diffusion length exceeds roughly 1.5–2 m, the core becomes effectively insulated and T_max,core saturates, which explains the quadratic form of Eq. (12). The core-to-surface temperature differential ΔT_cs,max, by contrast, is governed by the competition between the heat source intensity in the core and the convective heat loss at the surface: a larger cooling-surface density (higher M) accelerates surface cooling, while a longer heat-flow path (larger d_min) delays heat arrival at the surface, so both descriptors appear in Eq. (13). The crack index I_cr,min is a strength-to-stress ratio in which the tensile stress is driven primarily by the thermal incompatibility between the contracting surface and the expanding core; because this incompatibility is quantified directly by ΔT_cs,max, the inverse relation in Eq. (14) is not merely empirical but mechanistically grounded. Shape enters as a second-order effect through the restraint stiffness and the three-dimensional stress redistribution: a squat pile cap and a slender column can share the same M yet generate very different stress fields because their aspect ratios, free-face boundaries, and restraint configurations differ. This is why M alone cannot serve as a crack-control criterion for mixed-geometry design.

4.2. Practical Engineering Implications

Based on the parametric results obtained herein, the following practical recommendations are proposed for mass-concrete foundations cast under hot-humid tropical conditions. These recommendations should be regarded as screening-level guidance only, because the regression equations were calibrated on N = 8 geometry cases with a single baseline material system; project-specific or commercial application therefore requires recalibration using additional geometries, refined meshes, and locally measured material properties.

• (a) For the reference 3 × 3 × 3 m geometry, placement temperature should be limited to T₀ ≤ 25 °C. Above this level, the 70 °C DEF threshold is exceeded, and pre-cooling of aggregates or ice-water mixing becomes necessary. This operational limit is also consistent with Tasri (Tasri, 2025), whose V/A = 1.5 m⁻¹ curve gives a maximum allowable casting temperature of approximately 27.2 °C for ΔT ≤ 20 °C.

• (b) For thermally critical pours, cement content should preferably be limited to C ≤ 400 kg m⁻³ when placement temperature is at or above 15 °C. Beyond this level, I_cr,min falls below 1.5 in the present parametric study. Additional strength demand should therefore be satisfied, where possible, through supplementary cementitious materials such as GGBS, fly ash, or silica fume, which reduce the effective heat of hydration, as illustrated in (J.-H. Kim & Lee, 2004; Woo et al., 2010).

• (c) Insulated formwork, for example a steel form combined with a curing mat and an effective h_c of approximately kcal m⁻² h⁻¹ °C⁻¹, is one of the most effective passive measures for controlling ΔT_cs,max and maintaining I_cr above 1.5. Increasing surface convection may slightly reduce T_max, but it can simultaneously increase the core-to-surface temperature differential and reduce I_cr,min; aggressive surface cooling should therefore not be assumed to improve crack control automatically.

• (d) For cold-weather concreting, additional surface protection should be specified because lower ambient temperature increases the core-to-surface temperature differential even when T_max remains acceptable, reducing I_cr,min by as much as 0.42 over the 20–35 °C ambient range.

• (e) For nominally cubic foundation geometries, the local cubic-subset estimate of M_crit may be used as a preliminary screen under the baseline mixture and boundary-condition assumptions. For mixed or non-cubic geometries, however, Eqs. (12)–(14) should be used only as response-specific screening relations: Eq. (12) for peak temperature, Eq. (13) for ΔT_cs,max, and Eq. (14) for crack-index response. Final acceptance should still be based on explicit three-dimensional analysis following the framework described in Sections 2.1–2.6.

• (f) A 250 mm element size is suitable for comparative screening across the parametric campaign, but design checks near an acceptance boundary should be repeated with an element size of approximately 150 mm because Table 5 shows that I_cr,min remains mesh-sensitive between cases B4 and B6.

Pre-submission consistency check. The notation is used consistently as follows: T_max,core denotes peak core temperature; ΔT_cs,max denotes the maximum core-to-surface temperature differential used in Tables 4–10 and Eqs. (13)–(14); and ΔT_max,core in Table 12 denotes the range of T_max,core over each sensitivity series, not a core-to-surface differential. The response-specific regressions are reported as empirical screening relationships calibrated on the G-series dataset and should not be read as stand-alone code equations.

4.3. Limitations and Future Work

Several limitations constrain the generality of the proposed workflow. First, the regression models were calibrated on a limited G-series dataset (N = 8) and should therefore be interpreted as condition-specific screening relationships rather than universal design equations. Second, the reference 250 mm mesh is adequate for comparative parametric ranking but does not fully converge the localised corner-gradient response governing I_cr,min; cases close to I_cr = 1.2 or 1.5 should therefore be reassessed using refined meshes or a documented convergence correction. Third, the early-age mechanical properties were estimated from ACI 209 and ACI 318 (ACI 209.2R-08, 2008; ACI 318-19, 2022) relations rather than direct measurements of the actual Vietnamese mixture; project applications should therefore incorporate locally measured E_c(t), f_sp(t), and, where possible, creep- or relaxation-related properties. Fourth, the model assumes perfect mechanical bond at the concrete–subgrade interface and does not account for viscoelastic phenomena such as creep relaxation or interface slip. In field conditions, partial bond and early-age creep typically reduce peak tensile stress relative to the fully restrained elastic estimate; conversely, weak or soft subgrade can shift the restraint pattern and alter the critical stress location. The current idealised boundary conditions are appropriate for comparative screening across geometries, but design cases involving soft subgrade, intentional sliding layers, or long-duration loading should incorporate creep-modified constitutive models. Fifth, the model has not yet been validated against a full-scale instrumented foundation with measured temperature, strain, and cracking data. These limitations do not affect the internal comparative logic of the present parametric study, but they do limit external generalisation and should be kept in mind when interpreting the proposed screening equations.

The present surface-modulus study is also limited to block-, slab-, pile-cap-, and column-type geometries. More complex configurations, such as tapered WTG foundations with haunches, embedded ducts, or staged lifts, will require additional parametric analyses. Direct validation against field data from the BT1 Wind Farm foundation should therefore be regarded as a priority before the proposed screening equations are adopted for project-level design use.

A useful next step would be to define an explicit applicability envelope for Equations (12) – (14), for example in terms of geometry family, d_min range, convection level, mesh density, and source of material properties. Stating these bounds more formally in an expanded dataset would help distinguish screening use from design-level use and reduce the risk of interpreting the proposed regressions as universal code-type equations.

Three extensions are proposed for future work. First, the workflow should be applied to a full-scale WTG foundation at the BT1 Wind Farm project (V. Q. Nguyen et al., 2022) using thermocouple histories for quantitative validation of T_max, ΔT_cs,max, and I_cr,min. Second, Monte Carlo simulation can be incorporated to convert deterministic crack-index histories into thermal-cracking probability distributions, following the probabilistic framework of Ahorro and Estores (Ahorro & Estores, 2024). Third, the methodology may be extended to concrete containing coarse recycled aggregate, for which reduced thermal conductivity may lower peak core temperature while simultaneously increasing near-surface temperature gradients. Fourth, automated parametric scaling – for example via design-of-experiments sampling coupled with machine learning surrogate models such as Gaussian process regression or gradient-boosted trees – could substantially expand the geometry coverage beyond the present eight-case G-series and enable data-driven generalisation of the proposed screening relationships to a broader class of foundation shapes and boundary conditions.

A fourth avenue for future work is a systematic comparison with the European EN 1992-3:2006 restraint-factor framework (Bamforth, 2018), which addresses early-age cracking through an external-restraint formulation complementary to the internal-restraint crack-index method adopted herein. Such a comparison would clarify conditions under which each framework is more conservative and would enhance the applicability of the proposed workflow for structures designed to European standards.