Methodology for risk assessment of COVID-19 pandemic propagation.

This paper proposes a methodology to perform risk analysis of the virus spread. It is based on the coupling between CFD modelling of bioaerosol dispersion to the calculation of probability of contact events. CFD model of near-field sneeze droplets dispersion is developed to build the SARS-CoV-2 effect zones and to adequately capture the safe distance. The most shared classification of droplets size distribution of sneezes was used. Droplets were modeled through additive heating/evaporation/boiling laws and their impact on the continuous phase was examined. Larger droplets move behind the droplet nuclei front and exhibit greater vertical drop due to the effect of gravity. CFD simulations provided the iso-risk curves extension (i.e., the maximum distance as well as the angle) enclosed by the incident outcome effect zone. To calculate the risk indexes, a fault tree was developed and the probability of transmission assuming as of the top event "COVID-19 infection" was calculated starting from the virus spread curve, as main base case. Four phases of virus spread evolution were identified: initiation, propagation, generalised propagation and termination. For each phase, the maximum allowable close contact was computed, being fixed the values of the acceptable risk index. In particular, it was found that during the propagation case, the maximum allowable close contacts is two, suggesting that at this point lockdown should be activated. The here developed methodology could drive policy containment design to curb spread COVID-19 infection.

Introduction

COVID-19 outbreak demonstrates that respiratory infectious diseases can spread like wildfire in an intimately connected world. Unprecedented containment strategies have been applied to limit the diffusion of COVID-19, including lockdown, isolation, quarantine and social distancing (d ≥ 1–1.5 m). This distance is estimated by neglecting convection flow and atmospheric diffusivity that it is well known to significantly affect dispersion of exhaled aerosol clouds (Gorbunov, 2020).

Understanding the main virus transmission routes is crucial for policy containment design to curb spread COVID-19 infection. According to current evidence, COVID-19 virus is primarily transmitted between people through respiratory droplets and contact routes (Chan et al., 2020; Huang et al., 2020; Li et al., 2020; Liu et al., 2020).

Droplet transmission occurs when a person is close to a SARS-CoV-2 positive individual with respiratory symptoms (e.g., coughing and/or sneezing) as well as heavy breathing. In such condition, the person risks of having the mucosae exposed to potentially infective respiratory droplets (direct transmission) (World Health Organization, 2020a). Moreover, the droplet transmission can occur by indirect contact with surfaces in the immediate environment (fomite) around the infected individual (indirect transmission). The virus remains more stable on plastic and stainless steel than on copper and cardboard, and vital virus was detected up to 72 h after application on these surfaces, although the number of viral particles was largely reduced (van Doremalen et al., 2020). Surface cleaning and disinfection with different biocidal agents (e.g., sodium hypochlorite or ethanol) significantly reduces coronavirus infectivity on surfaces within 1 min exposure time (Kampf et al., 2020).

In the 1930's William F. Wells first introduced the classification of respiratory droplet emissions into large and small droplets (Wells, 1934). Larger droplets settle without appreciable loss by evaporation, contaminating the immediate surrounding area of the infected individual. These are responsible of the direct and/or indirect droplet transmission of virus (World Health Organization, 2020a). Small droplets evaporate faster than they settle, forming residual particulates of dried material from the original droplets. These residual particulates are referred to as droplet nuclei or aerosols (Wells, 1934). Airborne transmission refers to the presence of viruses within droplet nuclei, (particles < 5 μm in diameter) that can remain in the air for long periods of time (up to 3 h) and be transmitted to others over distances greater than 1 m (van Doremalen et al., 2020; World Health Organization, 2020a). Therefore, social distancing do not prevent infection by inhalation of pathogenic droplet nuclei exhaled by an infected individual that can travel distance of about 7–8 m (Bourouiba et al., 2014). According to the available results, it seems that SARS-CoV-2 infection is transmitted via large droplets (dimensions higher than 5–10 μm) (Bourouiba, 2020). However, in a recent work, Zhang et al. (2020) concluded that the airborne transmission is dominant for the spread of COVID-19 and that the simultaneous wearing of face masks, social distancing, quarantine and contact tracing is the most effective prevention against interhuman transmission (Zhang et al., 2020).

The most effective preventive measures include physical distance, the use of personal protective equipment (PPE, e.g., mask, gloves), handwashing, avoiding touching mucosae, routine cleaning of touched surfaces (World Health Organization, 2020b). As regards the protective action of PPE, Bae et al. (2020) reported that both surgical and cotton masks seem to be ineffective in preventing the dispersion of SARS–CoV-2 from the coughs of patients with Covid-19 to the environment and external mask surface, while Leung et al. (2020) stated that face masks significantly reduced detection of influenza virus RNA in aerosol particles. An investigation of effects of contaminations on the mask efficiency tests needs to be carried out (Bae et al., 2020; Leung et al., 2020).

The droplet concentration and sedimentation in indoor scenarios (e.g., workplace, school, universities, restaurants, hospitals, waiting rooms …) may be significantly affected by ventilation and air recirculation (Gao et al., 2016; Li et al., 2007; Morawska, 2006; Morawska and Cao, 2020; Qian and Zheng, 2018). Recently, Gorbunov (2020) showed that convection flow, atmospheric diffusivity and humidity on evolution and travel distances of exhaled aerosol clouds by an infected person play a significant role in affecting the transmission routes. The results obtained demonstrate that aerosol particles generated by coughing and sneezing can travel over 30 m (Gorbunov, 2020).

Several authors proposed CFD simulation of bioaerosol dispersion in indoor scenarios studying the effect of ventilation, sources position and air recirculation on far-field dispersion (distance above 1–2 m) (Aliabadi et al., 2010; Borro et al., 2020; Dhakar, 2019; King et al., 2015; La and Zhang, 2019; Thatiparti et al., 2017; Zhao et al., 2005).

CFD results show that ventilation play a major role in reducing the infection propagation. However, all the CFD simulations require to fix the spreader position, the spreader number, the other people positions, the chosen geometrical domain and the boundary conditions etc … As a consequence, the quantitative results obtained by hypothesising a specific scenario cannot be generalised.

In order to generalise the results, CFD simulations should be coupled to the probability analysis to eventually perform risk assessment. In this case, it is possible to derive risk indexes which contain the information about the effect zones (from CFD simulations) and the probability, thus providing generalised indications. It is worth noting that from the coupling of the CFD simulations to the pandemic curve evolution and to the frequencies/probability, the casualty due to the spreader position, number etc is intrinsically taken into account.

In this work, we performed the risk assessment of the contagion propagation. More precisely, we aimed at developing a methodology to keep the risk of contagion under control by identifying the maximum number of possible close contacts at each phase of the pandemic propagation, being fixed the values of the acceptable risk index.

The engineering methodology we here developed could be a guidance to minimise the risk of contagion propagation and then it could be used to organize work shifts and regulate the division between work in presence and smart working, as well as to public activities (e.g., restaurant, school, university etc.).

The risk assessment was performed starting from the evaluation of the SARS-CoV-2 effect zones to adequately capturing the safe distance. Effect zones and iso-risk curves extension were computed by means of a CFD model of near-field dispersion to fully capturing the aerosol dispersion at quiescent conditions and by taking into account the effect of ventilation air flow.

To evaluate the probability of transmission assuming as top event “COVID-19 infection”, we developed a fault tree.

Through a quantification of the impact of the number of possible close contacts and of active cases (e.g., a sensitivity analysis) on the top event occurrence frequency, we provided a simple as well as generalizable methodology to keep the risk of contagion under control.

Methodology

Isorisk maps

The Individual Risk (IR) is calculated to quantify the iso-risk curves which estimate risk at specific points by taking into account both the effect zone and the probability. The total individual risk at each point is evaluated as the sum of the individual risks, at that point, of all incident outcome cases associated with the hazard source. Individual Risk at location x, y (IRx,y,i) is given by the product between the frequency of the incident outcome case i (f ) and the probability (p ) the incident outcome case i will result in a fatality at location x,y. In this specific case, incident outcomes refer to the exposure to a transmission route (direct or indirect droplet transmission, airborne transmission) while the fatality represents the infection. For a more conservative assessment, we considered each incident outcome case with an equal impact (probability of fatality/contagion p  = 1) throughout its geographical impact zone: within the impact zone, the individual risk IRx,y,i is equal to the frequency of that incident outcome case, outside IRx,y,i is = 0. Moreover, the effect zone will include all the droplets coming from the release source. Considering the droplets as the vectors of infection, it is reasonable to assume the probability of contagion outside the effect zone equal to 0.

To assess the frequency of each transmission route exposure, fault tree probabilistic method was carried out. The approach starts with a well-defined top event and worked downward, through the intermediate events, towards the various basic events that can cause the accident. Basic events that must take place for the top event to occur were connected by an AND logic gate while those related in series were linked by an OR gate. In analytical methods for quantitative evaluation of fault trees, Boolean algebraic operations are used starting from basic events across the various logic gates to determine the probability of the top event (Crowl and Louvar, 2002; Lees, 1996). Probabilities of failure (or unreliabilities P) and failure rates μ of all the basic events have to be assigned. Failure probability and failure rate are connected as in the following formulathat is valid when the failure rate is reasonably constant (simplified approach) (Crowl and Louvar, 2002).

Probabilities of failure have to be set according to the situation being analysed thanks to the help of further professionals as experts of infection transmission, occupational medicine, biological risk in workplaces, risk assessment and management of indoor environment. Once probabilities are set, boolean algebraic operations are used across the various logic gates to determine the probability of the top event and in particular to quantify the frequency of each transmission route. For example, the contribution (i.e., probability) of aerosol transmission will be fundamental in indoor environments while it will be negligible outdoors.

CFD model

To quantify the iso-risk curves dimension, CFD simulation of a sneeze was carried out. Human exhalation such as coughing, sneezing and breathing can be considered as instantaneous airflows produced from a single source with a quite symmetrical and conical geometry. Although coughing and sneezing have gained much attention as potential, explosive sources of infectious aerosols, these are relatively rare events during daily life (high-consequence low-probability events) (Tang et al., 2013). Recent work has demonstrated that respiratory symptoms not only consist of droplets that quickly follow down by short-range deposition trajectories but, importantly, are made of a multiphase turbulent gas cloud instantaneously released (i.e., puff) that entrains ambient air and carries aerosol particles (Bourouiba et al., 2014; Scharfman et al., 2016). The computational domain and mesh were built and refined close to the surface of the body for CAD data in the shape of a person by means of the Design Modeler and Meshing packages of Ansys (Release 19). The head was then enclosed in a box, with the frontal wall and the lateral ones at 2 m and 0.5 m away, respectively. The top wall was put on the top of the head. The adequacy of the domain dimensions was preliminarily assessed using a larger domain (1.5 m high with side walls 0.8 m from the body, total cells 801227, minimum face area (m2): 1.7e-08, maximum face area (m2): 2.5e-03). From the simulations obtained with this enlarged domain, it was possible to observe that the boundary conditions set on the top and lateral surfaces do not influence the maximum extension of the cloud, always found in the x-direction. The reduction in the size of the domain allowed a reduction of the cells and of the computational cost. We obtained the mesh using Ansys Meshing. The full study of the adequacy of the domain dimensions was reported in detail in Supplementary Materials (Section S1). The element order was program controlled. The body surface was preliminary and randomly divided into 308 faces (total faces 314). We refined mesh in the proximity of the body by capturing curvature and setting a number of elements equal to 100 per body surface face. The maximum face size is maximum face area (m2): 6.72e-03 (maximum element size is 180 mm) for a total of 497118 elements. The target skewness is 0.9 and the smoothing is fast. The convergence was verified with three additional grids: 848601 elements (maximum size 50 mm), 1185228 elements (maximum size 25 mm) and a coarser grid 275955 elements (maximum size 360 mm). These grids were obtained by Ansys meshing varying the maximum size of the cells and the convergence was attested comparing the air velocity profiles along a horizontal line at 200 ms and the droplets cloud maximum extension at different times. The grid-independence analysis was detailed in Supplementary Materials (Section S2). In Fig. 1 and Fig. 2 show the computational domain considered (x-z view) and a section (x-z plane) of the used unstructured and non-uniform mesh, respectively.

Fig. 1

Sketch of the computational domain (x-z view).

Fig. 2

Geometry (x-z view) and Section (x-z plane) of the unstructured and non-uniform mesh used.

The model used for droplets dispersion simulation consists of the time-averaged Navier-Stokes equations (Eulerian approach) for the continuous phase (air), solved by using the standard k-ε model as turbulent sub-model. Indeed, sneezes enforce high velocities (up to 30 m/s) flows of air through narrow openings in a short time range (200–250 ms), resulting in turbulent flows (Re ≈ 104) (Scharfman et al., 2016). The k-ε model only solves the boundary layer in the log-layer, so it requires y + value between 20 and 300 for closure. The y + condition at the wall was ensured and was detailed in Supplementary Materials (Section S3).

The Navier-Stokes equations were discretized using a finite-volume formulation on the three-dimensional non-uniform unstructured grid. The semi-implicit method for pressure-linked equations (SIMPLE) was used to solve the pressure-velocity coupled equations. The spatial discretization of the model equations used first order schemes for convective terms and second order schemes for diffusion terms. First-order time integration was used to discretize temporal derivatives with a time step of 1•10−4 s. To verify that the choice was proper after the calculation is complete, we plotted contours of the Courant number within the domain. For a stable calculation, the Courant number should not exceed a value of 20–40 in most sensitive transient regions of the domain. The independence on the time step was ensured and was detailed in Supplementary Materials (Section S4). Discrete phase modelling (DPM) with the Lagrangian formulation was performed to simulate aerosol movement. Droplets can be modeled through additive heating/evaporation/boiling laws and their impact on the continuous phase can be examined (Fluent Inc, 2016). The momentum balance equation of the DPM reads as follows (Fluent Inc, 2016):where is the drag force per unit particle mass, function of the Reynolds dimensionless number according to the following equation

The Reynolds number is defined in the following

Here, u is the fluid phase averaged velocity, u is the particle velocity, μ is the molecular viscosity of the fluid, ρ is the fluid density, ρ p is the density of the particle, C is the drag coefficient and d p is the particle diameter. Equation (2) includes additional forces (f) per unit particle mass. The first of these is the “virtual mass” force, the force required to accelerate the fluid surrounding the particle. This force can be written as:

An additional force arises due to the pressure gradient in the fluid:

Several laws for drag coefficients, C , are available for the Euler-Lagrange model. The spherical drag law for smooth particles is defined in the followingwhere a , a and a are constants given by Morsi and Alexander (1972). Moreover, we included the turbulent dispersion activating the Discrete Random Walk (DRW) model. In this way, the interaction of a particle with a succession of discrete stylized fluid phase turbulent eddies is simulated (Fluent Inc, 2016). Fluent predicts the trajectory of a discrete phase particle by integrating the force balance on the particle. The unsteady particle tracking integration time step was taken equal to the fluid flow time step. Parallel calculations were performed by means of the segregated pressure-based solver of the code ANSYS Fluent (release 19). In order to achieve convergence, all residuals were set equal to 1•10−6. The simulation conditions are summarised in Table 1 . Fig. 3 shows the initial conditions of the aerosol cloud.

Table 1

Simulation conditions.

Parameter	Value
Type of injection	Surface injection (from nostrils)
Particle Type	Droplet
Material	Water-liquid
Density (kg m−3)	998.2
Drag-Law	Smooth spherical particles
Stochastic tracking	Discrete walk model
Dispersing phase in the room	80% RH
Ambient temperature (°C)	20
Droplets diameter (μm)	Size distribution by Duguid (Duguid, 1946)
Droplets number	500 k
Time step (s)	1•10−4
Velocity of the sneeze (m/s)	30
Sneeze air temperature (°C)	35
Sneeze duration (s)	0.2

Fig. 3

Initial conditions of cloud and quiescent background:initial air density,partial water pressure as function of ambient pressure and temperature,gravity acceleration,initial momentum of the cloud,initial cloud density,initial volume of the cloud,initial boundary of the cloud.

The injection of the aerosol was modeled as a 0.200 s continuous release from the nostrils with a fluid velocity of 30 m/s. As regards droplet size of sneezes, several researchers have focused on it, but high discrepancies on the size distribution can be found in the literature (Duguid, 1946; Han et al., 2013; Johnson et al., 2011; Morawska et al., 2009; Papineni and Rosenthal, 1997; Yang et al., 2007; Zayas et al., 2012). The most shared classification of droplets size distribution of sneezes is that showed by Duguid in the 1940s (Duguid, 1946). Using the Rosin-Rammler expression to represent the droplets size distribution, the mass fraction Y of droplets of diameter greater than d is given bywhere d is the size constant and n is the size distribution parameter. Fig. 4 shows the droplet size distribution of sneeze as found by Duguid through by microscopy measurements. Duguid reported droplet sizes in the range 1–360 μm with a peak of 4–8 μm for sneezes (Duguid, 1946). The value for d is obtained by noting that this is the value of d at which Y  = 0.368 (d  = 7.8 μm).

Fig. 4

Cumulative Size Distribution of Particles measured by Duguid (Duguid, 1946).

To verify the effect of ventilation on the near-field dispersion, two background were computed: quiescent background and wind at 1 m/s from the backwall. Typically, stationary air is not recommended since minimum air change rates are required to maintain indoor air quality. Conversely, wind speed equal to 1 m/s and above may be felt as a draught in a hot climate (ASHRAE Standard 55, 2010). However, the simulation of the near-field dispersion will allow to visualize the routes of transmission by following the paths of particles. Boundary conditions are shown in Table 2 .

Table 2

Boundary conditions.

Case	Left, top, right sides	Backwall	Front wall	Ground	Body	Nostrils
Quiescent air	Fluid: wallParticles: escape	Fluid: wallParticles: escape	Fluid: wallParticles: escape	Fluid: wallParticles: trap	Fluid: wallParticles: trap	Fluid: inlet (30 m/s, 0.2 s), after wallParticles: escape
Wind from the backwall	Fluid: wallParticles: escape	Fluid: inlet (1 m/s)Particles: escape	Fluid: pressure outflowParticles: escape	Fluid: wallParticles: trap	Fluid: wallParticles: trap	Fluid: inlet (30 m/s, 0.2 s), after wallParticles: escape

Results

CFD simulations

Fig. 5 shows the time sequences of the droplet diameter distribution ((x-z) view) at 5, 50, 100, 150 and 200 ms (a,b,c,d,e respectively) as simulated at quiescent conditions. The natural symmetric and conical shape of a sneeze can be appreciated. Considering the total number of injected droplets (500 k), at 200 ms, 0.2% of droplets evaporated while 1.3% is trapped by the surface at the bottom. At the end of the sneeze, larger droplets, responsible of the direct and/or indirect droplet transmission of the virus, tends to settle without appreciable loss by evaporation, contaminating the immediate surrounding area of the infected individual.

Fig. 5

Time sequence of the droplet diameter distribution as computed over (x–z) view at 5, 50, 100, 150 and 200 ms (a,b,c,d,e respectively).

At 340 ms (Fig. 6 ), the cloud mainly consists of small droplets (<200 μm) while the largest are trapped by the surface at the bottom. The maximum distance of the cloud (x direction) is equal to 70 cm, in agreement with that found by Scharfman et al. (2016) at 340 ms in their direct high-speed recordings of the droplets emitted at the exit of the mouth and nose during violent exhalation events (Scharfman et al., 2016).

Fig. 6

Droplets diameters distribution, (x–z) view at 340 ms

After 1 s from the beginning of the sneeze (Fig. 7 ), the cloud contains droplets with diameters < 100 μm. The maximum distance of the cloud (x direction) is equal to 1.20 m, greater than the minimum safety distance suggested by policymakers (World Health Organization, 2020a). This result confirms that even if the direct droplet transmission occurs only when a person is close (within 1 m) to a SARS-CoV-2 positive individual at the moment of respiratory symptoms, the large part of droplets (≪ 100 μm) contribute to airborne transmission while the rest with larger diameters settled contributing to the indirect droplet transmission.

Fig. 7

Droplet diameters distribution, (x–z) view at 1 s.

In addition to social distancing against direct droplet transmission and periodic surface cleaning against indirect droplet transmission, the airborne transmission has to be limited to effectively face the spread of COVID-19. There are strong evidences that demonstrate the association between ventilation, air recirculation in buildings and the transmission/spread of infectious diseases (Aliabadi et al., 2010; Borro et al., 2020; Dhakar, 2019; King et al., 2015; La and Zhang, 2019; Thatiparti et al., 2017; Zhao et al., 2005). Given the high number of droplet nuclei expelled through a sneeze, the design of the best and safest configuration and the position of the ventilation/aspiration is of crucial importance. The near-field dispersion is characterized by high turbulence, small time scales and rapid heat and mass transfer with background medium. Thus, the dispersion can be decoupled in the near-field (few meters, probably not influenced by ventilation systems) and the far-field (whole ventilation space) contributions. To verify the effect of ventilation on the near-field dispersion, two background were computed. Fig. 8 and Fig. 9 show the comparison between the aerosol dispersion and the velocity flow field at 1 s with quiescent background (a) and wind at 1 m/s from the backwall (b), respectively. In each case, after 1 s from the beginning of the sneeze, the cloud contains droplets with diameters <100 μm. The drag and gravity forces can be calculated as in the following for both 1 μm and 360 μm droplets where [kg/(m s)] is the viscosity of the fluid, [m] is the radius of the droplet, [m/s] is the mean velocity of the fluid, [kg/m3] is the droplet density, [kg/m3] is the fluid density, [m/s2] is the gravity acceleration. The drag forces are of the order of magnitude equal to 10−10 N and 10−7 N while the gravity forces are of the order of magnitude equal to 10−15 N and 10−6 N for 1 μm and 360 μm droplets respectively. Thus, for larger droplets the gravity force prevails over the drag force while for the smaller ones the opposite applies. The ventilation increases the maximum distance of the cloud (x direction) compared to the case of the quiescent background. Conversely, the air flow has negligible effects on dispersed droplets aliquot (18.2% and 17% of the initial droplets number (500 k) for quiescent background and backwall, respectively).

Fig. 8

Comparison between the aerosol dispersion at 1 s with quiescent background (a), air flow velocity = 1 m/s from the backwall (b), (x–z) view.

Fig. 9

Comparison between the velocity flow fields at 1 s with quiescent background (a), wind at 1 m/s from the backwall (b), (x–z) view.

In Fig. 10 the position of the aerosol with respect the infected person is shown as function of time in the case of quiescent conditions and in the presence of air flow (v = 1 m/s). The position of droplets with diameters larger than 10 μm in the case of quiescent conditions is also reported. It is noting that the effect of ventilation is not straightforward. It depends on the relative position of the ventilation source and the infected person. In this example, ventilation is backwall with respect to the infected person and then it increases the aerosol movement increasing the distance. As regards larger droplets, they move behind the aerosol front, probably reaching a plateau value of distance in a very short time (after 1 s from the beginning of the sneeze). As shown, the effect of ventilation is more pronounced as the droplets cloud moves away from the source.

Fig. 10

Aerosol position as function of the delay time with respect to the source.

From these results it can be concluded that the safe distance value depends on the route of transmission and in particular on the particle dimensions.

Particle nuclei: the distance of 1 m is not sufficient. Aerosol remains suspended in air and their dispersion depends on ventilation.

Particle of dimension higher than 5–10 μm: the distance of 1 m is sufficient since in a very short time the particle settle. Transmission may occur via the infected surface.

Fig. 11 shows the time sequence of iso-risk curves in quiescent conditions as function of time. The sneeze is directional thus it is represented by a pie-shaped section. Within the effect zone that includes all the suspended droplets, the individual risk IR is equal to the Top event frequency, that is the sum of direct droplet (f ), indirect droplet (f ) and airborne (f ) transmissions frequencies, outside IR = 0. As the time increases, the maximum distance as well as the angle enclosed by the incident outcome case effect zone (θ) increase. After 1 s, most of the largest droplets settle while the droplet nuclei continue to move with a fate dependent on ventilation and air recirculation. Consequently, once all the larger droplets are deposited, the effect zone will be supported by the droplet nuclei and the contagion risk index IR will coincide with the airborne transmission frequency. It is worth pointing out that in the calculations presented from now on, the highest IR will be considered, i.e. in the first moments from the start of the sneeze. The IR was computed by including all the three active contagion routes.

Fig. 11

Time sequence of the iso-risk curves ((x–y) view) at 5, 50, 100, 150, 200, 340 ms and 1 s (a,b,c,d,e,f,g respectively).

Fig. 12 shows the fault tree developed in this work to calculate the frequencies. It includes all the ways of transmission (direct and indirect droplet transmission and airborne transmission) with the relative basic events. A list of basic events is given in Table 3 . Basic events for each transmission route were set by following the literature indications showed in the Introduction (Section 1).

Fig. 12

Fault tree for “COVID-19 infection” top event.

Table 3

Basic events, corresponding description, failure rates and probabilities.

Name	Description	μ (failure/day)	P	Comments
Basic 1	PPE absence	0.05	0.05	1 out of 20 people does not use PPE
Basic 2	Close contacts with N people	Variable	Variable	From 1 to 100 close contacts/day
Basic 3	SARS-CoV- 2 positive people	Variable	Variable	Fig. 14
Basic 4	Surface contact	2.30	0.90	Probability: 90%
Basic 5	Hand-to-face contact	144	1.00	Once every 5 min in 12 h
Basic 6	Missed handwashing	0.05	0.05	1 out of 20 people does not wash his/her hands
Basic 7	Indoor environment	0.01	0.01	15 min/day
Basic 8	Low ventilation rate	0.7	0.50	Probability: 50%

It is worth noting that the probabilities of failure strongly depend on the situation being analysed and in particular on the number of infected people in the geographical area analysed with respect to total people. In order to quantify the frequency of each way of transmission, we set constant values for failure rates (and relative probability) based on considerations given in Table 3 in the Comments column.

To quantify the basic event “SARS-CoV- 2 positive people” (Basic 3), we considered the total COVID-19 trend in Italy as function of time (Fig. 13 a). The day zero corresponds to 25 February while the 170th day is 13 August with 325 and 252235 cases, respectively. The total COVID-19 cases data can be easily fit through a logistic function with 3 parameters (R2 = 0.9938) as reported in the followingwhere a is the horizontal asymptote (plateau) equal to 237855.06, t0 is the initial time equal to 40.54 and b is a parameter that affects the rate of exponential growth and is equal to 12.33. This curve can be used to visualize and monitor the progress of infections.

Fig. 13

Total cases (a) and active cases (b) profiles in Italy as function of time for the period from 25 February until August 16, 2020.

The infection propagation could be divided into four phases:

initial phase or initiation: slow (linear) increase of the case number.

propagation phase: exponential increase up to the inflection point that represents the maximum rate of cases number rise. This point is here defined as the turning point.

generalised spread phase: uncontrolled propagation of the virus.

termination phase (extinction or “freezing”): slowdown of the virus spread.

By removing deaths and recoveries from total cases curve, we get “currently infected cases” or “active cases” (cases still waiting for an outcome) curve (Fig. 13b). Dividing the active cases number by the total Italian population (60.36 million in 2019), we obtained the Basic 3 probability. In Fig. 14 the Basic 3 probability is shown as function of time. It is completely different in the 4 phases, reaching its maximum value in phase III.

Fig. 14

Basic 3 probability trend in Italy as function of time for the period from 25 February until August 16, 2020.

In Fig. 15 the Top event frequency (i.e. the IR index) is plotted as function of time, at different values of the Basic 2 frequency. Basic 2 represents the number of close contacts, its frequency is calculated as the close contacts per day (CC/day). Increasing the Basic 2 frequency by two orders of magnitude (from 1 to 100 CC/day), the Top event frequency increases from 1∙10−4 to 1∙10−2 CC/day.

Fig. 15

Top event frequency trends parametric in the Basic 2 frequency (black = 1 CC/day; light green = 2 CC/day; red = 5 CC/day; blue = 10 CC/day; light blue = 20 CC/day; light grey = 100 CC/day) in Italy as function of time for the period from 25 February until August 16, 2020. (For interpretation of the references to colour in this figure legend, the reader is referred to the Web version of this article.)

The evaluation of the acceptable risk index is very difficult in the risk analysis of chemical processes and it is well known that this value changes from country to country (American Petroleum Institute (API), 1995; Center for Chemical Process Safety (CCPS), 1999). In the field of public health, the problem persists. Inspired by the indication by the Health and Safety Executive (HSE) for the levels of risk, in terms of the probability of an individual dying in any one year, we set the IR value (Health and Safety Executive, 2021). By setting as maximum tolerable risk threshold value IRmax = 1∙10−4 events/day (Fig. 16 ), it was possible to derive the value of Basic 2 probability as the number of close contacts allowed to a single individual during each phase of the virus propagation so that the probability of contagion propagation is acceptable. It is worth noting that this high IR value can be tolerated for limited period for people equipped with appropriate PPE. In this work, the “maximum tolerable limit” was chosen to easily show the application of the methodology, finding the maximum number of allowable contacts. Surely, once the occurrence probabilities of basic events have been suitably established, also thanks to the help of experts, an appropriate acceptable level of risk can be properly used.

Fig. 16

Zoom of top event frequency trends parametric in the Basic 2 frequency (black = 1 CC/day; light green = 2 CC/day; red = 5 CC/day; blue = 10 CC/day; light blue = 20 CC/day; light grey = 100 CC/day) in Italy as function of time for the period from 25 February until August 16, 2020. The fixed risk threshold is also reported. (For interpretation of the references to colour in this figure legend, the reader is referred to the Web version of this article.)

From this analysis, we may identify the maximum number of close contacts in each phase of the virus propagation.

Initiation phase: each individual can have maximum 20 close contacts.

Propagation phase: only one close contact is allowed. This results in the need of lockdown.

Generalised spread phase: only one close contact is allowed. In this phase, the single meeting shows the greatest IR value, almost comparable to the risk threshold. Also in this phase, the lockdown seems to be of fundamental importance.

Termination phase: it starts with maximum 2 close contacts per day and it progressively increases to 8 close contacts per day on the 170th day.

To assess the effect of the each transmission route in the IR value, we considered the highest IR values for each phase (Table 4 ). From the results, we may observe that the highest frequency is found in the case of direct droplet transmission and it is the consequence of the minimal cut sets evaluation and more precisely of the basic event number required (i.e., cut set order) to reach the top event. A quantification of the impact of the single basic events on the top event was carried out through a sensitivity analysis. The sensitivity analysis was reported in detail in Supplementary Materials (Section S5) and as expected, since the Minimal cut set of the minimum order (3) consists of Basic 1+Basic 2+Basic 3, the top event frequency is very sensitive to the changing probabilities of these basic events.

Table 4

Frequency values for each transmission way and for top event.

Phase	t (days)	Allowed close contacts/day	fd,d(failure/day)	fd,i(failure/day)	fa(failure/day)	IR=fTop Event(failure/day)
I	13	20	1.00E-04	1.19E-07	1.47E-08	1.00E-04
II	41	1	7.93E-05	1.40E-06	1.72E-07	8.09E-05
III	53	1	9.01E-05	1.59E-06	1.96E-07	9.19E-05
IV	170	8	8.65E-05	1.90E-07	2.35E-08	8.67E-05

In the following, the step of the developed methodology are summarised.

The iso-risk curve size as well as the droplets dimensions strongly depend on the time elapsed since the beginning of the respiratory symptom and they can be computed by CFD simulations of the bio-aerosol near-field dispersion.

To calculate the value of the individual risk, a fault tree can be developed. The basic event frequencies as well as IR depend on the geographical location and on the phase of the spread.

Depending on the phase of the contagion propagation, through the calculation of the Top event frequency (IR is the same) and by varying the number of close contacts per day, it is possible to identify the point of the spreading curve at which activate measures to keep the risk of contagion under control. In particular, by fixing the acceptable risk index value, an estimate of the maximum allowed close contacts per day can help policymakers to make choices on lockdown activation or interruption and company managers to organize work shifts and smart working.

From the results of the fault tree analysis and the CFD simulations, the following criteria could be confirmed:

Direct transmission route. As confirmed by minimal cut sets and Table 4, it is the main route of virus transmission. It occurs via the emission of large droplets in open and indoor environment. As shown by CFD results, it may be prevented by keeping 1 m distance and by using PPE, reducing Basic 1 probability.

Airborne transmission route. It occurs via the emission of droplet nuclei in indoor environment. It is worth noting that the effect zones extension mainly depends on the distance covered by these droplets. It may be prevented by using PPE and by ensuring a correct ventilation, reducing Basic 1 and Basic 8 probabilities. It is worth noting this route is the only one that remains active after few seconds from the respiratory symptom.

Indirect transmission route. It occurs via the deposition of large droplets. It may be prevented by using PPE and handwashing. These prevention measures can guarantee the reduction of Basic 1 and Basic 6 probabilities.

Conclusions

Coupled CFD simulations and probabilistic analysis allows the quantification of the risk indexes for infection transmission.

CFD simulations show that after 1 s from the beginning of a sneeze, most of the particles with dimension higher than 5–10 μm is settled over the surface due to gravity. Conversely, aerosol nuclei (particles of dimension lower than 5–10 μm), remain suspended in air. From the CFD simulations, it has been possible to evaluate the effect zone dimension which is equal to 120 cm after 1 s and mainly consists of aerosol nuclei. In the case of air flow at velocity equal to 1 m/s, if it is backwall with respect to the infected person, it increases the distance. However, larger droplets move behind the aerosol front, reaching a plateau value of covered distance in a very short time.

The iso-risk zone dimension as well as the droplets dimensions strongly depend on the time elapsed since the beginning of the respiratory symptom and they can be computed from the CFD simulations of the bio-aerosol near-field dispersion. Within the impact zone, the individual risk IR is equal to the sum of direct droplet (f ), indirect droplet (f ) and airborne (f ) transmissions frequencies. It is worth noting that once all the larger droplets have been deposited, the effect zone will be supported by the droplet nuclei and the contagion risk index IR will coincide directly with the airborne transmission frequency. These frequency values have to be set thanks to the help of further professionals as experts of infection transmission, occupational medicine, biological risk in workplaces, risk assessment and management of indoor environment since they depend on the geographical location and on which phase of the virus spread is being analysed. Through the calculation of Top event frequency (or risk index, IR) varying the number of close contacts per day, we derived a simple methodology to keep the risk of contagion under control. In particular, an estimate of the maximum allowed close contacts per day can help policymakers to make choices such as a possible lockdown or the interruption of the same and company leaders to organize work shifts, manage smart working, etc.

Declaration of competing interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.