## Abstract

**Background** Natural disasters and infectious diseases are global issues, resulting in widespread disruption to human health and livelihood. At the scale of a global pandemic, the cooccurrence of natural disasters is inevitable. However, the impact of natural disasters on the spread of COVID-19 has not been extensively evaluated using agent based epidemiology models.

**Methods** We create an agent-based epidemiology model based on both COVID-19 clinical and epidemiological data and geographic data. We first model 35 scenarios with varying natural disaster timing and duration for a COVID-19 outbreak in a theoretical region. We then evaluate the potential effect of an eruption of Vesuvius volcano on the spread of COVID-19 in Campania, Italy. Our objective is to determine if the occurrence of a natural disaster during the COVID-19 pandemic is likely to increase infection cases and disease related fatalities.

**Results** In a majority of cases, the occurrence of a natural disaster increases the number of disease related fatalities. When the natural disaster occurs at the beginning of the outbreak within a given region, there is little to no increase in the progression of disease spread. However, the occurrence of a natural disaster close to the peak of infections may increase the number of fatalities by more than five-fold. In a theoretical test case, for a natural disaster that occurred fifty days after first infection case, the median increase in fatalities is 2%, 59%, and 180% for a 2, 14, and 31-day long natural disaster respectively, when compared to the no natural disaster scenario.

**Conclusion** We propose that the compound risk from natural disasters is greatest in the case of already widespread disease outbreak. The key risk factors for increase in spread of infection and disease related fatalities are the timing of the natural disaster relative to the peak in infections and the duration of the natural disaster.

## Introduction

As of September 2020, COVID-19 has spread to over 200 countries, infected more than 28.9 million, and killed over 920,000 people [1]. Natural disasters can threaten the measures in place to reduce disease transmission [2]. Understanding how disease spread can be worsened by natural disasters will aid preparedness and response planning. Although a small increase in risk of epidemic outbreak following natural disasters has been identified [2, 3], very little is known about the spread of infection following a natural disaster during a pandemic such as COVID-19.

Numerical modelling is an important tool for understanding and projecting the spread of infectious diseases [4], including COVID-19 [5, 6, 7]. A common approach is to divide a population up into susceptible, infected, recovered and deceased individuals [4, 8]. More sophisticated models incorporate a larger number of different states (e.g. asymptomatic or seriously ill individuals; [5]), and may account for other real world complications (e.g. age, vaccination campaigns, etc.).

The time varying proportion of each state may be calculated deterministically by solving a series of ordinary differential equations [4, 5]. This approach can accurately portray the spread of infectious diseases on a large scale [4, 5]. The deterministic approach is ill suited for simulation of small population sizes, heterogeneous populations and threshold behaviour [7]. A stochastic, agent based model structure may also be used [8, 9], in which case the variation between runs allows for quantification of uncertainty [6]. Agent based models are are well suited to represent heterogeneous or complex populations and spatial variability [8, 9].

## Methods

We use geoSIR, a geospatial, agent based model structure to account for the large scale movement of individuals within the model space from a natural disaster evacuation [8, 9]. Each cell may be in one of 7 states: susceptible (*S*), infected and mildly ill (*I _{m}*), infected and severely ill (

*I*), infected and asymptomatic (

_{s}*I*), recovered (

_{a}*R*), deceased (

*D*) or empty (

*X*). Empty space allows us to account for other geographic factors such as population density. Individuals may also self isolate or be quarantined. We base the definition of our states and model disease characteristics on COVID-19 data [10, 11, 12, 13, 14, 15].

_{e}In the initial state, all individuals are susceptible. The spread of the disease is initiated with the arrival of infected travellers. Two parameters may be adjusted to control the seeding of the disease: the number of incoming travellers (*Tr _{in}*) and the probability that each traveller is infected (

*P*). The number of imported cases per day

_{TrI}*C*is:

_{I}For the majority of the world, new COVID-19 cases were initially imported by a small number of travellers before spreading throughout the community [16]. As individuals are infected, the number of susceptible individuals evolves as:

With *∂I* the number of new infected individuals in day *t*. The total number of new infected individuals depends on the total number of infected individuals in the previous day , and a disease spread parameter *α*:

With *α* given by:

With *n _{ci}* the number of close encounters,

*P*the probability of infection for close encounters,

_{ci}*n*the number of distant encounters, and

_{di}*P*the probability of infection for close encounters. Close encounters can represent sustained contact with another individual within a household or among close friends. Distant encounters can represent an individual’s risk of catching an infection through indirect contact such as visits to supermarkets, fomites transmission or whilst commuting. The details of these two parameters can be modified to account for specific modes of transmission. We use data from two case studies of COVID-19 transmission in Wenzhou, China [17] and Bavaria, Germany [15] to estimate parameters. We also ensure that resulting reproductive number

_{di}*R*

_{0}(as calculated from equation 6) are consistent with the range estimated for COVID-19 [10, 14]. We use , and . The number of close and distant encounters are comparable to those used for an agent based simulation of influenza outbreaks in New York City [8].

The number of encounters and probability of infection transmission can be reduced through lockdowns and social distancing measures. Equally, they can be increased by large gatherings or movements of people. *ρ* is a measure of the local proportion of individuals that are susceptible:

*S _{I}, I_{I}, R_{I}, D_{I}*,and

*X*being respectively the number of susceptible, infected, recovered, dead and empty spaces in the relevant local travel region

_{I}*I*. The radius

*I*is larger for distant encounters than for close encounters in order to account for spread outside of an individual’s direct social circle. The mean reproductive number

*R*

_{0}(

*t*) of the disease at time

*t*is calculated as:

With the mean incubation time and the mean duration of symptoms. We use days, and days [12, 18, 13]. A bulk *R*_{0} of the disease is calculated at each time step, however internal variability within the simulated region can be high. This can result in periods of time in which the regional average reproductive number is less than 1, however the infection still spreads rapidly in local hotspots. Estimates of reproductive number are influenced by the local mitigation practices, environment and individual behaviour and comparison between regions is complex [19].

The daily change in total number of infected individuals is:

Where *∂R* is the new number of recovered individuals in a given day, and *∂D* the number of new deaths. includes asymptomatic, mild and severe cases in varying stages of the disease.

A proportion *P _{A}* of infections begin as asymptomatic, and the remaining 1

*- P*begin as mild. The proportion of asymptomatic COVID-19 patients has been estimated at around 20% through airport screenings [20]. We use

_{A}*P*= 0.2

_{A}*±*0.05. Asymptomatic infected individuals may infect others, but do not change state until they recover. Mild cases have a

*P*chance of transitioning to severe cases over the course of their illness (

_{S}*d*). We use a definition of severe cases that combines the ‘severe’ and ‘critical’ categories from Wu and McGoowan (2020), resulting in

_{sym}*P*= 0.19

_{S}*±*0.02 [11]. Each severely ill individual has a probability

*P*of dying. The expected probability of dying for any infected individual, termed infection fatality rate

_{D}*P*is thus:

_{IFR}Many of the regions worst affected by COVID-19 experienced their highest death rates at times when hospital capacity was reached or exceeded [21]. A hospital capacity term *H _{C}* allows for the simulation of scenarios in which hospital capacity becomes overwhelmed. This may represent number of intensive care unit (ICU) beds, number of ventilators or some other measure of how many severely ill patients can be treated at any given time. We also include a term that describes the probability of a seriously ill individual being hospitalised

*P*. The death rate is for seriously ill individuals is given by:

_{H}With *P _{D}*(0) the base death rate for seriously ill patients,

*γ*a scaling factor for increases in death rate and

A death rate of severely ill patients is capped at *P _{D}*(

*max*).

We use COVID-19 infection mortality rates from from a compilation of Chinese data [11,12], calculated to be 2.3% and 1.4% respectively. We define *P _{D}*(0) = 0.12

*±*0.02, giving an approximate infection mortality rate of 1.8%. We define hospital capacity according to local ICU bed capacities.

If an infected person does not die after the duration of their symptoms (*d _{sym}*), they recover. This model does not account for possible re-infection, and recovered individuals remain immune for the duration of the model run. The simulation of testing is possible in this model (see user manual), but we make the simplifying assumption that testing has not had a major impact on COVID-19 mitigation in the regions considered.

The model simulation grid may be built according to real-world geospatial data. We vary the initial proportion of susceptible individuals *S* and empty space *X _{e}* to account for differences in population density. High and medium natural hazard areas are also inputted. These determine the areas in which evacuation may be necessary. For Vesuvius, high and medium hazard zone masks are created based on the ‘red’ and ‘yellow’ zones defined in the most recent Vesuvius National Emergency Plan [22]. Individuals living in high or medium hazard zones are tagged, and retain more social interactions (e.g. a higher number of close or distant encounters) than other individuals following evacuation. Social distancing is challenging both during evacuation [23].

Information on changes in mobility following natural disasters is limited. Lu et al. (2012) analysed mobile phone data trends following the 2010 Haiti earthquake. They found that the num-ber of people travelling more than 20 km increased 86% in the day following the earthquake, and that increase in movement lasted two to three weeks[24]. These movements led to a 23% decrease in the population of capital city Port-au-Prince, equivalent to close to half a million outward journeys [24]. Lockdowns will be temporarily interrupted during evacuations, and many individuals will be unable to maintain social distancing or lockdown rules [2]. To account for this, we increase the number of encounters by a factor of 3 for individuals within the high and medium risk zones, and by a factor of 1.5 in the low risk zone. We model two different scenarios: an idealised region containing a small city and concentric geological hazard, and a theoretical eruption of Vesuvius during the COVID-19 pandemic in Campania, Italy.

### Theoretical region

In order to quantify disease spread during natural disasters in general, we create a non-specific experimental region. This region includes an idealised town or city with high population densities that radially decay away from the town centre. Initial effective population densities vary from 1 for the centre of town to 0.3 in rural areas. A concentric geological hazard is located within this study area, with the medium hazard region overlapping moderate population density outskirts of the city and rural areas. Hazard *N _{H}* is divided into high and medium hazard zones, with the high hazard zone entirely enveloped by the medium hazard zone. The theoretical region has a population of 100,000 individuals, and a total number of cells equal to . We find that a population of 100,000 is sufficiently large to accurately reproduce disease spread, yet low enough to remain computationally efficient. We conduct 100 different runs for each scenario, seeding the Mercenne Twister with system time to ensure independent runs.

### Campania

Campania, Italy has a population of 5.8 m inhabitants, includes Italy’s third largest city (Napoli), and to active volcano Vesuvius. With eight Plinian eruptions in the last 25,000 years and 32 confirmed eruption in the last 1000 years [25], further eruptions in the near future are likely. The potential damage associated with a large eruption of Vesuvius is severe. However, comprehensive risk mitigation strategies have been derived for the volcanic hazards [26]. The key mitigation strategy is a timely evacuation [26]. We investigate the effect of the evacuation of a large area around Vesuvius during an infectious disease outbreak similar to COVID-19, with the assumption that no one is directly harmed by the volcanic eruption. This disaster response results in widespread population displacement - a key risk factor in disease spread [2].

We use published data for the spread of COVID-19 to parametrise our model, and compare the outputs of the most realistic scenario (no eruption, lockdown initiated) to real world daily infection and fatality data. We use the number of ICU beds in Campania region (427) as a measure for hospital capacity [27]. Due to the large population size in this scenario (over 5 million individuals), we adapt the code to run on a 24 core Haswell E5-2680v3 processor node of the Minnesota Supercomputing Institute.

## Results

In the theoretical region, we perform 100 model runs for 35 different scenarios. These scenarios cover various permutations of lockdown, natural disaster timing and natural disaster duration. Both the number of infections and number of deaths are higher when a natural disaster occurs during the infection outbreak (Fig. 1). The scale of this increase depends on the relative timing of the natural disaster and peak in infection cases, and on the duration of the disaster (Figure 1). In the case with no natural disaster, the median number of infections is 4,900 (IQR 4,140-5,640) cases per 100,000 and the median number of deaths is 66 (IQR 53-76) per 100,000. The increase in deaths remains low when the natural disaster occurs during the onset of the outbreak (day 1, median 57, IQR 39-100 deaths per 100,000) and after the initial infection peak has subsided (day 200, median 73, IQR 53-92.5 deaths per 100,000 respectively). The number of infections and deaths is highest when the natural disaster occurs close to the peak of the outbreak. In this case, the number of infections are 153% higher and deaths are 602% higher when compared to the scenario with no natural disaster (NDt = day 20, median 12,390, IQR 7,100-16,820 cases per 100,000; median 463, IQR 126-847 deaths per 100,000).

For a given natural disaster timing (50 days), long disruptions cause more infections and deaths than short ones. A 3 day long natural disaster results in a similar number of deaths (median 66.5, IQR 55-80.5 per 100,000) as the case with no natural disaster. A 31 day long natural disaster, however, causes and 180% increase in median number of deaths (median 183.5 IQR 149.5-3172 deaths per 100,000). In the case where no lockdown is instigated, around 80% of the population is infected and death toll is extreme (median 3028, IQR 2674-3259 deaths per 100,000).

For Campania, we first model the COVID-19 outbreak in a scenario similar to reality (lock-down initiated, no Vesuvius eruption). Our model predicts a median of 580 (IQR 553-641) infections per 100,000, and a median of 8.5 (IQR 7.9-8.9) deaths per 100,000. Model infections are lower, but consistent with Vollmer et al., (2020)’s estimate of 590-950 cases per 100,000 [28] and model deaths are consistent with real world data (7.5 deaths per 100,000 as of July 2020) [1].

We then model five counterfactual scenarios, four of which include an eruption of Vesuvius (Figure 2) during Campania’s COVID-19 outbreak. The increase in total number of deaths remains low for an early (Day 2, mean 9.6, IQR 7.6-13.1 deaths per 100,000) or post infections peak (Day 100, median 10.9, IQR 10-12.1 deaths per 100,000) Vesuvius eruption. The median deaths rise to 68.5 (IQR 12.1-133.3) per 100,000 when a Vesuvius eruption occurs close to the infection peak (Day 50). In this case, the deaths per 100,000 are comparable to the states in Northern Italy worst hit by COVID-19 (168, 95 and 44 deaths per 100,000 in Lombardy, Piedmont and Veneto respectively as of September 2020) or other hard-hit regions such at New York State, USA (168 deaths per 100,000 as of September 2020) [1]. For an eruption occurring close to the peak in infections at day 50, there is a large variability in possible outcomes. The median death rate is 8 times higher than with no eruption (68.5 deaths per 100,000). However, in 32% of runs the death rate is no more than double that of the no eruption scenario, while in 12% of runs the death rate is more than 25 times higher.

Italy was one of the most hard hit countries during the COVID-19 pandemic [5], with over 33,000 deaths and 233,000 confirmed cases as of June 2020. Campania had more than 4800 confirmed infections and 400 deaths [1]. Although the probability of Vesuvius erupting during the COVID-19 pandemic is low, the wealth of volcanic evacuation plans and COVID-19 data make it a valuable case study. The exact nature of the geological hazard is not of primary importance, rather the parameters of the evacuation and increased human contact following influence the progression of the infectious disease outbreak. The Campania results may be used to understand the effect of other disasters requiring widespread evacuation during a pandemic.

Natural disasters have already occurred during the COVID-19 pandemic [29, 30], and will occur during future infectious disease outbreaks. Cyclone Amphan made landfall in Bangladesh and West Bengal, India on May 20, 2020, prompting widespread evacuations. The number of new COVID-19 infections in the first week of June was 3.5 times higher in West Bengal and 4.8 times higher in Bangladesh compared to the same period in May [1]. The number of COVID-19 related deaths rose by a similar percentage. Further research is required to determine whether Cyclone Amphan played a role in this increase in infections and deaths, which at this stage cannot be attributed to any particular cause.

The majority of previous studies on the relationship between infectious disease outbreak and natural disasters have focused on whether natural disasters initiate new outbreaks. In a small number of cases, a disease outbreak has followed a natural disaster [2]. However, a review of the topic concluded that “the risk for epidemics after a geophysical disaster is very low” [3]. Our results do not contradict this, but highlight the previously overlooked extreme risk in cases of already widespread infection.

The possibility of natural hazards interacting with COVID-19 has previously been raised [29, 30]. Quigley et al., 2020 use a phenomenological model to investigate the impact of natural disaster occurrence during the COVID-19 pandemic. They also find that this results in a larger total number of infections, although the simplicity of their model precludes detailed interpretation [29]. Phillips et al., 2020 propose that climate change is increasing the magnitude and frequency of climatic hazards, therefore raising the risk of natural disasters co-occurring with disease outbreaks [30]. Overall, our results show that the co-occurrence of a natural disaster increases COVID-19 spread. Furthermore, we demonstrate that the timing and duration of the natural disaster are key risk factors for increase in infections. Close links between the epidemiology response and natural disaster response communities are necessary for the formulation of timely risk assessments.

## Data Availability

All relevant data is available upon contact with the corresponding author.