Central American is a seismically active region where earthquakes have had an important impact in the social, commercial and financial activities of its inhabitants. Ash fallout in the other hand has not being that frequent, and in the rare instances they happened, were dealt on a rather empirical basis. In Costa Rica, two recent volcanic events: the 1973 Arenal Volcano eruption, consisting mostly of lava flows episodes, lasting until a few years ago, and the 1963-65 Irazu Volcano eruptions, basically of phreatomagmatic origin , joined with the recent influx of local geoscientists have lead to seek a more structured approach to deal with such possibility.

Irazu Volcano is located at the northern part of the Central Valley of Costa Rica, near the main population centers of the country. The increased vulnerability of telecommunication, medical equipment, air traffic, agriculture, etc. due to the ample use of miniaturized electronic devices turned it imperative to predict as accurate as possible where ash would fall should a new explosive event take place. This is no easy task considering that only two eruptions have been historically recorded: one in 1723, and another spanning from March 1963 till February 1965, a time when the country lacked neither prepared personnel nor scientific equipment to gather the required data to model it.

The 1963-65 events comprise of 3 major eruptions in April, July and December, the last one considered to be more explosive. In what follows, we will present an out of ordinary procedure we device to get a reasonable view of April 1963 episode, on account of the lack of information we faced. The main objective is to get as accurate as possible input parameters for TEPHRA2, to in turn improve the outputs resultant when different wind conditions are examined [1].

Historically, volcanic process and their impact have been dealt via a thorough description of the event, their main features, and lately by actual video records and satellite measurements. For the special case of the 1963-65 Irazu event, researchers have little to work with for isopach information was gathered 20 years after the event, after being weathered by tropical storms, [2], and in spite of having lasted 30 months newspaper coverage of the event is very scant for an editorial policy of the main newspaper, taken in June 1963, purposely suppressed all but the most important news related to Irazu activity. Therefore, indirect ways have to be found to find volcanological input parameters and realistically forecast future tephra deposition.

There exist several recent computer programs that have been used to calculate possible ash fallout patterns, namely HAZMAP, FALL3D, PUFF and TEPHRA2.

Being this our first attempt, we chose TEPHRA2 due to its model simplicity (based on Suzuki phenomenological approach), the small required set of initial parameters, its low demand on computer resources, and the availability of an open code source which allowed making the required adjustments to be run in the available computer equipment.

TEPHRA2 is a Linux based free software that can be downloaded from https://www.cas.usf.edu/~cconnor/vg@usf/tephra.html

Its Users Manual can be obtained from https://vhub.org/ resources/756/download/Tephra2_Users_Manual.pdf

Information about the use of the other mentioned programs can be found in Bonnadona [3] for HAZMAP, Folch [4] for FALL3D, and Kratzman [5] for PUFF and references there in.

TEPHRA2 has been the object of a thorough analysis in Bonadonna et al. [6] using a statistical view point. However, the applications of the ideas there proposed require knowledge and expertise not frequently available. We explored a simpler previous step that may help bridge the gap between two realities: one where qualified human resources in the geophysical and computational fields have the opportunity to productively interact in an already established ground; and that where such established ground is just in the making.

Using as a guide Scollo [7] observation that “results from our parametric study show that output of the models [mentioned above] can be strongly sensitive to the uncertainties and assumptions on input parameters, such as mainly mass eruption rate, column height, distribution of mass along the column, bulk grain size distribution”, supported by another TEPHRA2 user who repeatedly has used Suzuki Model “It is the experience of the writers that model results are strongly controlled by eruption volume, column height, and eruption velocity” [8] we explore the possibility of using first approximation initial values for TEPHRA2 to understand the short comings of the model, and of the initial values used.

Hence, we explored the possibility of using PLUMERIA as a tool to get a reliable value for the eruption velocity, and a plausible mass flux. PLUMERIA, software authored by Dr Mastin of USGS, allows the inclusion of key aspects as average weather conditions during the event, properties such as temperature and % of volatiles contained in the magma, and vent diameter. Although entrainment coefficients cannot be varied in the Windows version, and some of the assumptions may not fully apply to strombolian and vulcanian eruptions, it helped us include physical characteristics of the eruption process. Reader is warned that PLUMERIA does not include crosswind interaction. Windows based software PLUMERIA, can be downloadedfrom the site https://vulcan.wr.usgs.gov/Projects/Mastin/Publications/G3Plumeria/ framework.html

How the other required initial values for TEPHRA2 were acquired will be discussed below.

Suzuki model

The Suzuki model [8] is a convection-advection simple model that considers the eruptional column as composed by a gaseous phase that transports a Gaussian distributed tephra to a height determined by buoyancy. Horizontal uniform wind currents will then advent tephra to variable distances depending on its specific terminal fall velocity determined by the drag force suffered by each piece on account of its shape, texture, etc. It is based exclusively on mechanical considerations, so all the expected thermodynamics of the problem are lumped in two adjustable parameters: β and λ, that respectively control the average height of the particle Gaussian distribution maximum, and how sharp such maximum is. The complex mixing process caused by the interaction of the vertical volcanic column with the horizontal wind is not considered at all.

TEPHRA2 is based on this model, and requires the following input parameters: maximum φ, minimum φ, average φ, and it associated σ of the ash distribution; column height and exit velocity; mass rate of eruption, and of course, wind velocity and direction at the time of the event. In our case, the initial Suzuki model that included the β was substituted by Armienti formulation requiring the A coefficient. The maximum sharpness accounted by λ, taken since the beginning as λ=1 was later modified to take other positive integer values.

How Tephra2 initial values were calculated

Vent velocity calculation: To obtain an initial value as close to reality as possible, exit vent velocity was calculated with the help of PLUMERIA. This implied knowing vent diameter, magma temperature, and % of volatiles; the atmospheric conditions prevailing during the eruption, as well as geological parameters such as magma temperature and % of volatiles, as well as weather conditions at the time.

In our case, it allowed us to also include the fact that “the frequent coincidence of volcanic forcing with El Niño events disables the clear assignment of climate anomalies to either volcanic or El Niño forcing”, and that “the combined signal is different from a simple linear combination of the separate signals. It leads to a climate perturbation stronger than for forcing with El Niño or stratospheric aerosol alone and to a somewhat modified pattern”, as presented in Kirchner and Graf [9].

PLUMERIA offers the user two options to include the atmospheric conditions prevailing during the eruption:

1. A single layer atmosphere, which requires knowing the air temperature and relative humidity at the vent, the vent elevation as well as the lapse rate, tropopause elevation, tropopause thickness and the thermal lapse rate in the stratosphere. An educated guess can be easily done for the first three, not so for the last three, which requires the help of a knowledgeable expert for it varies with the geographical location of the eruption source as well as the month the vent took place.

2. A multilayer option, consisting of an independent text data file of relative humidity and temperature as functions of height. This alternative was used to increase the accuracy of the resulting initial vent velocity.

As satellite data at the time was basically non existence, we used NOAA reanalyzed diary record of atmospheric readings for year 1976. Original data was provided by Luis Alvarado, from his Master Thesis [10]. The great space-time scale of ENOS and the small size of the country lead us to consider it enough to use the profile data obtained from weather station at Juan Santamaría (10°N, 84°12´O, 920 a.s.l.) airport, with WMO code 78762. Besides, such a station is not far from Irazu Volcano.

To single out year 1976 as representative for 1963 weather conditions we paid attention to the 9 month medium intensity El Niño phenomena present from June 1963 till February 1964 and searched for a similar pattern in the following years. Although a very similar El Niño was also present during the June 1994 till March 1995 period, more than 30 years had elapsed since the eruptions, overshooting the official period for a meteorological cycle as proposed by the OMM. Instead, the 8 month medium intensity period of August 1976 till March 1977 was used as the closest one. That day of April 1976 that had the largest numbers of layers with acceptable readings was considered to be representative of the multilayer data for April 24, 1963 eruption. PLUMERIA also requires numerical input for the vent diameter magma temperature, and % of volatiles. We chose the following data:

Vent diameter: An estimated vent diameter of 500 meters was obtained from an in situ visual observation reported in Murata et al. [11]. As we are aware that this was the final diameter at June 1963, after the April eruptions, we tried a smaller value of 475 m.

Magma temperature and % of volatiles: a magma temperature of 1000°C [12], and 3% of volatiles [13] was input in PLUMERIA. A value of 0% water added was used for it was considered that March 1963 eruptions had evaporated the crateric lake present before the events took place.

One of the output results of PLUMERIA is a graph of column radius versus column height for a given vent velocity given by the user. An analysis of the corresponding graphs while the vent velocity is changed systematically shows that there exists a minimum vent velocity below which all the outputs show a collapsing column, reaching heights of no more than 500 m.

In situ observations revealed that during April 24, 1963 event tens of eruptions took place, all a few hundred meters in height, except two which reached heights of 6.8 Km and 8 Km. We propose that those two outliers were caused by mass ejected with the mentioned critical minimum vent velocity, which in our case was 89 m/s, a reasonable value [14].

A summary of values used in PLUMERIA:

Vent diameter (m): 475

Initial velocity (m/s): 89

Magma temperature (C): 000

Weight fraction gas: 0,030

Magma specific heat, (J/kg K): 1000

Magma density, (kg/m³): 2500

Mass fraction water added: 0,000

Mass flux, (kg/s): 6.00e+07

In summary, initial vent velocity as calculated with the help of PLUMERIA may be viewed as a way to include magma thermodynamics as well as atmospheric conditions in Suzuki model, giving a more realistic analysis. This approach also help update any hazard map if meteorological conditions change appreciable in time, as it may be expected with Climate Change predictions. This is done without changing the initial structure of the model, allowing comparison with other studies done using TEPHRA2.

In what follows, we discussed how the rest of the initial input required by TEPHRA2 was obtained.

Wind velocity: To be coherent, wind speed values were those associated to the April 1976 layers used in PLUMERIA for relative humidity and temperature mentioned above. Wind direction, however, was rotated some ten degrees to match that present at the time of the event, easing comparison with a hazard map published at the time.

Value of height: Although PLUMERIA also calculates a column height for the given conditions already discussed, the user is warned that such value does not include the effect of cross wind currents, required by Suzuki model. Rather, in situ observation of 8 Km was instead chosen.

Tephra gaussian distribution: Suzuki adopted a Gaussian tephra distribution to characterize the tephra transported by convection in the column, characterized by 4 terms: maximum φ, minimum φ, average φ, and its corresponding standard deviation σ. In Suzuki original paper as well as in those which followed, the corresponding values are derived from field observations. In Costa Rica, and possibly in other regions were eruptions are not frequent, it becomes a difficult task to find acceptable numbers for those 4 input values. Thus, we instead propose that as a first approximation we could use the curves presented in Woods and Bursi [15]. In our case we used the average values for Strombolian, Vulcanian, and Phretomagmatic eruptions Figure 1-3 [16]. Summary of values in Table 1.

Figure 1: Strombolianeruptions.

Figure 2: Vulcanian eruptions.

Figure 3: Phretomagmatic eruptions.

Table 1: Average Values for Strombolian, Vulcanian, Phreatomagmatic eruptions.

Suzuki β and Armienti a adjustable parameter: Suzuki original paper introduced β, an adjustable parameter designed to take into account the height of the maximum of the Gaussian distribution. Later the Armienti adjustable parameter A [17] was introduced, improving the model. Using the calculated relationship that the ratio of the height of the maximum concentration to the maximum column height equals (1-1/A), we found values of 4.0, 4.2, and 4.6 using 3 photos taken from afar, for eruptions at different times, quite near the A=4.2 quoted in Parra and Figueroa [18], and Pfeiffer et al. [19]. We settle for A=4, to ease comparison with other studies.

Suzuki λ adjustable parameter: Although Suzuki also introduced the adjustable parameter λ, to adjust how sharp the maximum is, a few lines after he argued in favor of λ=1. We propose that a more accurate model can be obtained if λ is allowed to vary as a positive integer [19]. A photo of Irazu eruption taken from afar shows a clear constant radius column for the first fourth of its height, more in line with a λ=3 case.

Mass eruption rate Q: Different relationships of mass eruption rate Q with height can be found in literature [14]. With the scant information at hand, we had no way to choose a given one. Besides, Clark et al. [2] observation that the largest historic eruption occurred in 1963–1965, and we estimate a minimum tephra volume of 3 × 10⁷ m³ for that eruption, based on an exponential linear regression of the log thickness isopachs, did not help much for the whole eruption lasted 30 months, an exceedingly long time, and included at least 4 major episodes lasting more than a day each. Using PLUMERIA values cited above, April 63 eruption is equivalent to 8.64 × 10⁷ m³, three times as much as the minimum estimated by Clark for the whole period. Reader is reminded that as Clark himself acknowledge, his raw data was acquired 26 years after the eruptions!

Thus, we walked a different path.

During Irazu eruption, a classified project known as Operation Ceniza-Arena [1] was set up to measure in situ rate of ash deposition in two plots of land selected from areas that had previously received moderate to heavy deposits of ceniza-arena from eruptions of Volcán Irazú. Plot 1 was near Ipis (9.4 miles at an azimuth 263º TN from the vent, an plot 2 near San Ramón (9, 4 miles at an azimuth 251º TN of the vent). The final report includes plenty of data of (gm/sq ft per hr) ash deposited during different periods, all of them measured before appreciable loss by weathering Table 2.

Table 2: Ash Deposited Rates.

Hence, we used the mass rate (Kg/s) calculated by PLUMERIA, into a mass (Kg/hr) and use it in TEPHRA2 as the eruption total mass. The idea was to simulate an eruption lasting exactly 1 hour, to be able to compare the actual data gathered in a given geographical location, with the resultant output in the correspondent point in the calculated map

For San José, the average deposition rates (gm/sq ft per hr) were:

An exponential fit reveals a decaying activity with a half-life of 0.74 month if maximum volcanic activity happened in mid-December. Although this value is in close agreement with the half-life of 0.87 month obtained from the two plots, in the June 1964 to February 1965, it cannot be extrapolated to April 1963. In addition, for Plot 1 sampled deposition rate amply varied during a given month: 0.8 for June 17 till 7.68 a day after (June 18); 0.0597 for July 20 till 2.78 for July 15, all in 1964.

In summary, the most we could do was to consider April 1963 simulated results subjectively reasonable well for they were roughly an order of magnitude lower than those measured in the same geographical points in 1964 Operation Ceniza-Arena [1].

Irazu April 24 1963 eruption

To test these ideas, three different scenarios were obtained by running TEPHRA2 with the three possible eruption modes: Strombolian, Vulcanian, and Phreatomagmatic, using the corresponding particle Gaussian distribution presented in the book The Physics of Explosive Volcanic Eruptions [16].

The PHREATOMAGMATIC result is the one that resemble the most the darkest part of the map published in La Nacion, the national newspaper, for April 1963 situation Figure 4.

Figure 4: Hazard Map.

The rapid increase in population density at the Central Valley of Costa Rica, joined to the corresponding health, government, telecommunication, recreational infrastructure needed has pushed local scientist to draw more accurate hazards maps that may consider different scenarios. To fulfill this requirement we should turn at least to available computer programs that run in powerful laptops. This forces to use the scant information available, and explore ways to combine two or more of such computer programs to calculate some of such initial conditions.

Using as a guide researches carried out in other latitudes, we here presented a way PLUMERIA was used to find a reasonable initial velocity, as well as to introduce thermo dynamical aspects in Suzuki model, a mechanically based model since the start.

Generally accepted particle Gaussian distributions were used as a first approximation to complete an acceptable set of initial values required to run TEPHRA2 that could then be used to generate hazard maps with different expected wind conditions.

This joint use of these two software helped finding acceptable initial volcanological conditions to be used in TEPHRA2, and eventually in other more complex programs.