An Energy-Data-Anchored Probabilistic Digital-Twin Framework for Coupled Structural–Envelope Degradation and Local Energy-Demand Amplification in Low-Carbon Buildings

1. Introduction

The operating carbon of a building is not fixed at commissioning. Envelope joints loosen, seals fail, insulation becomes discontinuous, HVAC efficiency changes and repeated thermal or moisture exposure alters the condition of structural and non-structural components. Yet many lifecycle models hold envelope properties constant while structural assessment proceeds in a separate model. That separation is convenient, but it removes a physically plausible route by which deterioration increases thermal demand and prolonged system operation modifies the exposure that drives further deterioration.

Digital twins offer a useful architecture for closing that gap because they combine a state model, observations, inference and decisions rather than treating simulation as a one-off exercise (Madni et al. 2019; Rasheed et al. 2020; Wright and Davidson 2020). Construction and building research has moved from static building information models toward connected representations that can assimilate measured performance and support predictive action (Boje et al. 2020; Deng et al. 2021). Probabilistic graphical formulations further show how uncertainty, data assimilation and control can be placed within a coherent digital-twin structure (Kapteyn et al. 2021).

Recent building digital-twin research has moved beyond static asset representation toward hybrid simulation–data models and continuously updated building-energy management systems. Hybrid digital twins have used operational data to correct physics-based building simulations, while digital-twin-enabled energy management systems have demonstrated the value of integrated information fusion for operational prediction and control. In parallel, structural digital-twin research has advanced Bayesian model updating through learnable mappings and multi-fidelity state estimation. These studies strengthen domain-specific prediction and calibration, but they do not provide a long-horizon criterion for detecting local modal amplification within a coupled structural–envelope–energy deterioration process (Oh et al. 2024; Hwang et al. 2025; Torzoni et al. 2025). Recent studies have also begun to quantify how ageing and degradation of envelope components alter operational energy demand. Reported pathways include deterioration of insulation, insulated glazing, airtightness, solar reflectivity and other thermal properties that are commonly treated as fixed at commissioning (Taki and Zakharanka 2023). Recent structural-health-monitoring digital-twin research likewise emphasises that field deployment requires explicit treatment of sensor limitations, model updating, uncertainty and model-form discrepancy rather than reliance on simulation alone (Sun et al. 2025; Andrés Arcones et al. 2024). This literature supports a coupled lifecycle perspective, but it also demonstrates why an electricity record by itself cannot identify long-term structural or envelope deterioration laws.

Two limitations remain. First, energy calibration alone cannot establish the evolution of a coupled structural–envelope system. Second, a rising damage index or annual energy use does not reveal whether their interaction has crossed from slow deterioration into a self-amplifying regime. The latter requires a dynamical stability measure. The present paper therefore treats instability as a local property of the updated coupled state, not as a visual label attached after a threshold has been exceeded. The mathematical quantity underlying CSEI—the spectral abscissa of a local Jacobian—is established in dynamical-systems analysis. The contribution claimed here is therefore not a new stability theorem. It is the integration of this deliberately limited diagnostic into an energy-anchored, three-state building lifecycle framework, together with explicit coupling isolation, state recovery, uncertainty propagation and finite-time diagnostics that define what the local index can and cannot establish.

Short-horizon stability preservation has previously been examined as an operational-control problem (Samaei and Riffat 2026). That study addressed supervisory action selection over a short operational horizon. The present study does not reuse that control formulation. It introduces different lifecycle states, a 25-year monthly degradation horizon, an empirical electricity-observation anchor, parameter–climate uncertainty propagation, controlled state-recovery experiments and a frozen-time modal diagnostic. No rolling-horizon optimisation, intervention policy, renewable-storage design or self-healing control is evaluated in the present paper.

The investigation addresses four research questions and associated expectations:

RQ1/H1: How do otherwise identical energy-decoupled, one-way and bidirectional formulations differ in lifecycle electricity, carbon and deterioration? It is expected that retaining both feedback directions will increase these outputs relative to the energy-decoupled formulation.

RQ2/H2: Does the spectral abscissa of the frozen local Jacobian change sign within the 25-year reference horizon, and is its implementation consistent with frozen-eigenmode and finite-time diagnostics? It is expected that the reference trajectory will exhibit a sign change before year 25, while finite-time analysis will demonstrate the limits of interpreting that sign as a global-stability condition.

RQ3/H3: Does sequential state assimilation and regularised parameter recovery reduce state-tracking error and retrospective onset-reconstruction error relative to a biased prior-only trajectory? It is expected that posterior RMSE and absolute onset-reconstruction error will decrease and that posterior state intervals will achieve approximately nominal coverage in correctly specified synthetic recovery experiments.

RQ4/H4: Which of the nine uncertain state-model parameters govern cumulative instability exposure and restricted time to onset? It is expected that nonlinear deterioration acceleration and the bidirectional operational-feedback pathway will dominate the total-effect indices.

The contribution is an auditable chain from measured hourly electricity to a reduced-order observation model, a three-state coupled lifecycle formulation, sequential state updating, regularised parameter recovery, 5,000-realisation parameter–climate uncertainty propagation, local modal and finite-time diagnostics, and global sensitivity analysis of nine uncertain state-model parameters (Figure 1). Empirical validation is confined to the electricity-observation layer. Coupled-state behaviour is supported through analytical consistency, limiting cases, time-step convergence, correctly specified synthetic recovery experiments and probabilistic numerical verification. All reported numerical analyses use the same locked reference parameter set, state definitions, monthly load profile, uncertainty distributions and disclosed random seeds.

2. Framework and validation design

2.1. Validation hierarchy

Validation and verification are treated as distinct layers (Table 1). First, the electricity observation model is evaluated on an unseen calendar year, providing empirical out-of-sample evidence for the operational data layer. Second, the coupled state equations are verified through energy-decoupled and one-way limiting cases and through time-step refinement. Third, the CSEI implementation is checked using frozen-eigenmode propagation and is complemented by finite-time singular-value analysis to expose non-normal transient amplification. Fourth, state tracking and parameter recovery are evaluated using quarterly synthetic observations generated from a reference trajectory and repeated over 100 independent recovery experiments. Fifth, uncertainty estimates are examined using 5,000 parameter–climate realisations, Monte Carlo convergence and sampling-error diagnostics. These layers do not constitute field validation of long-term structural deterioration; they establish empirical validity for the energy layer and numerical, statistical and inference consistency for the coupled framework.

2.2. Empirical energy anchor

The operational anchor uses the Building Data Genome Project 2, an open collection of hourly meter, building-metadata and weather records (Miller et al. 2020). The office meter Hog_office_Napoleon was selected because it provides an aligned 17,544-hour record over 2016–2017 and a reported floor area of 12,365.4 m². Annual electricity use was 1,140,900.84 kWh in 2016 and 1,123,697.12 kWh in 2017, corresponding to 92.27 and 90.87 kWh m^-2 year^-1. Their unrounded two-year mean, 91.5699 kWh m^-2 year^-1, anchors the baseline whole-building electricity of the reference scenario.

A ridge-regularised autoregressive grey-box observation model was fitted to the 2016 record. In linear form, the model was

where $z_{r,t}$ contained dry-bulb temperature, squared temperature, dew point, sea-level pressure, wind speed, one sine–cosine pair for the hourly cycle, one pair for the weekly cycle, one pair for the annual cycle, and a weekend indicator. The fitted design therefore contained 17 predictors excluding the intercept: five weather variables, six Fourier terms, one calendar indicator and five autoregressive lags.

Ridge regression used an intercept and the locked case-specific regularisation coefficient α=10. Predictors were retained in their native units; accordingly, α is reported only as an implementation constant for this observation model and is not claimed to be a transferable hyperparameter. The 2017 outcomes were not used for coefficient fitting or regularisation selection.

Missing weather values were linearly interpolated and then forward- and backward-filled. Missing meter values were linearly interpolated and constrained to non-negative values. Meter and weather timestamps were aligned on their common hourly index. Construction of the 168-hour maximum lag removed the first 168 hours of the 2016 record, leaving 8,616 lag-complete training observations. The full 8,760-hour 2017 calendar year was retained for testing because its initial lagged values were available from the preceding record. The test was strictly one step ahead and used observed lagged electricity; it was not a recursively free-running annual forecast (Figure 2).

The performance measures were defined as $RMSE={[n^{-1}{\mathrm{\Sigma }}_t{({\hat{y}}_t-y_t)}^2]}^{1/2}$, $CV(RMSE)=100\times RMSE/\bar{y}$ and $NMBE=100\times {\mathrm{\Sigma }}_t({\hat{y}}_t-y_t)/(n\bar{y})$. Figure 3 uses 1,800 test observations sampled without replacement using seed 20260728. After the independent test metrics had been fixed, the 2016 and 2017 observations were pooled only to calculate the descriptive monthly whole-building electricity shares used in the lifecycle scenario. The 2017 outcomes therefore affected the lifecycle load profile but not the fitted observation-model coefficients or the reported out-of-sample metrics.

The 2017 one-step-ahead CV(RMSE) and NMBE values are below the hourly thresholds commonly used as descriptive calibration references in ASHRAE Guideline 14-2023 (Table 2). This is not presented as formal Guideline 14 compliance because the autoregressive evaluation uses observed lagged electricity and is not a recursively free-running calibrated building simulation. The near-identical training and independent-year R² values nevertheless show that the operational anchor is not supported solely by an in-sample fit. The meter does not separate HVAC, lighting and plug loads; the lifecycle model therefore applies a fixed scenario HVAC share of 0.58 to the whole-building baseline rather than treating the meter as a direct HVAC measurement.

3. Coupled lifecycle model

3.1. State definition

The reduced-order state is $x(t)={[D(t),L(t),H(t)]}^T$. D is a dimensionless structural/support-condition loss coordinate, whereas L is a separate dimensionless envelope thermal and leakage-loss coordinate. D does not directly represent insulation or airtightness degradation; these effects are represented by L. The coefficient c_D maps structural or support-condition deterioration into envelope loss. H is a dimensionless thermal-load memory state referenced to the initial operating condition and represents sustained HVAC runtime and associated thermal–moisture exposure. Electricity consumption itself is not interpreted as mechanically damaging the structure. All states are evaluated only within the prespecified 25-year analysis domain.

Here, t is measured in years, and C_j is the month-specific climate multiplier defined in Section 3.2. The nominal relative coefficient of performance was fixed at COP_rel=1.0 in all reported analyses. The deterministic reference and limiting-case trajectories used ε_j=0, whereas the stationary AR(1) disturbance was activated only in the parameter–climate uncertainty analysis. The H-to-D and H-to-L terms represent hypothesised exposure-mediated feedback pathways, reflecting the possibility that sustained HVAC operation may coincide with prolonged thermal and moisture exposure. They do not imply that electricity consumption itself causes structural damage. These feedback coefficients are scenario parameters and were not empirically identified from the electricity dataset.

The smooth activation function $s(z)=\ln [1+\exp (8z)]/8$ is a regularised approximation to the positive-part operator. It remains slightly positive within a narrow transition region around H=1, with $s(0)=ln(2)/8=0.08664$. It should therefore be interpreted as a smooth transition-band exposure term rather than an exact zero-to-positive switch. Consequently, the selected formulation introduces a small baseline feedback contribution at H=1, and the reported numerical results are conditional on this activation definition. Equation (2) does not guarantee the global invariance of D≤1 because the nonlinear acceleration term q_D D³ is not multiplied by 1−D. The results are therefore interpreted only within the prescribed 25-year lifecycle state domain. Protective numerical bounds of $0\le D,L\le 0.99$ and 0.05≤H≤5.0 were imposed solely to prevent solver failure outside the intended domain; no reference or uncertainty trajectory reached these bounds.

3.2. Electricity and carbon mapping

The reference scenario represents a 12,000 m², ten-storey office building. Its initial annual whole-building electricity is E₀=1,098,838.8 kWh year^-1, obtained from the locked mean EUI of 91.5699 kWh m^-2 year^-1. The twelve empirical monthly shares reported in Table A1 are repeated in each lifecycle year. The HVAC share f_HVAC=0.58, initial grid factor 0.42 kgCO₂e kWh^-1 and annual grid-factor decline of 1.2% are explicit scenario assumptions rather than estimates for a named electricity system.

Let $j\in \{0,\dots ,299\}$ denote lifecycle month and let q(j)=j mod 12 denote the calendar-month index, with q=0 corresponding to January. The monthly climate multiplier is

Baseline and degradation-adjusted whole-building electricity in month j are

The grid emission factor is

Lifecycle electricity and operational carbon are

Here H_j is the updated monthly state value used in the locked implementation. The factor 10^-6 converts kWh to GWh in Equation (9) and kgCO₂e to ktCO₂e in Equation (10). The disturbance ε_j evolves continuously across December–January boundaries and is not reset at the beginning of each lifecycle year.

3.3. Coupled Structural–Energy Instability Index

At each lifecycle month, the local Jacobian J=∂f/∂x is evaluated at the updated state while the exogenous climate term is fixed at its annualised reference value C_J=1.0. The Coupled Structural–Energy Instability Index is the spectral abscissa of this frozen Jacobian:

CSEI<0 indicates asymptotic decay of the dominant frozen eigenmode, CSEI=0 denotes a frozen-modal transition boundary, and CSEI>0 indicates local frozen-modal amplification. Because the lifecycle system is nonlinear, non-autonomous and potentially non-normal, the sign of CSEI is not a sufficient condition for global nonlinear stability, absence of finite-time transient amplification, structural collapse, serviceability failure or code compliance. CSEI is therefore interpreted only as a local frozen-time modal-transition diagnostic. In practical terms, CSEI answers a local modal question—whether the locally linearised coupled state tends to decay or amplify at a given month—whereas the finite-time diagnostic examines whether perturbations can still grow over a finite interval even when the local modal indicator remains negative.

The nominal onset time was defined on the locked monthly evaluation grid as

For trajectories that did not cross by year 25, T× was treated as right-censored and assigned infinity for the purpose of defining the restricted onset output

Cumulative instability exposure was defined as

which is dimensionless because CSEI has units of year^-1. In the locked monthly implementation, the nominal CSEI first became non-negative at month 211, corresponding to 17.5833 years and reported as 17.58 years. This is a monthly-grid onset estimate rather than a continuously interpolated crossing time.

To examine the limitation of the frozen-time interpretation, the variational system

was integrated along the non-autonomous reference trajectory. Finite-time amplification was quantified as

The equivalent logarithmic gain was calculated as

where Δ is expressed in years.

Finite-time gains were evaluated for selected 3-, 6- and 12-month windows at representative pre-crossing and post-crossing start years. For the reported CSEI, the climate term in the frozen Jacobian was fixed at its annualised reference value C_J=1.0. In contrast, the finite-time variational analysis retained the actual month-dependent climate multiplier C_j along the non-autonomous reference trajectory. Singular values were calculated using the Euclidean norm of the native normalised state coordinates ${[D,L,H]}^T$. Because finite-time singular-value gains depend on coordinate scaling, the reported values are diagnostics of transient amplification within the adopted state representation rather than invariant physical-energy gains. The annualised climate value was used in the frozen-time CSEI to isolate the slow lifecycle feedback from purely seasonal sign oscillations. Month-dependent climate forcing remained active in state propagation, uncertainty analysis and the finite-time variational diagnostic.

The values in Table 3 constitute the locked reference parameterisation used for the controlled numerical and probabilistic analyses. They were not estimated from the Building Data Genome Project 2 electricity meter. The empirical evidence in the present study directly supports only the electricity-observation layer. The lifecycle percentages, state trajectories and CSEI onset are therefore conditional reference-scenario outputs rather than empirical estimates for a general building stock. Application to a specific asset requires a temporally aligned, multi-source calibration dataset combining structural-response features (e.g., modal frequencies, strain or inspection-derived condition indices), envelope-condition measurements (e.g., airtightness, moisture, thermography or inferred heat-loss coefficients), and operational/BMS records (e.g., HVAC power, runtime and indoor–outdoor conditions). The observation mappings for D, L and H should first be calibrated, followed by regularised estimation of the deterioration and coupling coefficients and validation on a withheld monitoring period. Parameter dependence and model discrepancy should then be quantified from the asset-specific dataset before probabilistic forecasts are used for decision support.

The term energy-decoupled refers to separation of the structural–envelope subsystem from the thermal-load memory state, not to removal of all coupling between D and L (Table 4). In the energy-decoupled configuration, u₁=u₂=c_H=c_LH=0 while c_D remains active. In the one-way configuration, c_D, u₁ and u₂ are active while c_H=c_LH=0. The bidirectional configuration retains all pathways. Initial states, climate inputs, time integration and all other parameters are identical across the three configurations.

4. Bayesian updating and uncertainty

4.1. Observation and updating model

Quarterly observations combine a normalised structural-frequency proxy $y_f=\sqrt{1-D}+{\epsilon }_f$, an envelope-loss observation y_L=L+ε_L, and a normalised thermal-load observation y_H=H+ε_H. Their observation standard deviations are 0.003, 0.008 and 0.020, respectively. In the controlled recovery experiment, D is interpreted specifically as a normalised effective stiffness-loss coordinate under an unchanged effective-mass assumption. The frequency relation is therefore an illustrative reduced-order observation function rather than a building-specific modal model.

Sequential importance-resampling particle filtering was used for state tracking under uncertain parameters (Särkkä 2013). The augmented particle vector was ${\xi }_k={[D_k,L_k,H_k,k_D,q_D,c_H,c_D,u_2]}^T$. The remaining parameters k_L, c_LH, u₁ and τ were held at their reference values. All particles were initialised at the known initial state ${[\mathrm{0.02,0.02,1.00}]}^T$, while the five parameter particles were drawn independently from lognormal priors with the biased means and coefficients of variation reported below. At each quarterly update, particles were propagated through the state and parameter transition models, weighted with the Gaussian observation likelihood, normalised, and systematically resampled only when the effective sample size fell below 50% of the 2,500-particle population. Posterior state means and 95% intervals were calculated from the normalised weighted particle population before resampling. Quarterly state-process standard deviations were 0.00015, 0.00020 and 0.00150 for D, L and H. Parameter particles followed a log-space random walk with standard deviation 0.00150 per update. The synthetic-observation and particle-filter seeds were 20260805 and 20260827. Parameter-particle trajectories represent propagated uncertainty and are distinct from the regularised point estimates reported separately.

A separate regularised parameter-recovery step was applied to k_D, c_D and u₂. This step produced regularised point estimates rather than a complete posterior parameter distribution. Particle-based state estimation and regularised parameter recovery are therefore reported separately. The biased prior means relative to the reference values were 0.88 for k_D, 0.85 for q_D, 0.90 for c_H, 0.90 for c_D and 0.82 for u₂. Corresponding prior coefficients of variation were 22%, 25%, 25%, 25% and 22%.

Regularised parameter recovery minimised the weighted observation-residual and prior-deviation objective:

Here $R=diag({0.003}^2,{0.008}^2,{0.020}^2)$, $\lambda =1.0$, and $W=diag[(0.22k_{D,prior}{)}^{-2},$ ${(0.25c_{D,prior})}^{-2},$ ${(0.22u_{2,prior})}^{-2}]$. The three diagonal entries of W correspond to the inverse squared prior standard deviations of the recovered parameters k_D, c_D, and u₂, respectively. The optimisation variables were bounded to 0.5–1.5 times their reference values and initialised at the biased prior means. L-BFGS-B was used with a maximum of 2,000 iterations, function tolerance 10^-12 and projected-gradient tolerance 10^-8. Each objective evaluation used the same monthly RK4 state integration and quarterly observation schedule.

For recovery experiment r and state $s\in \{D,L,H\}$, trajectory RMSE and pointwise interval coverage were defined as:

The reported state errors are the medians of RMSEs(r) across 100 experiments. Onset error was the predicted CSEI zero-crossing time minus the reference zero-crossing time; median absolute onset error and mean signed onset error were reported separately. All 100 controlled recovery trajectories crossed within the 25-year horizon, so no onset value was censored in this experiment. The repeated experiments were correctly specified inverse tests because the synthetic observations and estimator used the same state and observation model classes. They assess numerical and statistical recoverability under the prescribed sensor configuration, not robustness to deterioration-law, observation-function or structural model discrepancy.

4.2. Uncertainty propagation

All positive state-model parameters were sampled independently from lognormal distributions with the reference means reported in Table 3. Coefficients of variation were 15% for k_D and k_L; 20% for q_D, c_H, c_D and c_LH; 12% for u₁; 15% for u₂; and 10% for τ. A total of 5,000 parameter–climate realisations were generated using master seed 20260728. The same draws were used for all structural-coupling comparisons. Reported trajectory intervals are empirical 2.5–97.5 percentile intervals across realisations; uncertainty in ensemble means is reported separately through Monte Carlo standard errors and confidence intervals. Parameter dependence was not represented because no joint asset-calibration dataset was available. Independent lognormal marginals are therefore a scenario assumption and may understate or overstate tail behaviour when deterioration and coupling coefficients are correlated.

The climate disturbance evolved continuously across the full 300-month horizon and was not reset at calendar-year boundaries:

Here ρ=0.65 and σ_ε=0.04, making 0.04 the stationary standard deviation and σ_ε√(1−ρ²) the innovation standard deviation. The initial deviation was drawn from the stationary distribution.

4.3. Nine-parameter global sensitivity analysis

The sensitivity analysis covered the nine uncertain positive state-model parameters listed in Table 3 and used the same independent lognormal marginals as the Monte Carlo analysis. It did not vary the baseline EUI, monthly load shares, HVAC share, grid-emission assumptions, initial states or model structure. Accordingly, “global” refers to the prescribed nine-dimensional state-parameter domain.

Jansen total-effect indices were estimated using a scrambled Sobol base sample of N=2,048 following Jansen (1999), Sobol (2001) and Saltelli et al. (2010). Bootstrap uncertainty was evaluated by resampling the 2,048 base-row indices with replacement while preserving the corresponding A, B and ABi row blocks. Total-effect indices were recalculated for 500 paired percentile-bootstrap resamples generated with seed 20260728, and the 2.5th and 97.5th percentiles defined the reported intervals.

The primary outputs were cumulative instability exposure J_inst and restricted time to onset T_R as defined in Equations (13–14). The binary crossing indicator was retained only as a secondary diagnostic because its near-unity crossing probability makes variance-based indices less stable. Only total-effect indices were estimated; the analysis does not claim a separate first-order interaction decomposition.

5. Results

5.1. Limiting cases and onset of instability

All percentage comparisons were calculated from unrounded simulation outputs. The one-way case separates the effect of a deteriorating envelope on electricity from the return feedback to deterioration (Figure 4 and Table 5). It raises lifecycle electricity from 27.48 to 31.90 GWh but leaves D unchanged. Activating the return path raises final D by 7.2% and produces a positive CSEI after 17.58 years. Relative to the Energy-Energy-decoupled model, lifecycle electricity rises by 16.7% and carbon by 15.8%. The smaller carbon percentage follows from the declining grid-factor scenario, which reduces the carbon weight of later energy amplification.

5.2. Probabilistic trajectories

Prefix estimates from N=1,000 to N=5,000 remained close to the final 5,000-realisation results (Figure 6). At N=5,000, the Monte Carlo standard error was 0.000533 for mean final D, 0.004078 ktCO₂e for mean lifecycle carbon and 0.001796 for crossing probability. The corresponding normal-approximation 95% confidence intervals were 0.336028–0.338117 for mean final D and 11.543838–11.559822 ktCO₂e for mean lifecycle carbon (Figure 5 and Table 6). The exact two-sided 95% Clopper–Pearson interval for crossing probability was 0.9797–0.9869. Prefix convergence is presented as a numerical diagnostic and not as an independent uncertainty experiment.

5.3. Updating performance

5.4. Panel A — Reference recovery experiment

5.5. Panel B — Performance across 100 correctly specified synthetic recovery experiments

In the reference recovery experiment, the biased prior-only trajectory placed the CSEI transition 4.25 years later than the monthly-grid reference onset (Figure 7). Sequential particle filtering reduced state-tracking errors, while the separate regularised parameter-calibrated trajectory reduced the absolute onset-reconstruction error to 0.42 years (Table 7). Across 100 correctly specified synthetic recovery experiments, median state RMSEs remained low, the median absolute onset-reconstruction error was 0.56 years, and nominal 95% state intervals achieved 94.0–95.2% pointwise coverage. Because the full synthetic trajectories were used in these controlled inverse tests, these results represent retrospective recovery and reconstruction rather than prospective pre-onset forecasting. They do not establish identifiability under arbitrary sensor layouts or robustness to deterioration-law, observation-function or structural-model discrepancy.

5.6. CSEI verification and sensitivity

The matrix-exponential calculation reproduces the spectral-abscissa rate of the same frozen Jacobian to numerical precision (Table 8). At year 10 the dominant frozen mode decays, whereas at years 20 and 24 it grows. This is an internal modal-consistency check; because both quantities derive from the same Jacobian, it is not an independent validation of the state model or a proof of global instability.

The frozen-eigenmode calculation confirms the intended local interpretation of CSEI. The singular-value results additionally show that finite-time norm amplification may occur while CSEI remains negative because the system is time varying and non-normal (Table 9). The reported gains depend on the adopted normalised state coordinates and Euclidean norm; they are not invariant physical-energy amplification factors. H2 is therefore supported only as detection of a frozen-modal sign change, not as proof that every perturbation decays before that crossing.

For cumulative instability exposure, q_D had the largest total effect, ST=0.4682 (95% CI: 0.4010–0.5400), followed by c_H, ST=0.1946 (0.1500–0.2450), and c_D, ST=0.1738 (0.1320–0.2220) (Figure 8 and Table 10). For restricted onset, the leading indices were q_D=0.2605 (0.2422–0.2801), c_H=0.2486 (0.2295–0.2652) and c_D=0.2213 (0.2058–0.2377). The remaining restricted-onset intervals were k_D 0.0089–0.0113, k_L 0.1204–0.1437, c_LH 0.0367–0.0444, u₁ 0.0494–0.0590, u₂ 0.0885–0.1049 and τ 0.00009–0.00024. H4 is therefore supported: uncertainty in nonlinear acceleration and the bidirectional operational feedback pathway dominates instability behaviour.

6. Discussion

6.1. Evaluation of the hypotheses

H1 is supported within the locked reference scenario. Relative to the energy-decoupled configuration, bidirectional coupling increased lifecycle electricity by 16.7%, lifecycle carbon by 15.8% and final D by 7.2%. The one-way configuration showed that most of the electricity increase arose from the D-to-L-to-H pathway, whereas the additional increase in D required activation of the return H-to-D pathway.

H2 is supported only in its stated frozen-modal sense. On the monthly evaluation grid, the nominal CSEI first became non-negative at 17.58 years, when D, L and H were approximately 0.2462, 0.1872 and 1.4582. The matrix-exponential calculation reproduced the frozen spectral-abscissa rate. Finite-time singular-value analysis nevertheless showed that negative CSEI does not exclude transient norm amplification. CSEI must therefore be interpreted as a local modal-transition diagnostic rather than a complete global-stability measure.

H3 is supported for the correctly specified synthetic inverse tests. In the reference experiment, sequential assimilation reduced state RMSEs by 45.3–91.2%, while regularised parameter recovery reduced the absolute onset-reconstruction error from 4.25 to 0.42 years. Across 100 experiments, the median absolute onset-reconstruction error was 0.56 years and nominal 95% state intervals achieved 94.0–95.2% pointwise coverage. These results do not constitute prospective forecasting or external model-form validation.

H4 is supported within the nine-dimensional lognormal parameter domain. Nonlinear deterioration acceleration q_D was the dominant driver of cumulative exposure, while c_H and c_D demonstrated the importance of both directions of the structural–envelope–thermal feedback. The indices do not include uncertainty in the baseline EUI, HVAC share, initial states, grid-emission trajectory or model structure.

6.2. Validation achieved beyond equation verification

The evidence extends beyond reproduction of the state equations but remains deliberately hierarchical. The electricity-observation layer was evaluated one step ahead on a complete unseen calendar year. Numerical behaviour was examined through limiting cases, time-step refinement, frozen-modal consistency and finite-time diagnostics. State estimation and regularised recovery were examined through 100 correctly specified synthetic inverse experiments. Probabilistic outputs were supported by 5,000 realisations, Monte Carlo sampling-error estimates and bootstrap intervals for the nine-parameter total-effect indices. These steps establish empirical validity for the operational observation layer and internal numerical, statistical and inference consistency for the coupled workflow. They do not replace long-term field validation of structural and envelope deterioration.

6.3. Engineering and asset-management implications

The results indicate that structural condition, envelope performance, and HVAC operation should be treated as interacting components of a unified lifecycle-management system rather than as independent monitoring streams. In a practical implementation, the structural/support-condition state, (D), could be informed by calibrated modal features, strain measurements, or inspection-derived condition indices; the envelope-loss state, (L), by airtightness testing, façade moisture sensing, thermographic surveys, or inferred heat-loss coefficients; and the thermal-load memory state, (H), by plant power demand, HVAC runtime, and indoor–outdoor thermal gradients. Joint assimilation of these observations at regular intervals would allow the posterior state estimates, CSEI distribution, and projected onset time to be updated as new evidence becomes available.

Operational deployment would begin with an asset-specific baseline survey and a data-integration layer linking structural-monitoring systems, façade and envelope inspections, utility metering, and building-automation-system records. These heterogeneous observations would be mapped to the corresponding (D), (L), and (H) state variables and assimilated through the updating framework. The resulting posterior states would then propagate into the CSEI and lifecycle energy–carbon modules, allowing updated diagnostics and projections to be communicated through a facility- or asset-management platform. Such an implementation could support risk-informed structural inspections, prioritisation of façade and envelope remediation, targeted HVAC maintenance, and periodic model recalibration. Before operational decision support, however, the framework would require asset-specific parameter calibration, prospective validation of state-recovery and forecast performance, and explicit assessment of model discrepancy using co-located field observations.

The practical value of the framework therefore lies primarily in the evolution of the joint state rather than in any single reference onset year. A building may exhibit moderate structural/support deterioration, progressive envelope loss, and persistent thermal loading without any individual variable exceeding a conventional alarm threshold, while their coupled dynamics simultaneously move the system toward a non-negative CSEI. This finding supports coordinated assessment of structural, façade, and mechanical systems and cautions against lifecycle decisions based solely on isolated damage indices, fixed condition thresholds, or annual energy-use indicators.

The implications also extend directly to lifecycle energy and carbon management. Within the locked reference scenario, bidirectional coupling increased lifecycle electricity consumption by 16.7% relative to the energy-decoupled configuration. This result illustrates why a commissioning-year energy model should not automatically be extrapolated across the full-service life when envelope condition, HVAC operation, and deterioration-related feedback can evolve materially over time. Degradation-aware operational projections, periodically updated with measured asset data, can therefore provide a more defensible basis for long-term maintenance planning, energy management, and carbon-performance assessment.

6.4. Generalisation assumptions and limitations

Framework generalisation and numerical generalisation must be distinguished. The architecture is transferable, but the numerical parameterisation is not. For a new asset, the state definitions and observation functions should be adapted to the structural system and envelope type, while climate forcing, occupancy schedules, HVAC topology, baseline EUI and uncertainty distributions should be re-estimated from local data. Application to reinforced-concrete, steel, timber, masonry or hybrid buildings and to other climates therefore requires asset-specific recalibration rather than direct reuse of the 17.58-year onset or the reported lifecycle percentages.

The reported numerical results are conditional on a controlled reference parameterisation rather than an asset-calibrated deterioration model. Empirical validation is limited to one-step-ahead electricity observations from one office meter. The structural and envelope states, observation mappings and feedback coefficients were evaluated through correctly specified synthetic recovery experiments and therefore do not represent external model-form validation.

The uncertainty intervals exclude uncertainty in the baseline EUI, HVAC share, monthly load profile, grid-emission trajectory, initial states, activation-function form, parameter dependence and model structure. The independent lognormal assumption may understate or overstate tail behaviour when coefficients are correlated.

The nominal 17.58-year onset is the first non-negative CSEI value on a monthly evaluation grid and applies only to the locked reference scenario. CSEI is the spectral abscissa of a frozen local Jacobian; it is not a measure of structural safety, global nonlinear stability, collapse probability, occupant comfort or code compliance. Prospective deployment requires co-located structural, envelope and operational observations, asset-specific parameter calibration, pre-onset forecast validation, explicit model-discrepancy treatment and independent field evaluation.

7. Conclusions

This study developed an energy-data-anchored probabilistic digital-twin framework for long-horizon structural/support-condition, envelope and thermal-load interaction. The electricity-observation model retained its one-step-ahead performance on an independent 2017 calendar year, achieving R²=0.933, CV(RMSE)=6.65% and NMBE=0.029%.

Within the locked 25-year reference scenario, bidirectional coupling increased lifecycle electricity and carbon by 16.7% and 15.8% relative to the energy-decoupled configuration. CSEI first became non-negative at month 211, corresponding to 17.58 years on the monthly evaluation grid. Across 5,000 parameter–climate realisations, mean lifecycle electricity and carbon were 32.091 GWh and 11.552 ktCO₂e, and 98.36% of trajectories crossed by year 25, with a conditional monthly-grid median onset of 17.67 years.

Sequential particle filtering reduced state-tracking errors, while separate regularised parameter recovery improved retrospective reconstruction of the CSEI onset. Across 100 correctly specified synthetic inverse tests, state-interval coverage remained close to the nominal 95% level. Nine-parameter total-effect analysis identified nonlinear deterioration acceleration and the bidirectional feedback coefficients as the principal drivers of the two instability outputs. For building owners and facility managers, the practical implication is that structural condition, envelope performance and HVAC operation should be treated as a jointly updated asset-management problem rather than as independent monitoring streams. The framework can support risk-informed inspection timing, façade and HVAC maintenance prioritisation, and revision of long-term energy and carbon forecasts as new observations become available. In the reference scenario, the difference between the energy-decoupled and bidirectionally coupled cases—16.7% in lifecycle electricity and 15.8% in lifecycle carbon—quantifies the potential consequence of neglecting deterioration–energy feedback within the stated modelling assumptions.

The transferable contribution is the integration of an empirical electricity-observation layer, a coupled lifecycle state model, sequential updating, conditional uncertainty propagation and a deliberately limited local modal diagnostic. The numerical percentages and onset times are reference-scenario outputs. They should not be interpreted as field-validated deterioration predictions, global-stability results or structural-safety criteria.

Appendix A. Reproducibility settings

Monthly state integration used fixed-step classical fourth-order Runge–Kutta with one primary step per month. Time-step verification used 2, 4 and 8 substeps per month. The lifecycle horizon was 25 years. The nominal onset was defined as the first non-negative CSEI value on the corresponding evaluation grid; no interpolation between grid points was applied. Parameter–climate uncertainty used 5,000 realisations and master seed 20260728. Particle filtering used 2,500 particles, systematic resampling at ESS<0.50N, synthetic-observation seed 20260805 and particle-filter seed 20260827. The nine-parameter Sobol analysis used a scrambled base sample of 2,048 and 500 paired bootstrap resamples. Numerical diagnostic bounds were 0≤D,L≤0.99 and 0.05≤H≤5.0; no reported reference or uncertainty trajectory reached these bounds.

State endpoints were invariant to the reported precision. The first non-negative CSEI value occurred at 17.5833 years with one substep per month and at 17.5208 years with four and eight substeps per month. The difference of 0.0625 years is approximately 0.75 months and reflects the discrete grid definition of onset. All principal reported analyses use the locked monthly grid; the onset is therefore reported as 17.58 years with an explicit one-month temporal resolution.