Stochastic models of ventilation driven
by opposing wind and buoyancy

Veronica Andrian1 and John Craske1
1Department of Civil and Environmental Engineering,
Imperial College London, London SW7 2AZ, UK
(Updated )

## Abstract

> Stochastic versions of a classical model for natural ventilation are
> proposed and investigated to demonstrate the effect of random fluctuations
> on stability and predictability. In a stochastic context, the well-known
> deterministic result that ventilation driven by the competing effects of
> buoyancy and wind admits multiple steady states can be misleading, due to
> two distinct phenomena. First, with fluctuations in the buoyancy exchanged
> with an external environment, such systems eventually reside in the vicinity
> of global minima of their potential, rather than local minima associated
> with meta-stable equilibria. In the particular context of one heated space
> with a leeward low-level and windward high-level opening, a sustained
> buoyancy-driven flow opposing the wind direction is unlikely for wind
> strengths that exceed a statistically critical value, which is slightly
> larger than the critical value of the wind strength at which bifurcation in
> the deterministic system occurs. Second, fluctuations in the applied wind
> modify the topology of the system’s potential due to the nonlinear role that
> wind strength has in the equation for buoyancy conservation. Sufficiently
> large fluctuations in the wind rule out the possibility of ventilation
> opposing the wind direction at large base wind strengths. Although the
> phenomena described above might be perceived as making prediction easier,
> the results also highlight that certainty in the eventual state of the
> system goes hand in hand with uncertainty associated with longer transient
> effects. The work addresses growing interest in applying stochastic analysis
> to problems relating to building ventilation and urban fluid mechanics by
> describing a mathematically accessible example of the ‘stochasticisation’ of
> a canonical deterministic model, while highlighting the subtleties and
> challenges of developing stochastic models for ventilation in the future.

######

  1. 1 Introduction
  2. 2 Governing Equations
    1. 2.1 Deterministic Formulation
    2. 2.2 Stochastic (‘Brownian’) formulation
  3. 3 Analytical results
    1. 3.1 Coordinate transformation
    2. 3.2 Quadratic diffusion
    3. 3.3 Stationary density
    4. 3.4 The (statistically) critical wind strength
    5. 3.5 Kramer’s escape rate
  4. 4 Numerical simulations
  5. 5 An Ornstein-Uhlenbeck (‘Langevin’) model
  6. 6 Application
  7. 7 Conclusions
  8. 8 Acknowledgements

## 1 Introduction Introduction

By definition, buildings that are naturally ventilated are subjected to forces
from buoyancy and wind that fluctuate unpredictably in time [1, 2, 3].
Examples include the heat and motion from a variable number of occupants,
transient heat sources from equipment and solar gain, doors and windows that
are opened or closed intermittently, and the turbulence associated with
external wind loads [4, 5, 6, 7, 8, 9, 10, 11]. Indeed, the role played by
uncertainty has recently been highlighted more broadly in the context of the
succession events that lead to COVID-19 transmission indoors [12, 13]. From
the perspective of the deterministic macroscopic models that are typically
used to describe ventilation mathematically, the examples listed above are
random processes and must therefore be treated stochastically.

In spite of the uncertain environments to which they must be applied, most
existing ventilation models are deterministic and implicitly average over
fluctuations that would otherwise occur on relatively short time scales. Such
models include emptying filling box models, which have proven to be useful for
design and successfully capture the leading-order physics behind ventilation
[1]. Under controlled conditions, such models exhibit a good agreement with
experiments and field observations, and are a valuable design tool for
practitioners. Their elegance and simplicity belie the complex effects of
turbulence for which their parameters must necessarily account. Examples of
such parameters include the entrainment coefficient that is used to model
plumes in filling-box models [14, 15, 16] and discharge coefficients that are
used to account for dissipation across openings [17].

In contrast to the deterministic models described in the previous paragraph,
there is growing recognition that probabilistic methods are useful, if not
essential, in analysing and predicting the performance of buildings [18, 19].
For modelling occupancy [20], passive tracers [21] and wind speeds [11, 22],
such methods are already used. However, the incorporation of stochastic
processes in models for the fluid mechanics of natural ventilation is limited,
particularly in comparison with fields such as meteorology and climate
modelling [23].

Some stochastic effects in natural ventilation have previously been proposed
for an example involving buoyancy and wind [24], which motivates the present
study. The example consists of a single room with low and high level openings
subjected to internal heating that produces a uniform temperature and a wind
that opposes flow out of the high-level opening (see figure 1). By feeding
samples from probability distributions to account for occupancy, heating from
equipment and wind loads into equations describing a steady state, previous
work observed the effect of uncertainty in parameters on the structure of the
system’s solutions and their sensitivity [24]. However, the problem was not
formulated as a stochastic differential equation and analytical information
that can be obtained from the associated Fokker-Planck equation was not
utilised. More recently, stochastic fluctuations have been applied to the wind
opposing the displacement ventilation driven by a point source of buoyancy
[18]. Due to the nonlinear dependence of the ventilation on the opposing wind
strength, the incorporation of fluctuations in the wind modifies the interface
height associated with the resulting two-layer stratification.

The extent to which existing bulk models of ventilation implicitly account for
random fluctuations is dependent on the time scales under consideration and
therefore moot; while the effects of small scale turbulence and dissipation
are usually incorporated implicitly into discharge coefficients [17], slowly
increasing wind loads or changes in occupancy would typically appear
explicitly as a time-dependent forcing [25]. Being deterministic, they
describe a single trajectory or realisation, rather than a complete
description of the probabilities associated with a given observable, which
would correspond to infinitely more information about the process. Indeed,
deriving the probability densities associated with differential equations is
difficult in general, due to the ‘curse of dimensionality’ [26]. Each degree
of freedom of the system corresponds to an argument of the system’s
probability density, whose evolution is therefore governed by a partial
(rather than ordinary) differential equation [27]. Bulk models for
ventilation, however, are particularly amenable to probabilistic
interpretation because they typically involve a small number of dimensions.

A further motivation for adopting a probabilistic perspective comes from the
fact that ventilation driven by the combined effects of buoyancy and wind can
produce multiple equilibria (i.e. statistically steady states) [28, 29, 2]. In
contrast to systems with a single equilibrium, the eventual fate of systems
possessing multiple equilibria depends on their initial state. More generally,
in cases where a parameter can change, it is possible for a system to
transition between equilibria and exhibit hysteresis [30, 31]. In such cases
it is possible to identify the smallest perturbation in the wind that is
required to produce the transition for a single well-mixed space[32, 33] and
for the stratified environment created by a point source of buoyancy [34].
Systems subjected to random forces therefore have the propensity to flip
between meta-stable equilibria [35, 36, 37], as demonstrated by the classical
van der Pol oscillator subjected to white noise [38, 39] and ‘tipping points’
in climate models [40].

In a stochastic setting, the pictures of such bifurcations becomes blurred
[36]; the ‘skeleton’ provided by bifurcation diagrams is obscured by the
‘flesh’ of probability density functions that quantify the likelihood of
finding the system in a particular state and the parameter values at which
bifurcation occurs (in a sense that will be defined precisely in §2) are
modified [41]. Among the possible solutions identified in the previous
paragraph are some that are more probable than others, which means that
designers do not have to choose indiscriminately between possible outcomes.
Indeed, with density functions replacing deterministic values, the ‘flesh’ of
probability provides a rich and practical picture of the system’s behaviour,
which allows one to relate the amplification of input noise with the
underlying energetics of the system and to quantify the time one can expect to
wait for a system to transition between states [42, 43].

As a prototypical example for establishing a probabilistic picture of natural
ventilation driven by the competing effects of buoyancy and wind, this study
follows previous work[24, 44] and considers the single room with low and high
level openings described above and illustrated in figure 1. Previous
qualitative observations are made precise with analytical results and a more
rigorous formulation of the problem as a stochastic differential equation and
associated Fokker-Planck equation. The broader aim of the work is therefore to
demonstrate that useful information and insight can be obtained by
incorporating probabilistic effects into existing ventilation models using
accessible and well-established methods of stochastic analysis. At the same
time, we highlight various semantic challenges in developing a rigorous
foundation for stochastic models of ventilation.

The governing equations for the deterministic model are presented in §2.1 and
supplemented with a ‘Brownian’ stochastic term in §2.2, which is extended to a
‘Langevin’ formulation in §5 using an Ornstein-Uhlenbeck process. Analytical
results derived from the associated Fokker Planck equations are discussed in
§3 (and in §5 for the extension), and include predictions of the stationary
probability density. Numerical integration of the stochastic equations to
obtain realisations of the system and their comparison with the analytical
results are discussed in §4. In §6 the results are interpreted in the context
of a real application and conclusions are drawn in §7.

Figure 1: A space with lower and upper level openings, subjected to a
nondimensionalised distributed heat source which produces (by assumption) a
well-mixed interior of uniform buoyancy. The buoyancy-driven ventilation is
opposed by a pressure difference across each opening due to an opposing wind
of (nondimensionalised) strength $W$. The left-hand figure illustrates
‘forward flow’ and the right-hand figure illustrates ‘reverse flow’ (the
latter being possible when $W\geq 1$).

## 2 Governing Equations Governing Equations

### 2.1  _Deterministic Formulation_ Deterministic Formulation

Consider the confined space illustrated in figure 1 with an upper and lower
level openings subjected to a distributed heat source at floor level of
strength $\hat{F}$. The temperature (and, therefore, buoyancy $\hat{b}$) of
the air inside the space is assumed to be uniform [44, 25]. Using the
Boussinesq approximation, the system is governed by the following buoyancy
conservation equation:

$\hat{H}\hat{S}\dfrac{\mathrm{d}\hat{b}}{\mathrm{d}\hat{t}}=\hat{F}-\hat{Q}\hat{b}=\hat{F}-\hat{A}\bigg{\lvert}\hat{b}\hat{H}-\frac{\Delta\hat{p}}{\hat{\rho}_{0}}\bigg{\rvert}^{1/2}\hat{b},$
(1)

where $\hat{H}$ is the height of the space, $\hat{S}$ is the horizontal
(cross-sectional) area, $\Delta\hat{p}$ is the pressure difference induced by
the wind, $\hat{A}$ is an effective opening area [1] and $\hat{\rho}_{0}$ is a
reference density. The right-hand side of (1) represents the difference
between the buoyancy that is added to the space from heating and the buoyancy
that is lost from the space due to ventilation, which itself depends on the
pressure difference arising from both buoyancy and wind, and therefore makes
the equation nonlinear.

When the wind-induced pressure difference is below the critical value [44, 2,
28],

$\frac{\Delta\hat{p}_{c}}{\hat{\rho}_{0}}:=\left(\frac{27\hat{H}^{2}\hat{F}^{2}}{4\hat{A}^{2}}\right)^{1/3},$
(2)

equation (1) admits a single stable steady-state solution (for which
$\,\mathrm{d}\hat{b}/\,\mathrm{d}\hat{t}=0$) corresponding to forward flow.
When the pressure difference is greater than the critical value, (1) admits
two additional steady-state solutions: a stable and an unstable solution
corresponding to reverse flow in the direction of the wind (see figure 1). In
the parlance of dynamical systems theory, the steady-state solutions are known
as ‘fixed points’ of (1).

By introducing the following dimensionless quantities,

$b:=\left(\frac{4\hat{A}^{2}\hat{H}}{27\hat{F}^{2}}\right)^{1/3}\hat{b},\quad\quad
t:=\left(\frac{4\hat{A}^{2}\hat{F}}{27\hat{H}^{2}\hat{S}^{3}}\right)^{1/3}\hat{t},\quad\quad
W:=\frac{\Delta\hat{p}}{\Delta\hat{p}_{c}},$ (3)

equation (1) can be written as

$\dfrac{\mathrm{d}b}{\mathrm{d}t}=a(b;W):=1-c\lvert b-W\rvert^{1/2}\,b,$ (4)

where $c=\sqrt{27/4}$ and $W$, quantifying the relative strength of the wind
as a squared Froude number, is the only free parameter.

The fixed points of (4) correspond to either local minima (if they are stable)
or local maxima (if they are unstable) of the potential function

$V_{B}(b;W):=-b+\frac{2c}{15}(3b+2W)(b-W)\lvert{b-W}\rvert^{1/2},$ (5)

which is defined such that

$\dfrac{\mathrm{d}b}{\mathrm{d}t}=-\dfrac{\mathrm{d}V_{B}}{\mathrm{d}b}.$ (6)

Algebraically, the fixed points of (4) correspond to solutions to the cubic
equation $-\,\mathrm{d}V_{B}/\,\mathrm{d}b=a(b;W)=0$:

$b^{3}-Wb^{2}=\pm\frac{1}{c^{2}},$ (7)

where the right-hand side is positive for forward flow (one real solution,
$b_{3}$) and, when $W\geq 1$, negative for reverse flow (two real solutions,
$b_{1}$ and $b_{2}$, where $b_{1}<b_{2}<b_{3}$). Contours of $V_{B}(b;W)$ are
shown beneath the system’s bifurcation diagram in figure 2. The set of states
that eventually end up at a given fixed point is called the ‘basin of
attraction’ of that fixed point, or, ‘potential well’ in the context the
energetics implied by $V_{B}$.

Figure 2: Bifurcation diagram corresponding to solutions of $a(b;W)=0$ and
underlying potential function $V_{B}(b;W)$ for the deterministic system (4).
The thick solid lines denote stable equilibria and the dashed line denotes
unstable equilibria.

### 2.2  _Stochastic (‘Brownian’) formulation_ Stochastic (‘Brownian’)
formulation

There are many ways in which stochastic effects can be incorporated in (4).
However, in all cases it is necessary to reappraise the meaning of the
original deterministic equation. Does (4) represent an idealised system devoid
of disturbances or does it represent a system that has already been averaged
over sufficiently large time scales? Added to these semantic difficulties of
developing a stochastic model for ventilation are the choices that must be
made concerning the physical origin of the fluctuations. While large scale
fluctuations in the wind would be expected to change the bulk pressure
difference acting on the building, smaller-scale structures could produce
bidirectional flows through the openings that enhance buoyancy transfer
without necessarily changing the bulk ventilation. In this regard it should
also be noted that, in reality, the interior space is not of uniform buoyancy
and that temporal and spatial variations of buoyancy, and/or the heat applied
to the space, might also be treated as random variables.

If random noise is to be added to (4) to create a stochastic differential
equation, which would provide trajectories for individual realisations of the
system, choices concerning the rule used for time integration of the random
fluctuations present further technical and interpretative challenges [45].
Such ‘choices’ are discussed in §7 and disappear altogether with sufficient
information about the physical properties of the noise and/or a master
equation describing the evolution of the system’s probability density.

In view of the challenges outlined above and the current lack of fundamental
justification for adopting any particular mathematical form of the
fluctuations that could affect (4), the approach adopted here will be as
simple and transparent as possible. To account for the effects of rapid,
uncertain and external fluctuations applied to the system, a stochastic
buoyancy flux will be added to (4):

$\,\mathrm{d}B_{t}=a(B_{t};W)\,\mathrm{d}t+\sigma(B_{t};W)\circ\,\mathrm{d}\xi_{t},$
(8)

where $B_{t}$ is the random variable associated with buoyancy, $a(B_{t};W)$ is
defined in (4) and $\xi_{t}$ describes a Wiener process characterised by
normally-distributed jumps that are uncorrelated in time. More precisely,
$\xi_{t}$ has independent increments that are distributed, such that
$\xi_{t}-\xi_{s}\sim N(0,t-s)$ for $0\leq s<t$, and describes continuous
paths, with $\xi_{0}=0$ [46]. The standard deviation of the fluctuations
$\sigma(B_{t};W)$ appearing in (4) can be adapted in accordance with the
physical origin of the fluctuations.

The Wiener process in (8) is integrated using the Stratonovich interpretation
[45, 47], which, informally, means that $\sigma$ is evaluated at the mid-point
of each time step (see §4 for further details). The Stratonovich
interpretation used in (8) represents an external source of fluctuations in
the limit of its autocorrelation time tending to zero [45]. This, and other
subtleties associated with the incorporation of stochastic effects in (4) are
discussed further in §7.

Physically, the properties of $\xi_{t}$ correspond to the fluctuations
belonging to a normal distribution and occurring on time scales that are
significantly smaller than the integral timescales associated with the
deterministic version of the model (4). As demonstrated in §5, as an extension
of the more tractable model considered here, finite time correlations of
fluctuations in the wind can be incorporated by modelling $W$ in (4) as an
Ornstein-Uhlenbeck process with an autocorrelation that decays exponentially
in time [48, 49, 18].

Given (8), the probability density $f_{B}$ describes the likelihood of
observing a given buoyancy $B_{t}$, such that the probability

$\mathbbm{P}\\{0\leq B_{t}\leq
b\\}=\int\limits_{0}^{b}f_{B}(\beta;t)\,\mathrm{d}\beta.$ (9)

The probability density evolves according to the Fokker-Planck equation [50]:

$\dfrac{\partial f_{B}}{\partial t}=-\dfrac{\partial J_{B}}{\partial b},$ (10)

where

$J_{B}:=a(b;W)f_{B}-\frac{\sigma}{2}\dfrac{\partial}{\partial b}(\sigma
f_{B}),$ (11)

is a probability flux that satisfies $J_{B}(0;t)=J_{B}(\infty;t)=0\ \forall
t$, which ensures that probability is conserved in the domain, in the sense
that $\mathbbm{P}\\{0\leq B_{t}<\infty\\}=1\ \forall t$ for suitable choices
of the initial density. Indeed, equation (10) is a ‘mass conservation’
equation for probability and therefore accounts for advection along the
trajectory defined by (4) (typically referred to as ‘drift’), in addition to
diffusion due to the stochastic forcing in (8). The second term in the
probability flux (11) is diffusion that results from the Stratonovich
integration of the noise in (8). Integration of (10) against $b$ and commuting
$\partial_{b}$ and $\sigma(b;W)$ using the product rule in (11) shows that the
Stratonovich diffusion term modifies the evolution of the mean buoyancy.

The solution of (10) yields the entire probability density function for
$B_{t}$ and corresponds to Einstein’s original formulation of Brownian motion
[51]. ‘Brownian’ in this regard, and in the title of this section, refers to
the application of a stochastic force directly to the ‘configuration space’
corresponding to $B_{t}$. In §5, on the other hand, fluctuations are applied
to an additional ‘velocity space’, which is analogous to Langevin’s subsequent
formulation [52]. A further, independent, difference between the two
approaches is that while Einstein presented his model in terms of a Fokker-
Planck equation (cf. (10)), Langevin manipulated the stochastic differential
equations describing particle trajectories directly [53].

## 3 Analytical results Analytical results

Equation (10) is a one-dimensional advection diffusion equation and has
analytical solutions that can be expressed in terms of a potential function
that is scaled to account for the state-dependent variance
$\sigma(B_{t};W)^{2}$ of the destabilising random fluctuations. To derive the
solution it is first convenient to derive a potential function in terms of a
new random variable $X_{t}$, in terms of which the variance $\sigma^{2}$
becomes independent of state and, therefore, constant.

### 3.1  _Coordinate transformation_ Coordinate transformation

The transformation $g:b\mapsto x$ that produces a Fokker-Planck equation with
a constant diffusivity $\varepsilon^{2}/2$ in the new frame of reference is
[50, 31]:

$g(b):=\int_{0}\limits^{b}\frac{\varepsilon}{\sigma(\beta;W)}\,\mathrm{d}\beta,$
(12)

and the probability density or pushforward measure associated with the new
variable $x$, assuming that $g$ is differentiable and invertible, satisfies

$f_{X}(x)=f_{B}(g^{-1}(x))\dfrac{\mathrm{d}}{\mathrm{d}x}g^{-1}(x),$ (13)

or, in the other direction,

$f_{B}(b)=f_{X}(g(b))\dfrac{\mathrm{d}}{\mathrm{d}b}g(b).$ (14)

Unlike $f_{B}$, $f_{X}$ satisfies a Fokker-Planck equation with constant
diffusivity $\varepsilon^{2}/2$,

$\dfrac{\partial f_{X}}{\partial t}=-\dfrac{\partial J_{X}}{\partial x},$ (15)

where

$J_{X}=-\left(\dfrac{\mathrm{d}V_{X}}{\mathrm{d}x}+\frac{\varepsilon^{2}}{2}\dfrac{\partial}{\partial
x}\right)f_{X},$ (16)

for which the potential function (noting that $b$ must be obtained by
inverting (12)) is defined to within an arbitrary additive function $V_{0}$ of
$W$ according to

$V_{X}:=-\int\limits^{x}\frac{\varepsilon}{\sigma}\left[1-c\lvert{g^{-1}(x)-W}\rvert^{1/2}g^{-1}(x))\right]\,\mathrm{d}x+V_{0}(W).$
(17)

Equation (15) adopts a canonical frame of reference for obtaining analytical
solutions and is the setting of the classical stochastic problem of
determining the escape rate from a potential well [42], which will be examined
in §3.5. The formulation shows immediately that local minima and maxima of
$V_{X}$ correspond to zeros of the term in square brackets in (17) and
therefore coincide with the fixed points of the original system (4), which,
incidentally, would not have been the case if the stochastic term in (8) had
been integrated according to the Itô interpretation.

### 3.2  _Quadratic diffusion_ Quadratic diffusion

As a specific example of the model presented in the previous section, assume
that $\sigma(b;W)$ in (8) characterises an ‘external’ fluctuation in the
ventilation, such that the standard deviation of the forcing applied to (8)
(representing the strength of the resulting fluctuations in the buoyancy flux)
can be expressed as

$\sigma(b;W)=\varepsilon(W)b.$ (18)

For windward and leeward openings perpendicular to the wind direction, one
would expect $\varepsilon$ to increase with respect to $W$. However, for
oblique orientations it is conceivable that large fluctuations could occur at
relatively small values of $W$. Noting these possibilities, the results below
will be presented for constant $\varepsilon$ but can be adapted immediately
for particular relationships between $W$ and $\varepsilon$.

According to (17), a potential function associated with (18) can be expressed
as

$V_{X}(x):=e^{-x}+2cW^{1/2}\left(\lvert{e^{x}/W-1}\rvert^{1/2}-\tau^{-1}\lvert{e^{x}/W-1}\rvert^{1/2}\right)+V_{0}(W).$
(19)

where $e^{x}=:g^{-1}(x)=b$ and $\tau$ is a function that depends on whether
the flow is forward (first case) or reverse (second case):

$\tau:=\begin{cases}\mathrm{tan},\quad\quad&e^{x}=b\geq W,\\\
\mathrm{tanh},\quad\quad&e^{x}=b<W.\end{cases}$ (20)

Although the constant $V_{0}(W)$ is arbitrary and therefore does not affect
the dynamics induced by $V_{X}$, selecting $V_{0}(W)$ to ensure that

$\min\limits_{x}V_{x}(x;W)=0,$ (21)

is useful for the numerical conditioning of expressions involving $V_{X}$
below (particularly (23)). For plotting $V_{X}$, however, it will often be
convenient to adopt $V_{0}(W)\equiv 0$ (cf. (5) and figure 2). The potential
(19) is illustrated in figure 3 (with $V_{0}(W)\equiv 0$), along with
parabolic approximations in the vicinity of each local minima that will be
discussed in §3.3.

The local minima of $V_{X}$ correspond to the stable fixed points $b_{1}$ and
$b_{3}$ (denoted $x_{1}$ and $x_{3}$ in the $x$-domain, respectively) of (4).
Use of (18) implies that $V_{X}\rightarrow\infty$ as $b\rightarrow 0$
($x\rightarrow-\infty$), due to the first term on the right-hand side of (19)
and $V_{X}\rightarrow\infty$ as $b\rightarrow\infty$ ($x\rightarrow\infty$),
due to the remaining terms. Consequently, as will be seen in the next section,
the value and gradient of the probability density in these limits tends to
zero and therefore satisfies the boundary conditions for the buoyancy flux
discussed below (11).

Figure 3: Parabolic approximations corresponding to the first two terms of
(25) (dashed lines) to the exact potential (19) (solid line) at
$W=W_{s}=1.0887$, which is defined in §3.4 when $\sigma(b;W)=\varepsilon b$.

### 3.3  _Stationary density_ Stationary density

In a statistically stationary state $\partial_{t}f_{X}\equiv 0$, which, given
the boundary conditions $J_{X}(0)=J_{X}(\infty)=0$, implies a ‘detailed
balance’ [50] in which the probability flux vanishes everywhere:

$\dfrac{\partial f_{X}}{\partial
x}=-\frac{2}{\varepsilon^{2}}\dfrac{\mathrm{d}V_{X}}{\mathrm{d}x}f_{X},$ (22)

and, therefore,

$f_{X}(x;W)=N\,\mathrm{exp}\left(-\frac{2V_{X}(x;W)}{\varepsilon^{2}}\right),$
(23)

for which $N$ is the choice of integration constant that normalises $f_{X}$
(i.e. ensures that it is a probability density whose integral with respect to
$x\in(-\infty,\infty)$ is equal to $1$) and $f_{B}$ is obtained from (13).
Equation (23) is equivalent to the Boltzmann distribution in statistical
mechanics [54]. As can be seen from (23), since $\varepsilon$ is a constant,
local minima and maxima in $f_{X}$ correspond to local maxima and minima in
$V_{X}$, respectively, although the same cannot be said about the relationship
between $f_{B}$ and $V_{B}$ due to the pushforward transformation rule (13).

Figure 4 illustrates the dependence of the stationary probability density (23)
on the buoyancy $b$ and wind strength $W$ (in this regard, $f_{B}$ should be
viewed as a probability density conditional on a given value of $W$). As the
strength of the noise decreases, the probability density becomes concentrated
in narrow bands around the deterministic system’s bifurcation diagram. At a
particular value of the wind strength, which will be denoted $W_{s}$ and
defined precisely in §3.4, a transition in the most likely state of the system
occurs. When $W<W_{s}$ the system is predominantly found in forward flow,
whereas when $W>W_{s}$ the system is predominantly found in reverse flow.
Rather than ascribing equal probability to the forward and reverse branches of
the bifurcation diagram, the stationary probability density accounts for the
tendency of the system to minimise the underlying potential (19). In addition,
a crucial feature of the stationary probability density in the limit
$\varepsilon\rightarrow 0$ is that the point $W=W_{s}>1$, at which the system
transitions from forward to reverse flow, does not correspond to the point of
bifurcation in the underlying deterministic system.

Figure 4: The stationary probability density (shades of red) as a function of
the wind strength $W$ when $\sigma(b;W)=\varepsilon b$ for $\varepsilon=0.5\
(a)$, $\varepsilon=0.25\ (b)$ and $\varepsilon=0.125\ (c)$. As demonstrated in
§3.4, when $\varepsilon\rightarrow 0$ the system resides in a state of forward
flow for $W<W_{s}$ and reverse flow for $W_{s}<W$. The solid grey line
corresponds to $b=W$ and therefore separates points corresponding to forward
flow (above) and reverse flow (below).

In the vicinity of a local minima at $x=x_{m},\ m\in\\{1,3\\}$, of the
potential $V_{X}$, it is useful to scale the excursions $\Delta
x_{m}:=x-x_{m}$ by the magnitude of the noise and the restoring force provided
by the potential well:

$\Delta y_{m}:=\sqrt{V_{X}^{\prime\prime}(x_{m})}\frac{\Delta
x_{m}}{\varepsilon}.$ (24)

The definition of $\Delta y_{m}$ allows the potential function to be expressed
as

$V_{X}(x)\sim V_{X}(x_{m})+\frac{\varepsilon^{2}}{2}\Delta
y_{m}^{2}+\frac{\varepsilon^{3}\Delta
y_{m}^{3}}{3!}\frac{V_{X}^{\prime\prime\prime}(x_{m})}{V_{X}^{\prime\prime}(x_{m})^{3/2}}+\ldots,$
(25)

from which the term of order $\varepsilon$ is missing because the first
derivative $V_{X}^{\prime}$ is zero at the local minima $x_{m}$ by definition.

Under the assumption that

$\dfrac{\mathrm{d}^{n}V_{X}}{\mathrm{d}x^{n}}\left(\dfrac{\mathrm{d}^{2}V_{X}}{\mathrm{d}x^{2}}\right)^{-\frac{n}{2}}=o\left(\varepsilon^{2-n}\right)\quad\mathrm{as\
}\varepsilon\rightarrow 0\ \ \mathrm{at\ \ }x=x_{m}\ \ \mathrm{for\ \
}n=3,4,\ldots,$ (26)

all but the first two terms in the expansion (25) can be discarded as being
relatively small in the limit $\varepsilon\rightarrow 0$, leading to the
approximation of the potential function as a set of parabolas. Locally, the
expansion therefore produces Gaussian probability densities with respect to
$x$, conditional on the chosen minimum $m$:

$f_{X|M}(x;m)\sim\sqrt{\frac{V^{\prime\prime}(x_{m})}{\varepsilon^{2}\pi}}\mathrm{exp}\left(-\Delta
y_{m}(x)^{2}\right).$ (27)

The contribution that the local density (27) makes to $f_{B}$ can be
calculated by applying the transformation (14) and scaling the result by the
probability of observing $B_{t}\in\mathscr{B}_{m}$ for each $m$:

$f_{B}(b)\sim\sum\limits_{m=1,3}\mathbbm{P}\\{B_{t}\in\mathscr{B}_{m}\\}f_{X|M}(g(b);m)\dfrac{\mathrm{d}}{\mathrm{d}b}g(b),$
(28)

where $\mathscr{B}_{m}$ is the basin of attraction associated with the stable
fixed point at $b_{m}$ and $\mathbbm{P}\\{B_{t}\in\mathscr{B}_{m}\\}$ is the
corresponding approximation to the probability of finding the statistically
stationary buoyancy in the basin of attraction associated with the fixed point
$m$ as $\varepsilon\rightarrow 0$:

$\mathbbm{P}\\{B_{t}\in\mathscr{B}_{m}\\}\sim
N\sqrt{\frac{\varepsilon^{2}\pi}{V_{X}^{\prime\prime}(x_{m})}}\mathrm{exp}\left(-2\frac{V_{X}(x_{m})}{\varepsilon^{2}}\right).$
(29)

In this regard, the normalisation constant $N$ should be chosen to ensure that
$\mathbbm{P}\\{B_{t}\in\mathscr{B}_{1}\\}+\mathbbm{P}\\{B_{t}\in\mathscr{B}_{3}\\}=1$.
In order to satisfy the zero probability flux boundary condition at $b=0$
exactly, an image solution located at $b=-b_{m}=-g^{-1}(x_{m})$ could be added
to (28). However, since $\varepsilon$ is required to be small in order for
(28) to be valid, the error that would otherwise be induced by (28) at the
boundary $b=0$ will be exponentially small.

Figure 5 displays the stationary probability density $f_{B}$ and the
approximation (28) for different values of $\varepsilon$ at a value of the
wind strength just below $W_{s}$. When $\varepsilon=0.5$, figure 5 suggests
that the typical excursions of the system (characterised by the ratio of
$\varepsilon$ and the restoring force provided by the curvature of a potential
well) are at least comparable to the size of each basin of attraction. Indeed,
evaluation of $V_{X}^{\prime\prime}$ implies that when $\varepsilon=0.5$,
$\varepsilon/\sqrt{V_{X}^{\prime\prime}(x_{1})}\approx 0.5$, which is close to
the distance between the local minimum for reverse flow and the local maximum
(see figure 3). Consequently, the probability density in figure 5$(a)$,
although bimodal, suggests that the system flips regularly between forward and
reverse flow when $\varepsilon=0.5$. In such situations the approximation (28)
is inaccurate for certain values of the buoyancy ($b<0.3$, for example) but
does capture the density’s general shape. As $\varepsilon$ is reduced to
$0.125$, the probably density associated with the unstable fixed point (the
local maximum of $V_{X}$) approaches zero and the approximation (28) becomes
almost indistinguishable from the stationary probability density (23), which
consists of two disjoint densities centred on the potential wells.

Figure 5: The stationary probability density at $W=W_{s}-0.01$ (which means
that $\mathbbm{P}\\{b\in\mathscr{B}_{2}\\}\rightarrow 1$ as
$\varepsilon\rightarrow 0$) when $\sigma(b;W)=\varepsilon b$ for
$\varepsilon=0.5\ (a)$, $\varepsilon=0.25\ (b)$ and $\varepsilon=0.125\ (c)$.
The thick red line denotes the exact solution (23) and the blue dashed line
denotes the Gaussian approximation that becomes increasingly accurate as
$\varepsilon\rightarrow 0$. The vertical dashed lines separate reverse flow
($b<W$) from forward flow ($W<b$). The vertical black lines indicate the
location of the stable fixed points for forward flow and reverse flow.

### 3.4  _The (statistically) critical wind strength_ The (statistically)
critical wind strength

Figures 2 and 4 suggest that the asymmetry of the system’s underlying
potential biases the probability density that describes its eventual state. In
particular, as $\varepsilon\rightarrow 0$, the exponent on the right-hand side
of (23) is dominated by the smallest (most negative) value of $V_{X}$ and
converges weakly to a dirac measure. At small values of $W$, the deepest
potential well corresponds to forward flow, whereas for large values of $W$
the deepest potential well corresponds to reverse flow. Between these two
regimes, for some intermediate value $W_{s}>1$, the depths of the potential
wells are equal and the dirac measure, to which the distributions with
$\varepsilon\rightarrow 0$ converge, switches from forward flow to reverse
flow.

Stated differently, if one is prepared to wait long enough, the state of the
system subjected to asymptotically vanishing noise will be found in the
deepest potential well. It is therefore important to know the value $W_{s}$
that renders the two potential wells of equal depth, since such a value will
correspond to the transition in the expected state from forward flow to
reverse flow.

The required definition of $W_{s}$ is:

$\lim\limits_{\varepsilon\rightarrow
0}\mathbbm{P}\\{B_{t}\in\mathscr{B}_{1}(W)\\}=\begin{cases}0,&\quad\quad
W<W_{s},\\\ 1,&\quad\quad W>W_{s},\end{cases}$ (30)

which, regarding $\mathbbm{P}$ as referring to the stationary probability,
states that the statistically stationary system will almost always be in the
basin of attraction for forward or reverse flow as $\varepsilon\rightarrow 0$
when $W<W_{s}$ and $W>W_{s}$, respectively. The ‘statistically’ critical value
of the wind $W_{s}=1.0887$ makes the depth of the potential wells for reverse
flow and forward flow equal, and therefore satisfies

$V_{X}(x_{1}(W_{s}))=V_{X}(x_{3}(W_{s})),$ (31)

In spite of the fact that forward flow is a stable fixed point of the
underlying deterministic system, from the perspective of the system’s
energetics, and in the context of practical applications, one should expect
reverse flow to eventually dominate the system’s response for sufficiently
large opposing wind strengths. Globally, the fixed point associated with
forward flow for $W>W_{s}$ corresponds to a point of meta-stable equilibrium,
which, since the stochastic forcing does not have compact support, is not
stable enough to prevent an eventual transition to reverse flow. The next
section quantifies how long one would typically need to wait for this
‘eventual’ transition between forward flow and reverse flow to occur.

Figure 6: Steady-state probability $\mathbbm{P}\\{B_{t}\in\mathscr{B}_{1}\\}$
of finding the state in the basin of attraction for reverse flow as a function
of $W$ and $\varepsilon$ when $\sigma(b;W)=\varepsilon b$. At the top of the
figure are shown stationary probability densities for $\varepsilon=0.1$ at
various wind strengths.

For practical applications it is useful to identify combinations of $W$ and
$\varepsilon$ that lead to confidence in whether the statistically steady
state corresponds to forward or reverse flow. To this end, figure 6 plots
isolines of the stationary probability
$\mathbbm{P}\\{B_{t}\in\mathscr{B}_{1}\\}$, which is the probability of
finding the system in the basin of attraction for reverse flow. Because the
line $b=W$ demarcating forward and reverse flow (see, for example, figure 4)
does not coincide with the local maxima of $V_{B}$, reverse flow is
theoretically possible even when the system is in the basin of attraction for
forward flow. Consequently, probabilities corresponding to the likelihood of
being in the basin of attraction for reverse flow underestimate (by an amount
that diminishes with $\varepsilon$) the probability of forward flow, i.e.
$\mathbbm{P}\\{B_{t}\in\mathscr{B}_{1}\\}<\mathbbm{P}\\{B_{t}<W\\}$.

The contours in figure 6 are characterised by a fan shaped region of
uncertainty, in which the probability varies continuously between $0$ and $1$,
emanating from $(W,\varepsilon)=(W_{s},0)$. On either side of this region of
uncertainty are regions for which one can be close to certain that the system
will eventually be in either the basin of attraction for forward or reverse
flow. Plots such as (6) might provide a useful guide to engineers by
identifying conditions for which it would be difficult or inappropriate to
design for a single flow regime.

Figure 7: The logarithm of the inverse of Kramer’s escape rate $\log(1/r)$
from (33) as a function of the wind strength $W$ and fluctuation amplitude
$\varepsilon$ when $\sigma(b;W)=\varepsilon b$. The blue dashed line is the
contour for which $r=1$.

### 3.5  _Kramer’s escape rate_ Kramer’s escape rate

The stationary densities discussed in the previous section are solutions to
$\partial_{t}f_{B}=0$. How long one would have to wait in practice to see such
solutions, however, would depend on the initial state of the system. For
example, for small values of $\varepsilon$, wind strengths $W$ slightly larger
than $W_{s}$ and an initial density in the region $\mathscr{B}_{3}$ (forward
flow), one might have to wait a long time for the density to reach
$\mathscr{B}_{1}$ (reverse flow) in the way that figure 6 suggests.

The remarks above can be made precise by considering a potential barrier
$\Delta V_{X}$, quantified as the difference between the local maxima
$V_{X}(x_{2})$ and local minima $V_{X}(x_{3})<V_{X}(x_{2})$ associated with
forward flow,

$\Delta V_{X}:=V_{X}(x_{2}(W))-V_{X}(x_{3}(W))\geq 0.$ (32)

If the potential barrier is large relative to the diffusivity $\varepsilon$,
the rate $r(W,\varepsilon)$ at which states ‘escape’ from $\mathscr{B}_{3}$
(forward flow), quantified as the probability current divided by the
probability, can be approximated using Kramer’s escape rate formula [50]:

$r(W,\varepsilon)=\frac{1}{2\pi}\sqrt{V_{X}^{\prime\prime}(x_{3})\lvert
V_{X}^{\prime\prime}(x_{2})\rvert}\,\mathrm{exp}\left(-\frac{2\Delta
V_{X}}{\varepsilon^{2}}\right).$ (33)

The associated timescale, $1/r$, is plotted in figure 7 for the quadratic
diffusion model $\sigma(b;W)=\varepsilon b$. The figure shows that the time
one can expect to wait to reach a stationary state decreases with respect to
$W$, which controls the height of the potential barrier and increases with
respect to $\varepsilon$, which quantifies the ‘destabilising’ noise. Since
$V_{X}$ does not depend on $\varepsilon$, (33) suggests that the time $1/r$
that one would have to wait to reach a stationary density scales exponentially
with respect to $1/\varepsilon^{2}$.

Figure 8: Realisations of the system using (34) for $W=W_{s}-0.01$,
$\varepsilon=0.25\ (a)$ and $\varepsilon=0.5\ (b)$. The panels on the left
show the probability density $f_{B}$ collected from simulation data (blue
filled region) and compared with the analytical solution (23) (red line),
whereas the panels on the right show a single realisation of buoyancy $B_{t}$.
The dashed line corresponds to $b=W$ and therefore indicates the threshold
separating forward flow ($b>W$) from reverse flow $(b<W)$. The horizontal
black lines denote the location of the fixed points for forward flow and
reverse flow.

## 4 Numerical simulations Numerical simulations

To obtain individual realisations of the stochastic model, (8) can be
integrated directly. Using the Heun method [38], which accounts for the drift
induced by the Stratonovich interpretation of the noise without modification
of $a$ in (8), the discretised equations are

$B_{t_{n+1}}=B_{t_{n}}+\overline{a}(B_{t_{n}};W)\Delta
t_{n}+\overline{\sigma}(B_{t_{n}};W)\Delta\xi_{t_{n}},$ (34)

where $\Delta t_{n}:=t_{n+1}-t_{n}$,
$\Delta\xi_{t_{n}}:=\xi_{t_{n+1}}-\xi_{t_{n}}\sim N(0,\Delta t_{n})$ and

$\displaystyle\overline{a}(B_{t_{n}};W)$
$\displaystyle:=\frac{1}{2}\left[a(\overline{B}_{t_{n+1}};W)+a(B_{t_{n}};W)\right],$
(35) $\displaystyle\overline{\sigma}(B_{t_{n}};W)$
$\displaystyle:=\frac{1}{2}\left[\sigma(\overline{B}_{t_{n+1}};W)+\sigma(B_{t_{n}};W)\right],$
(36)

are intermediate values of the drift and standard deviation based on the
following estimate of the buoyancy at time $t_{n+1}$:

$\overline{B}_{t_{n+1}}:=B_{t_{n}}+a(B_{t_{n}};W)\Delta
t_{n}+\sigma(B_{t_{n}};W)\Delta\xi_{t_{n}}.$ (37)

The sample trajectories shown in figures 8 and 9 were obtained with $\Delta
t_{n}=0.05$ for $t_{n}\in[0,4000]$. A set of $100$ such realisations were
obtained for a given value of $W$ to construct the histograms shown on the
left-hand side of figures 8 and 9. The fixed point for forward flow was used
as an initial condition for all realisations and statistically transient data
for $t_{n}\in[0,100]$ were not included in the construction of the histograms.

Figure 8 shows realisations for $W=W_{s}-0.01$ (cf. figure 5, for which the
probability density would concentrate itself around the fixed point for
forward flow in the limit $\varepsilon\rightarrow 0$). In both cases
($\varepsilon=0.25$ and $\varepsilon=0.5$) the system’s trajectory crosses the
potential barrier frequently, which leads to large variance and significant
uncertainty regarding the direction of the ventilation. In figure 9 on the
other hand, where $W=1.2$, the system resides predominantly in reverse flow
for the reasons described in previous sections. Whereas $\varepsilon=0.5$
produces occasional excursions of the system into a state of forward flow,
$\varepsilon=0.25$ suggests that the flow is almost entirely in the reverse
direction.

The histograms shown in figures 8 and 9 exhibit a reasonably good agreement
with the analytical solution (23). The relatively small differences between
the analytical solution and the histograms obtained from realisations of (8)
in the vicinity $b=1$ in figure 8 illustrate the challenge of obtaining a
sufficient number of samples from direct simulation. Bearing in mind that the
governing equations are nondimensional, figure 8 shows that one needs to
collect data over extremely long durations/large ensembles in order to resolve
the underlying density, which makes the availability of analytical solutions
or approximations particularly attractive.

Figure 9: Realisations of the system using (34) for $W=1.2$, and
$\varepsilon=0.25\ (a)$, $\varepsilon=0.5\ (b)$. The panels on the left show
the probability density $f_{B}$ collected from simulation data (blue filled
region) and compared with the analytical solution (23) (red line), whereas the
panels on the right show a single realisation of buoyancy $B_{t}$. The dashed
line corresponds to $b=W$ and therefore indicates the threshold separating
forward flow ($b>W$) from reverse flow $(b<W)$. The horizontal black lines
denote the location of the fixed points for forward flow and reverse flow.
Figure 10: Contours of the potential $V_{B}^{*}$ modified by a stochastic wind
strength using an adiabatic approximation for $\gamma/\nu^{2}=1000$ $(a)$,
$\gamma/\nu^{2}=100$ $(b)$ and $\gamma/\nu^{2}=10$ $(c)$. The thick line
corresponds to (stable) local minima of the potential while the thick dashed
line corresponds to (unstable) local maxima of the potential. The dash-dot
line corresponds to $b=W$ and therefore demarcates forward flow from reverse
flow.

## 5 An Ornstein-Uhlenbeck (‘Langevin’) model An Ornstein-Uhlenbeck
(‘Langevin’) model

As discussed in §2.2, the effects of fluctuations can be incorporated in
several ways, which raises important questions concerning their physical
origin and the correct interpretation of the deterministic equations to which
they are added. One particular approach is to represent fluctuations in the
bulk pressure difference ($\Delta\hat{p}$ in §2.1) as an Ornstein-Uhlenbeck
process [18], which can account for temporal correlations. If the stochastic
buoyancy flux from §2.1 is retained, the governing equations become

$\displaystyle\,\mathrm{d}B_{t}$
$\displaystyle=a(B_{t},U_{t};W)\,\mathrm{d}t+\sigma(B_{t};W)\circ\,\mathrm{d}\xi_{t},$
(38) $\displaystyle\,\mathrm{d}U_{t}$ $\displaystyle=-\gamma
U_{t}+\nu\,\mathrm{d}\zeta_{t},$ (39)

where, following [18], $U_{t}$ corresponds to dimensionless fluctuations in
the wind velocity that are proportional to the square root of the resulting
pressure difference:

$a(B_{t},U_{t};W):=1-c\lvert B_{t}-(1+U_{t}^{2})W\rvert^{1/2}\,B_{t}.$ (40)

The parameter $\gamma$ determines the rate at which $U_{t}$ reverts to its
mean value of zero and $\nu\,\mathrm{d}\zeta_{t}$ denotes uncorrelated
increments of a zero-mean Gaussian process with time- and state-invariant
standard deviation $\nu$. One difference between (38)-(39) and the simpler
version studied in §2.2 is that the stochastic forcing appears in the
nonlinear drift term $a(B_{t},U_{t};W)$. A second difference is that the
probability density corresponding to (38)-(39) is a joint density $f_{BU}$ of
two (rather than one) variables.

The Fokker-Planck equation for the joint density $f_{BU}$ is

$\dfrac{\partial f_{BU}}{\partial t}=-\nabla\cdot\boldsymbol{J}_{BU},$ (41)

where

$\boldsymbol{J}_{BU}:=\begin{pmatrix}1-\sqrt{|b-(1+u)^{2}W|}b-\frac{1}{2}\sigma\partial_{b}\sigma\\\
-\gamma u-\frac{1}{2}\nu^{2}\partial_{u}\end{pmatrix}f_{BU},$ (42)

is the two-dimensional probability flux. Integrating (41) with respect to
$u\in(-\infty,\infty)$, gives an equation for the stationary marginal density
$f_{B}$, which is the quantity of interest:

$\dfrac{\partial}{\partial
b}\left(\mathbbm{E}[a(B_{t},U_{t}|W)|B_{t}=b]f_{B}\right)-\frac{1}{2}\dfrac{\partial}{\partial
b}\left(\sigma\left(\dfrac{\partial}{\partial b}\sigma f_{B}\right)\right)=0,$
(43)

where, making use of the conditional density $f_{U|B}:=f_{BU}/f_{B}$,

$\mathbbm{E}[a(B_{t},U_{t}|W)|B_{t}=b]:=\int\limits_{-\infty}^{\infty}a(b,u|W)f_{U|B}(u;b)\,\mathrm{d}u,$
(44)

is the unknown conditional expectation of the the drift $a(B_{t},U_{t}|W)$ for
a given value $b$ of $B_{t}$. Unless the full equation (41) is solved, an
assumption about the conditional expectation is required to make progress with
(43).

If $\gamma$ is sufficiently large, the modelled wind will evolve on a much
faster time scale than the buoyancy and can be regarded as being in a state of
quasi-equilibrium. Corresponding to such a separation of time scales is the
assumption that the conditional density $f_{U|B}$ is determined by a balance
of terms involving only relatively large derivatives with respect to $u$ in
(41) [31]:

$\dfrac{\partial}{\partial
u}\gamma(uf_{U|B})-\frac{\nu^{2}}{2}\frac{\partial^{2}}{\partial
u^{2}}f_{U|B}\sim 0,$ (45)

which has the well-known solution [55]

$f_{U|B}(u;b)\sim\sqrt{\frac{\gamma}{\pi\nu^{2}}}\exp\left(\frac{\gamma
u^{2}}{\nu^{2}}\right).$ (46)

In this particular case, for which $B_{t}$ does not feature in (39), $f_{U|B}$
in (46) corresponds to the marginal density $f_{U}$.

The subsequent use of $f_{U|B}$ obtained from (46) in (44) to account for the
effect that the fast variable has on the evolution of the slow variable, while
removing the former from the calculations, is known as ‘adiabatic elimination’
[56, 31]. With such an approximation, the drift term in (43) can be expressed
as the gradient of a potential that is convolved with the marginal density of
the fast variable:

$\dfrac{\partial}{\partial
b}\left(-\dfrac{\mathrm{d}V^{*}_{B}}{\mathrm{d}b}f_{B}-\frac{1}{2}\sigma\dfrac{\partial}{\partial
b}\left(\sigma f_{B}\right)\right)=0,$ (47)

where, noting that $f_{U|B}=f_{U}$ from (46) does not depend on $b$,

$V^{*}_{B}\left(b;\frac{\gamma}{\nu^{2}}\right):=\int\limits_{-\infty}^{\infty}V_{B}\left(b;(1+u)^{2}W\right)f_{U}(u)\,\mathrm{d}u.$
(48)

Figure 10 displays contours of $V^{*}_{B}$ along with ‘equilibrium’ paths
corresponding to its local extrema. Sufficiently large wind fluctuations
(small $\gamma/\nu^{2}$) modify the original potential (see cf. figure 2) via
the convolution (48) by removing the meta-stable equilibrium corresponding to
forward flow at large wind strengths $W$. When $\gamma/\nu^{2}=10$ the
potential is convex and has a unique minimum corresponding to forward or
reverse flow for small or large values of $W$, respectively. It is in this way
that the nonlinear effects of fluctuations in the wind strength in (38) modify
the underlying structure of the equation, regardless of noise in the buoyancy
flux that can be included separately in (38) with a particular
parameterisation for $\sigma$.

As illustrated in figure 11, which shows equilibria of the potential surface
$V_{B}^{*}$ as a function of both $W$ and $\gamma/\nu^{2}$, the cusp point of
the potential, above which a second local minima in $V_{B}^{*}$ emerges is
$\gamma/\nu^{2}=30.775$. The relevance of this number to practical
applications will be discussed at the end of the next section.

Figure 11: The modified potential surface $V_{B}^{*}$ derived in §5 as a
function of base wind strength $W$ and the parameters of the Ornstein-
Uhlenbeck process. The dotted line corresponds to $\gamma/\nu^{2}=30.775$,
which is the value above which the modified potential admits an unstable fixed
point and an additional stable fixed point corresponding to forward flow.
Figure 12: The stationary density for the parameters
$(W,\varepsilon)=(1.3,0.1)$ based on the application discussed in §6. The
thick red line denotes the exact solution and the blue dashed line denotes the
Gaussian approximation (28). The vertical dashed line separates reverse flow
($b<W$) from forward flow ($W<b$). The vertical black lines indicate the
location of the stable fixed points for reverse flow and forward flow.

## 6 Application Application

Although it was convenient to use the nondimensionalised equations (4), it is
useful to consider the wind strength $W$ and noise amplitude $\varepsilon$
that might be encountered in practice. Following the example used recently in
the context of ventilation driven by the competing effects of wind and
buoyancy [18], consider a conference room of cross-sectional area
$\hat{S}=50\,\mathrm{m}^{2}$, height $\hat{H}=4\,\mathrm{m}$ and effective
opening area $\hat{A}=0.63\,\mathrm{m}^{2}$. The room is subjected to a
buoyancy flux of $\hat{F}=0.055\,\mathrm{m}^{4}\mathrm{s}^{-3}$ (corresponding
to approximately 13 occupants) and the prevailing wind is assumed to create a
mean pressure difference of $\Delta\hat{p}=1.22\,\mathrm{Pa}$. Together, these
variables imply that, using (2) and (3),
$W:=\Delta\hat{p}/\Delta\hat{p}_{c}=1.30$, and

$\hat{t}:=\left(\frac{27\hat{H}^{2}\hat{S}^{3}}{4\hat{A}^{2}\hat{F}}\right)^{1/3}t\approx
852\,t\quad(\mathrm{seconds}).$ (49)

For the quadratic diffusion model applied to the buoyancy flux in §3.2,
$\varepsilon$ characterises the effects of fluctuations in a unidirectional
flow through an opening or various kinds of bidirectional flow that can occur.
While it is therefore difficult to relate wind speed statistics to
$\varepsilon$ in a precise way without a model developed from first
principles, it is instructive to consider the dimensional version of the
equations (cf. (1)), for which care is needed in deducing the correct
nondimensionalisation of $\varepsilon$:

$\hat{H}\hat{S}\,\mathrm{d}\hat{B}_{t}=(\hat{F}-\hat{Q}\hat{b})\,\mathrm{d}\hat{t}+\hat{\sigma}\,\mathrm{d}\hat{\xi}_{t}.$
(50)

In (50) the Wiener process has dimensions corresponding to the square root of
time and $\hat{\sigma}:=\hat{\varepsilon}\hat{B}_{t}$ is the dimensional
version of (18). Therefore, denoting the dimensions of length and time as $L$
and $T$, respectively, $[\hat{\varepsilon}]=L^{3}T^{-1/2}$.

If $\hat{v}$ corresponds to a velocity scale that characterises fluctuations
in the buoyancy flux through the opening, then it is reasonable to express
$\hat{\varepsilon}$ in terms of $\hat{v}$ and the opening area $\hat{A}$, such
that $\hat{\varepsilon}=\hat{v}^{1/2}\hat{A}^{5/4}$ on dimensional grounds;
hence, nondimensionalising $\hat{\xi}_{t}$ and $\hat{B}_{t}$ in (50),

$\,\mathrm{d}B_{t}=\ldots+\frac{\hat{\varepsilon}}{\hat{H}\hat{S}}B_{t}\,\mathrm{d}\hat{\xi}_{t}=\ldots+\frac{\hat{v}^{1/2}\hat{A}^{5/4}}{\hat{H}\hat{S}}\left({\frac{27\hat{H}^{2}\hat{S}^{3}}{4\hat{A}^{2}\hat{F}}}\right)^{1/6}B_{t}\,\mathrm{d}\xi_{t},$
(51)

which, when compared with (8), shows that

$\varepsilon=\left(\frac{27^{2}\hat{v}^{6}\hat{A}^{11}}{4^{2}\hat{H}^{8}\hat{S}^{6}\hat{F}^{2}}\right)^{1/12}.$
(52)

If, somewhat crudely, the velocity $\hat{v}$ characterising fluctuations in
the buoyancy flux is assumed to be the average wind velocity of
$1.4\,\mathrm{m}\mathrm{s}^{-1}$ from the recent related study[18] on which
the values adopted here are based, then $\varepsilon\approx 0.0969\approx
0.1$.

The stationary probability density corresponding to
$(W,\varepsilon)=(1.30,0.1)$ is plotted in figure 12 and agrees reasonably
well with the approximation (28). Despite the fact that a state of forward
flow is a permissible solution to the steady-state deterministic equations
(7), figure 12 shows that the eventual likelihood of such a system being found
in a state of forward flow is negligible for the parameters considered above.
From a design perspective, the stationary state is indistinguishable from a
Gaussian distribution centred on the fixed point for reverse flow. For
statistically transient states arising, for example, from changes in either
$W$ or $\varepsilon$, characterising the behaviour of the system
probabilistically is more challenging.

Multiplying the nondimensional time scale from figure 7 by (49) for
$(W,\varepsilon)=(1.30,0.1)$ shows that the time one would have to wait for
such a system to settle into the stationary state shown in figure 12, having
started in forward flow, is far in excess of the durations for which such a
room is likely to be used (i.e. minutes to hours). Considerations of the
transient evolution of the probability density would therefore play an
important role in such applications. While relatively small fluctuations
correspond to certainty regarding the eventual state of the flow, they
correspond to uncertainty in producing a longer transient period than
relatively large fluctuations.

For the Ornstein-Uhlenbeck process discussed in §5, the random variable
$U_{t}$ corresponds to fluctuations in the wind speed divided by the base
velocity scale $\sqrt{2\Delta\hat{p}/\hat{\rho}_{0}}$ used to define $W$ in
(3). The statistics reported in a recent study on stochastic models for
wind[22, 18] suggest that plausible values of $\gamma/\nu^{2}$ in (46) lie
between $10$ and $1000$ and therefore encompass the bifurcation illustrated in
figure 11. Relatively small values of $\gamma/\nu^{2}$ within this range
(relatively long, large gusts of wind) will lead to reverse flow for large
base wind strengths, whereas relatively large values of $\gamma/\nu^{2}$
(relatively short, small gusts of wind) do not necessarily rule out the
possibility of forward flow at large base wind strengths.

## 7 Conclusions Conclusions

This work has described stochastic models for the ventilation driven by
opposing wind and buoyancy of well-mixed spaces. The results expose several
properties concerning the system’s predictability and stability for the first
time and highlight challenges associated with the incorporation of stochastic
effects into existing deterministic models for ventilation.

In the absence of rigorous justification for a particular form of stochastic
forcing from either theory or experiments, the models chosen for this study
were simple and of general applicability. In §2.2 a stochastic buoyancy flux
was added to the deterministic evolution in the form of a Stratonovich
integral of a Wiener process (8) as the limit of a stochastic process with
relatively short autocorrelation time. In general, the approach can account
for fluctuations of the buoyancy flux that depend on the system’s state and
the base wind strength $W$. The evolution of the resulting probability density
is determined by a Fokker-Planck equation which accounts for the stabilising
effects of the underlying potential and the destabilising effects of noise. An
analytically-defined Boltzmann-like distribution describes the corresponding
stationary density. In §5 the wind itself was modelled as an Ornstein-
Uhlenbeck process at the expense of introducing an additional dimension, which
was subsequently approximated using adiabatic elimination.

Although the deterministic system for this problem admits multiple steady
states, for sufficiently large opposing wind strengths and small amounts of
noise, the system’s stationary state is unlikely to be found in the meta-
stable region associated with forward flow. One reason for this is that, as
the wind strength increases, fluctuations in the buoyancy flux are, with
increasing liklihood, able to overcome the potential barrier separating
forward and reverse flow, whose height decreases. Another reason is that when
the wind itself is modelled as a stochastic process (see §5), fluctuations in
the wind convolve with the nonlinear function describing the system’s
underlying potential, which can lead to the removal of the meta-stable
equilibrium altogether, as shown in figure 10. In these respects, the state of
the stochastic system, which accounts for disturbances that are likely to be
present in practical applications, is easier to ‘determine’ in a design
context than it would be for the idealised ‘deterministic’ system. Although
the particularities of such results depend on the choice of stochastic model,
their qualitative properties are generic and expected to be found in a broad
range of possible stochastic models.

An important caveat to the conclusions above is that greater certainty in the
eventual state of the system is typically consistent with an increase in the
duration of the transient evolution of the probability density. It is
therefore difficult to be prescriptive about whether the system is likely to
be in a particular state without careful consideration of the base wind
strength and relative size of the fluctuations. Base wind strengths that are
close to critical are the most difficult in this regard, because they are
associated with bimodal probability densities with transients of the largest
duration.

A preeminent challenge to the future incorporation of stochastic effects in
existing deterministic models of ventilation is finding a physically faithful
interpretation and representation of the fluctuations, and reconciling them
with the established deterministic part of such models. In bulk models of
turbulence it is not always clear whether the established deterministic
equations represent the limit of a system with zero external noise or as the
macroscopic average of a system with finite (internal) noise. The former
interpretation is somewhat easier to approach theoretically and, if regarded
as the limit of a process with relatively short autocorrelation times, is
consistent with our use of the Stratonovich integral in (8) [45].

When systems involve internal noise (which is arguably always the case where
turbulence is concerned) the situation is more complicated because the
fluctuations can’t be removed. In that case it would be necessary to return to
fundamentals and develop a master equation for the system’s probability
density from first principles, rather than trying to ascribe meaning to an
otherwise ill-defined deterministic equation (4) from a Langevin perspective
[45]. With a suitable master equation determining every aspect of the system’s
evolution, the frequently debated ‘choice’ between Itô or Stratonovich
integration, would be irrelevant because the corresponding Langevin equation
for either approach would be defined to ensure consistency with the underlying
master equation.

Hopefully the considerations above will motivate dedicated experiments and
direct simulations of turbulence to underpin a rigorous formulation of the
required stochastic equations, which are likely to play an important role in
future applications.

## 8 Acknowledgements Acknowledgements

This work was supported by the Engineering and Physical Sciences Research
Council [grant number EP/V033883/1] as part of the [D∗]stratify project.

## References

  * [1] Linden PF. 1999 The fluid mechanics of natural ventilation. Annual Review of Fluid Mechanics 31, 201–238. (10.1146/annurev.fluid.31.1.201)
  * [2] Hunt GR, Linden PF. 2005 Displacement and mixing ventilation driven by opposing wind and buoyancy. Journal of Fluid Mechanics 527, 27–55.
  * [3] Hunt GR, Linden PP. 1999 The fluid mechanics of natural ventilation–displacement ventilation by buoyancy-driven flows assisted by wind. Building and Environment 34, 707 – 720.
  * [4] Ouf MM, O’Brien W, Gunay B. 2019 On quantifying building performance adaptability to variable occupancy. Building and Environment 155, 257–267. (https://doi.org/10.1016/j.buildenv.2019.03.048)
  * [5] Hong T, Taylor-Lange SC, D’Oca S, Yan D, Corgnati SP. 2016 Advances in research and applications of energy-related occupant behavior in buildings. Energy and Buildings 116, 694–702. (https://doi.org/10.1016/j.enbuild.2015.11.052)
  * [6] Delzendeh E, Wu S, Lee A, Zhou Y. 2017 The impact of occupants’ behaviours on building energy analysis: A research review. Renewable and Sustainable Energy Reviews 80, 1061–1071. (https://doi.org/10.1016/j.rser.2017.05.264)
  * [7] Sun K, Hong T. 2017 A framework for quantifying the impact of occupant behavior on energy savings of energy conservation measures. Energy and Buildings 146, 383–396. (https://doi.org/10.1016/j.enbuild.2017.04.065)
  * [8] Laaroussi Y, Bahrar M, Elmankibi M, Draoui A, Si-Larbi A. 2019 Occupant behaviour: a major issue for building energy performance. In IOP Conference Series: Materials Science and Engineering vol. 609 p. 072050. IOP Publishing.
  * [9] Boettcher F, Renner C, Waldl HP, Peinke J. 2003 On the statistics of wind gusts. Boundary-Layer Meteorology 108, 163–173.
  * [10] Coceal O, Thomas TG, Castro IP, Belcher SE. 2006 Mean flow and turbulence statistics over groups of urban-like cubical obstacles. Boundary-Layer Meteorology 121, 491–519.
  * [11] Efthimiou G, Hertwig D, Andronopoulos S, Bartzis J, Coceal O. 2017 A statistical model for the prediction of wind-speed probabilities in the atmospheric surface layer. Boundary-Layer Meteorology 163, 179–201.
  * [12] Burridge HC, Bhagat RK, Stettler MEJ, Kumar P, De Mel I, Demis P, Hart A, Johnson-Llambias Y, King MF, Klymenko O, McMillan A, Morawiecki P, Pennington T, Short M, Sykes D, Trinh PH, Wilson SK, Wong C, Wragg H, Davies Wykes MS, Iddon C, Woods AW, Mingotti N, Bhamidipati N, Woodward H, Beggs C, Davies H, Fitzgerald S, Pain C, Linden PF. 2021 The ventilation of buildings and other mitigating measures for COVID-19: a focus on wintertime. Proceedings of the Royal Society A: Mathematical, Physical and Engineering Sciences 477, 20200855. (10.1098/rspa.2020.0855)
  * [13] Bhagat RK, Davies Wykes MS, Dalziel SB, Linden PF. 2020 Effects of ventilation on the indoor spread of COVID-19. Journal of Fluid Mechanics 903, F1. (10.1017/jfm.2020.720)
  * [14] Batchelor GK. 1954 Heat convection and buoyancy effects in fluids. Quarterly Journal of the Royal Meteorological Society 80, 339–358.
  * [15] Linden PF, Lane-Serff GF, Smeed DA. 1990 Emptying filling boxes: the fluid mechanics of natural ventilation. Journal of Fluid Mechanics 212, 309–335.
  * [16] Baines WD, Turner JS. 1969 Turbulent buoyant convection from a source in a confined region. Journal of Fluid Mechanics 37, 51–80. (10.1017/S0022112069000413)
  * [17] Karava P, Stathopoulos T, Athienitis A. 2004 Wind Driven Flow through Openings – A Review of Discharge Coefficients. International Journal of Ventilation 3, 255–266. (10.1080/14733315.2004.11683920)
  * [18] Vesipa R, Ridolfi L, Salizzoni P. 2023 Wind fluctuations affect the mean behaviour of naturally ventilated systems. Building and Environment 229, 109928. (https://doi.org/10.1016/j.buildenv.2022.109928)
  * [19] Leprince J, Madsen H, Miller C, Real JP, van der Vlist R, Basu K, Zeiler W. 2022 Fifty shades of grey: Automated stochastic model identification of building heat dynamics. Energy and Buildings 266, 112095. (https://doi.org/10.1016/j.enbuild.2022.112095)
  * [20] Wolf S, Calí D, Krogstie J, Madsen H. 2019 Carbon dioxide-based occupancy estimation using stochastic differential equations. Applied Energy 236, 32–41. (https://doi.org/10.1016/j.apenergy.2018.11.078)
  * [21] Macarulla M, Casals M, Forcada N, Gangolells M, Giretti A. 2018 Estimation of a room ventilation air change rate using a stochastic grey-box modelling approach. Measurement 124, 539–548. (https://doi.org/10.1016/j.measurement.2018.04.029)
  * [22] Ma J, Fouladirad M, Grall A. 2018 Flexible wind speed generation model: Markov chain with an embedded diffusion process. Energy 164, 316–328. (https://doi.org/10.1016/j.energy.2018.08.212)
  * [23] Palmer T, Williams P. 2008 Introduction. Stochastic physics and climate modelling. Philosophical Transactions of the Royal Society A: Mathematical, Physical and Engineering Sciences 366, 2419–2425. (10.1098/rsta.2008.0059)
  * [24] Fontanini A, Vaidya U, Ganapathysubramanian B. 2013 A stochastic approach to modeling the dynamics of natural ventilation systems. Energy and Buildings 63, 87–97. (https://doi.org/10.1016/j.enbuild.2013.03.053)
  * [25] Coomaraswamy I, Caulfield C. 2011 Time-dependent ventilation flows driven by opposing wind and buoyancy. Journal of fluid mechanics 672, 33–59.
  * [26] Bellman R, Kalaba R. 1959 A mathematical theory of adaptive control processes. Proceedings of the National Academy of Sciences 45, 1288–1290.
  * [27] Van Kampen N. 2011 Stochastic Processes in Physics and Chemistry. North-Holland Personal Library. Elsevier Science.
  * [28] Li Y, Delsante A. 2001 Natural ventilation induced by combined wind and thermal forces. Building and Environment 36, 59–71.
  * [29] Chenvidyakarn T, Woods A. 2005 Multiple steady states in stack ventilation. Building and Environment 40, 399 – 410. (https://doi.org/10.1016/j.buildenv.2004.06.020)
  * [30] Strogatz S. 1994 Nonlinear Dynamics and Chaos: With Applications to Physics, Biology, Chemistry, and Engineering. Advanced book program. Westview Press.
  * [31] Haken H. 2013 Synergetics: An Introduction. Springer Series in Synergetics. Springer Berlin Heidelberg.
  * [32] Yuan J, Glicksman LR. 2008 Multiple steady states in combined buoyancy and wind driven natural ventilation: The conditions for multiple solutions and the critical point for initial conditions. Building and Environment 43, 62–69. (https://doi.org/10.1016/j.buildenv.2006.11.035)
  * [33] Lishman B, Woods AW. 2009 On transitions in natural ventilation flow driven by changes in the wind. Building and Environment 44, 666–673.
  * [34] Craske J, Hughes G. 2019 On the robustness of emptying filling boxes to sudden changes in the wind. Journal of Fluid Mechanics 868, R3. (10.1017/jfm.2019.199)
  * [35] Grafke T, Vanden-Eijnden E. 2017 Non-equilibrium transitions in multiscale systems with a bifurcating slow manifold. Journal of Statistical Mechanics: Theory and Experiment 2017, 093208. (10.1088/1742-5468/aa85cb)
  * [36] Meunier C, Verga A. 1988 Noise and bifurcations. Journal of statistical physics 50, 345–375.
  * [37] Kuehn C. 2011 A mathematical framework for critical transitions: Bifurcations, fast–slow systems and stochastic dynamics. Physica D: Nonlinear Phenomena 240, 1020–1035. (https://doi.org/10.1016/j.physd.2011.02.012)
  * [38] Burrage K, Burrage PM, Tian T. 2004 Numerical methods for strong solutions of stochastic differential equations: an overview. Proceedings of the Royal Society of London. Series A: Mathematical, Physical and Engineering Sciences 460, 373–402. (10.1098/rspa.2003.1247)
  * [39] Schenk-Hoppé KR. 1996 Bifurcation scenarios of the noisy Duffing-van der Pol oscillator. Nonlinear Dynamics 11, 255–274.
  * [40] Mendez A, Farazmand M. 2021 Investigating climate tipping points under various emission reduction and carbon capture scenarios with a stochastic climate model. Proceedings of the Royal Society A: Mathematical, Physical and Engineering Sciences 477, 20210697. (10.1098/rspa.2021.0697)
  * [41] Namachchivaya NS. 1990 Stochastic bifurcation. Applied Mathematics and Computation 38, 101–159.
  * [42] Kramers H. 1940 Brownian motion in a field of force and the diffusion model of chemical reactions. Physica 7, 284–304. (https://doi.org/10.1016/S0031-8914(40)90098-2)
  * [43] Virgin LN, Plaut RH, Cheng CC. 1992 Prediction of escape from a potential well under harmonic excitation. International Journal of Non-Linear Mechanics 27, 357–365. (https://doi.org/10.1016/0020-7462(92)90005-R)
  * [44] Gladstone C, Woods AW. 2001 On buoyancy-driven natural ventilation of a room with a heated floor. Journal of Fluid Mechanics 441, 293–314. (10.1017/S0022112001004876)
  * [45] van Kämpen NG. 1981 Itô versus Stratonovich. Journal of Statistical Physics 24, 175–187. (10.1007/BF01007642)
  * [46] Resnick S. 1992 Adventures in Stochastic Processes. Adventures in Stochastic Processes. Birkhäuser Boston.
  * [47] Moon W, Wettlaufer JS. 2014 On the interpretation of Stratonovich calculus. New Journal of Physics 16, 055017. (10.1088/1367-2630/16/5/055017)
  * [48] Doob JL. 1942 The Brownian Movement and Stochastic Equations. Annals of Mathematics 43, 351–369.
  * [49] Uhlenbeck GE, Ornstein LS. 1930 On the Theory of the Brownian Motion. Phys. Rev. 36, 823–841. (10.1103/PhysRev.36.823)
  * [50] Risken H, Frank T. 1996 The Fokker-Planck Equation: Methods of Solution and Applications. Springer Series in Synergetics. Springer Berlin Heidelberg.
  * [51] Einstein AB On the movement of small particles suspended in a stationary liquid demanded by the molecular-kinetic theory of heart. Annalen der Physik 17, 549–560.
  * [52] Langevin P. 1908 Sur la théorie du mouvement brownien. Compt. Rendus 146, 530–533.
  * [53] Lemons DS, Gythiel A. 1997 Paul Langevin’s 1908 paper “On the Theory of Brownian Motion” [“Sur la théorie du mouvement brownien,” C. R. Acad. Sci. (Paris) 146, 530–533 (1908)]. American Journal of Physics 65, 1079–1081. (10.1119/1.18725)
  * [54] Landau L, Lifshitz E. 2013 Statistical Physics: Volume 5. Number v. 5. Elsevier Science.
  * [55] Gardiner C. 1990 Handbook of Stochastic Methods: For Physics, Chemistry and the Natural Sciences. Lecture Notes in Mathematics. Springer Berlin Heidelberg.
  * [56] Theiss W, Titulaer U. 1985 The systematic adiabatic elimination of fast variables from a many-dimensional Fokker-Planck equation. Physica A: Statistical Mechanics and its Applications 130, 123–142. (https://doi.org/10.1016/0378-4371(85)90100-1)