Evidence of explosive hydromagmatic eruptions during the emplacement of the North Atlantic Igneous Province

Early Eocene sediments in northwest Denmark contain over 180 well-preserved volcanic ash layers, likely sourced from the North Atlantic Igneous Province (NAIP) between 56.0 and 54.6 Ma. Most of these ashes are basaltic, widespread, and represent a phase of unusually large and explosive eruptions that is coincident with the opening of the northeast Atlantic Ocean. Explosive basaltic eruptions of this magnitude are unheard of in historical times and in the current geological record. Here, we combine analyses of glass sulfur concentrations and variations in the morphology and vesicularity of pristine volcanic glass grains to explore the possible eruptive processes promoting such widespread basaltic ash dispersal. We suggest that these ashes formed in shallow subaqueous environments (<200 m water depth) where they fragmented and rapidly quenched during explosive hydromagmatic activity. We speculate that magma-water interaction during the opening of the northeast Atlantic Ocean was the main cause of this unusual explosive basaltic activity.


Introduction
The North Atlantic Igneous Province (NAIP) is a vast complex composed of massive continental flood basalts, widespread sill intrusions, and large volcanic centres surrounding the margins of the present-day northeast Atlantic Ocean and Labrador Sea [ Figure 1; Á Horni et al. 2017;Saunders et al. 1997;Storey, Duncan, and Tegner 2007]. The province began to form at~63 Ma, but the bulk of the igneous material (5-10ˆ10 6 km 3 ) was emplaced between 56 and 54 Ma, coincident with the breakup of Greenland from Eurasia and the formation of the northeast Atlantic Ocean [Larsen et al. 2016;Saunders 2016;Saunders et al. 1997;Storey, Duncan, and Tegner 2007;Wilkinson et al. 2017]. Most of the continental flood basalts in East Greenland, the sill intrusions in the Vøring and Møre basins (Figure 1), and the seaward dipping reflectors along the North Atlantic margins were emplaced during this time [Á Horni et al. 2017;Planke et al. 2005;Storey, Duncan, and Tegner 2007;Wilkinson et al. 2017]. The formation of the NAIP is of significant scientific interest because it coincides with both the breakup of the northeast Atlantic Ocean and the extreme climate perturbations of the late Paleocene and early Eocene, including the Paleocene-Eocene Thermal Maximum (PETM; [Saunders 2016;Schmitz et al. 2004;Storey, Duncan, and Swisher 2007;Storey, Duncan, Corresponding author: e.w.stokke@geo.uio.no and Tegner 2007; Svensen et al. 2004]). Explosive volcanism played an integral part of the NAIP activity, as evidenced by the preservation of hundreds of tephra layers in the North Sea region [Knox and Morton 1988;Larsen et al. 2003; Morton and Knox 1990]. Tephra is especially abundant in the earliest Eocene stratigraphy, suggesting increased explosive activity during the later stages of NAIP emplacement and the start of seafloor spreading [e.g. Larsen et al. 2003]. To our knowledge, this is the first tephra sequence that traces the activity of a Large Igneous Province (LIP) during continental breakup that has been recognised in the geological record. The unique preservation of the NAIP compared to other LIPs enables us to recognize processes that are normally not observable, improving our understanding of the explosive components of LIPs.
Particularly well-preserved tephra layers from the NAIP are exposed subaerially in northwest Denmark (Figures 1 and 2) due to recent glaciotectonic activity [Pedersen 2008]. These layers are all within the ash fraction with grain sizes <2 mm, and will henceforth be referred to as ash. More than 180 individual ash layers are preserved in the earliest Eocene stratigraphy, comprising organic rich marine clays of the Stolleklint Clay (syn-PETM) gradually overlain by the diatomite-rich Fur Formation (post-PETM). The ashes are stratigraphically sorted into a negative (´39 to´1) and positive (`1 to`140) ash series [Bøggild 1918]. The negative series is composed of a heterogeneous range of ash compositions (from basalt to rhyolite) that are often heavily G r e e n l a n d E u r a s i a N o r t h A m e r i c a Fur V ø r in g b a s in M ø r e b a s in Figure 1: Map of the known extent of the North Atlantic Igneous Province (NAIP) in a regional palaeogeographic reconstruction from 56 Ma. The striped area indicates the most likely source area for the Danish ash layers based on Larsen et al. [2003]. The yellow points show selected known localities of NAIP-derived ash (ash locations from Haaland et al. [2000], Jolley and Widdowson [2005], Morton and Keene [1984], and Morton and Knox [1990]), although NAIP-derived ash is also distributed far outside this map frame such as Goban Spur and in the Austrian Alps [e.g. Egger et al. 2005]. Dark red: volcanic centres. Red areas: the known extent of subaerial and submarine extrusive NAIP volcanism. Purple areas: the known extent of only submarine volcanism. Dark grey areas: known extent of NAIP sill intrusions. Blue lines: plate boundaries. Black lines: present-day coastlines. Light blue: shelf areas. Dark blue: ocean basins. Figure modified from Abdelmalak et al. [2016]; Á Horni et al. [2017]; Jones et al. [2019], and references therein.
altered. In contrast, the positive series is dominated by near-pristine tholeiitic basalts ]. All the ashes represent primary fall deposits with no evidence of post-depositional transport or redeposition [Pedersen et al. 1975;Pedersen and Surlyk 1983]. Normal grading in the ash layers indicates that each layer is the product of a single explosive eruption ( Figure 3).
It is generally accepted that the ashes are sourced from the NAIP, based on correlative chemistry and the lack of any other appropriate volcanic centres of the same age nearby ; Morton and Knox 1990]. However, the NAIP source volcanoes were located between 700 and 1500 km from the depositional environment of the eastern North Sea at that time [Figure 1;Abdelmalak et al. 2016]. These are exceptionally large transport distances for basaltic eruptions, during which the bulk of erupted material is more commonly deposited close to source [e.g. Parfitt 1998]. If the source area of the volcanoes is correct (Figure 1; ]), then each of these 180 ash layers represents an extremely large magnitude explosive eruption. For example, the Grímsvötn volcanic system in Iceland frequently produces explosive eruptions aided by hydromagmatic interactions with the overlying Vatnajökull ice cap. The early Holocene G10ka tephra series from Grímsvötn record several basaltic eruptions that are among the largest known Icelandic eruptions in the Neogene and covered an area of 2 million km 2 [Óladóttir et al. 2020]. However, only millimetre thick ash layers from the G10ka tephra series can be found at similar distances from the source as the Danish ash layers [Óladóttir et al. 2020]. By comparison, the largest basaltic ash in the Danish series is 16 cm thick [Bøggild 1918, 14 cm at Fur], and more than 20 of these tholeiitic basalts can be traced all the way to the Austrian Alps [Egger et al. 2005]. A key interval (Ash`19) has an estimated volume of 1200 km 3 (dense rock equivalent; DRE) and an areal extent of 16 million km 2 [Egger and Brückl 2006]. This study focuses on the early Eocene tholeiitic ash layers preserved in the diatomitic Fur Formation, and explores the eruptive processes that promoted such widespread ash dispersal. The large ash grain sizes (up to 500 µm), the thicknesses of the basaltic ash layers (2-160 mm), and the large distances from possible source volcanoes (700-1500 km) indicate that these were extremely explosive eruptions with high magma fluxes. Figure 4 shows the plume heights and wind speeds necessary to transport a 250 µm ash particle 750 km, based on the model of Stevenson et al. [2015]. This suggests that wind speeds need to be in the range of 35-60 m s -1 and plume heights in the range of 15-30 km, or possibly higher ( Figure 4). Therefore, it is only possible to transport basaltic ash grains of this size >750 km if the eruption was Plinian or similar, and if the stratospheric winds were in the right direction and at the upper end of present-day velocities. Note that the model assumes constant wind speeds at all altitudes, so given the lower average wind speeds in the troposphere, these values are likely to be minimum estimates. This strongly suggests that each of these ash layers represents a highly explosive (VEI 7-8) and large volume (M7-8; [Mason et al. 2004]) eruption. It is plausible that lateral stratospheric velocities were augmented by the high magnitudes and high intensities of the eruptions. Dynamical modelling of the formation of giant umbrella clouds (>600 km diameter) suggests that radial spreading in the umbrella region enhances lateral transport in the stratosphere and leads to a more circular shape than for less violent eruptions [Baines and Sparks 2005;Mastin et al. 2014]. At mid-to high-latitudes, these giant ash clouds are inherently unstable due to internal differences in Coriolis forces, meaning that at least part of the cloud is drawn towards the equator [Baines et al. 2008]. In the case of the NAIP-derived eruptions, the North Sea region to the southeast would have been ideally positioned to receive additional ash fall through this effect. This would reduce the wind velocities needed to transport ash particles such great distances in the early Eocene, but requires that the eruptions be both high and sustained in intensity.
A review of bore holes in the Pacific Ocean shows that the magnitude of eruptions from volcanoes in Japan correlate with both tephra layer thicknesses and distance from the source volcano [Mahony et al. 2016]. The only tephra layers thicker than 1 cm and over 1000 km from the source in their dataset were M8 eruptions [Mahony et al. 2016] defined as eruptions with >400 km 3 DRE magma volume [Mason et al. 2004]. Importantly, all known M8 eruptions have involved silicarich magmas. The high viscosities of these magmas inhibit the escape of volatiles, enabling the generation of overpressure and eventually highly energetic fragmentation [Cashman and Scheu 2015;Dingwell 1996;Papale 1999;Wallace and Edmonds 2011]. In contrast, the low viscosities of silica-poor basaltic melts do not sufficiently restrict bubble growth and instead enable gas-melt separation and non-explosive degassing of exsolved volatiles [Cashman and Scheu 2015;Mangan and Cashman 1996]. Sub-Plinian and Plinian eruptions of basaltic magmas are therefore rare, and those that are known have largely been attributed to rapid ascent and syn-eruptive crystallisation of hydrous basalts, and the associated increase in magma viscosity [Arzilli et al. 2019;Costantini et al. 2010;Houghton et al. 2004;Sable et al. 2006;Walker et al. 1984;Williams 1983]. Critically, these rheological constraints imply that the tholeiitic basaltic ash layers in the positive ash serieswhich we show below are largely crystal-free-likely reflect a type and scale of explosive basaltic eruption not observed in the historical record.
Without historical examples as a reference, we need to consider the possible mechanisms for producing widespread basaltic tephra. Fundamentally, a high magma flux is required. However, enhanced explosivity and fine fragmentation of basaltic magma would also require either (1) increased viscosity due to rapid syn-eruptive crystallisation, enabling brittle fragmentation mechanisms more similar to that of silicic magmas, or (2) magma-water interactions amplifying brittle fragmentation and generating the heat transfer necessary to drive high plumes for widespread ash dispersal. The second option is conceivable, as the . The model uses conservative estimates of 250 µm basaltic ash particle size and a 750 km transport distance (shown as a black line). The input parameters use a basaltic glass density of 2800 kg m -3 , a grain sphericity of 0.7, a size-dependent density function [Bonadonna and Phillips 2003], and a Ganser fall velocity model. The model makes the unrealistic assumption of constant wind velocities at all altitudes and in a straight line from source to deposition, meaning that the output is a minimum estimate of the required plume heights/wind speeds to distribute ash particles these distances.
rupturing of Greenland from Eurasia led to seafloor spreading through the middle of the Greenland-Faroe continental flood basalt provinces in the early Eocene Storey, Duncan, and Swisher 2007;Storey, Duncan, and Tegner 2007]. This rifting would have allowed a seawater incursion at the heart of the elevated magmatic activity, which could have produced the explosive basaltic ashes preserved in Denmark and the North Sea. A hydromagmatic origin for the ashes has been proposed previously by several authors, based on textural observations and the widespread ash distribution [Haaland et al. 2000; Morton and Evans 1988;Pedersen and Jørgensen 1981], but has yet to be rigorously tested. Explosive hydromagmatic eruptions occur when erupting magma interacts with a body of external water, resulting in rapid quenching of the melt and extensive melt fragmentation [Büttner et al. 1999;Morrisey et al. 2000;Németh and Kósik 2020;Thórdarson et al. 1996]. Rapid quenching arrests the degassing process, leading to elevated residual volatile contents compared to magmatic tephra degassed to atmospheric pressure [Mastin et al. 2004;Németh and Kósik 2020;Óladóttir et al. 2007]. Therefore, one method that can potentially identify hydromagmatic eruptions is measuring the volatile contents of quenched glass, particularly late-stage degassing species such as sulfur and chlorine [Liu et al. 2018;Mastin et al. 2004;Óladóttir et al. 2007;Thórdarson et al. 2003]. Recent studies of basaltic eruptions in Iceland have used this approach to identify changes in glacial thickness above Katla volcano [Óladóttir et al. 2007] and to differentiate between coeval magmatic and hydromagmatic eruptions along a single fissure [Liu et al. 2017;Liu et al. 2018]. Coupled with volcanic glass morphology, Liu et al. [2018] demonstrated that glass sulfur concentrations were elevated in the hydromagmatic products due to quenching at pressures greater than atmospheric (1 atm). Although the Fur Formation ashes are considerably older (56.0-54.6 Ma) than these Holocene tephras, they are remarkably well-preserved. We therefore measure glass volatile contents (S, Cl), in combination with component analysis of ash particle morphology, throughout the ash series to test the hypothesis that the NAIP tephras were a product of hydromagmatic activity. We use these data to infer relative changes in the Volcanica 3(2): 227 -250. doi: 1 .3 9 9/vol. 3. 2.22725 pressure they were erupted under (i.e. depth below sea level) over time. Our results advance our understanding of the evolution of the NAIP and the physical volcanology of hydromagmatic eruptions more generally; an eruptive style that is currently under-represented in the literature, particularly for eruptions in deep time.

The island of Fur
Ash samples used for this study were collected on the island of Fur, northwest Denmark (Figures 1  and 2). During the Paleocene-Eocene transition and early Eocene, this area was part of a marginal basin within the eastern part of the larger semi-enclosed epicontinental North Sea Basin. The stratigraphy exposed at Fur reflects both small and large scale climatic and tectonic changes. A marine regression due to thermal uplift of the NAIP and the onset of PETM global warming marks the Paleocene-Eocene transition in Denmark. This caused a lithological shift from the latest Paleocene condensed, bioturbated clays of the Holmehus Formation, into the earliest Eocene dark, laminated Stolleklint Clay [ Figure 5; Heilmann-Clausen 1995;Heilmann-Clausen et al. 1985]. The PETM onset has been identified at the base of the Stolleklint Clay by a 4.5 % negative shift in organic δ 13 C values [Jones et al. 2019]. Towards the end of the PETM the Stolleklint Clay grades upward into the 60 m thick Fur Formation ( Figure 5), composed of diatomite with variable clay contents. The sediments are thermally immature and have experienced very little lithification and consolidation ]. However, Quaternary glaciotectonic activity has generated large scale deformation in northern Denmark, creating abundant small scale folding and thrusting throughout the Fur stratigraphy [Pedersen 2008]. Despite this, ash layers are mostly laterally continuous and easy to trace, without any obvious overthickening.

The Danish ash series
More than 180 ash layers are interbedded in the stratigraphy, with the majority (~140) found within the diatomitic Fur Formation ( Figure 3A; Figure 5). The volcanic ashes are grouped into a negative and positive ash series [Bøggild 1918], with additional ash layers [termed SK1, SK2, and SK3; Jones et al. 2019] within the base of the Stolleklint Clay ( Figure 5). Numerous efforts have been made to constrain the timing of explosive activity preserved in Denmark. The lowermost ash layers (SK1 and SK2) are found stratigraphically just below the base of the PETM carbon isotope excursion (CIE), which is estimated to have begun between 56.0 and 55.9 Ma [ Figure 5; Westerhold et al. 2017 Figure 2). Pa. = Paleocene; Th. = Thanetian; H./Ø. = Holmehus/Østerrende Formation. For additional lithological information, see Stokke et al. [2020]. Extent of PETM carbon isotope excursion (CIE) taken from Stokke et al. [2020]. Grey lines indicate ash layers, with main ash layers labelled and maximum ash layer thickness measured at Fur indicated on the top. Note that there are some thickness variations within Denmark [Bøggild 1918]. 1 Charles et al. [2011], assuming the timings of the Svalbard and Fur CIEs are coeval; 2 Storey, Duncan, and Tegner [2007]; 3 Westerhold et al. [2009]; 4 King [2016]. recovery of the PETM CIE has a U-Pb radioisotopic age of 55.785˘0.086 Ma [Charles et al. 2011]. If the Svalbard and Danish CIEs are coeval, this would be located stratigraphically about 1.5 m above Ash´33 in the Fur Formation [ Figure 5; Stokke et al. 2020]. Approximately 10 m up the section, 40 Ar/ 39 Ar dating of Ash layer´17 gives a corrected absolute age of 55.60˘0.12 Ma [ Figure 5; Jones et al. 2019;Storey, Duncan, and Swisher 2007]. Westerhold et al. [2009] estimated a 200 kyr duration between Ash´17 and`19, based on correlation to the orbital cyclostratigraphy of the DSDP 550 core at Goban Spur. The uppermost basaltic Ash`140 marks the top of the~60 m thick Fur Formation [Bøggild 1918;Heilmann-Clausen et al. 1985]. There is no evidence of the magnetic reversal from C24r to C24n within the Fur Formation nor the second Eocene thermal maximum (ETM2;~54.1 Ma; [King 2016;Stokke et al. 2020]). The positive series basalts correlates with subphase 2b of pyroclastic deposition in the lower Balder Formation [Knox and Morton 1988], which suggest that the top of the Fur Formation is no younger than 54.6 Ma ( Figure 5; Heilmann-Clausen [pers. comm.]; [King 2016;Watson et al. 2017]). Therefore, the evidence of explosive volcanism in Denmark likely spans~1.4 Myr from 56.0 to 54.6 Ma, with deposition of the basaltic positive series initiating sometime after 55.6 Ma ( Figure 5).
The provenance of the Danish ash layers was a source of some contention. Earlier work assumed that the ashes were derived from a nearby volcanic source in the Skagerrak or North Sea [Nielsen and Heilmann-Clausen 1988;Pedersen et al. 1975]. However, subsequent geophysical and seismic surveys failed to find any volcanic centres in these areas, suggesting that the ashes originated from further afield. A strong geochemical correlation between North Sea ashes and NAIP volcanic centres ], coupled with widespread ash sequences across the North Sea and northern Europe [King 2016;Morton and Evans 1988;Morton and Knox 1990], strongly support the hypothesis that these ash layers originated from NAIP sources. A comprehensive attempt to correlate the ashes with different NAIP volcanic centres divided the ashes into four stages of volcanic activity ]. The first three stages involve ashes SK1-3 and the negative ash series. These include a heterogeneous mix of variably-preserved ashes, likely sourced from different volcanic centres on the British Isles, East Greenland, and around the northeast Atlantic margin ]. Stage 4 comprises the entire positive ash series, which temporally and chemically correlates to the offshore Balder Formation ashes and likely originated from the main rift zone during opening of the North Atlantic [ Figure 1; Larsen et al. 2003; Morton and Knox 1990]. The positive ash series is composed of well-preserved, glassy tholeiitic basalts, with the exception of the rhyolitic layers Ash`13 and Ash 19 Pedersen and Jørgensen 1981].

Sample description
A total of 21 ash layers are included in this study, all from the post-PETM Fur Formation and collected from two different outcrop exposures on Fur Island (Figure 2). Ash´19 was sampled at the Stolleklint beach (56°50'29"N, 8°59'33"E), while Ash´13 and all of the positive series ashes from the Knudeklint beach locality (56°50'14"N 8°57'26"E). Additional ashes from the negative series were also sampled, but these turned out to be too altered to yield useful analyses and are therefore excluded from this study. The main part of the samples are therefore from the positive series, covering the interval from Ash`1 to Ash`118. The ash series does continue up to Ash`140, but these are thinner and difficult to sample due to poor exposure and preservation. The ash layers are relatively distinct and easy to trace laterally. They vary in thickness from <1 to 20 cm (Figure 5) with grain sizes up to 500 µm. The ash layers typically have a sharp lower boundary followed by normally-graded bedding consistent with primary deposition ( Figure 3); importantly, these properties suggest that no significant post-depositional overthickening has occurred. Except for bioturbation at the top of some of the ash beds and occasional fluid escape structures, there are no sedimentary structures within the ash, indicating there has been no later redeposition.
The ash layers are variably consolidated and cemented. Hard, cemented layers were sampled as solid blocks using a hammer and prepared as 30 µm thin sections at the University of Oslo. Samples that were consolidated but not cemented were carved into suitable pieces with a knife, then dried at 80°C and glued with epoxy, before being prepared as 30 µm thin sections. Unconsolidated samples were oven dried at 80°C, then dry-sieved into six individual size fractions: >2 mm, 1-2 mm, 500 µm-1 mm, 250-500 µm, 125-250 µm, and <125 µm. The four smallest fractions were arranged in quadrants and mounted in epoxy grain mounts, which were then polished to intersect the ash grains.

Electron Probe Microanalysis (EPMA)
Polished and carbon coated grain mounts were analysed on a Cameca SX100 microprobe equipped with 5 wavelength dispersive spectrometers (WDS) at the University of Cambridge, UK. Analyses were conducted with 15 kV accelerating voltage and a 10 nA beam current using a defocused beam size of 10 µm. Counting times were 5 s for K, Na, and Si; 10 s for Ca and Fe; 15 s for P, Al, and Cl; 20 s for Ti, Mn, and Mg; and 60 s for S. A combination of mineral and glass standards were used for primary calibration. Repeat analyses of secondary standards were used to monitor instrumental drift during and between analytical sessions. Analyses were conducted in two sessions, each For each ash layer, we analysed the matrix glass compositions of between 20 and 24 separate particles, distributed between the four size fractions. Most analyses were from ash particles within the 125-250 µm fraction; the 10-µm beam size made it difficult to ensure a sufficient area of clear glass in the smallest fraction and the larger size fractions were mostly composed of partly consolidated ash aggregates and smectite clay.

Image acquisition and textural analysis
Polished and carbon-coated grain mounts were imaged using a Quanta-650F Field Emission Gun Scanning Electron Microscope (FEG-SEM) at the University of Cambridge, UK, operating in backscattered electron (BSE) mode. We acquired images of individual ash particles at high vacuum and a 15kV accelerating voltage. Composite images of the 125-250 µm quadrant of the grain mounts (8ˆ8 tiles) were acquired using MAPS 2.0 software (FEI, Thermo-Fischer), with a three-point focus interpolation and a zoom of~350ˆ.
We analysed the morphology and texture of the ash particles using backscattered SEM image maps, by manually counting and classifying 400 particles within a defined area of each sample image. The degree of post-depositional alteration of the Fur ashes prevented reliable automated shape analysis [Dellino and Volpe 1996;Liu et al. 2015]. During manual classification, the ash particles were first categorised into four groups: altered particles, microcrystalline grains, lithic particles, and silicate glass ('glassy') grains ( Figure 6). Altered particles include grains were the primary textures are no longer recognisable due to alteration and devitrification. Microcrystalline grains contain plagioclase microlites within a silicate glass matrix. Lithic particles are holocrystalline, and likely sourced from syneruptive erosion of conduit walls and surrounding rock material. Silicate glass particles are crystal-free, and the main focus of this study. These were further classified into three different categories based on their vesicle texture and external morphology: dense, vesicular, and shards ( Figure 6). Dense, glassy grains are characterised by a blocky morphology with angular fracture surfaces and poor internal vesicularity (ď20 area %; hereafter we simply use % for area %). Shards are characterised by concave edges that comprise >50 % of the perimeter of the grain (interpreted as portions of bubble walls) with few internal vesicles. Vesicular grains are characterised by >20 % internal vesicularity. Vesicles are generally circular (in cross-section) with relatively thick bubble walls.
The two ash layers from the negative series, Asheś 19 and´13, have slightly different glass chemistries than the positive series. Ash´19 has lower SiO 2 and higher TiO 2 than the positive series (Figure 8), while Ash´13 is slightly outside the main trend with a marginally more alkaline composition (Table S1, Supplementary Material) and a large variance in glass oxide concentrations ( Figure 8). However, Ash´13 is also the most affected by alteration of the analysed ashes, as indicated by the overall physical appearance ( Figure 6A).
The positive ash series shows a stronger trend for both CaO vs MgO and Al 2 O 3 vs MgO (Figure 8). The correlation is particularly evident between MgO and CaO, with the linear regression characterised by an r 2 of 0.90. This relationship suggests that the stratigraphically highest layers (Ash`94,`114, and`118) are the most evolved, with low MgO, CaO, and Al 2 O 3 , and high K 2 O. However, the linear trend is not a simple function of stratigraphic depth, as the higher MgO concentrations of Ashes`54 and`102 suggest they are more primitive than the stratigraphically lower layers 18 and`46 ( Figure 8). The linear trend in matrix glass chemistry involving plagioclase-compatible elements suggests that plagioclase fractionation is at least partly controlling the range of glass (melt) compositions. In contrast, the correlation is very weak between S and K 2 O (r 2 <0.05), even for the more chemically homogenous positive ash series (Figure 8). This suggests sulfur concentrations are controlled by a process other than crystal fractionation.

Mineralogy and texture
Although the partly consolidated nature of the layers inhibited a detailed grain size analysis, we find that the distribution is skewed towards finer grain sizes. Most of the grains within size fractions >500 µm are composed of aggregates of smaller glass grains and fragments, while the <125 µm fraction was difficult to analyse due to the small size, and abundance of clay and devitrified grains. The 125-250 µm fraction was therefore selected for morphological analysis, as this provided most information. Ash samples are predominantly composed of glass grains (altered or pristine; Figure 9). Altered and devitrified grains are common throughout; up to 75% of grains are altered in some layers, including Ashes`18 and`31 ( Figure 9). The lowermost Ash´19 is the least altered sample, while the positive series shows an upward decrease in altered and devitrified glass grains (Figure 9). Microcrystalline and lithic particles are generally sparse and make up a negligible portion of the total particles, with only minor portions (<5 % combined) of lithics, microcrystalline glass, and crystals ( Figure 9). The mineral assemblage is composed mostly of plagioclase and Fe-Ti oxides (mainly ilmenite), while pyroxenes and olivines are only observed occasionally. Figure 10 shows how the texture and morphology of pristine glass grains (i.e. excluding altered grains) varies stratigraphically. The dense glassy fraction varies between 17 % and 52 %, with a minor overall increase from Ash`18 upwards. Ash`1 deviates from this stratigraphic trend with a high fraction of dense glassy grains of about 50 %. The vesicular fraction varies between 12 % and 54 %. Although less distinct, it follows the dense fraction with an overall decrease upsection. The fraction of shards is more stable and shows no obvious stratigraphic variation, fluctuating between 18 % and 39 % (Figure 10).

Dissolved volatile contents
Average sulfur concentration for all analysed grains is 592˘381 (2σ) ppm S, while all individual analyses vary between 171 and 1131 ppm S. Each individual analysed glass grain was also classified as dense, vesicular or shard, in order to investigate the relationship between morphology and residual volatiles in individual grains ( Figure 11). Although the dense and vesicular  Baragar (1971). Both plots include all analyses, data in Table S1, Supplementary Material. fractions mostly overlap, there seems to be a slight tendency for higher sulfur concentrations within the dense fraction. The fraction of shards is under-represented as their narrow dimensions made them difficult to analyse.
Sulfur concentrations are plotted stratigraphically together with the textural analysis in Figure 10, in order to evaluate the changes in residual sulfur over time. While the distribution of analyses within some samples are skewed towards higher or lower sulfur concentrations, most samples show a relatively normal distribution ( Figure 10). Average sulfur concentrations (where averages are calculated from all analyses within an individual ash layer) show abundant stratigraphic variation between a minimum of 431˘191 (2σ) ppm (Ash`90) and a maximum of 916˘120 (2σ) ppm (Ash`118). There appears to be a slight overall stratigraphic increase in the mean sulfur concentrations from 476˘327 (2σ) ppm S at the base (Ash´19) to 922˘239 (2σ) ppm S at the top (Ash`118; Figure 10), although given the limitations of this dataset this observation is not statistically significant (r 2 median = 0.12; Figure 10).
The average chlorine concentration for all analyses is 139˘88 (2σ) ppm Cl, with individual analyses varying between 43-371 ppm Cl. Mean chlorine concentrations for each ash sample vary between 88˘50 (2σ) ppm (Ash`102) and 200˘77 (2σ) ppm (Ash´19). Again, there is a strong division between the negative and positive ashes, with the negative Ashes´19 and´13 having comparatively higher chlorine concentrations than the positive series ( Figure 12). There also seems to be a stronger linear correlation within the positive series, than when the negative series is included (Figure 12).
The opposite is true when chlorine is plotted as a function of K 2 O (Figure 12), suggesting that the chlorine contents in the positive and negative ashes were initially different or controlled by different processes.

Morphological and textural evidence
Textural evidence shows that the proportions of dense glassy grains are generally high for a basaltic eruption throughout the stratigraphy. Most samples contain >30 % dense glassy grains, and in the uppermost 10 m of the stratigraphy (Ashes`46 to`118) there is consistently >40 %. Vesicles are macroscopic evidence of degassing, and typical vesicularities for basaltic tephra erupted in Hawaiian fountains generally range between 45 and 95 % [Parcheta et al. 2013;Porritt et al. 2012;Stovall et al. 2012]. Therefore, the high abundance of dense blocky glass morphologies with poor vesicularity suggests that degassing was arrested prematurely [Cashman and Scheu 2015;Graettinger et al. 2013;Houghton and Gonnermann 2008]. The blocky morphologies with brittle, angular surfaces and poor vesicularity are consistent with brittle fragmentation due to thermal fractures, rather than inertial fragmentation due to bubble expansion [Houghton and Gonnermann 2008;Mangan and Cashman 1996;Morrisey et al. 2000]. Pyroclastic textures and morphologies are often used to determine the fragmentation process during an eruption, primarily differentiating between "wet" hydromagmatic and "dry" magmatic eruptions, although     Table S1, Supplementary Material. use of this morphological approach alone can be ambiguous [Mastin et al. 2004;White and Valentine 2016]. Volcanic particles resulting from brittle fragmentation due to magma-water interaction are generally characterised by dense and blocky fracture-bounded glass grains, with low but variable vesicularity [Büttner et al. 2002;Büttner et al. 1999;Dellino and Volpe 1996;Graettinger et al. 2013;Liu et al. 2017;Liu et al. 2015;Mastin et al. 2004;Morrisey et al. 2000]. Hydromagmatic explosions are associated with highly efficient fragmentation, and therefore the deposits are typically more fine-grained compared to basaltic "dry" magmatic deposits [Zimanowski et al. 2003]. Liu et al. [2015] found that the amount of dense glassy grains increased with decreasing grain size, and decreased with increasing distance to the source. The high abundance of dense glassy grains in the Fur ashes despite the long transport distance, coupled with the dominance of fine ash-sized material, strongly supports a hydromagmatic origin for the ashes.
Highly explosive eruptions are generally associated with rapid magma ascent, such as during Plinian to ultra-Plinian eruptions [Cashman and Scheu 2015], but Plinian-type eruptions are rare for low viscosity basaltic magmas [Houghton and Gonnermann 2008].
The few examples in the geological record have been attributed to extensive syn-eruptive microlite crystallisation, such as the 122 BCE Mount Etna eruption [e.g. Sable et al. 2006] or the 1886 Tarawera eruption [e.g. Houghton et al. 2004]. The associated increase in melt viscosity will inhibit bubble expansion and promote rapid ascent and explosive fragmentation [Arzilli et al. 2019;Costantini et al. 2010;Szramek 2016]. However, the Fur ashes show relatively low crystallinities, with small fractions of both phenocrysts and microlites (Figure 9). This suggests that syn-eruptive crystallisation processes cannot be primarily responsible for the enhanced explosivity required to generate the Danish  Figure 6) of pristine glass, altered grains, and microcrystalline and lithics in each sample. Data in Table S4, Supplementary Material.
ash layers. While there is a high proportion (up to 50 %) of dense (probable) hydromagmatic glass grains, vesicular grains and shards are prevalent up-section indicating at least partial exsolution of a vapour phase (Figure 10). Importantly, however, the presence of vesicular glass grains does not exclude a hydromagmatic origin. Instead, it shows that magma-water interaction oc-curred after the onset of decompression-driven exsolution. Several studies have demonstrated that solely using pyroclastic texture and morphology to differentiate between hydromagmatic and magmatic fragmentation, largely on the basis of vesicularity, may be too simplistic [Liu et al. 2017;White and Valentine 2016]. If the magma is quenched during rather than prior to vesiculation, then we would expect a spectrum of vesicularities and thus morphologies. In addition, the entire mass of erupted material will not be in direct contact with water during a hydromagmatic eruption, leading to variations in the resulting fragmented material [Houghton et al. 1996;Schipper and White 2016;White 1996]. The inner part of an eruption column can also be thermally insulated [e.g. Cole et al. 2001], leading to continued vesiculation and degassing during eruption. Overall, it is very likely, and even expected, that hydromagmatic erupted material will exhibit considerable variation in texture and morphology, even within a single eruption.

Evidence from volatile analyses
Additional diagnostic properties can be used to further elucidate the eruptive mechanisms and degassing histories of eruptions. Volatile analysis is an efficient tool to uncover the degassing history [Dixon et al. 1991;Edmonds 2008;Thórdarson et al. 2003;Thórdarson et al. 1996;Wallace and Edmonds 2011]. Elevated volatile concentrations in glass-above that expected in equilibrium with atmospheric pressures-indicate quenching before the melt is fully degassed. This would be most likely attributable to quenching under elevated pressures, such as during subaqueous or subglacial magma-water interactions [Dixon et al. 2002;Liu et al. 2018;Mastin et al. 2004;Métrich et al. 1991;Óladóttir et al. 2007]. Although most quantitative assessments of degassing histories have focussed on CO 2 and H 2 O exsolution [e.g. Newman and Lowenstern 2002], sulfur in glass is more sensitive to small variations in quench pressure at shallow depths due to its greater solubility compared to CO 2 and H 2 O in tholeiitic basaltic melts. However, there are many possible factors affecting the residual sulfur concentrations of volcanic glass [Mavrogenes and O'Neill 1999;Wallace and Carmichael 1992;Wallace and Edmonds 2011], and it does not necessarily reflect only variations in hydromagmatic activity. Here, we consider the various processes that could produce the observed trends in sulfur content.
(a) Fractionation: In the absence of degassing or formation of a sulfur-rich phase (e.g. sulfide minerals), sulfur should behave like an incompatible element and increase with increasing magmatic fractionation [Wallace and Edmonds 2011]. A positive correlation between sulfur and K 2 O would therefore indicate that the sulfur concentrations preserved in matrix glasses are controlled by fractionation rather than degassing. For the Fur ash samples there is a poor correlation between S and  (b) Post-fragmentation degassing: Residual sulfur concentrations could also have been affected by post fragmentation degassing, depending on the time between fragmentation and quenching. However, this depends strongly on the cooling rate of individual glass grains [Lloyd et al. 2012]. Cooling duration is dependent on clast size and distance from the exterior surface [e.g. Mastin 2007]. Fast cooling rates of 10 3.9 to 10 5.1 K s -1 is determined for 0.5 mm angular hyaloclastite particles [Helo et al. 2013], similar to that calculated for cooling of outer surfaces under water spray impact [Mastin 2007]. Glass grains analysed from the Fur ashes are all <500 µm, and most within the 125-250 µm fraction, suggesting they quenched rapidly and that significant post-fragmentation degassing of sulfur is unlikely.
Assuming that the range of residual sulfur concentrations in the Fur ash layers is indeed indicative of differences in quench pressure, we now consider what threshold separates hydromagmatic from magmatic erupted glass. A study of Hawaiian basalts found that the average dissolved sulfur of partially degassed hydromagmatic tephra was about 330 ppm S, compared to 100-150 ppm S expected for fully degassed Hawaiian tephra [Mastin et al. 2004]. However, Hawaiian basalts are generally lower in iron than typical Icelandic basalts, resulting in lower sulfur solubilities and subsequently lower initial sulfur concentrations [Liu et al. 2018;Wallace and Carmichael 1992]. The sulfur content expected for fully degassed Icelandic melts is therefore slightly higher, varying between 300-500 ppm depending on iron content [Thórdarson et al. 2001]. In order to interpret the Fur ash layers, some assumptions have to be made regarding their similarity to modern volcanic systems. Several previous studies have remarked on the similarity between the Fur ashes and Palaeogene and present-day Icelandic basalts Morton and Knox 1990;Pedersen et al. 1975]. Morton and Knox [1990] noted how the Danish ashes resemble Icelandic tholeiites in the Th-Hf-Ta diagram of Wood [1980], both plotting within the E-type MORB and tholeiitic withinplate basalt field. Several studies have presented a summary of Icelandic glass volatile contents in the form of a S vs TiO 2 /FeO scatter plot (Figure 13). This clearly shows (a) the progressive decrease in glass sulfur content from melt inclusions, to hydromagmatic glasses, to magmatic glasses, and finally to lava selvages, and (b) the effect of Fe on the initial sulfur content [Liu et al. 2018;Óladóttir et al. 2007;Thórdarson et al. 2001;Thórdarson et al. 2003]. The Fur ashes plot well within the typical TiO 2 /FeO values for modern Icelandic basalts (Figure 13), with intermediate TiO 2 concentrations close to those of Laki [Thórdarson et al. 1996], Grímsvötn [Sigmarsson et al. 2013], and Surtsey [Schipper et al. 2015]. The median sulfur concentrations of the Fur ashes vary between~445-753 ppm S, which overlap with both the hydromagmatic tephra of Surtsey and Grímsvötn, and the magmatic eruptions of Laki ( Figure 13).
Chlorine concentrations are less diagnostic of degassing processes, as Cl is soluble until low pressures in basaltic magma [e.g. Sigvaldason and Óskarsson 1976]. Our analyses show that chlorine concentrations are quite variable ( Figure 12); which may reflect some shallow post-fragmentation degassing. Chlorine is also an incompatible element, and has been shown to correlate strongly with K 2 O [Davis et al. 2003;Sigvaldason and Óskarsson 1976], which we also see in our data (Figure 12). Previous studies on Icelandic tephras have shown that the Cl contents of magmatic and hydromagmatic samples overlap, with the magmatic glasses extending to slightly lower Cl concentrations indicating partial degassing [Liu et al. 2018]. Liu et al. [2018] found that Icelandic hydromagmatic glasses generally had chlorine content above 140 ppm Cl, while a study of Hawaiian basalts found that subaerial basalts varied between 64 and 119 ppm Cl and submarine from 70 to 560 ppm Cl [Davis et al. 2003]. The positive series ashes have median values from about 80-170 ppm Cl, and correlate moderately well with sulfur contents particularly at higher sulfur concentrations (Figure 12). This may reflect some partial degassing potentially due to hydromagmatic activity. Ashes´19 and 13 have higher chlorine contents than any of the positive ashes, yet below average sulfur contents indicate that this is not due to arrested degassing. We show in Figure 7 that these two ashes have a different chemistry to the positive ash series. Indeed, they are older and likely erupted prior to the main rift phase, and therefore likely have different equilibrium Cl concentrations, though the high chlorine contents could also be due to some crustal contamination [e.g. Davis et al. 2003].  Surts ey (Schipper et al., 2015) Katla (Óladóttir et al., 2007) Fur (this study) gether (Figure 10), there seems to be a clear positive co-variation between the fraction of dense glassy grains and the median residual sulfur concentrations in glass. This suggests that residual sulfur concentrations are a function of the extent of degassing, where the least degassed melts are also the least vesicular, and vice versa. This supports the conclusions drawn from the morphological evidence that the Fur ash sequence quenched under varying pressures higher than 1 atm.

Relative degree of sulfur degassing
The residual sulfur in quenched silicate glasses can be used to indicate the total amount of sulfur degassing, if the initial sulfur concentration of the melt is known [Devine et al. 1984]. Melt inclusions in mineral phenocrysts are an ideal source to measure the initial sulfur concentration [Davis et al. 2017;Métrich and Wallace 2008;Thórdarson et al. 2003]. Unfortunately, no viable melt inclusions have been identified in the Fur ashes. Instead, we estimate the un-degassed sulfur content of the melt using theoretical models of the sulfur concentration at sulfide saturation [SCSS; Blake et al. 2010;Smythe et al. 2017;Wallace and Carmichael 1992]. A major assumption when applying SCSS models is that the magma was sulfide-saturated. Both MORB and Icelandic melts are thought to be sulfide-saturated for much of their ascent path [Edmonds and Wallace 2017].
Considering the geochemical similarity between Icelandic melts and the Fur ashes (Figure 13), we assume that the Fur ashes were similarly sulfide-saturated.
There are several available models for calculating the SCSS, most based on a linear function of sulfur and iron concentrations of the silicate melt [Blake et al. 2010;Wallace and Carmichael 1992]. This is largely due to the important role of iron in enhancing sulfur solubility by combining with S 2to form an immiscible Fe-S-O sulfide liquid phase [Mavrogenes and O'Neill 1999;Wykes et al. 2015]. However, recent studies demonstrate that at fixed temperature, pressure, and silicate melt compositions, the SCSS also depends largely on the composition of the immiscible FeS-NiS-CuS 0.5 sulfide liquid itself [Patten et al. 2013;Smythe et al. 2017]. Consequently, the relative input of Fe-Ni-Cu needs to be taken into account, and we therefore apply the SCSS model of Smythe et al. [2017]. In the absence of a known sulfide melt composition for the Fur ashes, we apply the average sulfide globule chemistry from the Northern Mid-Atlantic Ridge [Keith et al. 2017] as the closest available analogue for the NAIP. At temperatures >1100°C and at pressures about 0.1 GPa, the sulfide melt is typically still a liquid phase "precrystallisation" [Patten et al. 2013]. We subsequently apply a pressure of 0.1 GPa, and a temperature of 1130°C calculated with the glass/liquid thermometer (Equation 1) of Montierth et al. [1995]: T p˝Cq " 23.0MgO liq`1 012˝C (1) where MgO liq is the average MgO content for Fur glasses of 5.14 wt % (Tables S1 and S3, Supplementary Material). Figure 14 shows the modelled SCSS for each Fur ash sample plotted according to its average iron concentration. Using these values for the SCSS, which represent the theoretical un-degassed sulfur content of the melt, the total amount of sulfur degassed can be estimated by calculating the difference between the initial and the residual glass sulfur concentrations for each Fur ash layer ( Figure 10). The ashes record total amounts of degassing between 55-80 % of the initial sulfur content, with the least sulfur degassing in the uppermost Ash`118 ( Figure 10). Changes in pressure have a limited effect on the modelled SCSS, and within the tested temperature range of 1000-1500°C the total degassing vary with no more than about˘20 % (Table S3, Supplementary Material). Therefore, even if the NAIP source had anomalously high magma temperatures, the total degassing is unlikely to be more than 85 % of the initial sulfur content for any given sample (Table S3, Supplementary Material).
Comparing to other basaltic systems, Mastin et al. [2004] documented a similar 70-80 % degassing(relative tothe initial sulfur content) in hydromagmatic tephra from the Kīlauea Volcano (Hawai'i), attributed to partial degassing due to quenching at higher than atmospheric pressures. Partial degassing was also observed in several Icelandic hydromagmatic systems, such as a 17-70 % degassing of the preeruptive sulfur content from the Hverfjall Fires [Liu et al. 2018] and the 18-75 % degassing from the subglacial 1783 Grímsvötn eruption [Métrich et al. 1991]. Highly variable volatile loss appears to be a general feature of hydromagmatic glasses, and may reflect variable degrees of magma-water interaction and quenching rates [e.g. White and Valentine 2016]. Overall, the consistent partial degassing of the Fur ashes is a strong indicator that the entire section is the result of hydromagmatic activity.

Estimated total mass of emitted SO 2
Using the modelled values for initial sulfur concentrations, we have also estimated the total mass of sulfur emitted into the atmosphere from each eruption. We do not consider the issue of excess sulfur, as this problem is negligible or absent in Hawaiian and Icelandic basaltic magmas [Sharma et al. 2004]. We calculate the total mass by taking the difference between initial and residual mass of sulfur and multiplying by the erupted volume, following the approach of Thórdarson et al. [1996] (Equation 2):   [Smythe et al. 2017]. Modelled using an average Northern Mid-Atlantic Ridge (NMAR) sulfide liquid composition [Keith et al. 2017], average silicate composition for each Fur ash layer, temperature of 1130°C, and pressure 0.1 GPa. Vertical lines indicate error (1σ) for each modelled SCSS. Data in Table S3, Supplementary Material. where m t is the total emitted mass of SO 2 , m i the initial mass of SO 2 , m r the residual mass of SO 2 , V the eruptive volume (DRE) of NAIP tephra, ρ the magma density (2750 kg m´3), and υ SO2 the initial (i) and residual (r) mass fraction of SO 2 . We use the average content of sulfur for the whole dataset; converted to SO 2 , which is the species assumed to be present in the magma. This amounts to a modelled initial sulfur content of 0.38˘0.10 wt% SO 2 , and a measured residual sulfur content of 0.12˘0.04 wt% SO 2 . The rhyolitic Ash`19 is one of the thickest (20 cm) and most widely distributed ashes, and has an estimated total volume of 1200 km 3 [DRE; Egger and Brückl 2006]. Although this ash was rhyolitic and potentially more explosive, basaltic layers have been traced as far as the rhyolitic suggesting a similar distribution. We define Ash`19 as the upper limit for erupted volume, and apply a potential range of eruptive volumes (V) for each eruption of 100-1000 km 3 (DRE). However, calculating DRE equivalents is complicated by the competing effects of vesicularity and post-depositional compaction, so these values needs to be treated carefully. Keeping in mind the many uncertainties with these calcula-tions, this gives us a potential range of 0.72˘0.18 Gt to 7.2˘1.8 Gt SO 2 for each eruption. Using a similar approach, Blake et al. [2010] estimated a total yield of the Grande Ronde basalts of the Columbia River Basalts (a smaller and more recent LIP than the NAIP) of 1000 Gt delivered in intermittent bursts, each between <1 and 30 Gt SO 2 . Considering that these estimates are from fully-degassed flood basalts, in contrast to the only partially degassed ashes of the NAIP, our results seem comparable to previous estimates.
On average over the past decade, the volcanic SO 2 sources consistently detected from space have discharged a total of~63 kt d´1 SO 2 during passive degassing, or~23˘2 Mt yr´1 [Carn et al. 2017]. The estimated sulfur emissions from the NAIP would therefore significantly increase the total atmospheric sulfur content, with potentially large environmental consequences. A recent study also from the Fur Formation in Denmark has documented a potential sea surface temperature cooling using the organic palaeothermometer TEX 86 ]. Post-PETM cool conditions were initiated just below Ash´19 and persisted to varying degrees throughout the eruptions of the main phase of explosive basaltic volcanism, raising the question whether there was a causal relationship. A modelling study of the Roza member of the Columbia River Basalts found that a decade long eruption with annual SO 2 emissions of 1.2 Gt could lead to a´3.0 K temperature decrease [Schmidt et al. 2016]. Considering the magnitude of the Fur Formation ashes and the estimated sulfur emissions up to 1.8 Gt per eruption transported to the stratosphere, it is likely that these eruptions would have instigated a powerful surface cooling [Jones et al. 2005]. However, it is important to stress that there are a lot of assumptions and uncertainties regarding these estimates, and they should therefore be interpreted with caution.

Eruption depth and implications for the NAIP
We have shown that the sulfur content preserved in Fur glasses is controlled by variable amounts of degassing, which can be explained by arrested volatile loss due to magma-water interaction and quenching. Most of the ash layers are therefore likely erupted during hydromagmatic activity. In addition, we have shown that there is stratigraphic variation with a potential minor overall upward increase in the proportion of dense glassy grains and in the total residual sulfur content ( Figure 10). Assuming equilibrium degassing, and if the extent of degassing varies as a function of quench pressure, then these variations may reflect changes in the eruption depth [e.g. Moore and Schilling 1973]. However, the potential for disequilibrium degassing due to rapid ascent of clasts mean that other processes associated with thermal gradients in the eruptive jet could also explain these variations. During the 1963 eruption of Surtsey, uprush jets maintained hot magmatic temperatures in the core while cooling inwards from the margins [Moore 1985;Þórarinsson 1967]. Differential clast quenching rates under these conditions would produce a range of residual glass S contents, even if the depth of magma-water interaction remained constant [Mastin 2007;Mastin et al. 2004]. Changes in the size or speed of eruptive jets through time, and thus changes in the extent of equilibrium versus disequilibrium degassing, could therefore also account for the stratigraphic variations in residual sulfur content observed in the Danish ash series. Though, keep in mind it is difficult to assess the importance of this effect for the Fur ashes, given the extreme differences in scale and intensity between these eruptions and Surtseyan events.
Several studies have shown a link between the residual volatile content and depth of magma-water interaction [Davis et al. 2003;Liu et al. 2018;Moore and Schilling 1973;Unni and Schilling 1978]. Even accounting for the likelihood of some degree of disequilibrium degassing, a quantitative assessment of volatile saturation pressures can provide an upper bound on eruption depth when combined with an appropriate pressure gradient. Hydromagmatic activity can occur deep in the vent, in which case the residual volatiles would reflect both hydrostatic and lithostatic overpressure. If magma-water interaction involves a groundwater incursions into the volcanic conduit, then lithic particles are typically a large proportion of the erupted material due to explosive disruption and evacuation of subsurface country rock [Lorenz 1986;Morrisey et al. 2000]. The very low fraction of lithic clasts within the Fur ashes (<1 %; Figure 9) supports a shallow magmawater interaction involving surface water, rather than deep groundwater [e.g. Graettinger et al. 2013].
Quantitative assessments of eruption depth relies on thermodynamic modelling of residual H 2 O and CO 2 contents, for which pressure-and temperaturedependent solubilities are well-constrained [Newman and Lowenstern 2002;Witham et al. 2012]. This multispecies approach was applied to the Hverfjall lavas in Iceland, where hydrostatic pressures corresponding to an average fragmentation depth of 210˘30 m (water depth) were estimated based on a combination of H 2 O, CO 2 , and S [Liu et al. 2018]. Due to a lack of H 2 O of CO 2 measurements, we instead present a qualitative comparison to similar volcanic systems to infer an upper bound fragmentation (quench) depth. The Fur ashes have slightly lower sulfur concentrations than seen at Hverfjall (Figure 13), suggesting a potentially shallower water depth. Dredged basalts along the Reykjanes Ridge show that at a depths <200 m the sulfur concentrations decrease from 843 ppm S down to a minimum of 425 ppm S at 43 m depth [Moore and Schilling 1973]. This suggests that sulfur exsolution takes place very rapidly at low pressures in tholeiitic Icelandic basalts. Vesiculation is also increasingly variable with decreasing pressure, with some samples showing up to 50 % vesicularity at depths <200 m [Moore and Schilling 1973]. A similar shallow-marine subaqueous environment (water depths <200 m)with variable vesiculation at the point of magmawater interaction and variation in quench rate following fragmentation-would together explain the range of elevated glass S contents and ash textures that we observe.
A shallow-marine subaqueous environment suggest that the ashes erupted during the opening of the northeast Atlantic Ocean. Gradual submergence of the NAIP through time is documented on both sides of the opening northeast Atlantic rift [Á Horni et al. 2017]. Massive hyaloclastites within seaward dipping reflectors on the Vøring margin [Abdelmalak et al. 2016] (Figure 1) and offshore of the Faroe Islands [Jerram et al. 2009] reflect a gradual transition from subaerial to submarine basalt flows. The age of the positive ash series deposits (<55.6 to >54.6 Ma) also suggest that the Danish ashes coincide with the start of seafloor spreading, and post-date much of the voluminous flood basalts in East Greenland and the Faroes [Saunders 2016;Wilkinson et al. 2017]. The G10ka tephra series from Grímsvötn Evidence of explosive hydromagmatic eruptions Stokke et al., 2020 in Iceland demonstrates that wide distribution of hydromagmatic ash (in this case subglacial) is plausible [Óladóttir et al. 2020]. However, while Larsen et al. [2003] suggest that the Danish ashes are sourced from somewhere within the opening rift (Figure 1), it is not possible to locate the exact source within this rift system. We still speculate that both the initiation and the eventual termination of explosive volcanic activity could at least partly be a result of a local marine transgression during seafloor spreading and opening of the northeast Atlantic Ocean. Indeed, the subtle increase in residual volatile contents (reduced degassing) and high fraction of dense glass grains toward the top of the stratigraphy (Figure 9) are consistent, at least conceptually, with gradually increasing overlying pressure arresting volatile degassing and thus suppressing explosive magma-water interaction. Evolving eruptive style with changing confining pressure has been reported elsewhere. Although on a much smaller scale, a reduction in overlying hydrostatic pressure was invoked to explain transitions from pillow basalts into voluminous basaltic tephras at Askja volcano in Iceland [Graettinger et al. 2013] and during the early stages of the 1963 Surtsey eruption [Moore 1985], with no evidence for changing geochemistry or magma flux.
Previous studies have argued that the early Eocene shift to explosive volcanic activity during the emplacement of the NAIP was caused primarily by an increased magma flux [e.g. Saunders 2016]. Substantial magma fluxes would no doubt be necessary to produce the observed volumes and large distribution of the ashes, but magmatic fluxes were already exceptionally high during the late Paleocene emplacement of the continental flood basalts in East Greenland [Larsen and Tegner 2006;Wotzlaw et al. 2012]. We argue that the shift to an overall explosive volcanic environment could have been caused by a change to a subaqueous hydromagmatic environment with no need for a sudden increase in magma flux. We hypothesize that large-scale magmawater interactions would in fact be essential to achieve the unusually wide distribution of tholeiitic basaltic tephra observed in our study.

Conclusion
In this study we analysed the matrix glass in basaltic ash layers deposited in early Eocene sedimentary strata in Denmark. These ashes were sourced from the North Atlantic Igneous Province during the opening of the northeast Atlantic Ocean. The thickness and form of the ash layers, coupled with the large distances from possible source volcanoes suggest they erupted during explosive eruptions of unusually large magnitude and intensity, especially for basaltic magma compositions. The exceptional preservation of these layers provides a unique insight into volcanic processes that occur during plate tectonic breakup and LIP activity. Image anal-ysis of internal texture and external morphology show that the ashes are composed of 25-65 % well-preserved glass grains, with a minor amount (<5 %) of microcrystalline and lithic grains. The pristine glass fraction is dominated by dense, low vesicularity glass grains, relative to vesicular glass and shards. Geochemical analyses of pristine matrix glasses show an overall correlation between elevated residual sulfur concentrations and high fraction of dense glassy grains. Using initial sulfur concentrations calculated using theoretical models of sulfur content at sulfide saturation (SCSS), we find that these glasses are only partly degassed in sulfur (55-80 %). We also find that sulfur degassing during eruption of these ashes could potentially have led to emission of 0.72˘0.18 to 7.2˘1.8 Gt SO 2 for each eruption, with potentially large climatic consequences.
These data indicate that the glass quenched at pressures greater than expected for equilibration under atmospheric conditions, such as under hydrostatic pressure. Together, the fine grain size, dominance of dense angular ash morphologies, and elevated glass volatile contents suggest a hydromagmatic origin for the Danish ashes. Explosive hydromagmatic activity enhanced fine fragmentation and provided the high rates of heat transfer required to promote widespread distribution of large ash volumes. By comparing the residual sulfur concentrations with other basaltic hydromagmatic eruptions, we suggest that the eruptions likely occurred in a shallow marine subaqueous environment (<200 m water depth) during the opening of the northeast Atlantic Ocean. We propose that explosive hydromagmatic activity was the main driver of these unusually large explosive basaltic eruptions.