2.1. Heating ventilation and air conditioning (HVAC)

The laboratory is mechanically ventilated via a typical variable air volume (VAV) system controlled by a building management system (BMS). Temperature-controlled fresh air is pumped into a plenum above the room (see the blue pipe in Figure 1). Four fan coil units (FCUs) (Quartz Sapphire model SPR9) draw air from the plenum and heat/cool the air using a hot/cold water system according to the temperature setpoint for the room. The conditioned air is ducted from the FCUs to supply the room via linear slot diffusers placed at the north and south sides of the room (blue pipes and arrows in Figure 1). The rate of supply to the room is usually different from the rate of extraction, which is constant and occurs through five square grilles positioned along the longitudinal axis of the ceiling (red pipes in Figure 1).

To account for the possible mismatch between the rate at which air is supplied by the FCUs and extracted according to the VAV system, the plenum and the room exchange air through two additional grilles on the ceiling (not shown). In addition, sections of the linear slot diffusers are not connected to the FCUs, and allow air exchange between the room and the plenum. These air exchanges can occur in either direction (from the room towards the plenum and vice versa) depending on the air handling unit (AHU) regime.

2.2. Occupancy

Occupancy estimates are inferred from Imperial’s deployment of the HubStar (formerly LoneRooftop) Building Insights Dashboard, which records the number of Wi-Fi connections (hourly average and maximum values are available). Wi-Fi estimated occupancy has several advantages, such as protecting the privacy of occupants and achieving a high accuracy ( $96\,\%$ according to the study by Simma et al. (Reference Simma, Mammoli and Bogus2019)). These advantages and the wide availability of Wi-Fi signals in buildings have made it a convenient means of occupancy estimation in many studies in the past ten years (Zhao et al. Reference Zhao, Li, Liang and Wang2022). Nevertheless, estimates have uncertainties associated with the varying number of Wi-Fi devices each person carries. In our case, the time resolution and the anonymisation of the recording of Wi-Fi connections lead to the estimated uncertainty of the order of 30 % (see Appendix B.1 in the supplementary material available at https://doi.org/10.1017/flo.2026.10045 for more details).

2.3. Measurements

The laboratory is monitored continuously with 24 temperature sensors (Trend thermistors model T/TFR 4) and 8 combined CO $_2$ /temperature/humidity sensors (TREND space sensors model RS-WMB-THC) connected to the BMS (via 4 Trend IQ4 controllers). The temperature sensors are installed on 4 vertical risers (i.e. 6 per riser) at different positions in the room. The combined sensors are situated on vertical walls at different horizontal positions in the room (see Figure 1 for more details). The CO $_2$ and temperature measurements were further supplemented by sensors in the supply and extract ducting of the BMS system.

Figure 2. Monthly distributions of vertical heterogeneities (expressed as $T-\overline {T}$ , with $\overline {T}$ the room averaged temperature at a given time) from the temperature sensors for June–November 2022.

The dataset presented in this paper spans from June 2022 to November 2022, with measurements sampled and recorded every 5 minutes. An overview of the temperature data collected from the sensors on the risers is shown in Figure 2. Each panel shows the distribution of the temperature measured by the sensors at the given height, after subtracting the mean room temperature at each measurement time. The largest departures from the mean occur at the lowest two measurement levels. The sensor at 0.5 m is located beneath the desks of the room and the one at 1.0 m is close to the computers; both are therefore influenced by local heat sources and partial shielding, which explains their larger deviations from the rest of the vertical profiles. For the other heights, the distributions lie within 1 $^{\circ }$ C, with the exception of a small number of outliers. This narrow distribution for most of the data suggests that the room temperature is rather uniform along the vertical.

For the rest of the present study, we focus on the data provided by the 8 combined CO $_2$ /temperature/humidity sensors. We assume that the occupants are the only source of CO $_2$ in the room.

We supplement the data acquired in the room with a range of standard measurements recorded by the BMS at 5-min intervals. These measurements include the supply and extraction volume flow rates, $Q_{\text{in}}$ and $Q_{\text{out}}$ , respectively, and the supply and return temperatures from the FCU. The flow rates $Q_{\text{in}}$ and $Q_{\text{out}}$ are measured by sensors placed in the air ducts. The data are presented in § 4.

3.1. CO $_2$ budget

Let $C_{0}$ and $C_{1}$ be the bulk excess CO $_2$ concentration (relative to outdoors) in the room and plenum, respectively, and let $C_{0}'$ and $C_{1}'$ be the quasi-steady-state excess concentration associated with the ceiling zone and FCU, respectively (see Figure 3 for a schematic representation of the spaces in a vertical cross-section of the room). We use $q$ and $q'$ to denote the rates at which fresh air is supplied to the plenum and at which air is drawn into the plenum from the ceiling zone, respectively, and define $Q:=q+q'$ (see Figure 3). The variable $Q$ , therefore, corresponds to the volume flux driven by the FCUs.

Figure 3. (a) Schematic diagram (not to scale) of a vertical section through the ceiling plenum and computer laboratory below, with excess concentrations $C_1$ and $C_0$ , respectively; the corresponding control volumes are highlighted with a red dashed line. $C'_1$ and $C'_0$ indicate the excess CO $_2$ concentration in the FCU and ceiling zone (highlighted with a darker shade of grey). As illustrated by the white arrows, conditioned air is supplied to the plenum with a volume flux $q$ . A fraction $\gamma _{1}\in [0,1]$ of that supply air mixes with the air within the plenum; the rest is drawn directly into the fan coil unit (FCU), whose fan drives a total volume flux $Q=q+q'$ . After being fed into the computer laboratory, a fraction $\gamma _{0}\in [0,1]$ of air from the FCU either mixes with the air in the room, while the rest comprises a ‘short circuit’ and remains in the ceiling zone. Due to the fact that $Q\neq q$ , the FCU drives a secondary volume flux $q'=Q-q$ between the computer laboratory and plenum. (b) Graph corresponding to the governing equations described in (3.2). Each node corresponds to a control volume (labelled with CO $_2$ concentration) and each branch corresponds to a flux of ${CO}_{2}$ between control volumes (labelled with the corresponding volume flux).

First, we establish an equation for the CO $_2$ concentration in the FCUs, the combined volume of which is labelled $V_{1}'$ . Unless stated otherwise, all CO $_2$ concentrations given refer to excess values.

The dimensionless parameters $\gamma _{0}$ and $\gamma _{1}$ encode how the FCU flow is partitioned between different pathways within the plenum–room system. Since the supply air duct is not connected directly to the FCUs, let $\gamma _{1}\in [0,1]$ represent the fraction of the supply rate $q$ that first mixes with the air in the plenum before entering the FCU:

(3.1) \begin{equation} V_{1}'\frac {\mathrm{d} C_{1}'}{\mathrm{d} t}=\underbrace {(\gamma _{1}q+q')C_{1}}_{\mathrm{from\ plenum}} -\underbrace {QC_{1}'}_{\mathrm{to\ room}}, \end{equation}

for which we have assumed that $\gamma _{1}q+q'\geq 0$ , i.e. the FCU always extracts air from the plenum.

For the control volume just below the room’s ceiling, we adopt a similar approach to the FCU, using $\gamma _{0}\in [0,1]$ to denote the fraction of FCU outlet air that first mixes with the rest of the room (before entering the ceiling zone):

(3.2) \begin{equation} V_{0}'\frac {\mathrm{d} C_{0}'}{\mathrm{d} t}= \underbrace {\gamma _{0}QC_{0}}_{\mathrm{from\ room}}+\underbrace {(1-\gamma _{0})QC_{1}'}_{\mathrm{from\ FCUs}}-\underbrace {q'\phi }_{\mathrm{to\ plenum}} -\underbrace {qC_{0}'}_{\mathrm{to\ extract}}, \end{equation}

where

(3.3) \begin{equation} \phi =\begin{cases} C_{0}' &\quad q'\geq 0, \\ C_{1} &\quad q'\lt 0, \end{cases} \end{equation}

allows the flux of CO $_{2}$ to be represented as either from the ceiling zone to the plenum ( $q'\geq 0$ ) or from the plenum to the ceiling zone ( $q'\lt 0$ ).

Let us now consider the CO $_2$ concentration in the room and plenum. We assume that the rate of CO $_2$ generation per person is $F = 0.012$ g s⁻¹person⁻¹, consistent with the range $F = 0.009{-}0.012$ g s⁻¹person⁻¹ for occupants aged 21–30 years undertaking typical office work (Persily & de Jonge Reference Persily and de Jonge2017). Indeed, our data (to be discussed in Figure 8) will indicate that $F = 0.012 \pm 0.001$ g s⁻¹ person⁻¹ provides the best fit for the model, as evaluated by the highest $R^2$ value. Furthermore, we assume that the infiltration and exfiltration of air are negligible (the validity of such an assumption is discussed in Appendix B.2). The system of linear ordinary differential equations (ODEs) describing the time-variation of concentration in the room and the plenum respectively is therefore

(3.4) \begin{align} V_{0}\frac {\mathrm{d} C_0}{\mathrm{d} t} & = \overbrace {\gamma _{0}QC_{1}'}^{\mathrm{from\ FCU}}-\overbrace {\gamma _{0}Q C_{0}}^{\mathrm{to\ ceiling\ zone}}+NF, \\[2pt] \nonumber \end{align}

(3.5) \begin{align} V_{1}\frac {\mathrm{d} C_1}{\mathrm{d} t} & = \underbrace {q'\phi }_{\mathrm{from\ ceiling\ zone}} - \underbrace {(\gamma _{1}q+q')C_{1}}_{\mathrm{to\ FCU}}, \\[10pt] \nonumber \end{align}

where $N$ is the number of occupants. The system above is similar to the two-zone transient model proposed by Lawrence and Braun (Reference Lawrence and Braun2006), with the difference that we have relaxed the condition of complete mixing in the plenum ( $V_1$ ) and occupied zone ( $V_0' + V_0$ ) by adding the $\gamma$ parameters. The governing equations can be graphically represented in the network diagram in Figure 3b, where the nodes represent control volumes labelled with concentrations. The volume fluxes between control volumes are represented by the network’s branches.

Let us now re-write the governing equations (3.1), (3.2), (3.4) and (3.5) in non-dimensional form by introducing the following variables:

(3.6) \begin{equation} \tau := \frac {q}{V_1}t, \quad \varepsilon := \frac {V_1}{V_0}, \quad \zeta := \frac {q'}{q}, \quad f_{0} := \frac {NFV_1}{q V_0}. \end{equation}

Physically, $\tau$ is a dimensionless time based on the filling time scale associated with the plenum (in terms of its volume $V_{1}$ and the AHU supply volume flux $q$ ), $\varepsilon$ represents the ratio of the plenum and room volumes, $\zeta$ is the ratio of the rates of room-plenum re-circulation and ambient fresh air introduction, and $f_0$ is the dimensionless forcing.

We further assume that the FCUs and ceiling buffer zone volume are much smaller than the plenum volume, i.e. $V_1'/V_1 \ll 1$ and $V_0'/V_1 \ll 1$ . In physical terms, storage of CO $_2$ in the FCU and ceiling buffer zone, represented by the time derivative in (3.1)–(3.2), is therefore assumed to be insignificant, such that (3.1) and (3.2) give explicit algebraic expressions relating $C_0'$ and $C_1'$ to $\gamma _0$ , $\gamma _1$ , $q$ and $q'$ .

Under these assumptions, the resulting non-dimensionalised equations are

(3.7) \begin{align} \frac {\mathrm{d} C_0}{\mathrm{d} \tau } & = - \varepsilon \gamma _0(\zeta +1) C_0 + \varepsilon \gamma _0(\zeta + \gamma _1)C_1 + {f_0} , \\[-10pt] \nonumber\end{align}

(3.8) \begin{align} \frac {\mathrm{d} C_1}{\mathrm{d} \tau } & = \zeta \phi - (\zeta + \gamma _1)C_1. \\[8pt] \nonumber \end{align}

We can express the system (3.7) and (3.8) in the form

(3.9) \begin{equation} \frac {\mathrm{d}\boldsymbol{C}}{\mathrm{d}\tau } = \boldsymbol{A}\boldsymbol{C}+\boldsymbol{f}, \end{equation}

where $\boldsymbol{A}$ is a coefficient matrix, $\boldsymbol{C} = (C_0, C_1)^{\top }$ is the state vector and $\boldsymbol{f}=(f_0, 0)^{\top }$ is the forcing vector.

Assuming that $\boldsymbol{A}$ is diagonalisable (see discussion in § 5.2 for the general case in which this assumption is not valid), we can derive the solution for (3.9),

(3.10) \begin{equation} \boldsymbol{C} = \boldsymbol{R}\exp (\boldsymbol{\Lambda } \tau )\boldsymbol{R}^{-1}(\boldsymbol{C}(0) +\boldsymbol{A}^{-1}\boldsymbol{f})-\boldsymbol{A}^{-1}\boldsymbol{f}, \end{equation}

where $\boldsymbol{C}(0)$ are the initial conditions, $\boldsymbol{\Lambda }$ is the eigenvalue matrix and $\boldsymbol{R}$ is the corresponding (right) eigenvector matrix. The derivation of (3.10) and explicit expressions for the eigenvalues and eigenvectors are provided in Appendix A.

3.2. Phase space diagrams

Each possible state of the system is represented by a point in phase space, and a sequence of consecutive states forms a trajectory. Therefore, phase space diagrams are useful for representing systems geometrically, and investigating their equilibria, stability and evolution.

Figure 4. Phase space diagram corresponding to the variable $\boldsymbol{y}$ (3.12, which accounts for the stationary equilibrium $\boldsymbol{A}^{-1}\boldsymbol{f}$ ) for fixed $\varepsilon = $ 0.17, $\zeta =$ 1.8, and several combinations of $\gamma _0$ and $\gamma _1$ values. The grey arrows show the flow evolution trajectories for different initial conditions and the black dots show one particular solution of the dynamical system at uniform time intervals. The red and blue arrows show the two eigenvectors.

If the forcing $\boldsymbol{f}$ is regarded as constant, then (3.9) can be expressed as the autonomous system:

(3.11) \begin{equation} \frac {\mathrm{d}\boldsymbol{y}}{\mathrm{d}\tau } = \boldsymbol{A}\boldsymbol{y}, \end{equation}

where

(3.12) \begin{equation} \boldsymbol{y} := \boldsymbol{C}+\boldsymbol{A}^{-1}\boldsymbol{f}. \end{equation}

Figure 4 illustrates the phase space associated with (3.11). We choose $\varepsilon = 0.17$ and $\zeta =1.8$ in (3.7) and (3.8) because they represent the measurements in the living laboratory and are characteristic of a ventilated office space. The values of $\gamma _0$ and $\gamma _1$ are varied to understand how the (parametrised) mixing affects the rate of convergence towards one manifold of the system. The grey arrows in Figure 4 mark flow evolution trajectories for different initial conditions. The red and blue arrows mark the directions of the two eigenvectors and correspond to the fast (large negative eigenvalue) and slow (small negative eigenvalue) dynamics of the system, respectively.

The black dots show one particular solution of the dynamical system at uniform time intervals. They show that the solutions converge rapidly towards the slow eigenmode, which dominates the adjustment to the equilibrium state $\boldsymbol{y}=\boldsymbol{0}$ .

In terms of the parameters, we can easily see that $\gamma _1$ has far less influence on the dynamics than $\gamma _0$ . This behaviour comes from the condition $\varepsilon \ll 1$ , which corresponds to the volume of the plenum being much smaller than the volume of the room below ( $V_1 \ll V_0$ ). When $\gamma _0 \to 1$ and $\gamma _1 \to 0$ , the CO $_2$ concentrations in the room and the plenum become comparable, as can be seen by the blue vector approaching a slope of one in the bottom right panel in Figure 4. When $\gamma _0 = 1$ and $\gamma _1 = 0$ , the slow manifold maps onto $y_0=y_1$ , while the fast manifold corresponds to the difference between the CO $_2$ concentration in the room and plenum. In all other cases shown in Figure 4, the concentration is higher in the room than in the plenum. More generally, for other combinations of $\gamma _{0}$ and $\gamma _{1}$ , the physical interpretation of the slow and fast manifolds is less straightforward, as they involve other linear combinations of the concentrations in the rooms, which depend on the mixing properties of the system.

5.1. Model reduction

We revisit the model presented in § 3.1 in the light of the behaviour highlighted by the comparison with the measured data. Numerous situations have emerged in which there is a sufficiently large separation of time scales that, from the perspective of the relatively slow evolution of the overall system, the transient effects associated with the fast dynamics can be neglected. Exploiting a separation of time scales in this way is used widely to derive a simplified low-dimensional description of a physical system and, in broad terms, corresponds to the ‘adiabatic elimination’ of fast variables in stochastic systems (Haken Reference Haken1983).

The equation governing the second (fast) eigenmode can be written (Appendix A; (35))

(5.1) \begin{equation} \frac {\mathrm{d} z_2}{\mathrm{d}\tau } = \lambda _2 z_2+ \frac {1}{\lambda _2 R_{21}}\frac {\mathrm{d} f_0}{\mathrm{d} \tau } , \end{equation}

where $|1/\lambda _2|$ is the response time scale of the eigenmode, $z_2$ is the amplitude of the eigenmode and $\mathbf{R}$ is the eigenvector matrix. This reduction is valid only in parameter regimes where $\rho :=\max (\lambda _1, \lambda _2)/\min (\lambda _1, \lambda _2)\ll 1$ ; in the low- $\zeta$ limit with $\gamma _0 \approx 1$ or small $\gamma _1$ , where $\lambda _1 \approx \lambda _2$ , the fast–slow splitting becomes ill-conditioned and the full two-dimensional system must be used (see § 6).

Since $\lambda _{2}\lt 0$ , we can deduce that an upper bound for finite changes $\Delta z_{2}$ is obtained by ignoring the first term on the right-hand side of (5.1): $\Delta z_{2}\leq \Delta f_{0}/(\lambda _{2} R_{21})$ . The upper bound becomes tighter in the limit that the corresponding $\lambda _{2}\Delta t\rightarrow 0$ . For larger values of $\Delta t$ , the upper bound becomes increasingly conservative. In other words, it changes in $f_0$ over time scales shorter than $1/\lambda _2$ that are likely to lead to deviations from $z_{2}\approx 0$ , which would remain a valid approximation, and therefore enable a reduction of the dimension of the system if $\Delta f_{0}\ll \lambda _{2}R_{12}$ . In practice, this condition might not be fulfilled when there are rapid changes in occupancy, such as at the beginning and end of a class.

By assuming that the equation for $z_{2}$ plays no role, the (one-dimensional) reduced model can be written (Appendix A; (35)) as

(5.2) \begin{equation} \frac {\mathrm{d} z_1}{\mathrm{d}\tau } = \lambda _1 z_1+ \frac {1}{\lambda _1}R^{-1}_{11}\frac {\mathrm{d} f_0}{\mathrm{d} \tau } , \end{equation}

and the solution (3.10) reduces to

(5.3) \begin{equation} \boldsymbol{C} = \boldsymbol{R}_{1}\exp (\lambda _1 \tau )\boldsymbol{R}_{1}^{-1}(\boldsymbol{C}_0 +\boldsymbol{A}^{-1}\boldsymbol{f})-\boldsymbol{A}^{-1}\boldsymbol{f}, \end{equation}

where $\boldsymbol{R}_{1}$ is the column vector of the eigenvector matrix (corresponding to the two components of one eigenvector) and $\boldsymbol{R}_{1}^{-1}:=\boldsymbol{R}_{1}^{\top }/(\boldsymbol{R}_{1}^{\top }\boldsymbol{R}_{1})$ is the left inverse. Note that (5.3) depends only on the eigenvalue and eigenvector associated with the slow dynamics.

Figure 9. (a) Comparison of the reduced (green line) and full model (blue line) with the observational data (black dots) for the month of August. The parameters for the two models are the same as the one used in Figure 6. (b) Occupancy data.

The full and reduced models are compared in Figure 9a (blue and green lines, respectively), where the time evolution of $C_0$ is plotted for a chosen day in August (same data as in Figure 6). The reduced model prediction (given by (5.3)) is consistent with the full model prediction (given by the integration of (3.4)). We can notice that the rapid changes in occupancy are approximated in the reduced order system by an instantaneous adjustment of the fast mode and are responsible for the differences between the outputs from the full and reduced order models (panel a).

In our study, the pronounced separation between fast and slow time scales arises from the specific geometry and flow configuration of the living laboratory. In particular, the plenum volume is much smaller than the occupied-room volume ( $\varepsilon \ll 1$ ) and when the FCUs are operating (ON regime), the data are well described by $\gamma _0 \approx 1$ , indicating that most of the supply air is well mixed in the occupied room. Under these conditions and assuming that the effective ceiling layer is thin, the excess CO $_2$ concentration in that layer is close to the bulk room concentration, so that additional storage in the ceiling zone does not introduce a distinct dynamical time scale. Consequently, the dynamics collapse rapidly onto a one-dimensional slow manifold associated with the room–plenum exchange. Although these characteristics might arise in other spaces, we want to emphasise that this behaviour is not universal. In configurations where the ceiling zone has a larger effective volume or where persistent stratification leads to significant departures of the ceiling concentration from the bulk room concentration, the ceiling layer can support an additional slow mode. In such cases, the projection onto a single slow eigenmode may be inappropriate and a higher-dimensional reduced model (or even the full multi-zone description) may be required.

5.2. Generalisation to higher dimensional models

The model presented in § 3.1 can be generalised to analyse the dynamics of an arbitrary $n$ -zone configuration. In that case, the dynamics of the system would still be expressed as (3.9), but $\boldsymbol{C}$ would be a vector of size $n$ and $C_i$ would represent the CO $_2$ concentration in the $i$ th zone. Prior to any reduction of the system, an element $A_{ij}$ for $i\neq j$ of the $n\times n$ matrix $\boldsymbol{A}$ corresponds to a normalised ventilation into zone $i$ from zone $j$ , whereas the diagonal elements $A_{ii}\leq 0$ would represent the total normalised volume flux directed out of zone $i$ , as described by Parker and Bowman (Reference Parker and Bowman2011). Although $\boldsymbol{A}$ cannot, in general, be diagonalised because it might not have $n$ independent eigenvectors, it can be transformed into Jordan normal form. A trivial example of a configuration for which $\boldsymbol{A}$ cannot be diagonalised comes from $n$ rooms connected in series, such that fresh air is supplied to the first room in the series, while air from the last room discharges to a surrounding environment. In that case, the matrix $\boldsymbol{A}$ , when written in Jordan normal form, contains a single ‘Jordan block’, with non-zero entries along the diagonal (corresponding to the ventilation rate) and ones along the superdiagonal. The ventilation rate corresponds to the only eigenvalue of the system with geometric multiplicity equal to one (i.e. the $n$ eigenvectors of the system point in the same direction). More generally, systems consisting of multiple interacting flow loops will contain several Jordan blocks. In such cases, it would, in principle, be possible to apply the approach described in this paper by examining the eigenvalue and, therefore, time scale associated with each block and its corresponding eigenvector. Model reduction would then correspond to projecting the system onto the product of eigenspaces associated with slow dynamics.