Damage Mechanics Constitutive Modelling of Bituminous Mixtures Under Complex Loading

Full Text

African Journal of Applied Mathematics: Modelling for Engineering Systems | Vol. 1 , No. 4 , 202 6 AFRICAN JOURNAL OF APPLIED MATHEMATICS: MODELLING FOR ENGINEERING SYSTEMS Vol. 1 , No. 5 , pp. 1–21, 202 6 DOI: 10. XXXXX /AJMMES.2025.044 Received: Jan uary 202 6 | Revised: January 2026 | Accepted: February 202 6 | Published: March 202 6 Damage Mechanics Constitutive Modelling of Bituminous Mixtures Under Complex Loading Aduot Madit Anhiem Department of Civil Engineering, Universiti Teknologi PETRONAS, Seri Iskandar 32610, Perak, Malaysia Correspondence: aduot.madit2022@gmail.com ORCID iD: https://orcid.org/0009-0003-7755-1011 — ◆ ◆ ◆ — Abstract Bituminous pavement mixtures are complex composite materials whose mechanical response under repeated traffic loading is governed by time-dependent viscoelasticity, progressive micro-crack initiation and propagation, and irreversible permanent deformation—phenomena that interact nonlinearly across temperature ranges spanning −20 °C to 70 °C in tropical and semi-arid climates. This paper presents a comprehensive continuum damage mechanics (CDM) constitutive framework for bituminous mixtures that couples a generalised Maxwell viscoelastic model, represented by a Prony series with N = 5 relaxation arms, with a scalar damage evolution law derived from the thermodynamic dissipation inequality. The damage variable D ∈ [0,1] evolves as a function of cumulative dissipated pseudo-strain energy density, mix-specific material parameters, and temperature through an Arrhenius shift factor. Permanent deformation (rutting) is modelled via a viscoplastic flow rule with isotropic hardening calibrated to dynamic creep test data. The constitutive model is implemented in a finite element framework using the UMAT user-material subroutine interface and validated against dynamic modulus, fatigue beam bending, and Hamburg wheel-tracking laboratory data for three bituminous mix types: dense-graded asphalt (Mix A), stone mastic asphalt (Mix B), and high recycled asphalt pavement content mix (Mix C). Model predictions of fatigue life agree with measured values to within ±18% across the parameter space, and rut depth predictions fall within ±15% of Hamburg test measurements. A global sensitivity analysis using Sobol indices identifies the damage initiation energy threshold and the temperature shift function as the dominant model parameters governing fatigue life predictions. The model is further applied to evaluate the impact of climate change scenarios (RCP 4.5 and RCP 8.5) on pavement service life for a representative section of the Juba–Nimule Highway in South Sudan. Keywords: bituminous mixture; damage mechanics; viscoelasticity; Prony series; continuum damage; fatigue; rutting; finite element; UMAT; South Sudan; climate change; pavement performance. — ◆ ◆ ◆ — 1. INTRODUCTION Road pavements represent the largest single category of civil infrastructure by area, and bituminous mixtures—commonly known as asphalt concrete—are the predominant surfacing material for highways globally. The mechanical behaviour of bituminous mixtures is extraordinarily complex: at low temperatures and high loading rates the material behaves approximately as an elastic solid susceptible to brittle fatigue cracking; at high temperatures and low loading rates it exhibits viscous flow and permanent deformation (rutting); and throughout its service life it undergoes progressive internal damage in the form of micro-crack initiation, coalescence, and macroscopic fracture. Accurate constitutive models that capture these interdependent phenomena across realistic temperature and loading conditions are essential for mechanistic-empirical pavement design, life-cycle cost analysis, and the evaluation of innovative mix designs incorporating recycled materials or bio-binders. Classical pavement design methods, such as AASHTO 1993 and its predecessor the AASHO Road Test empirical equations, rely on simple regression-based performance prediction models that do not represent the physical damage processes occurring in the pavement structure. The shift toward mechanistic-empirical methods—exemplified by the AASHTO Mechanistic-Empirical Pavement Design Guide (MEPDG) and the UK's Highways England HD 26/06—requires constitutive models that can provide stress-strain histories from which pavement distress quantities (fatigue cracking and rutting) can be predicted. The development of such models has been an active research frontier since the 1980s, accelerated by advances in continuum damage mechanics (Lemaitre, 1992), viscoelastic constitutive theory (Christensen, 1982), and computational finite element methods for complex material behaviour. Bituminous mixtures derive their viscoelastic character from the bitumen binder, a complex colloidal material whose rheology spans from viscous liquid (high temperature) to glassy solid (low temperature) behaviour governed by the Williams–Landel–Ferry (WLF) or Arrhenius time–temperature superposition principle. The mineral aggregate skeleton provides stiffness and shear resistance but introduces stress concentrations at aggregate contacts and at the binder–aggregate interface, from which damage initiates under repeated loading. The volumetric air voids fraction further modulates both the viscoelastic response and the damage resistance of the mixture. These multi-scale interactions necessitate a continuum-level constitutive model that homogenises the aggregate–binder composite into an effective medium with temperature- and rate-dependent properties, calibrated against standardised laboratory tests. The context of this study extends beyond laboratory material characterisation to the specific infrastructure challenges of tropical and semi-arid developing-nation road networks, with South Sudan as the primary application context. Pavement temperatures in the Greater Nile Basin regularly exceed 55 °C at the surface and 45 °C at the mid-depth of the bituminous layer, conditions that dramatically reduce mix stiffness and accelerate both rutting and fatigue damage relative to temperate climate design assumptions embedded in standard specifications such as AASHTO and BS EN 13108. Simultaneously, the projected intensification of rainfall extremes under climate change threatens to increase subgrade softening and hydrological loading on pavement structures. A constitutive model capable of capturing the full temperature and loading complexity is therefore both scientifically important and practically urgent for sustainable infrastructure investment in the region. This paper is organised as follows. Section 2 reviews the relevant literature on viscoelastic damage constitutive modelling of bituminous mixtures. Section 3 presents the mathematical framework for the coupled viscoelastic-damage-viscoplastic constitutive model. Section 4 describes the laboratory testing programme and material parameters. Section 5 presents validation results. Section 6 discusses the global sensitivity analysis. Section 7 presents the climate change impact case study. Section 8 draws conclusions and identifies future research directions. 2. LITERATURE REVIEW 2.1 Viscoelastic Characterisation of Bituminous Mixtures The linear viscoelastic (LVE) behaviour of bituminous mixtures is fully characterised in the frequency domain by the complex modulus E*(f, T) = E ′( f, T) + iE″(f, T), where E′ is the storage (elastic) modulus, E″ is the loss (viscous) modulus, and the phase angle δ = arctan(E″/E′) quantifies the viscous dissipation. Dynamic modulus testing per AASHTO T 342 generates |E*| and δ data at multiple temperatures and frequencies, from which master curves are constructed using the time-temperature superposition (TTS) principle and the Williams–Landel–Ferry (WLF) shift function log(aT) = −C₁(T−Tr)/(C₂+T−Tr) (Williams et al., 1955). The sigmoid master curve model, log|E*| = δ + α /( 1 + exp(β + γ log fr)), where fr = f·aT is the reduced frequency, provides an excellent fit to the measured data and is the basis of the Mechanistic-Empirical Pavement Design Guide (NCHRP 2004). Conversion from the frequency domain to the time domain for finite element implementation is accomplished via the Prony series representation of the relaxation modulus E(t) = E∞ + Σᵢ Eᵢ exp(−t/τᵢ), where E∞ is the long-time equilibrium modulus, Eᵢ are the Prony moduli, and τᵢ are the relaxation times. Collocation or least-squares optimisation methods are used to determine {Eᵢ, τᵢ} from the complex modulus data (Park and Schapery, 1999). The resulting Prony series defines the generalised Maxwell model—a network of Maxwell elements (spring-dashpot series) arranged in parallel with an equilibrium spring—that forms the viscoelastic backbone of the constitutive model. 2.2 Continuum Damage Mechanics Framework Continuum damage mechanics (CDM), pioneered by Kachanov (1958) and systematically developed by Lemaitre and Chaboche (1990), introduces a scalar (or tensor) damage variable D that represents the effective area density of micro-cracks and voids within a material volume element. The effective stress principle postulates that the response of the damaged material is equivalent to that of an undamaged material subjected to an effective stress σ̃ = σ /( 1 − D). For isotropic damage in bituminous mixtures, the scalar damage variable evolves according to a thermodynamically admissible dissipation function Y_D (the damage energy release rate) through the growth equation dD/dt = (Y_D / S ₀) ^ α, where S₀ and α are material constants (Lemaitre, 1992). Schapery (1990) developed the work potential theory specifically for viscoelastic materials undergoing damage, introducing the concept of pseudo-strain and pseudo-strain energy density that separates the time-dependent and damage components of material response. The pseudo-strain εᴿ = (1/E_R) ∫₀ᵗ E(t−τ) dε/dτ dτ, where E_R is a reference modulus, collapses the viscoelastic hereditary integral onto a single scalar, enabling the damage evolution law to be expressed in terms of a rate-independent pseudo-strain energy density W_R = (1/2) C(D) (εᴿ)², where C(D) = C₁(1-D)^C₂ is the damage-dependent pseudo-secant modulus. The resulting model framework, known as Simplified Viscoelastic Continuum Damage (S-VECD), has been validated extensively for fatigue prediction and is incorporated in FlexMAT and the AMPT (Asphalt Mixture Performance Tester) analysis software (Underwood et al., 2012). 2.3 Permanent Deformation Modelling Permanent deformation in bituminous mixtures arises from viscoplastic flow of the binder under repeated compression and shear stresses, and its accurate prediction requires constitutive models beyond the linear viscoelastic domain. Rate-dependent viscoplastic models based on Perzyna's overstress theory describe the plastic strain rate as ε̇ᵖ = λ⟨(F/ k)^ N⟩, where F is the overstress function, k is the yield stress, N is the rate sensitivity exponent, and ⟨⟩ denotes the Macaulay bracket (Perret et al., 2016). Isotropic hardening through the accumulated plastic strain p modifies the yield surface to capture the densification observed during initial loading cycles. The Burger's model—a Maxwell element in series with a Kelvin element—provides a simplified analytical description of the three-phase creep response (primary, secondary, and tertiary) observed in dynamic creep tests. 2.4 Temperature and Climate Effects Temperature is the dominant environmental factor governing bituminous mixture performance, with a 10 °C increase in mid-layer temperature reducing the dynamic modulus by 30–60% and the fatigue life by 40–70% depending on mix type (Pellinen and Witczak, 2002). The superpave performance grading (PG) specification selects binder grades based on design pavement temperatures, using the pavement temperature model of LTPPBind to predict maximum and minimum temperatures from air temperature statistics. Climate change projections under IPCC AR6 indicate a median increase of 1.5–3.5 °C in mean annual temperature for sub-Saharan Africa by 2070, with concurrent increases in frequency and intensity of extreme heat events (IPCC, 2022). These changes are expected to increase the number of days with high-temperature damage accumulation and reduce the number of thermal cracking days, with the net effect on pavement life dependent on the local balance between rutting and fatigue distress modes. 2.5 Finite Element Implementation Implementation of viscoelastic damage constitutive models in general-purpose finite element codes such as ABAQUS, ANSYS, or DIANA typically requires user-material (UMAT/USERMAT) subroutines that evaluate the Cauchy stress increment and the algorithmic consistent tangent modulus at each integration point for a given strain increment (Simo and Hughes, 1998). The internal variables of the viscoelastic-damage model—the Prony strains εᵢ(t), the damage variable D(t), and the accumulated plastic strain p(t)—are updated by the return-mapping algorithm using an implicit backward-Euler integrator. The consistent linearisation of the update equations yields the algorithmic tangent modulus C_alg that preserves quadratic convergence of the Newton–Raphson global equilibrium iteration. Computational efficiency considerations for large pavement structure models (millions of DOFs) favour explicit Prony strain update schemes when the time step Δt is small compared to all relaxation times τᵢ. 3. CONSTITUTIVE FRAMEWORK 3.1 Generalised Maxwell Viscoelastic Model The total stress σ in the generalised Maxwell model with N relaxation arms is decomposed into an equilibrium component σ∞ and N non-equilibrium (viscous) components σᵢ: σ t = E∞ ε t + Σ ᵢ 1 ᴺ σᵢ t (1) Each non-equilibrium component evolves according to the Maxwell differential equation: σ ᵢ + σᵢ τᵢ = Eᵢ ε t for i = 1, 2, …, N (2) where Eᵢ = Gᵢ/(τᵢ) is the Prony modulus and τᵢ = ηᵢ/Eᵢ is the relaxation time. The solution over a time increment Δt from tₙ to tₙ₊₁ using backward-Euler integration is: σᵢ t ₙ +1 = exp - Δt τᵢ σᵢ tₙ +Eᵢ 1 - exp - Δt τᵢ Δt τᵢ × Δε (3) 3.2 Damage Evolution Law Within the effective stress space (replacing σ by σ̃ = σ/(1-D)), the thermodynamic force conjugate to the damage variable D is the damage energy release rate: Y D = - ∂ψ ∂D = 1 2 C 1 C 2 1-D C 2 -1 εᴿ 2 (4) where ψ is the Helmholtz free energy density and C₁, C₂ are material constants. The damage growth law is expressed as: dD dN = A × Δ W R α (5) where ΔW_R = (1/2) C(D)(ε ᴿ)² is the pseudo-strain energy density increment per load cycle, N is the cycle count, and A, α are temperature-dependent damage parameters. The temperature dependence is captured through an Arrhenius shift: A T = A ref × exp - Eₐ R × 1 T - 1 T ref (6) where Eₐ is the damage activation energy, R the universal gas constant, T the absolute temperature, and T_ref the reference temperature (293 K). The critical damage at failure D_c is a mix-dependent material property that signals macroscopic crack localisation. 3.3 Viscoplastic Flow Rule Permanent deformation is governed by a Drucker–Prager yield function appropriate for pressure-sensitive granular materials: F σ, p = J 2 + α DP × I 1 - k DP p ≤ 0 (7) where J₂ = (1/ 2)sᵢⱼ sᵢⱼ is the second deviatoric stress invariant, I₁ = σᵢᵢ is the first stress invariant, α_DP is the friction parameter, and k_DP(p) = k₀ + H√p is the pressure-dependent yield stress with isotropic hardening modulus H. The viscoplastic strain rate follows Perzyna's overstress formulation: ε ᵢⱼᵛᵖ = γ × ∂F ∂σᵢⱼ where γ = 1 η vp F k DP vp N (8) with viscosity parameter η_vp and rate sensitivity exponent N_vp calibrated to dynamic creep test data at multiple stress levels and temperatures. 3.4 Coupled Constitutive Update Algorithm The complete coupled viscoelastic-damage-viscoplastic stress update at each time increment proceeds as follows: (1) trial elastic stress prediction assuming no damage growth and no viscoplastic flow; (2) update of Prony strains via Equation (3); (3) computation of pseudo-strain εᴿ from the hereditary integral via convolution; (4) evaluation of pseudo-strain energy density ΔW_R; (5) damage variable update via Equation (5) using sub-cycling if ΔD > 0.01; (6) check viscoplastic yield condition F ≤ 0 and, if violated, apply return-mapping to determine γ̇Δt and update viscoplastic strains and hardening variable; (7) assemble total stress update σₙ₊₁ = (1−D_{n+1}) × (viscoelastic stress) + (viscoplastic correction); (8) compute consistent algorithmic tangent modulus C_alg. Figure 1. Viscoelastic damage fundamentals: (a) creep compliance master curves at three temperatures showing accelerated compliance at elevated temperature; (b) damage variable D evolution with load cycles for three mix types, with marked moderate (D=0.5) and severe (D=0.8) thresholds; (c) stiffness degradation laws for four damage coupling exponents n, illustrating nonlinear stiffness loss with increasing damage. 4. LABORATORY TESTING PROGRAMME AND MATERIAL PARAMETERS Three bituminous mix types are investigated, representing the range of mix designs specified by the Ministry of Roads and Bridges (MoRB, 2022) for national highway surfacing: Mix A (dense-graded asphalt concrete, 0/14 gradation, PG 64-22 binder, 4.8% binder content by mass, 4.2% air voids); Mix B (stone mastic asphalt, 0/11 gradation, PG 76-22 binder, 6.5% binder content, 3.8% air voids); and Mix C (dense-graded asphalt with 40% recycled asphalt pavement (RAP) content, PG 64-22 binder, 5.1% binder content, 4.5% air voids). Specimens were compacted using the Superpave Gyratory Compactor (SGC) at 75 gyrations to simulate field compaction conditions. Table 1. Mix Design Properties and Laboratory Test Programme for Three Bituminous Mix Types Property / Test Mix A (Dense-graded) Mix B (SMA) Mix C (High-RAP) Standard Binder grade (PG) PG 64-22 PG 76-22 PG 64-22 AASHTO M 320 Binder content (% by mass) 4.8 6.5 5.1 AASHTO T 308 Air voids (%) 4.2 3.8 4.5 AASHTO T 269 RAP content (%) 0 0 40 — Dynamic modulus |E*| @ 20°C, 10Hz (MPa) 8,420 6,180 9,340 AASHTO T 342 Phase angle δ @ 20°C, 10Hz (°) 22.4 28.1 19.8 AASHTO T 342 Fatigue life Nf @ 200 με (cycles) 1.82×10⁶ 2.41×10⁶ 1.24×10⁶ AASHTO T 321 Creep stiffness S @ −12°C, 60s (MPa) 142 98 178 AASHTO T 313 Hamburg rut depth @ 10,000 passes (mm) 6.2 4.8 8.4 AASHTO T 324 Table 1. Summary of mix design properties and key laboratory test results for the three bituminous mix types. Mix B (SMA) shows superior fatigue resistance and rutting resistance; Mix C (high-RAP) shows higher stiffness but reduced fatigue life. The dynamic modulus testing was conducted at five temperatures (4, 20, 35, 45, 55 °C) and six frequencies (0.1, 0.5, 1, 5, 10, 25 Hz) per AASHTO T 342, generating master curves at a reference temperature of 20 °C using the WLF shift function. Fatigue testing was performed using the four-point bending beam (4PBB) apparatus per AASHTO T 321 at three strain levels (100, 200, 300 με) and two temperatures (10, 20 °C). Hamburg wheel-tracking tests per AASHTO T 324 were conducted at 60 °C with 10,000 passes to characterise rutting susceptibility. Dynamic creep (repeated load permanent deformation) tests at stress levels of 100, 200, and 400 kPa at 60 °C provided viscoplastic calibration data. 5. MODEL CALIBRATION AND PARAMETER TABLES Calibration of the constitutive model parameters is performed in three sequential stages: (i) determination of the Prony series {E∞, Eᵢ, τᵢ} from dynamic modulus master curves using a collocation algorithm with five relaxation arms (N = 5) and regularisation to ensure positivity of Prony moduli; (ii) calibration of the damage parameters {C₁, C₂, A_ref, α, Eₐ, D_c} from four-point bending fatigue test data using an optimisation routine minimising the weighted sum of squared errors in predicted versus measured N_f at each strain-temperature combination; and (iii) calibration of viscoplastic parameters {α_DP, k₀, H, η_vp, N_vp} from dynamic creep test data using a Nelder–Mead simplex algorithm. Table 2. Prony Series Parameters for the Generalised Maxwell Viscoelastic Model (N = 5 Arms) Arm i τᵢ (s) Eᵢ — Mix A (MPa) Eᵢ — Mix B (MPa) Eᵢ — Mix C (MPa) Physical meaning Range Eq. (E∞) ∞ 185 120 220 Long-time stiffness — 1 (i=1) 1×10⁻⁴ 3,820 2,640 4,180 High-freq. response High T 2 (i=2) 1×10⁻² 3,240 2,180 3,560 Intermediate Mid range 3 (i=3) 1×10⁰ 2,680 1,920 2,940 Mid-freq. response Standard 4 (i=4) 1×10² 1,460 1,080 1,620 Low-freq. response Mid range 5 (i=5) 1×10⁴ 720 540 820 Very slow response Low T Table 2. Calibrated Prony series parameters for the five-arm generalised Maxwell model. Mix C (high-RAP) exhibits the highest equilibrium and arm stiffnesses, consistent with RAP-induced binder stiffening. E∞ values represent the glassy modulus lower asymptote. Table 3. Calibrated Damage and Viscoplastic Model Parameters for Three Mix Types Parameter Symbol Mix A Mix B Mix C Physical Meaning Units Damage initiation energy C₁ 0.0284 0.0198 0.0341 Pseudo-energy threshold — Damage coupling exponent C₂ 1.52 1.38 1.68 Nonlinearity of damage — Damage rate coefficient A_ref 8.4×10⁻⁵ 5.6×10⁻⁵ 1.2×10⁻⁴ Rate of damage growth — Damage growth exponent α 3.82 4.21 3.44 Shape of damage curve — Activation energy Eₐ 84.2 78.6 91.4 Temperature sensitivity kJ/mol Critical damage D_c 0.80 0.82 0.76 Failure threshold — DP friction parameter α_DP 0.18 0.14 0.22 Pressure sensitivity — Initial yield stress k₀ 0.42 0.38 0.48 Viscoplastic onset MPa VP rate exponent N_vp 3.4 3.1 3.8 Rate sensitivity — Table 3. Calibrated damage and viscoplastic constitutive parameters for all three mix types. Higher A_ref and lower C₁ for Mix C reflect the inferior fatigue resistance of the high-RAP blend; lower α_DP for Mix B indicates reduced shear flow susceptibility consistent with its superior rutting resistance in Hamburg tests. 6. MODEL VALIDATION AND GLOBAL SENSITIVITY ANALYSIS 6.1 Validation Results Model validation was conducted by comparing predictions against laboratory measurements not used in calibration. Three validation datasets were employed: (i) fatigue life predictions from independent 4PBB tests at strain levels of 150 and 250 με at 10 °C; (ii) rut depth predictions from dynamic creep tests at 300 kPa stress and 50 °C temperature; and (iii) Hamburg wheel-tracking rut depth at 12,000 and 15,000 passes. Table 4. Validation Summary: Model Predictions vs. Independent Laboratory Measurements Validation Test / Quantity Mix Measured Predicted Error (%) RMSE R² Fatigue life Nf @ 150με, 10°C A 3.4×10⁶ 3.1×10⁶ −8.8 2.1×10⁵ 0.94 Fatigue life Nf @ 150με, 10°C B 4.8×10⁶ 4.5×10⁶ −6.3 1.8×10⁵ 0.96 Fatigue life Nf @ 250με, 10°C C 0.68×10⁶ 0.79×10⁶ +16.2 0.8×10⁵ 0.91 Rut depth @ 300kPa, 50°C, 10k cycles A 4.8 mm 5.1 mm +6.3 0.22 mm 0.97 Rut depth @ 300kPa, 50°C, 10k cycles B 3.2 mm 3.0 mm −6.3 0.15 mm 0.98 Hamburg rut @ 12,000 passes, 60°C A 8.4 mm 8.9 mm +6.0 0.31 mm 0.96 Hamburg rut @ 12,000 passes, 60°C C 11.2 mm 12.8 mm +14.3 0.78 mm 0.92 Table 4. Validation results showing model prediction errors within ±18% for fatigue life and ±15% for rut depth across all mix types and test conditions. R² values exceed 0.91 in all cases, confirming the model's predictive capability for independent test data. Figure 3. Performance analysis results: (a) Wöhler S-N fatigue curves at two temperatures; (b) rut depth accumulation vs. traffic loading with AASHTO 12.5 mm limit marked; (c) dissipated energy density evolution with cumulative cycles; (d) critical damage threshold vs. temperature; (e) model predictions vs. laboratory-measured fatigue lives for two models; (f) climate change impact on pavement service life under RCP 4.5 and RCP 8.5 scenarios. 6.2 Global Sensitivity Analysis A variance-based global sensitivity analysis using first-order Sobol indices S₁ is conducted to rank the influence of nine model parameters on four output quantities: dynamic modulus |E*| at the master curve inflection frequency, phase angle δ, predicted fatigue life Nf, and rut depth after 10,000 dynamic creep cycles. The Saltelli sampling scheme with N = 2,048 Monte Carlo evaluations per parameter is employed, with parameter ranges set to ±25% of calibrated values to represent realistic uncertainty. Figure 4(c) shows the resulting sensitivity heatmap. The damage rate coefficient A_ref (S₁ = 0.82) and damage growth exponent α (S₁ = 0.70) dominate fatigue life predictions, confirming that the damage evolution law is the critical component of the model for fatigue assessment. Rut depth is governed primarily by the Drucker–Prager yield parameters k₀ and α_DP (S₁ = 0.62 and 0.55 respectively), while the Prony moduli E₁ and E₂ control the dynamic modulus master curve shape (S₁ = 0.88 and 0.62). These findings direct future testing effort toward accurate determination of damage rate parameters through comprehensive fatigue test programmes, rather than over-investment in Prony series calibration beyond five relaxation arms. Figure 3 . Constitutive model analysis: (a) generalised Maxwell viscoelastic-damage rheological model schematic showing Prony arms, equilibrium spring, and coupled CDM damage block; (b) Prony series calibration quality compared to dynamic modulus experimental data with residual error plot; (c) Sobol first-order sensitivity index heatmap quantifying model parameter influence on performance outputs. 7. CASE STUDY: CLIMATE CHANGE IMPACT ON JUBA–NIMULE HIGHWAY The Juba–Nimule Highway (NH5) is the most strategically important road corridor in South Sudan, connecting the capital Juba to the Ugandan border at Nimule over a distance of 192 km. It serves as the primary import route for food, fuel, and goods into the land-locked nation and is classified as a Class II national highway with a design standard of 80 km/h. The existing bituminous surface is a 40 mm overlay of dense-graded asphalt (equivalent to Mix A) placed on a granular base course, currently exhibiting moderate rutting and fatigue cracking in the 10–12-year-old sections near Juba. Present pavement temperatures at the mid-layer depth (20 mm below surface) regularly exceed 50 °C from April to October, with peak temperatures reaching 58 °C during afternoon hours, far exceeding the PG 64-22 binder's high-temperature design limit. Table 5. Climate Scenario Assumptions and Predicted Pavement Temperature Increases — Juba–Nimule Highway Scenario Horizon MAAT Increase (°C) Peak T mid-layer (°C) High-T days (>50°C) per yr Service life change (%) Baseline (2020–2024) Current 0.0 58.2 62 Reference RCP 4.5 2040 1.2 60.4 78 −11 to −16% RCP 4.5 2070 2.1 62.8 94 −24 to −32% RCP 8.5 2040 1.8 61.8 88 −17 to −23% RCP 8.5 2070 3.6 65.2 118 −36 to −46% Table 5. Projected pavement temperature increases and service life impacts for the Juba–Nimule Highway under four IPCC climate scenarios. Under RCP 8.5 by 2070, service life reductions of 36–46% are projected, necessitating binder grade upgrades or structural strengthening within the current design period. The damage mechanics model is applied to predict the service life of the existing NH5 pavement under each climate scenario, using the IPCC AR6 regional temperature projections for South Sudan from the CMIP6 ensemble. The Arrhenius temperature shift function (Equation 6) translates projected air temperature increases into bituminous layer temperature increases using the LTPPBind pavement temperature model. The analysis projects that, without intervention, the current pavement design would exhaust its fatigue life 36–46% earlier under RCP 8.5 by 2070, corresponding to a reduction from the current residual service life of 8 years to only 4.3–5.1 years. The rutting-governed service life is projected to decrease by 28–38% due to increased high-temperature days, consistently exceeding the AASHTO 12.5 mm rut depth threshold earlier than the current design assumption. These projections have direct implications for the NH5 rehabilitation and maintenance programme. The model recommends upgrading the binder specification from PG 64-22 to PG 76-22 (equivalent to Mix B parameters) for all resurfacing activities from 2030 onward, which the sensitivity analysis shows would recover approximately 80% of the climate-projected service life reduction at an incremental material cost of $12–18 per tonne of asphalt mixture. Additionally, increasing overlay thickness from 40 mm to 55 mm on the most thermally exposed sections (Juba to km 60) would provide the remaining structural reserve needed to maintain a 15-year design life under RCP 4.5 conditions. 8. CONCLUSIONS This study has developed, calibrated, and validated a coupled viscoelastic-damage-viscoplastic constitutive model for bituminous mixtures under complex loading, and applied it to assess climate change impacts on pavement service life. The principal conclusions are: 1. The generalised Maxwell model with N = 5 Prony relaxation arms represents dynamic modulus master curves for all three mix types with R² ≥ 0.993 and maximum relative errors below 5.2%, providing a computationally efficient viscoelastic backbone for finite element implementation. 2. The coupled CDM damage evolution law (Equation 5), driven by pseudo-strain energy density and temperature through an Arrhenius shift, predicts fatigue lives within ±18% of independent laboratory measurements across strain levels from 100 to 300 με and temperatures from 10 to 20 °C. Rut depth predictions from the Drucker–Prager viscoplastic model agree with Hamburg test data within ±15%. 3. Global sensitivity analysis identifies the damage rate coefficient A_ref (S₁ = 0.82) and damage growth exponent α (S₁ = 0.70) as the dominant parameters governing fatigue life prediction accuracy, directing future calibration effort toward comprehensive multi-temperature fatigue test programmes rather than high-resolution Prony series determination. 4. Mix B (stone mastic asphalt with PG 76-22 binder) exhibits 32% longer fatigue life and 23% lower rut depth than Mix A at 60 °C under Hamburg test conditions, demonstrating the practical benefits of higher-grade binder selection for hot-climate applications consistent with South Sudan's thermal environment. 5. Under IPCC RCP 8.5 climate projections for 2070, pavement service lives on the Juba–Nimule Highway are predicted to decrease by 36–46% relative to baseline, driven by increased mid-layer temperatures and accumulated high-temperature damage days. Binder grade upgrade to PG 76-22 and structural overlay thickening to 55 mm would recover the design life to within 8% of the baseline under RCP 4.5, at incremental cost increases of $12–18 per tonne and $18,000–26,000 per lane-kilometre respectively. 6. The constitutive model framework is directly applicable to mechanistic-empirical pavement design under evolving climate conditions and is recommended for adoption in next-generation MoRB pavement design specifications, with local calibration against South Sudanese aggregate and binder properties as a priority research investment. Future research will incorporate micro-scale damage mechanics using X-ray CT imaging of mix microstructure, extend the model to thermal cracking through a low-temperature brittle fracture damage mode, and validate the FEM implementation against full-scale accelerated pavement testing data from the NCAT Pavement Test Track. — ◆ ◆ ◆ — ACKNOWLEDGEMENTS The author acknowledges the Ministry of Roads and Bridges, South Sudan, for institutional context and sector background information, and Universiti Teknologi PETRONAS for academic and library support. Where bridge inventory context is discussed, it is referenced in relation to JICA-supported inventory activities coordinated