High-resolution modelling of atmospheric dispersion of dense gas using TWODEE-2.1: application to the 1986 Lake Nyos limnic eruption

Atmospheric dispersal of a gas denser than air can threat the environment and surrounding communities if the terrain and meteorological conditions favour its accumulation in topographic depressions, thereby reaching toxic concentration levels. Numerical modelling of atmospheric gas dispersion constitutes a useful tool for gas hazard assessment studies, essential for planning risk mitigation actions. In complex terrains, microscale winds and local orographic features can have a strong influence on the gas cloud behaviour, potentially leading to inaccurate results if not captured by coarser-scale modelling. We introduce a methodology for microscale wind field characterisation based on transfer functions that couple a mesoscale numerical weather prediction model with a microscale computational fluid dynamics (CFD) model for the atmospheric boundary layer. The resulting time-dependent high-resolution microscale wind field is used as input for a shallow-layer gas dispersal model (TWODEE-2.1) to simulate the time evolution of CO2 gas concentration at different heights above the terrain. The strategy is applied to review simulations of the 1986 Lake Nyos event in Cameroon, where a huge CO2 cloud released by a limnic eruption spread downslopes from the lake, suffocating thousands of people and animals across the Nyos and adjacent secondary valleys. Besides several new features introduced in the new version of the gas dispersal code (TWODEE-2.1), we have also implemented a novel impact criterion based on the percentage of human fatalities depending on CO2 concentration and exposure time. New model results are quantitatively validated using the reported percentage of fatalities at several locations. The comparison with previous simulations that assumed coarser-scale steady winds and topography illustrates the importance of highresolution modelling in complex terrains.


Introduction
The atmospheric dispersion of gases (of natural, accidental or intentional origins) can be very hazardous to life and the environment.In industry, historic examples of tragic accidents include the dioxin release in Seveso (Italy, 1976), the methyl isocyanate in Bophal (India, 1984), or the petroleum gas explosions in Mexico City (Mexico, 1984), among several others (e.g.Britter, 1989).From the point of view of natural hazards, cold CO 2 gas released from natural Earth degassing can be of concern because, being denser than air at standard temperature and pressure, CO 2 gas clouds hug the ground and spread downslopes governed by local topography and near-surface winds.Under particular conditions (e.g.atmospheric stability), gas can accumulate in topographic depressions reaching concentration levels toxic for humans and animals.For example, many gas manifestations in Central and Southern Italy, characterised by persistent CO 2 emissions, have caused several periodic accidents (Chiodini et al., 2004).Another case of natural origin, much rarer than diffuse emissions but potentially much more hazardous, are limnic eruptions triggered during overturning of CO 2 -rich volcanic lakes (e.g.Zhang, 1996), such as the eruptions that occurred at lakes Monoun andNyos (Cameroon) in 1984 and1986, respectively.The most tragic event occurred on 21 August 1986 at Lake Nyos, when a CO 2 gas cloud spread across the surrounding valleys, suffocating inhabitants and animals.On the 30th anniversary of this disaster, the Commission on Volcanic Lakes of the International Association of Volcanology and Chemistry of the Earth's Interior (IAVCEI) organised the workshop "CVL9-Cameroon" (Yaounde, 14-20 March 2016) to gather scientific experts from several disciplines with the objective of reviewing the groundbreaking research advances in that field and discuss future road maps for research and risk assessment and mitigation (Rouwet et al., 2016).One of the outcomes from this workshop was that, despite the successful mitigation measures taken at these two lakes (e.g.installation of degassing pipes, relocation of local communities to safer settlements) reducing the level of hazard notably, there is not yet a quantitative hazard assessment from volcanic lakes in Cameroon or elsewhere in Africa (e.g.Lake Kivu in Democratic Republic of Congo and Rwanda).In this regard, dense gas dispersal simulations and characterisation of eruption/gas emission scenarios constitute the backbone of a probabilistic quantitative hazard assessment.
Few studies have been conducted to simulate atmospheric dispersion of dense CO 2 clouds from natural Earth degassing (e.g.Pierret et al., 1992;Costa et al., 2005Costa et al., , 2008;;Chiodini et al., 2010;Costa and Chiodini, 2015;Burton et al., 2017).In particular, Costa and Chiodini (2015) have recently used the TWODEE-2.0model (Hankin and Britter, 1999;Folch et al., 2009) to simulate the CO 2 cloud from the 21 August 1986 Lake Nyos limnic eruption by considering various scenarios for cloud volume and eruption duration estimated from previous studies (Evans et al., 1994;Kanari, 1989;Tuttle et al., 1987).For their simulations, Costa and Chiodini (2015) considered a digital terrain elevation model with a resolution of 90 m.Results evidenced the strong effect of the topography and local wind field, leading to different gas flow patterns across the different valleys of the complex topography of the area.Despite some scenarios could reproduce lethal CO 2 concentrations in many of the sites where fatalities did actually occur, these simulations were unable to capture the full dispersal pattern given some model limitations including, but not limited to, steady winds and gas emission rates or insufficient accuracy of near-surface winds.This later aspect is particularly critical on very complex terrains, where local-scale (tens of metres) wind patterns and the resolution of the digital elevation model (DEM) can have a strong influence on model results.
Based on all these previous considerations, the objective of this paper is to address high-resolution (tens of metres) numerical modelling of CO 2 atmospheric dispersal in complex terrains.To this purpose, we introduce a methodology for local wind field characterisation based on transfer functions that couple mesoscale numerical weather prediction (NWP) models with a microscale Reynolds-averaged Navier-Stokes (RANS) computational fluid dynamics (CFD) model for the atmospheric boundary layer.The resulting high-resolution time-dependent wind field is given as input to the TWODEE shallow-layer gas dispersal model to simulate the evolution of gas concentration with time at different heights above the terrain.In addition to coupling with mesoscale/microscale 4-D meteorological models, the TWODEE-2.0code has been improved (updated to version 2.1) in order to overcome a few limitations, including the option to easily describe timedependent gas sources, heterogeneous terrain roughness, and assessment of the impact through a novel probabilistic approach that estimates the percentage of fatalities depending on gas concentration and exposure time (dose) at nearground levels.The new version of the code is applied to reconstruct gas dispersion from the 1986 Lake Nyos event compared with previous simulations (Costa and Chiodini, 2015), and results are used to illustrate the model gain on the time evolution of the gas source.
In this paper, Sect. 2 overviews the main characteristics of the TWODEE model and the modifications introduced to account for heterogeneous high-resolution 4-D wind fields, time-dependent source term, and terrain variable roughness.The new probabilistic impact criterion for TWODEE model validation is also presented.Section 3 summarises the events occurred at Lake Nyos and surroundings during 21 and 22 August 1986 and the results (and limitations) from previous simulation studies (Costa and Chiodini, 2015).A novel strategy for high-resolution local wind field characterisation on very complex terrains is presented on Sect. 4. The methodology uses the concept of transfer functions to couple the mesoscale Weather Research and Forecasting (WRF) (Skamarock et al., 2008) model with the CFD code Alya (Houzeaux et al., 2009) adapted to atmospheric boundary layer flows (Avila et al., 2013).Section 5 shows the TWODEE-2.1 model results for the 1986 Lake Nyos case of study and discusses how and why the model new features increase the accuracy of simulations.In addition, we also perform a model parametric study varying the source term to constrain its intensity, evolution, and duration.Finally, Sect.6 summarises and discusses the main results and future developments.
2 The TWODEE dense gas dispersal model TWODEE-2 (Hankin and Britter, 1999;Costa et al., 2008;Folch et al., 2009) is a FORTRAN90 code for the atmospheric dispersal of dense gases based on the shallow-layer approach.Under the assumption that h/L 1 (h being the gas cloud depth and L a characteristic length), the 2-D shallow-layer approach allows a compromise between more realistic but computationally demanding 3-D CFD models and simpler 1-D integral models.The TWODEE family models build on the depth-averaged equations for a gas cloud resulting from mixing a gas of density ρ g with an ambient fluid (air) of density ρ a (ρ g > ρ a ).The integration of volume, mass, and momentum balance equations over the mixed cloud depth from the ground to the top of the cloud yields to (e.g.Hankin and Britter, 1999) where h is the cloud depth (defined as the height below which 95 % of the buoyancy is located), ρ is the depth-averaged cloud density, u = (u x , u y ) is the depth-averaged cloud velocity, u a ) = (v a , w a ) is the ambient fluid (air) velocity vector, u e is the ambient fluid entrainment velocity modulus, u s is the gas inflow velocity modulus (source term), e = e(x, y) is the terrain elevation, C D is a drag coefficient, F is the turbulent shear stress force (per unit area), and S 1 ≈ 0.5 and κ are semi-empirical parameters.The terms in the momentum equation (Eq. 3) include the local time derivative, the convective term, the pressure gradient (assumed hydrostatic although the density profile can be non-uniform), the effect of terrain slope, the surface shear stress (depending on the terrain roughness and characterised by the drag coefficient), the force per unit area exerted by turbulent shear stress, and, finally, the leading edge terms that account for interaction among dense and ambient fluids.Given closure equations for the drag coefficient C D , shear stress force F and entrainment velocity u e (Hankin and Britter, 1999;Folch et al., 2009), the set of equations above can be resolved numerically to obtain cloud height and vertically averaged density and velocity depending on terrain, source term (definition of u s ), and ambient fluid velocity (wind field).Although TWODEE is a shallow water model, it can also estimate the vertical density profile from the depth-averaged density ρ assuming an empirical exponential decay (Hankin and Britter, 1999): from which the vertical concentration profile c(z) (expressed in ppm) and the dosage Do(t, z) during a time interval (0, t) can be computed as where c b is the background concentration (in ppm) and n is the so-called toxicity exponent.

TWODEE-2.1
The TWODEE code version 2.0 (Folch et al., 2009)  In the case of WRF, a meteorological pre-processor reads the WRF model outputs for a user-defined time interval (n t WRF output time steps) and performs an interpolation from the WRF native grid (Arakawa staggered C grid in the horizontal, pressure levels for the vertical) to a series of userdefined z cuts on a terrain-following regular grid.In the case of microscale simulations, the meteorological pre-processor reads outputs from the ALYA-CFDWind model (see Sect. 4.2) and, for each WRF time slice, performs the mesoscale-to-microscale downscaling using transfer functions as explained in Sect.4.3.This results in n t downscaled wind fields on the same set of terrain-following z cuts.Whichever the approach considered for the wind field (i.e.mesoscale WRF or WRF downscaled with ALYA-CFDWind), TWODEE-2.1 reads winds at the gridded z cuts and then interpolates to its computational mesh at a user-defined reference height.
4. It can assess the impact of dense CO 2 gas dispersal.According to the Occupational Safety and Health Administration (OSHA), the long-term exposure limit (LTEL) of CO 2 is set at air concentrations of 0.5 % (5000 ppm) for exposures up to 8 h per day, whereas the short-term exposure limit (STEL) is set at air concentrations of 3 % (30 000 ppm) for exposures up to 15 min.However, at higher concentrations/dosage, CO 2 causes several adverse health effects when inhaled.Experience shows that CO 2 air concentrations of around 5 % (50 000 ppm) produce heavy breathing, sweating, quicker pulse, weak narcotic effects, and headaches.Under these concentrations, the exposure time to avoid the development of adverse health symptoms is only a few minutes.For example, a 30 min exposure to 5 % concentration (50 000 ppm) produces intoxication manifested as headaches, dizziness, restlessness, breathlessness, increased heart rate and blood pressure, and visual distortion.At around 10 % concentration (100 000 ppm) humans are affected by respiratory distress, impaired hearing, nausea, vomiting, and loss of consciousness in only 10 min.Finally, CO 2 air concentrations > 15 % (150 000 ppm) are considered lethal causing coma, convulsions, and rapid death.The UK Health and Safety Executive (HSE; http://www.hse.gov.uk)developed an assessment of dangerous toxic substances, including CO 2 , defining the specified level of toxicity (SLOT) and the significant likelihood of death (SLOD) depending on concentration and duration of exposure (Harper, 2011).Considering these values (Harper, 2011), we assume a cumulative normal distribution for the percentage of human fatalities: where P is the probability of death, c is the CO 2 concentration (expressed in %), d is the exposure duration (expressed in minutes), and a i and b i are empirical constants.After calibration of Eq. ( 7) with the HSE tabulated values (assuming SLOT at 3 %), we obtained the following values for the constants: a 0 = 5.056, b 0 = 17.885, c 0 = 0.357, a 1 = 0.662, b 1 = 2.421, and c 1 = 0.354.Results are shown in Fig. 1, plotting the percentage of fatalities (probability of death in %) as a function of CO 2 concentration for different exposure durations.
An impact criterion based on these empirical curves was added in TWODEE-2.1 to compute, at each point of the computational domain, the predicted percentage of fatalities at user-defined heights.
3 The 1986 Lake Nyos event and previous modelling results Lake Nyos (Fig. 2) is one of the ∼ 40 volcanic lakes scattered along the 1600 km long Cameroon Volcanic Line (e.g.Lockwood and Rubin, 1989).This lake became famous worldwide on Thursday 21 August 1986 after the occurrence of the most tragic limnic eruption ever registered.During few (< 5) hours, a huge (0.1-1 km 3 ) CO 2 gas cloud (Tuttle et al., 1987;Kanari, 1989;Evans et al., 1994) released during the Lake Nyos overturning spread downslopes from the lake (1100 m a.s.l.), filling up the underlying Nyos valley and suffocating around 1700 people and 3000 cattle (e.g.Kling et al., 1986).Evidence from eyewitness reports indicated that the cloud was directed primarily W-NW at around 21:00 LT (20:00 UTC) affecting the bottom valleys of Cha and Fang (see Fig. 3), but without causing reported deaths at the latter location (Baxter and Kapila, 1989).Following this initial dispersal phase, the gas cloud direction shifted towards NE, probably as a consequence of a sudden wind veer, filling up the Nyos valley down to the Subum village (∼ 10 km line of sight from the lake), where the largest number of casualties occurred.Deaths in humans and animals (including birds) occurred at a distance of up to 20 km across the main and adjacent secondary valleys (Le Guern et al., 1992).
Le Guern et al. (1992) collected multiple testimonies of survivors that allowed the authors to reconstruct the timeline of the event, locating where casualties occurred, and constructing a map of the areas impacted by the disaster (see Table 1 and Fig. 1 in Le Guern et al., 1992).It should be stressed that the percentages of fatalities reported by Le Guern et al. (1992) resulted from posterior interviews by anthropologists and are subject to large uncertainties related to translation and interpretation, approximate location of sites, and actual number of casualties.Notwithstanding these limitations, data inferred from eye witnesses constitute the only source of information available for indirect model validation given the lack of any wind or gas concentration measurement at Nyos on that time.At present, the level of risk at lakes Nyos and Monoun has decreased notably after the successful deployment of degassing pipes in March 2001 and February 2003, respectively.The progressive gas removal resulted in considerable deepening of the level of gas-rich water in a short period of time (Kling et al., 2005;Kusakabe et al., 2008).However, there is still the recognised need to perform a quantitative CO 2 hazard assessment for several reasons, including the possibility of future gas build-up or a breakthrough of the dam built at the northern shore of the lake.Thus, for example, Aka and Yokoyama (2013) estimated that dam break could cause a sudden drop in lake level by 40 m, followed by decompression and inevitably a new gas burst.Costa and Chiodini (2015) have recently used the TWODEE-2.0model to simulate four different scenarios (released gas mass ranging from 0.29 to 1.95 Tg, wrongly reported in their table as Gg) using a 90 m resolution DEM.Surface wind data resulted from applying the DWM massconsistent pre-processor to a constant wind profile extracted from the closest point of the NCEP/NCAR Reanalysis 1 (Kalnay et al., 1996).As a result of this limitation, none of their simulations was able to capture properly the cloud dispersion pattern, strongly influenced by a sudden wind veering.Figure 4 shows, for illustrative purposes, results for their scenario II, the one that better reproduced the observations (see Fig. 5 in Costa and Chiodini, 2015, for other scenario results).In the following sections, these previous simulations are revisited considering a higher DEM accuracy (30 m instead of 90 m; from ASTER G-DEM, a product of METI and NASA) and high-resolution transient microscale surface winds derived from downscaling the Advanced Research WRF (WRF-ARW) mesoscale winds with the ALYA-CFDWind model.The WRF is a fully compressible, nonhydrostatic mesoscale NWP model and atmospheric simulation system designed to serve both operational forecasting and atmospheric research needs (Skamarock et al., 2008).The model uses finite differences schemes on a staggered horizontal Arakawa C grid and a terrain-following vertical coordinate system to solve the atmospheric flow.Here, the version 3.4.1 of the dynamical solver WRF-ARW was configured with the physical parameterisations and schemes summarised in Table 2. Dudhia (Dudhia, 1989) For the Nyos application, the WRF-ARW simulation starts on Thursday 21 August 1986 at 00:00 UTC lasting 48 h (around 18 h are allowed for model spin-up).Initial and daily (four times) boundary conditions driving WRF come from the NCEP/DOE Reanalysis 2. Figure 5 shows the five domains used, consisting of one parent grid at 81 km horizontal resolution, covering the African continent, and four nested domains at 27, 9, 3, and 1 km horizontal resolutions, all centred around Lake Nyos.The regional synoptic situation on Thursday, 21 August 1986 is illustrated in Fig. 6, showing WRF-ARW results for the first nest domain (27 km resolution) at different times.It can be observed how a lowpressure (< 800 hPa at 2000 m) region and its cyclonic circulation located NE of Cameroon coexists with strong winds from both south and north and a region of weak winds over NW Cameroon at around 16:00 UTC (Fig. 6a).The simulations indicate that this situation lasted for about 6 h until 22:00 UTC, when winds over Nyos became stronger and pointed SE.In contrast, near-surface winds followed a very different pattern (Fig. 7) reflecting the topographyinduced forcing at lower atmospheric levels.However, given the orographic complexity, the WRF-ARW results at 1 km resolution (inner nest) are still insufficient for driving highresolution gas transport simulations, indicating the need for an ulterior downscaling.

Microscale wind modelling using ALYA-CFDWind
ALYA (e.g.Houzeaux et al., 2009;Vázquez et al., 2016) is a high performance computing (HPC) multi-physics parallel solver based on a finite element method able to run with thousands of processors with an optimal scalability.Within this multi-physics general framework, Avila et al. (2013) implemented a solver for the atmospheric boundary layer based on the RANS equations and a κ − turbulence model adapted to atmospheric flows in complex terrains.This model, called ALYA-CFDWind, can handle thermal coupling assuming the Boussinesq approximation although here we constrain to neutral atmospheric stability for simplicity.The ALYA-CFDWind module, originally developed in the context of wind energy, includes Coriolis effects, a consistent limitation of the mixing length, a wall law for atmospheric boundary layers (logarithmic profile depending on terrain roughness and wind friction velocity), automatic meshing, and generation of boundary conditions for atmospheric boundary layer wind profiles over a flat terrain.In order to have consistent inflow boundary conditions (i.e.flat terrain inflow profile) and also to prevent flow recirculation at the outflow boundary, the ALYA-CFDWind computational domain is made of an external flat buffer designed to accommodate the flow, an adjacent transition zone, and an inner higher-resolution zone with the real topography and where the flow is actually computed.
For the Nyos case, the ALYA-CFDWind computational domain consists of an inner zone of 20 × 20 km 2 at 50 m horizontal resolution, a transition zone of 15 km, and a flat buffer zone of 10 km to accommodate the flow.Along the vertical direction, the structured hexahedral grid extends up to 5 km above the terrain with 64 vertical layers growing geometrically in size from 0.5 m at surface to 250 m at the top of the computational domain.The resulting computational mesh (see Fig. 8) has a total number of grid points of about 30 million.The Coriolis force was set to that of a latitude of 6 • N and the maximum mixing length was calculated depending on the wind at top as in Apsley and Castro (1997).As boundary conditions, we prescribed the wind at top (geostrophic wind) to a reference value of 10 ms −1 considering different geostrophic wind directions (sectors)   at 15 • intervals.One ALYA-CFDWind simulation was performed for each reference direction of the geostrophic wind (i.e.24 different runs are necessary to scan all possible geostrophic wind directions).Assuming self-similarity, these reference runs can be scaled to actual velocity and interpolated in direction during the meso-to-micro downscaling strategy (see Sect. 4.3) depending on the mesoscale (WRF) time-dependent geostrophic wind direction and intensity.

Meso-to-micro downscaling strategy
Mesoscale NWP models like WRF-ARW can be used to forecast winds at horizontal resolutions down to around 1 km.This grid resolution may be insufficient to drive subsequent gas dispersion simulations over complex terrains, where sub-grid-scale topographic features can alter nearsurface winds and therefore the resulting gas dispersal pattern.For this reason, a meso-to-micro downscaling strategy may be necessary in order to capture wind-forcing effects caused by the local sub-grid topography.
At present, model downscaling is the subject of active research within the atmospheric and wind engineering communities, with two well-differentiated strategies dominating the scene.On the one hand, statistical downscaling methodologies (e.g.Ranaboldo et al., 2013;Devis et al., 2013;Kirchmeier et al., 2014) build on finding correlations between global/regional model simulations and observations during long periods of time in order to identify patterns used to forecast in time and space by extrapolation.This has proven to be effective for some applications (e.g.wind farm power production forecast) but requires long series of wind observations that, for the application considered here, do not exist.On the other hand, dynamical or physical downscaling (e.g.Castro et al., 2014) couples different nested models so that the outer mesoscale model furnishes boundary conditions to the inner microscale solver at each model time integration step.The inner microscale model can range in complexity from simpler mass-consistent diagnostic models to a CFD solver.This option is clearly more attractive but, apart from the higher computational cost, still presents some challenges related to not well-resolved inconsistencies in the physics of models across scales.On top of this, a pure dynamical downscaling approach becomes computationally prohibitive for hazard assessment purposes, where climatically representative wind series have to be considered (note that this implies thousands of coupled simulations in order to statistically cover all meteorological situations with its associated probability).Given these constrains, we adopt an intermediate physical-statistical strategy based on transfer functions, a concept inspired on methodologies used for microscale wind resource assessment over regional scales (e.g. Sanz Rodrigo et al., 2010).
The idea behind transfer functions is simple.Given a mesoscale wind field (in our case WRF-ARW at 1 or 3 km grid resolution, see Sect.4.1) and a set of microscale wind fields, each characterised by a reference wind direction φ and intensity (in our case 24 ALYA-CFDWind runs at 50 m grid resolution and 15 • geostrophic wind binning, see Sect.4.2), the transfer functions determine a new (downscaled) wind field as where u down (x, y, z, t) is the resulting downscaled wind velocity, f = f (θ, t) is the point-dependent transfer function at time t, u WRF→CFD is the WRF wind velocity interpolated to the (finer) microscale mesh, u θ CFD is the ALYA-CFDWind velocity for a reference wind direction θ , and < u θ CFD > R is the CFDWind velocity for the same direction θ averaged over a radius of influence R (of size similar to that of the WRF cells).Note that, by construction, at each point, the CFD field u θ CFD averaged over its circumference of influence R equals the WRF velocity interpolated at that point.In other words, the resulting downscaled wind field coincides on average (over a WRF cell) with that of the mesoscale model but, at the same time, it has all the local wind fluctuations around this mean (caused by the microscale topographic forcing) which cannot be captured by the (coarser) mesoscale grid.Moreover, the methodology is thought to preserve flow characteristics due to the thermal stratification present in the WRF model.
A central point in this hybrid downscaling methodology is how to determine the reference wind direction θ and then how to link it consistently with the mesoscale flow.For this, we adopt the following strategy to obtain downscaled 2-D fields at given heights z t above the terrain: 1.The first step consists of extracting from the microscale computational domain (ALYA-CFDWind) a series of 2-D plains 2-D at user-defined heights above the terrain (e.g.z t = 2, 10, and 50 m).Each of these 2-D planes is then decomposed in a series of structured patches or segments S ij , allowing for some overlap at the borders of the sub-domains (i.e.∪S ij = 2-D ; ∩S ij = ∅).
In particular, we consider here one squared patch around each WRF mass point with sizes 2dx WRF × 2dy WRF , dx WRF × dy WRF being the WRF cell area (see Fig. 9).2. For each sub-domain S ij , the reference wind direction θ at time t is computed as the averaged WRF velocity over the segment.Note that, for small computational domains or flat terrains, few variations in θ are expected across different segments S ij .However, over large areas or in very complex terrains, synoptic-scale effects may result on variations of θ at different segments.
3. For each sub-domain S ij and time t, the microscale (ALYA-CFDWind) reference solution u θ CFD is determined by performing a linear interpolation between the two reference runs φ 1 and φ 2 that limit the bin direction containing θ (i.e.θ ∈ (φ 1 , φ 2 )).For example, if for a given segment and time one has θ = 5 • , then the binbounding solutions u 0 CFD and u 15 CFD (i.e.precomputed solutions for θ = 0 • and θ = 15 • , respectively) are combined to obtain u 5 CFD .
4. A smoothing operation in the overlap regions between neighbouring sub-domains S ij is performed.This is necessary because different values of θ in adjacent segments can result on discontinuous values of u θ CFD across segments.In particular, we consider a simple weighted interpolation in the regions with overlap. 5. Finally, the transfer functions is applied using Eq. ( 10) to scale the microscale wind modulus and obtain u down .
Figure 10 compares the 3 km resolution WRF results at 10 m above the terrain with the downscaled field at the area of interest around the Lake Nyos.This figure highlights the local information added by the microscale model over mountainous areas and valleys, where WRF shows a smoother behaviour.Unfortunately, no surface wind data existed at that time to validate these results.However, some consistency with reports exists.For example, a strong microscale wind rotation from NE to SW (wind origin direction) is clearly visible along the Nyos valley from 18:00 to 21:00 UTC, in agreement with cloud dispersal reports indicating two differentiated dispersal phases.In order to illustrate this phenomenon, Fig. 11 plots time series of 10 m wind velocity and direction from 16:00 to 24:00 UTC at two locations of special interest, the lake itself and a point at the bottom of the Nyos valley (see Fig. 3).Note how near-surface wind direction changes sharply at both locations between 17:30 (wind from NE) and 19:30 (wind from SW) UTC.A systematic WRF model error phase of about 3 h seems to exist because dispersal reports suggest such a strong wind veering occurring after 22:00 LT (21:00 UTC).In any case, these wind field variations strongly indicate the need for including time-dependent heterogeneous winds for the gas dispersal simulations.Finally, it is worth looking at the atmospheric stability during the time of the event since this parameter also favours the accumulation of dense gas at depressions.Figure 12 plots the static stability versus time, clearly reflecting a unstable to stable transition related to the diurnal cycle occurring at around 16:00 UTC.

High-resolution dispersal modelling results
Once the high-resolution winds (i.e.WRF-ARW winds downscaled with ALYA-CFDWind using transfer functions) have been obtained for the period and area of interest, 10 m wind values every 20 min (i.e. three times hourly) were given to the TWODEE-2.1 model together with the 30 m resolution DEM and the definition of the source term to simulate the CO 2 cloud dispersal.Our simulations assume the eruption started at 17:30 UTC (18:30 LT), i.e. 2-3 h advanced with respect to eyewitness reports.This source term shift in the dispersal simulations was necessary because of the WRF model phase error discussed in Sect.4.3.However, given the large uncertainties in the source term concerning eruption duration, intensity, and evolution of CO 2 mass flow rate with time, we performed a source term characterisation considering different sets of simulations, each set with different source term characteristics (see Table 3): -Eruptions lasting 3 h with a constant CO 2 emission but varying source intensity.This is similar to the scenarios considered by Costa and Chiodini (2015), who ranged the source intensity between 1.4 × 10 5 kg m −2 d −1 (based on Kanari, 1989) and 4.3 × 10 5 kg m −2 d −1 (based on Evans et al., 1994) (Group 1 in Table 3).
-Eruptions lasting 3 h with a CO 2 emission decreasing linearly and considering different time rates and source intensities (Group 2).
-Eruptions lasting 3 h with a CO 2 emission increasing linearly and considering different time rates and source intensities (Group 3).
-Eruptions lasting 3 h with a CO 2 emission decreasing exponentially considering different decay rates and source intensities (Group 4).
-Eruptions lasting 3 h with a CO 2 emission increasing exponentially considering different growth rates and source intensities (Group ).
-Eruptions with an initial phase with low constant CO 2 emission followed by a burst and exponential decay considering different durations for each phase, decay rates and source intensities (groups 6 and 7).In addition, eruption durations for these groups were varied between 2 and 4 h as in Costa and Chiodini (2015).
Results for all these simulations (a total of 62 source term scenarios) were validated on a point-by-point basis comparing the TWODEE-2.1 predicted probability of fatalities with the actual percentage of fatalities reported at 53 locations (Table 1).When defining a metric for the skill score of a simulation, it is important to consider whether observations are distributed well across the possible range of values or not.In our case, for example, 29 observations of 53 (i.e.55 %) and 8 observations of 53 (i.e. 15 %) have 0 and 100 % of observed fatalities, respectively.It means that around 70 % of the values are at the tails of the distribution and, consequently, the evaluation of scores without binning would lead to biassed skills, favouring those scenarios in which the source intensity (emitted mass) is largely under-or overpredicted.In order to prevent this, we adopt binning strategy and evaluate scores first for each discrete bins and then globally by averaging over all bins.Ideally, the number of bins should scale as √ n so that, in our case (n = 53), around seven bins would be recommendable.However, in order to avoid having bins with little or even no observations, we considered the following five-class binning for the percentage of fatalities: the first bin 0-5 % (no significant mortality), the second bin 5-35 %, the third 35-65 %, the fourth 65-95 %, and the fifth 95-100 % (total mortality).This allows us to compute the Pearson cumulative test statistic χ 2 as where n b = 5 is the number of bins, O i is the number of observations (localities) within the ith bin, and M i is the number of model localities laying in the same bin.However, we also compute the total mean absolute error (MAE) as with    K km −1 ) versus time at the same two points than in Fig. 11 (lake in blue and valley in red) according to WRF simulations at 1 km resolution.The potential temperature gradient has been computed using WRF σ levels 1 and 6, at roughly 12 and 100 m above terrain, respectively.Positive values indicate stable atmosphere.Note the transition from unstable to stable stratification occurring after 16:00 UTC.where (MAE) i is the bin absolute error, m the number of observations (localities) in the j th bin, PO is the observed percentage of fatalities, and PM is the modelled probability at the same locality.
As a best-fit criterion we considered the lowest value of χ 2 but trying also to minimise MAE.We found that the higher-score simulations belong to the same group of source term runs -a source with two phases: an initial phase with a constant CO 2 emission followed by a second phase with exponential source strength decay (group 6).This is a consistent result.Table 4 and Fig. 13 summarise the characteristics of the source term for the highest ranked run.The total CO 2 emitted mass is 0.86 Tg, released over 3 h.This value is in good agreement with previous independent estimations, ranging from 0.29 Tg (Evans et al., 1994) to 1.95 Tg (Tuttle et al., 1987).In contrast, the values of χ 2 and MAE are 4.98 and 34.7 %, respectively, giving a very good fit across all bins (see Fig. 14).For comparison, the scenario II in Costa and Chiodini (2015) (see Fig. 4) gives a χ 2 of 7.61 and a MAE of 36.6 %, substantially higher even if the characteristics of the source term are similar to our best-fit runs both in terms of duration (4 h) and total mass (1.33 Tg).These differences can be explained because of the high-resolution time-dependent winds, which are able to capture the wind veering leading to the differentiated cloud dispersal branches (see Fig. 13).This can also be observed by looking at the evolution of concentration with time at different locations.As observed in Fig. 15, the TWODEE-2.1 simulations can reproduce a first W-NW branch affecting the bottom valleys of Cha and Fang (location L34) followed by a gas cloud direction shift towards NE affecting Subum village (location L36).According to simulations, this occurred around 2 h after the eruption onset (i.e.20:30 LT).Considering the 2-3 h shift necessary to correct the meteorological (WRF) phase error, this is in excellent agreement with eyewitness reports.

Summary and discussion
We have developed a high-resolution numerical model for CO 2 atmospheric dispersal in complex terrains by introwww.nat-hazards-earth-syst-sci.net/17/861/2017/Nat.Hazards Earth Syst.Sci., 17, 861-879, 2017 ducing a methodology for time-dependent microscale wind characterisation based on transfer functions that couple a mesoscale numerical weather prediction model with a microscale computational fluid dynamics model for the atmospheric boundary layer.The model was applied to reconstruct the source conditions and the catastrophic CO 2 dispersal at Lake Nyos in 1986.Simulation results were compared with the observed fatalities through a novel probabilistic approach.We found that the new model with high-resolution time-dependent winds shows better agreement with the ob-servations compared with previous simulations, indicating that the model is capable of performing gas dispersal hazard assessment on very complex terrains in the future.The optimal runs (higher scores for percentage of fatalities) shared same source evolution pattern: an initial phase with a low constant CO 2 emission and a second phase with a burst followed by exponential decay in CO 2 mass flux.This suggests additional information about physical processes at Lake Nyos during the 1986 limnic eruption.The scenario is compatible with a gas eruption that started locally in some region of the lake, continuously emitting a relatively low amount of CO 2 near the equilibrium state; then a perturbation, likely due to the gas ascent itself or to other mechanism(s), grew that destabilised the whole lake and triggered the exponential phase by overturn of the different layers in the lake in a supersaturation state.One simple model explaining the exponential decay during the second phase is that the driving force for CO 2 emission is an overpressure in the lake, such as caused by CO 2 supersaturation, and the temporal change in the overpressure is controlled by the emitted CO 2 flux, like it occurs in a volatile supersaturated magma chamber.This model can be simply formulated as (e.g.Huppert and Woods, 2002) and where Q is the emitted CO 2 flux, p is the overpressure, and C 1 and C 2 are two constants.These equations lead to the solution of exponential decay form for Q: where Q 0 is the initial value of Q.Although the precise mechanism in the lake during the limnic eruption is still under discussion, we can infer from the above model that the second phase is simply explained by a relaxation process of supersaturation state in the lake.Although further investigations for the dynamics of lake water are necessary for reproducing the physical process of the limnic eruption in detail, the features of the CO 2 emission obtained in this study may provide strong constraints on physical modelling in the lake.

Figure 1 .
Figure 1.Percentage of fatalities (probability of death in %) depending on exposure duration (in min) according to Eq. (7).Values are shown for different values of CO 2 concentration ranging from 7 to 14 % vol.

Figure 2 .
Figure 2. Map of Cameroon showing the Lake Nyos area (red square).

Figure 3 .
Figure 3. Topographic map of the area around Lake Nyos showing elevation contours (in m a.s.l.).The locations listed in Table 1 are shown by halved circles coloured according the percentage of fatalities (in %) reported by Le Guern et al. (1992).

Figure 4 .
Figure 4. Previous modelling results from Costa and Chiodini (2015) scenario II (4 h of gas emission assuming a constant mass flux of 1.4 × 10 5 kg m −2 d −1 from a diffuse source of 235 × 235 m 2 ).Left: maximum CO 2 concentration (%vol.)achieved at 1 m height.Right: percentage of fatalities (in %) predicted by the model applying Eq. (7) at 1 m height.Halved circles show the actual reported percentages at locations using the same colour scale.

Figure 5 .
Figure 5. WRF-ARW domains for the Nyos run.The model configuration consists of one parent domain (d01) at 81 km horizontal resolution and four nests (d02 to d05) at 27, 9, 3, and 1 km resolution centred at the Lake Nyos (red triangle, not on scale).Colour contours indicate the WRF-ARW model topography at each resolution.

Figure 6 .
Figure 6.Synoptic meteorological situation according to the WRF-ARW simulation showing pressure contours (hPa) and wind vectors (m s −1 ) for domain d02 (first nest, 27 km resolution domain) at 2000 m height above sea level.Results for 21 August 1986 at 16:00 (left), 20:00 (centre), and 24:00 (right) UTC.The red triangle shows the location of the Nyos Lake.For clarity, only a few model wind vectors are shown.

Figure 7 .
Figure 7. Local meteorological situation according to the WRF-ARW simulation showing terrain height contours (m a.s.l.) and wind vectors (m s −1 ) for domain d05 (last nest, 1 km resolution domain) at 10 m above the surface.Results for 21 August 1986 at 16:00 (left), 20:00 (centre), and 24:00 (right) UTC.The red triangle shows the location of the Nyos Lake.For clarity, only few model wind vectors are shown.

Figure 8 .
Figure 8. Detail of the ALYA-CFDWind computational mesh (50 m horizontal resolution) around the Nyos Lake.The bottom-left inset shows the extent of the computational domain composed of three differentiated zones, flat buffer (red), transition (pale blue), and inner domain (green) at 50 m resolution containing the detailed terrain information.The arrows indicate the approximate gas flow path according to observations.

Figure 9 .
Figure 9. Segmentation strategy adopted for downscaling on a 2-D plane 2-D .A patch or segment S ij is defined for each WRF-ARW grid mass point (black squares) containing many ALYA-CFDWind grid points (red small squares).The last plot shows the overlapped region.

Figure 10 .
Figure10.Wind vectors at 10 m above terrain as given by WRF-ARW at 3 km resolution (left column) and the downscaling (right column) for 21 August 1986 at 18:00 (top), 19:00 (middle), and 21:00 (bottom) UTC.For visualisation and point-to-point comparison purposes, WRF results have been interpolated to the CFD mesh and, when necessary, extrapolated below its lower level.The red triangle shows the location of the Lake Nyos.

Figure 11 .
Figure11.Time series of 10 m wind speed and direction at two points located at the Nyos Lake (a, c) and at the bottom of the Nyos valley (b, d).Results from WRF simulations at 1 km resolution (blue lines) and downscaling using ALYA-CFDWind at 50 m (red lines).Wind direction criterion is that of the coming direction; i.e. 0 • indicates wind coming from the N, 90 • coming from the E. Note the sudden short-lived wind veering from NE to SW starting at around 18:00 UTC.

Figure 13 .Figure 14 .
Figure 13.Best-fit simulation results.Left: maximum CO 2 concentration (%vol.)achieved at 1 m height.Right: percentage of fatalities (in %) predicted by the model applying Eq. (7) at 1 m height.Halved circles show the actual reported percentages at locations using the same colour scale.

Figure 15 .
Figure 15.Evolution with time (in h) of CO 2 concentration (%vol.) at different locations for the best-fit simulation.Results at 1 m height.Location L05 is near the Lake, L24 at the Nyos valley, L34 near Cha (W branch), and L36 near Subum (E branch).See Fig. 3 for details.

Table 1 .
List of affected points with coordinates and percentage of observed fatalities.(x,y)coordinates are in UTM zone 32N (datum WGS84).Points starting with an L have been digitalised from LeGuern et al. (1992); points starting with an N were obtained by the authors after interviewing survivors in situ.Comments as reported by LeGuern et al. (1992)in their Table1.

Table 2 .
WRF-ARW model configuration and physical parameterisations used for the 1986 Lake Nyos simulations.

Table 3 .
Source term settings used in the different groups of simulations.Table shows the considered ranges for eruption starting time, total CO 2 emitted mass, CO 2 mass flux, and constant of decay (or growth) for the exponential phase.Static stability (gradient of potential temperature in

Table 4 .
Source term characteristics of the best-fit simulation showing total CO 2 emitted mass, mass flux variation, χ 2 , and mean absolute error (MAE) considering the five-class binning criterion.