Eruptive history and ⁴⁰Ar∕³⁹Ar geochronology of the Milos volcanic field, Greece

High-resolution geochronology is essential for determining the growth rate of volcanoes, which is one of the key factors for establishing the periodicity of volcanic eruptions. However, there are less high-resolution eruptive histories (> 10⁶ years) determined for long-lived submarine arc volcanic complexes than for subaerial complexes, since submarine volcanoes are far more difficult to observe than subaerial ones. In this study, high-resolution geochronology and major-element data are presented for the Milos volcanic field (VF) in the South Aegean Volcanic Arc, Greece. The Milos VF has been active for over 3 Myr, and the first 2 × 10⁶ years of its eruptive history occurred in a submarine setting that has been emerged above sea level. The long submarine volcanic history of the Milos VF makes it an excellent natural laboratory to study the growth rate of a long-lived submarine arc volcanic complex. This study reports 21 new high-precision ${}^{\mathrm{40}}\mathrm{Ar}{/}^{\mathrm{39}}\mathrm{Ar}$ ages and major-element compositions for 11 volcanic units of the Milos VF. This allows us to divide the Milos volcanic history into at least three periods of different long-term volumetric volcanic output rate (Q_e). Periods I (submarine, ∼ 3.3–2.13 Ma) and III (subaerial, 1.48 Ma–present) have a low Q_e of 0.9 ± 0.5 × 10⁻⁵ and 0.25 ± 0.05 × 10⁻⁵ km³ yr⁻¹, respectively. Period II (submarine, 2.13–1.48 Ma) has a 3–12 times higher Q_e of 3.0 ± 1.7 × 10⁻⁵ km³ yr⁻¹. The Q_e of the Milos VF is 2–3 orders of magnitude lower than the average for rhyolitic systems and continental arcs.

1 Introduction

Short-term eruptive histories and compositional variations in lavas and pyroclastic deposits of many arc volcanic fields are well established. However, high-resolution eruptive histories that extend back > 10⁵–10⁶ years have been determined only for a handful of long-lived subaerial arc volcanic complexes. Some examples are Mount Adams (Hildreth and Lanphere, 1994), Tatara–San Pedro (Singer et al., 1997), Santorini (Druitt et al., 1999), Montserrat (Cole et al., 2002), Mount Baker (Hildreth et al., 2003a), Katmai (Hildreth et al., 2003b) and Ceboruco–San Pedro (Frey et al., 2004). To establish the growth rate of volcanic complexes and disentangle the processes responsible for the eruption, fractionation, storage and transport of magmas over time, comprehensive geological studies are required. These include detailed field mapping, sampling, high-resolution geochronology and geochemical analysis. Based on these integrated studies, the growth rate of volcanoes can be determined to establish the periodicity of effusive and explosive volcanism.

The Milos volcanic field (VF) is a long-lived volcanic complex that has been active for over 3 Myr. The Milos VF erupted for a significant part of its life below sea level, similar to the other well-studied volcanic structures in the eastern Mediterranean (e.g. Vougioukalakis et al., 2019). The eruptive history of the Milos VF has been examined with a broad range of chronostratigraphic techniques such as K–Ar, U–Pb, fission track ¹⁴C and biostratigraphy (e.g. Angelier et al., 1977; Fytikas et al., 1976, 1986; Traineau and Dalabakis, 1989; Matsuda et al., 1999; Stewart and McPhie, 2006; Van Hinsbergen et al., 2004; Calvo et al., 2012). However, most of the published ages have been measured using the less precise K–Ar or fission track methods, and modern, high-precision ${}^{\mathrm{40}}\mathrm{Ar}{/}^{\mathrm{39}}\mathrm{Ar}$ ages for the Milos VF have not been published so far. In this study, (1) we provide high-precision ${}^{\mathrm{40}}\mathrm{Ar}{/}^{\mathrm{39}}\mathrm{Ar}$ geochronology of key volcanic units of the Milos VF and (2) refine the stratigraphic framework of the Milos VF with the new high-precision ${}^{\mathrm{40}}\mathrm{Ar}{/}^{\mathrm{39}}\mathrm{Ar}$ ages and major-element composition. (3) We also quantify and constrain the compositional and volumetric temporal evolution of volcanic products of the Milos VF.

1.1 Geological setting

The Milos VF is part of the South Aegean Volcanic Arc (SAVA), an arc which was formed in the eastern Mediterranean by subduction of the African plate beneath the Aegean microplate (Fig. 1; Nicholls, 1971; Spakman et al., 1988; Duermeijer et al., 2000; Pe-Piper and Piper, 2007; Rontogianni et al., 2011). The present-day Benioff zone is located approximately 90 km underneath Milos (Hayes et al., 2018). The upper plate is influenced by extensional tectonics (e.g. McKenzie, 1978; Pe-Piper and Piper, 2013), which is evident on the island of Milos as horst and graben structures (Fig. 2).

Figure 1 Map of the South Aegean Volcanic Arc (SAVA). Red triangles indicate volcanic fields (VFs): the Susaki, Methana and Milos VFs in the western SAVA, Santorini VF in the centre, and Nisyros VF in the eastern SAVA. Red contour lines show the depth to the Benioff zone (Hayes et al., 2018). The white arrow represents the GPS-determined plate velocity of the Aegean microplate relative to the African plate from Doglioni et al. (2002).

Figure 2 Distribution of the proximal and medial facies of the submarine pumice-cone/crypto-dome volcanoes, submarine, submarine–subaerial and subaerial domes, and rhyolitic complexes (tuff cone and associated lava) of Milos, modified after Fytikas et al. (1986) and Stewart and McPhie (2006). The distal facies of Stewart and McPhie (2006) is not shown.

The Milos VF is exposed on the islands of the Milos archipelago: Milos, Antimilos, Kimolos and Polyegos. The focus of this study is Milos, which has a surface area of 151 km². The geology and volcanology of Milos have been extensively studied in the last 100 years. The first geological map was produced by Sonder (1924). This work was extended by Fytikas et al. (1976) and Angelier et al. (1977) and the subsequent publications of Fytikas et al. (1986) and Fytikas (1989). Interpretations based on volcanic facies of the complete stratigraphy were made by Stewart and McPhie (2003, 2006). More detailed studies of single volcanic centres (e.g. the Bombarda volcano and Fyriplaka complex) were published by Campos Venuti and Rossi (1996) and Rinaldi and Venuti (2003). Milos has also been extensively studied for its epithermal gold mineralization, summarized by Alfieris et al. (2013). Milos was known during the Neolithic period for its export of high-quality obsidian. Today the main export product is kaolinite mined from hydrothermally altered felsic volcanic units in the centre of the island (e.g. Alfieris et al., 2013).

Figure 3 Simplified geological map of Milos with ${}^{\mathrm{40}}\mathrm{Ar}{/}^{\mathrm{39}}\mathrm{Ar}$ ages of this study and sample locations of key volcanic deposits, modified after Stewart and McPhie (2006) and Grasemann et al. (2018). The stratigraphic units of Milos are from Fytikas et al. (1986). Age data from this study are in black; published ages are shown in red (Angelier et al., 1977; Fytikas et al., 1986; Traineau and Dalabakis, 1989; Stewart and McPhie, 2006). The “green lahar” (Fytikas, 1977) consists of deposits from multiple phreatic explosions and contains fragments of metamorphic, sedimentary and volcanic rocks.

The geology of Milos can be divided into four main units: (1) metamorphic basement, (2) Neogene sedimentary rocks, (3) volcanic sequences and (4) the alluvial cover. The metamorphic basement crops out at the south-west, south and south-east of Milos (Fig. 3) and is also found as clasts in many volcanic units. The metamorphic rocks include lawsonite-free jadeite eclogite, lawsonite eclogite, glaucophane schist, quartz–muscovite–chlorite and chlorite–amphibole schist (Fytikas et al., 1976, 1986; Grasemann et al., 2018; Kornprobst et al., 1979). The exposed units belong to the Cycladic Blueschist Unit (Lower Cycladic nappe), whereas eclogite pebbles in the phreatic eruption products called “green lahar” by Fytikas (1977) are derived from the Upper Cycladic nappe (Grasemann et al., 2018).

On top of this metamorphic basement, Neogene fossiliferous marine sedimentary rocks were deposited (e.g. Van Hinsbergen et al. 2004). This sedimentary sequence can be divided into a lower unit A and upper unit B that is unconformably overlain by volcaniclastic sediments (Van Hinsbergen et al., 2004). Unit A is 80 m thick and consists of fluviatile–lacustrine, brackish and shallow marine conglomerate, sandstone, dolomite and limestone. Unit B is 25–60 m thick and consists of sandstone overlain by a succession of alternating marls and sapropels, suggesting a deeper marine setting (Van Hinsbergen et al., 2004). Five volcanic ash layers that contain biotite are found in this Neogene sedimentary sequence, either suggesting that volcanic eruptions in small volume already occurred in the Milos area or that these ash layers are derived from larger eruptions of volcanic centres further away from Milos (van Hinsbergen et al., 2004). Age determinations by bio-magneto- and cyclo-stratigraphy suggested that deposition of Unit A started at approximately 5 Ma, and that Milos subsided 900 m in 0.6 Myr (Van Hinsbergen et al. 2004) due to extension. This subsidence happened ca. 1.0–1.5 Myr before the onset of the main phase of Pliocene–recent volcanism on Milos.

The Pliocene–recent volcanic sequence of Milos has been subdivided into different units by Angelier et al. (1977) and Fytikas et al. (1986). In addition, Stewart and McPhie (2006) provided a detailed facies analysis of the different volcanic units. The subdivision by Angelier et al. (1977) is not constrained well due to their limited amount of age data. The subdivision of volcanic units by Fytikas et al. (1986) and facies descriptions of Stewart and McPhie (2006) are summarized below. It is important to note that according to Stewart and McPhie (2006), the five volcanic cycles described by Fytikas et al. (1986) are difficult to match with existing age data and the continuous progression in volcanic construction (Fig. 4). For example, the first phase of Fytikas et al. (1986), the Basal Pyroclastic Series, contains the large pumice-cone/crypto-dome volcanoes according to Stewart and McPhie (2006). Two of these pumice-cone/crypto-dome volcanoes are much younger and intercalated between the Complex of Domes and Lava Flows (CDLF) of Fytikas et al. (1986).

Figure 4 Previous proposed stratigraphic frameworks for Milos by Angelier et al. (1977), Fytikas et al. (1986), and Stewart and McPhie (2006). Volcanic unit II of Angelier et al. (1977) contains unit I. Stewart and McPhie (2006) described the volcanic facies of Milos mainly based on the geochronological studies of Angelier et al. (1977) and Fytikas et al. (1986). Abbreviation: SFCPCV – submarine felsic pumice-cone/crypto-dome volcanoes.

The first volcanic unit deposited in the Milos area is the Basal Pyroclastic Series (BPS) (Fytikas et al., 1986) or submarine felsic pumice-cone/crypto-dome volcanoes (Stewart and McPhie, 2006, Figs. 2–4). This unit consists of thickly bedded pumice breccia with a rhyolitic–dacitic composition. These rhyolites–dacites are aphyric or contain quartz–feldspar ± biotite phenocrysts. Graded sandstone and bioturbated and fossil-rich (in situ bivalve shells) mudstone are intercalated, indicating a marine environment and a water depth of several hundreds of metres (e.g. Stewart, 2003; Stewart and McPhie, 2006), whereas later degassed magmas with a similar composition intruded as sills and crypto-domes. The BPS has been strongly affected by hydrothermal fluids, especially the proximal deposits (e.g. Kilias et al., 2001).

The second volcanic unit was named the Complex of Domes and Lava Flows (Fytikas et al., 1986), and the volcanic facies of this unit are described as submarine dacitic and andesitic domes by Stewart and McPhie (2006). This phase of effusive submarine volcanism was predominantly andesitic and dacitic in composition and produced microcrystalline rocks with phenocrysts of pyroxene, amphibole, biotite and plagioclase. The eruption centres were mainly located along NNE faults and formed up to 300 m thick deposits extending over areas of 2.5 to 10 km² around the eruption centres. In the north-eastern part of Milos, an andesitic scoria cone provided scoria lapilli and bombs to deeper water settings. Sandstone intercalated in the CDLF contains both igneous and metamorphic minerals suggesting input from the basement. Rounded pebbles of rhyolite and dacite indicate that some of the volcanic deposits were above sea level or in very shallow, near-shore environments (e.g. Stewart and McPhie, 2006).

The third volcanic unit is called the Pyroclastic Series and Lava Domes (PSLD) by Fytikas et al. (1986) and belongs to the submarine-to-subaerial dacitic and andesitic lava domes of Stewart and McPhie (2006). This highly variable group is dominated by rhyolitic, dacitic and andesitic lavas, domes, pyroclastic deposits and felsic pumiceous sediments (Stewart and McPhie, 2006). Thickness varies between 50–200 m, and the deposits are located in the eastern and northern parts of Milos (Figs. 2 and 3). The initial pyroclastic layers were subaqueously deposited and the extrusion of a dome resulted in the deposition of talus around the margins by mass flow. On top of the dome sand- and siltstone with fossils (Ostrea fossil assemblage) and traction–current structures suggest that the top of the dome was above wave base. The youngest deposits of this unit are dacitic and andesitic lavas and domes. These domes generated subaerial block-and-ash flow and surge deposits. Paleosols within these deposits are a clear indicator that some areas were above sea level. The last unit of the PSLD is represented by large subaerial rhyolitic lava that contains quartz and biotite phenocrysts and is found near Halepa in the southern central part of Milos.

The fourth unit consists of the subaerially constructed rhyolitic complexes of Trachilas and Fyriplaka (CTF) (Fytikas et al., 1986), which Stewart and McPhie (2006) interpreted as subaerial rhyolitic lava–pumice cones. These two volcanic complexes are built from rhyolitic pumice deposits and lavas that contain quartz and biotite phenocrysts (10 modal %–20 modal %). The deposits have a maximum thickness of 120 m and decrease to several metres' thickness in the distal parts. Basement-derived schist is found as lithic clasts (Fytikas et al., 1986). In addition, the Kalamos rhyolitic lava dome, which outcrops on the southern coast of Milos, produced lava that spread westwards to the Fyriplaka beach (Fig. 2). This lava belongs to this fourth phase and is probably derived from an older volcano and not the Fyriplaka complex (Campos Venuti and Rossi, 1996).

The fifth volcanic unit comprises deposits from phreatic activity, especially in the northern part of the Zefiria Graben and near Agia Kiriaki (Fig. 2 of Stewart and McPhie, 2006). Many overlapping craters are surrounded by lithic breccias that are composed of variably altered metamorphic basement clasts and volcanic clasts. This phreatic activity has continued into historic times (Traineau and Dalabakis, 1989). Fytikas et al. (1986) referred to this unit as “green lahar”, although it is indicated that this deposit is not a lahar but the product of phreatic eruptions in the last 0.2 Myr.

1.2 Previous geochronological studies

Previous geochronological work is summarized in Table 1. Angelier et al. (1977) reported six K–Ar ages (0.95–2.50 Ma). These ages were used in combination with field observations to divide the Milos volcanic succession into four units. However, the samples from Fyriplaka, the fourth unit, were too young to be dated by Angelier et al. (1977). Fytikas et al. (1976, 1986) published 16 K–Ar ages for Milos (0.09–3.50 Ma) including an age of 0.09–0.14 Ma for the Fyriplaka complex. Fytikas et al. (1986) also obtained three K–Ar ages for Antimilos (0.32 ± 0.05 Ma), Kimolos (3.34 ± 0.06 Ma) and Polyegos (2.34 ± 0.17 Ma). Traineau and Dalabakis (1989) dated the very young phreatic deposits by ¹⁴C dating and found ages between 200 BCE and 200 CE. Matsuda et al. (1999) published two K–Ar ages of 0.8 ± 0.1 (MI-1) and 1.2 ± 0.1 Ma (MI-4) for the Plakes dome that was also studied by Fytikas et al. (1986). Bigazzi and Radi (1981) published two fission track ages of 1.54 ± 0.18 and 1.57 ± 0.15 Ma for obsidians of Bombarda–Adamas and Dhemeneghaki, respectively. Later fission track studies by Arias et al. (2006) (1.57 ± 0.12 and 1.60 ± 0.06 Ma) confirmed these ages. The fission track ages are younger than the K–Ar ages given by Angelier et al. (1977; 1.84 ± 0.08 Ma for Dhemeneghaki) and Fytikas et al. (1986; 1.71 ± 0.05 Ma for Bombarda). In the most recent geochronological study of the Milos VF, Stewart and McPhie (2006) published four SHRIMP U–Pb zircon ages: Triades dacite facies (1.44 ± 0.08 and 2.18 ± 0.09 Ma), Kalogeros crypto-dome (2.70 ± 0.04 Ma) and the Filakopi Pumice Breccia (2.66 ± 0.07 Ma). All uncertainties reported here are 1 standard deviation uncertainties as reported in the original publications, except for the ¹⁴C ages for which uncertainties were not specified.

Table 1 Published eruption ages of stratigraphic units of the island of Milos.

Angelier et al. (1977) do not provide sample names but only numbers for the sample locations. Here the location is given after “Angelier_” (Angelier et al. 1977, their Fig. 3).
Abbreviations are as follows: BPS – Basal Pyroclastic Series; CDLF – Complex of Domes and Lava Flows; PSLD – pyroclastic series and lava domes; CTF – complexes of
Trachilas and Fyriplaka. See more details in Fig. 4.

The previous geochronological work for the Milos VF is mainly based on K–Ar ages. However, K–Ar ages may show undesirable and unresolvable scatter due to various problems: (1) inaccurate determination of radiogenic argon due to either incorporation of excess argon or incomplete degassing of argon during the experiments; (2) inclusion of cumulate or wall rock phenocrysts in bulk analyses; (3) disturbance of a variety of geological processes such as slow cooling or thermal reheating; (4) unrecognized heterogeneities due to separate measurements of potassium and argon content by different methods; (5) requirement of relatively large quantities (milligrams) of pure sample (e.g. Lee, 2015). In addition to these methodological issues, in the case of Milos we observe that hydrothermal alteration caused substantial kaolinitization, in particular of the felsic volcanic samples, that most likely has affected the K–Ar systematics. Some of these issues are also valid for the ${}^{\mathrm{40}}\mathrm{Ar}{/}^{\mathrm{39}}\mathrm{Ar}$ method. However, the K–Ar method does not allow testing whether ages are compromised.

${}^{\mathrm{40}}\mathrm{Ar}{/}^{\mathrm{39}}\mathrm{Ar}$ ages only need isotopes of argon to be measured from a single aliquot of sample with the same equipment that can eliminate some of the problems with sample inhomogeneity. Furthermore, step heating and multiple single-fusion experiments can shed light on sample inhomogeneity due to partial alteration effects. The high sensitivity of modern noble gas mass spectrometers for ${}^{\mathrm{40}}\mathrm{Ar}{/}^{\mathrm{39}}\mathrm{Ar}$ measurements results in very small sample amounts needed for analysis, which can yield more information on the thermal or alteration histories than larger samples. Moreover, other argon isotopes (³⁶Ar, ³⁷Ar and ³⁸Ar) can be used to infer some information about the chemical compositions (i.e. Ca and Cl) of samples. A high-resolution laser incremental heating method of ${}^{\mathrm{40}}\mathrm{Ar}{/}^{\mathrm{39}}\mathrm{Ar}$ dating allows us to resolve the admixture of phenocryst-hosted inherited ⁴⁰Ar in the final temperature steps of the incremental-step heating experiments.

2 Methods

2.1 Mineral separation and sample preparation

Samples were collected from all major volcanic units on Milos island based on the studies of Fytikas et al. (1986), Stewart and McPhie (2006), and our own observations in the field. Photos of the sample locations and thin sections can be found in Supplement file I. Approximately 2 kg of fresh pumice clasts or lava was sampled from each unit. Samples were cut into ∼ 5 cm³ cubes using a diamond saw to remove potentially altered surfaces and obtain the fresh interior parts. These cubes were ultra-sonicated for 30 min in demi-water to remove dust and seawater and dried in an oven overnight at 50 ^∘C. Dry sample cubes were crushed in a steel jaw crusher, and this fraction was split into two portions of roughly equal size. One of them was powdered in an agate shatter box and agate ball mill to a grain size of less than 2 µm for the major-element analysis. The second fraction was sieved to obtain a grain size of 250–500 µm for ${}^{\mathrm{40}}\mathrm{Ar}{/}^{\mathrm{39}}\mathrm{Ar}$ dating.

Heavy-liquid density separation techniques (IJlst, 1973) were used to purify mineral separates (groundmass, biotite, amphibole) required for the ${}^{\mathrm{40}}\mathrm{Ar}{/}^{\mathrm{39}}\mathrm{Ar}$ dating. Different densities of heavy liquids were used to obtain groundmass (2700 $\le \mathit{\rho }\le$ 3000 kg m⁻³), biotite (2900 $\le \mathit{\rho }\le$ 3100 kg m⁻³) and amphibole (∼ 3100 $\le \mathit{\rho }\le$ 3200 kg m⁻³). A Frantz isodynamic magnetic separator was used to remove the magnetic minerals from the non-magnetic minerals and groundmass. The samples for ${}^{\mathrm{40}}\mathrm{Ar}{/}^{\mathrm{39}}\mathrm{Ar}$ analysis were purified by handpicking under a binocular optical microscope to select mineral grains without visible alteration and inclusions.

2.2 ${}^{\mathrm{40}}\mathrm{Ar}{/}^{\mathrm{39}}\mathrm{Ar}$ dating

The mineral and groundmass samples were wrapped in either 6 or 9 mm aluminium foil and packed in 20 mm aluminium cups, which were vertically stacked. Based on stratigraphy and previous geochronological constraints > 1 Ma samples and the < 1 Ma samples were irradiated for 7 and 1 h, respectively, in irradiation batches VU108 and VU110 in the Cadmium-Lined In-Core Irradiation Tube (CLICIT) facility of the Oregon State University Training Research, Isotopes, General Atomics (TRIGA) reactor. The neutron flux for all irradiations was monitored by standard bracketing using the Drachenfels sanidine (DRA; 25.52 ± 0.08 Ma, modified from Wijbrans et al., 1995, and calibrated relative to Kuiper et al., 2008) and Fish Canyon Tuff sanidine (FCs; 28.201 ± 0.023 Ma, Kuiper et al., 2008) with Min et al. (2000) decay constants.

In total, 24 samples (8 groundmasses, 15 biotites and 2 amphiboles; for sample G15M0026 both biotite and amphibole were analysed) were measured by either ${}^{\mathrm{40}}\mathrm{Ar}{/}^{\mathrm{39}}\mathrm{Ar}$ fusion and/or incremental heating techniques. For incremental heating experiments, 80–100 grains per sample were loaded into a 25-hole (surface per hole ∼ 36 mm²) copper tray together with single-grain standards in ∼ 12 mm² holes. The tray was prebaked in a vacuum (10⁻⁵–10⁻⁶ mbar) at 250 ^∘C overnight to remove atmospheric argon and subsequently baked overnight at 120 ^∘C in the ultra-high vacuum sample chamber (< 5 × 10⁻⁹ mbar) and purification system connected to a Thermo Scientific Helix multi-collector (MC) mass spectrometer.

Samples and standards were heated with a focused laser beam at 8 % power using a 50 W continuous-wave (CW) CO₂ laser. The released gas was cleaned by exposure to a cold trap cooled by a Lauda cooler at −70 ^∘C, a SAES NP10 at 400 ^∘C, a Ti sponge at 500 ^∘C and cold SAES ST172 Fe–V–Zr sintered metal. The five isotopes of argon were measured simultaneously on five different collectors: ⁴⁰Ar on the H2-Faraday, ³⁹Ar on the H1-Faraday or the H1-CDD, ³⁸Ar on the AX-CDD, ³⁷Ar on the L1-CDD, and ³⁶Ar on the L2-CDD for 15 cycles with 33 s integration time (CDD: compact discrete dynodes; AX: axial; H: high-mass side; L: low-mass side). The Faraday cups on H2 and H1 were equipped with amplifiers with 10¹³ Ω resistors in their feedback loop. Procedural blanks were measured every two or three analyses in different sequences, and air shots were measured every 8–12 h to correct the instrumental mass discrimination. The gain between different collectors was monitored by measuring CO₂ on mass 44 in dynamic mode on all collectors. Gain was generally stable over periods of weeks. Note that because samples, standards and air calibration runs are measured during the same period, gain correction does not substantially change the final age results. The raw mass spectrometer data output was converted by an Excel macro script designed in-house to be compatible with the ArArCalc 2.5 data reduction software (Koppers, 2002). The ${}^{\mathrm{40}}\mathrm{Ar}{/}^{\mathrm{36}}\mathrm{Ar}$ atmospheric air value of 298.56 from Lee et al. (2006) is used in the calculations. The correction factors for neutron interference reactions are (2.64 ± 0.02) × 10⁻⁴ for (${}^{\mathrm{36}}\mathrm{Ar}{/}^{\mathrm{37}}\mathrm{Ar}{)}_{\mathrm{Ca}}$, (6.73 ± 0.04) × 10⁻⁴ for (${}^{\mathrm{39}}\mathrm{Ar}{/}^{\mathrm{37}}\mathrm{Ar}{)}_{\mathrm{Ca}}$, (1.21 ± 0.003) × 10⁻² for (${}^{\mathrm{38}}\mathrm{Ar}{/}^{\mathrm{39}}\mathrm{Ar}{)}_{\mathrm{K}}$ and (8.6 ± 0.7) × x10⁻⁴ for (${}^{\mathrm{40}}\mathrm{Ar}{/}^{\mathrm{39}}\mathrm{Ar}{)}_{\mathrm{K}}$. All uncertainties are quoted at the 1σ level and include all analytical errors (i.e. blank, mass discrimination with neutron interference correction and analytical error in J-factor, the parameter associated with the irradiation process).

A reliable plateau age is defined as experiments with at least three consecutive steps overlapping at 2σ, containing > 50 % of the ³⁹Ar_K, with a mean square weighted deviate (MSWD) value < 2.5 and with an ${}^{\mathrm{40}}\mathrm{Ar}{/}^{\mathrm{36}}\mathrm{Ar}$ inverse isochron intercept that does not deviate from atmospheric argon at 2σ. All the inverse isochron ages used the same steps as used in the weighted mean ages, and all relevant analytical data for the age calculations following standard practices (Schaen et al., 2020) can be found in Supplement file II.

2.3 Whole-rock major-element analysis by XRF

Major-element concentrations were measured by X-ray fluorescence spectroscopy (XRF) on a Panalytical AxiosMax. A Panalytical Eagon2 was used to create 40 mm fused glass beads of ${\mathrm{Li}}_{\mathrm{2}}{\mathrm{B}}_{\mathrm{4}}{\mathrm{O}}_{\mathrm{7}}/{\mathrm{LiBO}}_{\mathrm{2}}$ (65.5:33.5 %, Johnson & Johnson Spectroflux 110) with a 1:6 sample–flux ratio that were melted at 1150 ^∘C. Sample powders were ignited at 1000 ^∘C for 2 h to determine loss on ignition (LOI) before being mixed with the ${\mathrm{Li}}_{\mathrm{2}}{\mathrm{B}}_{\mathrm{4}}{\mathrm{O}}_{\mathrm{7}}/{\mathrm{LiBO}}_{\mathrm{2}}$ flux. Interference-corrected spectra intensities were converted to oxide-concentrations against a calibration curve consisting of 30 international standards. The precision, expressed as the coefficient of variation (CV), is better than 0.5 %. The accuracy, as measured on the international standards AGV-2, BHVO-2, BCR-2 and GSP-2, was better than 0.7 % (1 relative standard deviation, RSD) (Supplement file III).

2.4 Eruption volume calculation

The minimum and/or maximum eruption volume of each volcano during each eruption period is derived from the ranges of thickness and surface areas that are reported in Campos Venuti and Rossi (1996) and Stewart and McPhie (2006). We converted these volumes to dense rock equivalent (DRE) based on the magma type of different deposits. This analysis only includes the onshore deposits and results in a smaller estimate for larger pyroclastic volumes. The DRE volume is calculated using the equation of Crosweller et al. (2012):

$$\begin{array}{}\text{(1)}& \mathrm{\text{DRE}}\left({\mathrm{km}}^{\mathrm{3}}\right)={\displaystyle \frac{\text{tephra vol}\left({\mathrm{km}}^{\mathrm{3}}\right)\text{tephra density}\left(\mathrm{kg}{\mathrm{m}}^{-\mathrm{3}}\right)}{\text{magma density}\left(\mathrm{kg}{\mathrm{m}}^{-\mathrm{3}}\right)}}.\end{array}$$

Tephra density is assumed to be 1000 kg m⁻³ (Crosweller et al., 2012). Magma density varies depending on the magma type. Here we used 2300 kg m⁻³ for rocks with a SiO₂ range of 65 wt %–77 wt % and 2500 kg m⁻³ for all samples with SiO₂ < 65 wt %. DRE corresponds to the unvesiculated erupted magma volume. Therefore, we did not convert the volume of some crypto-dome and lavas from Profitis Illias (G15M0017), Triades (G15M0021-24), Dhemeneghaki (G15M0032B) and Halepa (G15M0013) to the DRE since they contain less than 5 % vesicles.

3 Results

3.1 ${}^{\mathrm{40}}\mathrm{Ar}{/}^{\mathrm{39}}\mathrm{Ar}$ age results

In this section, we present our groundmass, biotite and amphibole ${}^{\mathrm{40}}\mathrm{Ar}{/}^{\mathrm{39}}\mathrm{Ar}$ results for 11 volcanic units of Milos. The ${}^{\mathrm{40}}\mathrm{Ar}{/}^{\mathrm{39}}\mathrm{Ar}$ ages range from 0.06 to 4.10 Ma and cover most of the major volcanic units of Milos. Tables 2 and 3 show the ${}^{\mathrm{40}}\mathrm{Ar}{/}^{\mathrm{39}}\mathrm{Ar}$ results of incremental heating steps and single-grain fusion analyses, respectively. Note that the Irr-ID column in these two tables represents the irradiation ID of the analytical experiment (e.g. VU108-, VU110-) and the top right superscripts (G, B, A, O) in the sample IDs (e.g. G15M0029^G, G15M0021^B) refer to groundmass, biotite, amphibole and obsidian.

Table 2 Incremental heating ${}^{\mathrm{40}}\mathrm{Ar}{/}^{\mathrm{39}}\mathrm{Ar}$ results of the Milos volcanic field.

The age in bold is regarded as the best estimate of the eruptive age. The ⁴⁰Ar^∗ (%) is the average radiogenic ⁴⁰Ar of the analyses included in the weighted mean. The
experiment was analysed on biotite^B, obsidian^O, amphibole^A and groundmass^G of a sample. The same steps were used for the calculation of isochron ages as used in the
weighted mean ages.

Table 3 ${}^{\mathrm{40}}\mathrm{Ar}{/}^{\mathrm{39}}\mathrm{Ar}$ results of single-grain fusion analyses on the Milos volcanic field.

The age in bold is regarded as the best estimate of the eruptive age. The ⁴⁰Ar^∗ (%) is the average radiogenic ⁴⁰Ar of the analyses included in the weighted mean. ^* The $\mathrm{K}/\mathrm{Ca}$
ratio is calibrated by removing the total fusion with excess ³⁷Ar (Ca) (fA > 1). ^B The experiment was analysed on biotite of the sample. The same steps were used for the
calculation of isochron ages as used in the weighted mean ages.

3.1.1 Groundmass ${}^{\mathrm{40}}\mathrm{Ar}{/}^{\mathrm{39}}\mathrm{Ar}$ plateau and/or isochron ages

All groundmass samples yielding ${}^{\mathrm{40}}\mathrm{Ar}{/}^{\mathrm{39}}\mathrm{Ar}$ plateau and isochron ages with more than 50 % ³⁹Ar_K and less than 2.5 MSWD included in their age spectrum are shown in Fig. 4 and reported in Table 2. The ${}^{\mathrm{40}}\mathrm{Ar}{/}^{\mathrm{36}}\mathrm{Ar}$ isochron intercepts do not deviate from atmospheric argon at the 2σ level, unless stated otherwise (Table 3). Sample G15M0016 was collected from a dyke at Kleftiko in the south-west of Milos (Fig. 2). Three incremental heating experiments were performed on the groundmass of this sample (Fig. 5a). The first experiment (VU108-Z8a) produced a weighted mean age of 2.71 ± 0.02 Ma (MSWD 2.31; ³⁹Ar_K 79.6 %; inverse isochron age 2.65 ± 0.10 Ma). The other two, VU108-Z8a_4 and VU108-Z8b_1, have plateau ages of 2.61 ± 0.03 Ma (MSWD 0.93; ³⁹Ar_K 57.4 %; inverse isochron age 2.69 ± 0.10 Ma) and 2.67 ± 0.01 Ma (MSWD 1.50; ³⁹Ar_K 65.57 %; inverse isochron age 2.55 ± 0.05 Ma), respectively. The three experiments are remarkably similar. Although the amount of radiogenic ⁴⁰Ar is low (< 20 %), a combined age of 2.66 ± 0.01 Ma is considered to be the best estimate with a relatively high MSWD value (2.51).

Figure 5 Groundmass ${}^{\mathrm{40}}\mathrm{Ar}{/}^{\mathrm{39}}\mathrm{Ar}$ plateau ages for samples G15M0016 (a), G15M0032B (b), G15M0019 (c) and G15M0020 (d). The Mavro Vouni dome (a), Dhemeneghaki volcano (b) and Kontaro dacitic dome (c, d) are located in, respectively, the south-western, north-eastern and eastern parts of the Milos VF (see Fig. 2). Final age calculation is reported with 1σ errors. See the individual steps of samples G15M0016, G15M0019 and G15M0029 in Supplement file II.

Two lava samples, G15M0019 and G15M0020, were collected from Kontaro in north-eastern Milos (Fig. 2). Three replicate incremental heating step experiments on groundmass from sample G15M0019 (VU108-Z6a_4; VU108-Z6a_5 and VU108-Z6b_1; Fig. 5b) were performed that are not reproducible. Their plateau ages range from 1.55 Ma to 1.62 Ma with relatively high MSWD (3.8–4.5), 56 %–95 % of the total ³⁹Ar_K, 34 %–53 % of radiogenic ⁴⁰Ar, 0.88–1.02 of $\mathrm{K}/\mathrm{Ca}$ and an atmospheric isochron intercept of 297–315. We regard the isochron age from the last experiment (VU108-Z6b_1) as the reliable age (1.48 ± 0.02 Ma, MSWD 0.44) because its MSWD value is the only one smaller than 2.5 in this experiment and therefore the best estimate for the eruption age. Three replicate incremental heating step experiments on groundmass from sample G15M0020 (VU108-Z5a_5; VU108-Z5b_1 and VU108-Z5b_2; Fig. 5c) were analysed. These experiments are similar at the lower-temperature heating steps. They produced statistically meaningful plateau ages ranging from 1.52–1.56 Ma with 41 %–62 % of the total ³⁹Ar_K, 18 %–48 % of radiogenic ⁴⁰Ar, 1.51–1.73 of $\mathrm{K}/\mathrm{Ca}$ and an atmospheric isochron intercept of 295–300. Their combined weighted mean age is 1.54 ± 0.01 Ma (MSWD 3.06; ³⁹Ar_K 57.32 %) with 25.31 % of ⁴⁰Ar^∗.

Sample G15M0032B (obsidian) was collected from a pumice-cone volcano at Dhemeneghaki (Fig. 2). One incremental heating experiment on this sample (VU108-Z18, Fig. 5d) yielded a plateau age of 1.825 ± 0.002 Ma (MSWD 0.91; ³⁹Ar_K 98.6 %). The ⁴⁰Ar^∗ is 93.86 %. The inverse isochron age is identical to the weighted mean plateau age of 1.825 ± 0.002 Ma. The age of 1.825 ± 0.002 Ma is considered the best estimate for the eruption age of the Dhemeneghaki obsidian.

3.1.2 Groundmass ${}^{\mathrm{40}}\mathrm{Ar}{/}^{\mathrm{39}}\mathrm{Ar}$ plateau and/or isochron ages (25 %–40 % ³⁹Ar_K released)

The results shown in Fig. 5 did not yield weighted mean plateau ages according to standard criteria including ³⁹Ar_K > 50 % but still provide some useful age information. Sample G15M0017 was collected from a crypto-dome of the Profitis Illias volcano in south-western Milos (Fig. 2). Three replicate incremental heating experiments – VU108-Z7a, VU108-Z7a_4 and VU108-Z7b_1 – were performed on this sample, which resulted in disturbed age spectra (Fig. 6a). The consecutive lower-temperature steps of all experiments define ages of < 2.5 Ma, which is much younger than the ages of the submarine pyroclastic products of the lower series at Kleftiko and/or Profitis Illias (3.0–3.5 Ma, Fytikas et al., 1986, and Stewart and McPhie, 2006). At the consecutive higher-temperature heating steps, these experiments yielded 3.64 ± 0.08 Ma (${}^{\mathrm{40}}\mathrm{Ar}{/}^{\mathrm{36}}\mathrm{Ar}$ 293.87 ± 4.77; VU108-Z7a), 4.10 ± 0.06 Ma (${}^{\mathrm{40}}\mathrm{Ar}{/}^{\mathrm{36}}\mathrm{Ar}$ 298.44 ± 15.51; VU108-Z7a_4) and 3.41 ± 0.05 Ma (${}^{\mathrm{40}}\mathrm{Ar}{/}^{\mathrm{36}}\mathrm{Ar}$ 295.97 ± 7.34; VU108-Z7b_1). The total fusion and inverse isochron ages of the three experiments gave large ranges of 2.25–3.23 and 3.68–4.14 Ma, respectively, and none of these high-temperature heating steps produced a statistical plateau (all MSWD > 2.0). The amount of radiogenic ⁴⁰Ar of both the ${}^{\mathrm{40}}\mathrm{Ar}{/}^{\mathrm{39}}\mathrm{Ar}$ result from our sample and the K–Ar age data from previous studies (Fytikas et al., 1986) is rather low (< 15 %) for a sample of this age based on our laboratory experience. Therefore, the estimated age range for the oldest volcanic products of the Milos VF should be confirmed by other dating techniques.

Figure 6 Groundmass ${}^{\mathrm{40}}\mathrm{Ar}{/}^{\mathrm{39}}\mathrm{Ar}$ plateau or inverse isochron ages for samples G15M0017 (a), G15M0015 (b) and G15M0029 (c). Individual steps and final age calculation are reported with 1σ errors. The Profitis Illias volcano (a, b) and dacitic Korakia dome (c) are located in the south-western and north-eastern parts of the Milos VF, respectively (Fig. 2). See the individual steps of samples G15M0015 and G15M0029 in Supplement file II.

Sample G15M0015 is also a crypto-dome breccia from Profitis Illias (Fig. 2). Two replicate incremental-step heating experiments were performed on the groundmass of this sample (VU108-Z9a and VU108-Z9b_1, Fig. 6b). The experiment VU108-Z9a groundmass shows a disturbed age spectrum and ages increase from ∼ 3 Ma in the initial heating steps to ∼ 3.2 Ma, followed by a decrease to ∼ 3 Ma in the high-temperature heating steps. The consecutive heating steps only exist at the lower-temperature steps, yielding a “plateau” of 3.12 ± 0.02 Ma (MSWD 9.07). Due to the excess argon (${}^{\mathrm{40}}\mathrm{Ar}{/}^{\mathrm{36}}\mathrm{Ar}$ 304.19 ± 1.25 comprising 43.07 % of the released ³⁹Ar_K), the inverse isochron of 3.06 ± 0.02 Ma (MSWD 0.01) is more reliable for this analysis. The inverse isochron age of the second groundmass (VU108-Z9b_1) is identical at 3.04 ± 0.02 Ma (MSWD 1.14; ³⁹Ar_K 27.00 %) and ${}^{\mathrm{40}}\mathrm{Ar}{/}^{\mathrm{36}}\mathrm{Ar}$ of 293.83 ± 1.38 obtained at high-temperature steps. The two experiments are remarkably similar. Although the sample does not formally fulfil the definition of a plateau age comprising > 50 % ³⁹Ar_K released, a combined age of 3.06 ± 0.02 Ma (MSWD 1.14; ³⁹Ar_K 22.79 %, ⁴⁰Ar^∗ 41.77 %) most likely represents the eruption age. This ${}^{\mathrm{40}}\mathrm{Ar}{/}^{\mathrm{36}}\mathrm{Ar}$ age is consistent with the K–Ar age from the same lithology of 3.08 ± 0.08 Ma (Fytikas et al. 1986).

Sample G15M0029 is an andesite collected from Korakia in the north-east of Milos (Fig. 2). Two incremental heating experiments (VU108-Z16a and VU108-Z16b_1, Fig. 6c) were performed on this sample. The two experiments are remarkably similar and show a decreasing age from ∼ 2.85 Ma at the lower-temperature heating steps to 2.65 Ma at the higher temperatures. The higher-temperature heating steps of both experiments yielded weighted mean plateau ages of 2.67 ± 0.01 Ma (MSWD 0.96; ³⁹Ar_K 23.61 %, ⁴⁰Ar^∗ 56.34 %; inverse isochron age 2.68 ± 0.02 Ma) and 2.69 ± 0.01 Ma (MSWD 1.32; ³⁹Ar_K 27.08 %, ⁴⁰Ar^∗ 55.78 %; inverse isochron age 2.67 ± 0.03 Ma). The isochron intercepts for both experiments are atmospheric. The combined age of 2.68 ± 0.01 Ma should be considered with caution due to the rather low amount of released ³⁹Ar (23 %–28 %).

3.1.3 Single biotite grain ${}^{\mathrm{40}}\mathrm{Ar}{/}^{\mathrm{39}}\mathrm{Ar}$ fusion and/or isochron ages

Results of nine single-fusion experiments are given in Fig. 7. Nine or ten replicate single-fusion experiments were conducted on 5–10 grains biotite per fusion. Sample G15M0006 is from dacite with columnar joints from the Kalogeros crypto-dome in the north-east of Milos (VU108-Z11, Fig. 7a). The sample shows a weighted mean age of 2.72 ± 0.01 Ma for 9 out of 10 total fusion experiments (MSWD 1.95; $\mathrm{9}/\mathrm{10}$) with an average 47.9 % of radiogenic ⁴⁰Ar. The inverse isochron age is 2.62 ± 0.04 Ma (MSWD 0.99). Note that excess argon (${}^{\mathrm{40}}\mathrm{Ar}{/}^{\mathrm{36}}\mathrm{Ar}$ 310.2 ± 4.0) is present. Hence the inverse isochron age is younger compared to the weighted mean age. The isochron age of 2.62 ± 0.04 Ma is regarded as the best estimate for the emplacement age.

Figure 7 Biotite ${}^{\mathrm{40}}\mathrm{Ar}{/}^{\mathrm{39}}\mathrm{Ar}$ total fusion ages for samples G15M0006 (a) and G15M0025-26 (b, c), G15M0022-24 (d–f), G15M0013 (g) and G15M0033-35 (h–j). Data outside the shaded area are not included in the weighted mean. Individual steps and final age calculation are reported with 1σ errors. The Kalogeros crypto-dome and Mavros Kavos lava dome are located in the north-eastern and south-western parts of the Milos VF, respectively, and the Triades lava dome, Halepa lava dome, Trachilas complex and Kalamos lava are situated in the southern, northern and south-eastern parts of the Milos VF (see Fig. 2).

Sample G15M0025 was collected from the Mavros Kavos lava dome located in the west of Milos (Fig. 2). The biotite of this sample (VU108-Z2, Fig. 7b) shows a weighted mean age of 2.36 ± 0.01 Ma (MSWD 0.70; $\mathrm{9}/\mathrm{10}$; ⁴⁰Ar^∗ 37.60 %, inverse isochron age 2.34 ± 0.04 Ma) with an ${}^{\mathrm{40}}\mathrm{Ar}{/}^{\mathrm{36}}\mathrm{Ar}$ intercept of 300.6 ± 3.5. The age of 2.36 ± 0.01 Ma is considered the best eruption age estimate for this sample.

Samples G15M0023 and G15M0024 are from the Triades lava dome north-east of Milos (Fig. 2). A mafic enclave G15M0022 (host rock G15M0021) was collected from a lava near Cape Vani (Fig. 2). The total fusion experiments of the biotites show that their initial ${}^{\mathrm{40}}\mathrm{Ar}{/}^{\mathrm{36}}\mathrm{Ar}$ estimates overlap with air (296–300). The total fusion ages gave the best estimates for their eruption ages of 2.10–2.13 Ma using 22 out of 31 fusions with a range of radiogenic ⁴⁰Ar between 30 %–36 % (Fig. 7b).

Sample G15M0013 is from the rhyolitic Halepa lava dome in the south of Milos (Fig. 2). The total fusion experiment (VU108-Z13, Fig. 7c) on biotite of this sample produced a weighted mean age of 1.04 ± 0.01 Ma (MSWD 1.62; $\mathrm{9}/\mathrm{10}$, ⁴⁰Ar^∗ 26.3 %; inverse isochron age 1.02 ± 0.04 Ma) with an initial ${}^{\mathrm{40}}\mathrm{Ar}{/}^{\mathrm{36}}\mathrm{Ar}$ estimate of 299. 8 ± 4.1. The best estimate for the eruption age of the Halepa rhyolite is 1.04 ± 0.01 Ma.

Samples G15M0034 and G15M0035 were collected from a lava dome located south-east of the Trachilas cone (Fig. 2). Nine total fusion experiments (VU108-Z21, Fig. 7c) were performed on biotite of sample G15M0035 and yielded the age of 0.63 ± 0.02 Ma (MSWD 1.26; $\mathrm{6}/\mathrm{9}$; ⁴⁰Ar^∗ 4.9 %; inverse isochron age 0.77 ± 0.13 Ma). The atmospheric isochron intercept overlaps with air at 2σ (296.4 ± 1.7). The 4.9 % of radiogenic ⁴⁰Ar is so low that we should consider the age of 0.63 ± 0.02 Ma with caution. For biotite of sample G15M0034 (VU108-Z20, Fig. 7c) one total fusion experiment produced a weighted mean age of 0.51 ± 0.02 Ma (MSWD 0.95; $\mathrm{6}/\mathrm{10}$; ⁴⁰Ar^∗ 3.5 %; inverse isochron age 0.61 ± 0.08 Ma) with an atmospheric isochron intercept. The age of 0.51 ± 0.02 Ma also needs to be regarded as possibly suspect due to the low amount of radiogenic ⁴⁰Ar.

Sample G15M0033 was collected from the Kalamos lava along the coast of the south-west of the Fyriplaka rhyolitic complex (Fig. 2). Biotite of this sample (VU108-Z19, Fig. 7c) yielded 0.412 ± 0.004 Ma (MSWD 1.10; $\mathrm{8}/\mathrm{10}$; inverse isochron age 0.39 ± 0.02 Ma) with ∼ 22.2 % of radiogenic ⁴⁰Ar, which is regarded as the eruption age for the Kalamos lava.

3.1.4 Multiple biotite grain ${}^{\mathrm{40}}\mathrm{Ar}{/}^{\mathrm{39}}\mathrm{Ar}$ incremental heating plateau and/or isochron ages

Figure 8 displays the biotite ${}^{\mathrm{40}}\mathrm{Ar}{/}^{\mathrm{39}}\mathrm{Ar}$ ages measured by the incremental heating steps method. Sample G15M0021 is the host lava of mafic enclave G15M0022. Twelve replicate total fusion experiments on its biotite (VU110-Z4, Table 3) produced an age of 2.48 ± 0.04 Ma (MSWD 1.49; $\mathrm{4}/\mathrm{12}$, ⁴⁰Ar^∗ 36.09 %; inverse isochron age 3.44 ± 0.46 Ma). Although this suggests a correct age, the large analytical error of each fusion (> 0.3 Ma on average) and poor reproducibility ($\mathrm{4}/\mathrm{12}$) of this experiment probably result in an unreliable age. Therefore, two more incremental heating experiments were performed on this sample (VU110-Z4_2 and VU110-Z4_2b, Fig. 8a), which gave an age of 1.97 ± 0.01 Ma (MSWD 1.66; ³⁹Ar_K 63.8 %, ⁴⁰Ar^∗ 54.7 %; inverse isochron age 1.97 ± 0.03 Ma) and 2.01 ± 0.01 Ma (MSWD 6.76; ³⁹Ar_K 75.39 %, ⁴⁰Ar^∗ 57.84 %; inverse isochron age 2.04 ± 0.05 Ma), respectively. The scatter in the latter is too high to define a reliable plateau age and the first incremental heating experiment is regarded as the best estimate of the eruption age of this sample.

Figure 8 Biotite ${}^{\mathrm{40}}\mathrm{Ar}{/}^{\mathrm{39}}\mathrm{Ar}$ plateau ages for samples G15M0021 (a), G15M0007 (b), and G15M0009 (VU110-Z23_combined), G15M0012 (VU110-Z24_combined) and G15M0008 (VU110-Z22_combined) (c). The numbers in red represent negative ages. Individual steps and final age calculation are reported with 1σ errors. The Triades lava dome, Trachilas and Fyriplaka complexes are located in the north-western, northern and south-eastern parts of the Milos VF, respectively (see Fig. 2). See the individual steps of samples G15M0021, G15M0007, G15M0009, G15M0012 and G15M0008 in Supplement file II.

Sample G15M0007 was collected from the rhyolitic Trachilas complex in the north of Milos (Fig. 2). Twenty-two total fusion (VU110-Z12, Table 3) and two incremental heating experiments (VU110-Z12a and 12b, Fig. 8b) were performed on biotite of this sample. The total fusion experiments did not result in a reliable age due to the large errors of single steps (± 0.19 Ma on average) and the rather low amount of radiogenic ⁴⁰Ar (9.1 %). On the other hand, the first incremental heating experiment produced a plateau age of 0.30 ± 0.01 Ma (MSWD 4.61; ³⁹Ar_K 56.60 %; inverse isochron age 0.28 ± 0.05 Ma) including 14.51 % of radiogenic ⁴⁰Ar. The second incremental heating experiment yielded a plateau of 0.317 ± 0.004 Ma (MSWD 1.29; ³⁹Ar_K 74.05 %; inverse isochron age 0.31 ± 0.03 Ma) with a higher amount of radiogenic ⁴⁰Ar (18.30 %). The isochron intercepts of both incremental heating experiments are atmospheric. The second experiment is the best estimate for the eruption age, since it contained the largest amount of radiogenic ⁴⁰Ar and has a better reproducibility of single heating steps.

Three pumice clasts (G15M0008-9 and G15M0012) were sampled from different layers of the Fyriplaka complex (Fig. 2). The first incremental-step heating experiment on biotite from sample G15M0009 (VU110-Z23a, Fig. 8c) gave negative ages at the lower-temperature heating steps. Four consecutive higher-temperature heating steps seem to define a plateau of 0.11 ± 0.02 Ma (MSWD 1.37) only using 18.33 % of the total ³⁹Ar_K with 1.65 % of radiogenic ⁴⁰Ar. The second experiment (VU110-Z23b) also yielded a plateau of 0.11 ± 0.03 Ma (MSWD 6.77) at higher-temperature heating steps including 41.05 % of the total ³⁹Ar_K and 3.13 % of radiogenic ⁴⁰Ar. The significantly larger error of the isochron age may be due to the clustering of data close to 0 on the y axis. The two experiments (VU110-Z23a and Z23b) are comparable. The combined age of 0.11 ± 0.02 (MSWD 3.5) is consistent with the age of 0.09–0.14 Ma from Fytikas et al. (1986). Although only 29.50 % of the released ³⁹Ar_K was used for this sample, this age still probably represents the eruption age of this layer in the Fyriplaka complex.

For biotite of sample G15M0012, both incremental-step heating experiments are comparable. Both of them yielded plateau ages of 0.05 ± 0.01 Ma (VU110-Z24a; MSWD 3.09; ³⁹Ar_K 38.89 %, ⁴⁰Ar^∗ 2.89 %; inverse isochron age 0.14 ± 0.03 Ma) and 0.09 ± 0.02 Ma (VU110-Z24b; MSWD 8.16; ³⁹Ar_K 48.04 %, ⁴⁰Ar^∗ 4.59 %; inverse isochron age 0.09 ± 0.05 Ma) at higher-temperature heating steps (Fig. 8c). The clustering of data points of experiment VU110-Z24a could result in the lower initial estimate of ${}^{\mathrm{40}}\mathrm{Ar}{/}^{\mathrm{36}}\mathrm{Ar}$ (285.98 ± 4.76). However, the combined age of 0.07 ± 0.01 Ma, using 43.53 % of the total ³⁹Ar_K with an atmospheric isochron intercept (295.67 ± 7.39), could be the representative age of eruption.

Biotite of sample G15M0008 did not result in a reliable plateau in the first incremental-step heating experiment (VU110-Z22a, Fig. 8c) but shows a very disturbed age spectrum. The second experiment (VU110-Z22b) yielded 0.062 ± 0.003 Ma (MSWD 0.91) using 71.81 % of the total ³⁹Ar_K with 2.69 % of radiogenic ⁴⁰Ar as the best estimate of the eruption age.

3.1.5 Multiple amphibole grain ${}^{\mathrm{40}}\mathrm{Ar}{/}^{\mathrm{39}}\mathrm{Ar}$ multi-grain incremental heating plateau and/or isochron ages

There are only two amphibole samples that yielded ${}^{\mathrm{40}}\mathrm{Ar}{/}^{\mathrm{36}}\mathrm{Ar}$ plateau and/or isochron ages (Fig. 9a and b). Sample G15M0004 was collected from the pyroclastic series of Adamas from the PSLD (Fytikas et al., 1986), to the north of Bombarda (Fig. 2). Two replicate heating experiments of G15M0004 amphibole (VU108-Z10_1 and VU108-Z10_2) were performed yielding 2.99 ± 0.11 Ma (MSWD 1.00; ³⁹Ar_K 87.31 %, ⁴⁰Ar^∗ 16.36 %; inverse isochron age 7.89 ± 2.46 Ma) and 2.86 ± 0.09 Ma (MSWD 1.50; ³⁹Ar_K 86.18 %, ⁴⁰Ar^∗ 17.58 %; inverse isochron age 0.70 ± 0.29 Ma). The variable atmospheric isochron intercept of both experiments (${}^{\mathrm{40}}\mathrm{Ar}{/}^{\mathrm{36}}\mathrm{Ar}$ 202.39 ± 48.47 and 348.91 ± 27.33) is due to the clustering of the data points. Note that also the amount of radiogenic ⁴⁰Ar is rather low (∼ 17 %). The two experiments are remarkably similar. A combined inverse isochron age of 1.95 ± 0.45 Ma (MSWD 1.17; ${}^{\mathrm{40}}\mathrm{Ar}{/}^{\mathrm{36}}\mathrm{Ar}$ 319.51 ± 14.70) is considered the best estimate, but ideally this age should be checked by other techniques.

Figure 9 Amphibole ${}^{\mathrm{40}}\mathrm{Ar}{/}^{\mathrm{39}}\mathrm{Ar}$ plateau or inverse isochron ages for samples G15M0004 (a) and G15M0026 (b). Final age calculation is reported with 1σ errors. The Adamas and Mavros Kavos lava domes are located in the northern and south-western parts of the Milos VF, respectively (see Fig. 2). See the individual steps of samples G15M0004 and G15M0026 in Supplement file II.

Sample G15M0026 is from the same location as sample G15M0025, which gives us the opportunity to compare the biotite age with the amphibole age. One total fusion experiment on biotite (VU108-Z1b) yielded a weighted mean age of 2.35 ± 0.01 Ma (MSWD 1.36; ⁴⁰Ar^∗ 38.6 %). The atmospheric isochron intercept is low (${}^{\mathrm{40}}\mathrm{Ar}{/}^{\mathrm{36}}\mathrm{Ar}$ 292.01 ± 2.92), and the inverse isochron age of 2.42 ± 0.04 Ma (MSWD 0.93) is considered the best result from the biotite. Two incremental heating experiments for amphibole (VU108-Z1b_1 and VU108-Z1b_2) gave plateau ages of 2.67–2.70 Ma, which are much higher values than the biotite inverse isochron ages (2.28–2.31 Ma). This result could be caused by the high ${}^{\mathrm{40}}\mathrm{Ar}{/}^{\mathrm{36}}\mathrm{Ar}$ isochron intercepts (> 320) with large uncertainties of ∼ 29. Therefore, on the basis of the remarkable similarity of the two experiments, the combined inverse isochron age of 2.31 ± 0.28 Ma (MSWD 0.93, ³⁹Ar_K 71.36 %, ⁴⁰Ar^∗ 34.97 %) is regarded as the best estimate from amphibole, which overlaps with the biotite age of 2.42 ± 0.03 Ma. This biotite age of 2.42 ± 0.03 Ma is considered to be the best approximation of the eruption age.

3.2 Major-element results

Major-element results are given in Table 4. The SiO₂ compositions range from 54 to 78 wt % (basaltic andesite to rhyolite, see Fig. 10a). The most felsic samples (SiO₂ > 75 wt %) belong to the Fyriplaka and Trachilas complexes. Our data overlap with those of previous studies and display a similar range in SiO₂-K₂O (Francalanci and Zellmer, 2019, and references therein). The samples of Polyegos are similar to the Fyriplaka and Trachilas complexes, whereas the older Milos samples overlap with Kimolos and Antimilos (Fytikas et al., 1986; Francalanci et al., 2007).

Table 4 Major-element composition of volcanic samples from the Milos volcanic field.

The classification of rock type for each sample is on the basis of field observation and the SiO₂ vs. K₂O plot of Bas et al. (1986). All iron expressed as Fe₂O₃T(otal). LOI: loss on ignition.

Figure 10 SiO₂ vs. K₂O (a) and AFM (b) diagrams for the Milos volcanic field with data of this study as solid circles. Published data are represented by shaded fields (Francalanci and Zelmer, 2019, and references therein). Fields for the tholeiite, calc-alkaline, high-K calc-alkaline and shoshonitic series are from Peccerillo and Taylor (1976). Vertical lines defining fields for basalt, basaltic andesite, andesite, dacite and rhyolite are from Bas et al. (1986). The solid line dividing tholeiitic and calc-alkaline fields is from Irvine and Baragar (1971).

Although some samples of Antimilos are tholeiitic, all of the Milos volcanic units belong to the calc-alkaline and medium- to high-K series (Fig. 10b). A mafic inclusion, sample G15M0022, has high K₂O (6 wt %), similar to sample G15M0021 (7.2 wt %). Both of them were collected from the Cape Vani area (Fig. 2). The SiO₂ wt % vs. our ${}^{\mathrm{40}}\mathrm{Ar}{/}^{\mathrm{39}}\mathrm{Ar}$ ages diagram (Fig. 11b) shows that there is a tendency for the volcanic units to become more felsic over time. In the diagram with ${\mathrm{K}}_{\mathrm{2}}\mathrm{O}/{\mathrm{SiO}}_{\mathrm{2}}$ vs. age there is no significant change (Fig. 11c).

Figure 11 Eruption age vs. (a) cumulative eruption volume for the volcanic deposits of Milos, (b) SiO₂wt % and (c) ${\mathrm{K}}_{\mathrm{2}}\mathrm{O}/{\mathrm{SiO}}_{\mathrm{2}}$, of the Milos volcanic units of this study and previous studies. The maximum (red line) and minimum (dashed red line) cumulative eruption volume curves were estimated from Campos Venuti and Rossi (1996) and Stewart and McPhie (2006). Q_e is the long-term volumetric volcanic output rate (see discussion). The exact volume of volcanic products between 4.1 and 3.08 Ma is not well constraint and indicated with a question mark. The major-element data of the old pumices of Filakopi volcanoes (2.66 Ma) are from Stewart (2003). The major-element data of the Plakes lava dome is from Fytikas et al. (1986). Geochemical data of the old pumices of Profitis Illias (∼ 3.08 Ma) is lacking due to the severe alteration.

3.3 Variations in eruption volume with ages

Figure 11a shows the cumulative volcanic output volume of the Milos VF over time. This diagram shows that the Milos VF can be separated into three periods: periods I (∼ 3.3–2.13 Ma) and III (1.48–0.00 Ma) are characterized by low volcanic output volumes, whereas Period II (2.13–1.48 Ma) shows a rapid increase in volcanic output volume. Periods I and II are built up in a submarine setting, whereas Period III is in a subaerial setting. The Milos VF was largely (∼ 85 % by volume) constructed in a submarine setting before ∼ 1.48 Ma (periods I and II) (Fig. 11a). During Period III (1.48 Ma–present), only a small volume (∼ 15 %) of rhyolitic magma was added from different eruption vents. See the details of periods I–III in Sect. 4.3.2.

4 Discussion

4.1 Comparison with the previous geochronological studies on the Milos VF

More than half of our ${}^{\mathrm{40}}\mathrm{Ar}{/}^{\mathrm{39}}\mathrm{Ar}$ ages derived for this study are based on high-resolution laser incremental heating method. All incremental-step heating experiments are reproducible, except for the sample G15M0017 which gave the oldest age. The total fusion experiments of this study gave an at least 5 times smaller analytical uncertainty (1 SE on average ≤ 0.01 Ma) than the previous studies using conventional K–Ar (Angelier et al., 1977; Fytikas et al., 1976, 1986; Matsuda et al., 1999) and SHRIMP U–Pb zircon methods (Stewart and McPhie, 2006). Fission track dating on obsidians of the Milos VF produced two ages (Bigazzi and Radi, 1981; Arias et al., 2006), which seem to overlap with the K–Ar and ${}^{\mathrm{40}}\mathrm{Ar}{/}^{\mathrm{39}}\mathrm{Ar}$ ages but with larger uncertainty. U–Pb zircon ages could indicate the timing of zircon formation at high temperature (> 1000 ^∘C) in magma chambers significantly prior to volcanic eruption (e.g. Flowers et al., 2005). On the other hand, the lower closure temperature of K-rich minerals (< 700 ^∘C) makes the K–Ar and ${}^{\mathrm{40}}\mathrm{Ar}{/}^{\mathrm{39}}\mathrm{Ar}$ ages better suited to determine the timing of the extrusion of volcanic products (e.g. Grove and Harrison, 1996; Cassata and Renne, 2013).

The MSWD value, as a measure of the scatter of the individual step ages, is based on the error enveloping around the data point. The decrease in error will automatically cause an increase in MSWD (e.g. York, 1968; Wendt and Carl, 1991). The MSWD values reported in this study are relatively high. In part this is caused by the fact that modern multi-collector mass spectrometers used for ${}^{\mathrm{40}}\mathrm{Ar}{/}^{\mathrm{39}}\mathrm{Ar}$ dating can measure the isotope ratios very precisely, which in turn would increase the MSWD. It will be more valuable and challenging to find a plateau or isochron age which meets the MSWD criteria (< 2.5) by modern multi-collector ${}^{\mathrm{40}}\mathrm{Ar}{/}^{\mathrm{39}}\mathrm{Ar}$ dating than by K–Ar or ${}^{\mathrm{40}}\mathrm{Ar}{/}^{\mathrm{39}}\mathrm{Ar}$ dating using a single detector instrument (e.g. Mark et al., 2009).

Potential drawbacks of the ${}^{\mathrm{40}}\mathrm{Ar}{/}^{\mathrm{39}}\mathrm{Ar}$ method are its dependence on neutron irradiation causing the production of interfering argon isotopes that need to be corrected for. The uncertainty in the ages of standards that are required to quantify the neutron flux also needs to be incorporated in the final ages as there are uncertainties related to decay constants (Supplement file II). Finally, recoil can occur during irradiation. Minerals such as biotite can be prone to recoil, yielding slightly older ages (e.g. Hora et al., 2010).

In this section, our ${}^{\mathrm{40}}\mathrm{Ar}{/}^{\mathrm{39}}\mathrm{Ar}$ results are compared with previously published geochronological data and subsequently used to refine the stratigraphy of the Milos VF. In the last part, we will discuss the temporal variations in major elements and the volumetric volcanic output rate of the Milos VF.

Figure 12 The ${}^{\mathrm{40}}\mathrm{Ar}{/}^{\mathrm{39}}\mathrm{Ar}$ ages of this study (x axis) compared to the $\mathrm{K}--\mathrm{Ar}$ ages (Angelier et al., 1977; Fytikas et al., 1986), U–Pb zircon ages (Stewart and McPhie, 2006) and fission track ages (Bigazzi and Radi, 1981; Arias et al., 2006) (y axis) for the same volcanic units. Ages which deviate from the 1:1 correlation line are discussed in Sect. 4.1.

Figure 12 compares previous published K–Ar, U–Pb zircon and fission track ages from the same volcanic units with the new ${}^{\mathrm{40}}\mathrm{Ar}{/}^{\mathrm{39}}\mathrm{Ar}$ data of this study. In general, there is a good agreement; however, 6 ages out of 23 differ significantly from previous studies and will be discussed below.

The obsidian fission track ages (Bigazzi and Radi, 1981; Arias et al., 2006) for the Dhemeneghaki volcano are 0.25 Myr younger than the K–Ar ages (1.84 Ma, Angelier et al., 1977) and the ${}^{\mathrm{40}}\mathrm{Ar}{/}^{\mathrm{39}}\mathrm{Ar}$ age of this study (1.825 Ma, G15M0032B). The good agreement between the K–Ar and ${}^{\mathrm{40}}\mathrm{Ar}{/}^{\mathrm{39}}\mathrm{Ar}$ ages suggests that the fission track ages record a different, lower-temperature event than the K–Ar and ${}^{\mathrm{40}}\mathrm{Ar}{/}^{\mathrm{39}}\mathrm{Ar}$ ages. In addition, the larger uncertainty of fission track ages (> 0.05 Ma) also overlaps with the ${}^{\mathrm{40}}\mathrm{Ar}{/}^{\mathrm{39}}\mathrm{Ar}$ age at 2σ. We assume that the ${}^{\mathrm{40}}\mathrm{Ar}{/}^{\mathrm{39}}\mathrm{Ar}$ age is the correct extrusion age for the obsidian of the Dhemeneghaki volcano.

Angelier et al. (1977) reported one dacite sample in the north-west of Milos with an age of 1.71 Ma (Angelier_3, location 3 in Fig. 3 of Angelier et al., 1977). Argon loss could result in these ages (Angelier_3–5 in Fig. 12) being younger than our ${}^{\mathrm{40}}\mathrm{Ar}{/}^{\mathrm{39}}\mathrm{Ar}$ groundmass ages of 1.97 ± 0.01 Ma (dacite sample G15M0021 and −22).

The amphibole of sample G15M0004 of the Adamas dacitic lava dome, located ∼ 1 km north of the rhyolitic Bombarda volcano, gave an inverse isochron age of 1.95 Ma ± 0.45 Ma. This age overlaps with the K–Ar age for the Adamas lava dome of 2.03 ± 0.06 Ma (dacite M 66) of Fytikas et al. (1986). The large analytical uncertainty of our sample G15M0004 is caused by a combination of low ⁴⁰Ar^∗ yields and clustering of data points that define the inverse isochron showing the excess argon was identified by the ${}^{\mathrm{40}}\mathrm{Ar}{/}^{\mathrm{39}}\mathrm{Ar}$ method (${}^{\mathrm{40}}\mathrm{Ar}{/}^{\mathrm{36}}\mathrm{Ar}$ 319.51 ± 14.70), whereas the presence of excess argon cannot be tested by the K–Ar technique, implying that the Fytikas et al. (1986) might be slightly old.

The Korakia andesite has an age of 1.59 ± 0.25 Ma (M 103, Fytikas et al., 1986) and was deposited in a submarine–subaerial environment on top of the Sarakiniko Formation, which was dated based on paleomagnetic polarity in combination with a K–Ar age (1.80–1.85 Ma, Stewart and McPhie, 2003, and references therein). The much older ${}^{\mathrm{40}}\mathrm{Ar}{/}^{\mathrm{39}}\mathrm{Ar}$ groundmass age (2.68 ± 0.01 Ma) of Korakia andesite sample G15M0029 is unreliable, and it could indicate the emplacement age of the Kalogeros crypto-dome (2.70 ± 0.04 Ma, Stewart and McPhie, 2006) or represent a geologically meaningless age with only 23 %–27 % of the total ³⁹Ar released in the plateau. In this case, the K–Ar age of 1.59 ± 0.25 Ma is regarded as the likely eruption age for the Korakia andesite, although its argon loss or excess Ar component is unknown.

Figure 13 Nine selected stratigraphic columns covering the (a) young (< 1.4 Ma) and (b) old (> 1.4 Ma) volcanic deposits of Milos modified after Stewart and McPhie (2006), except for (7) Dhemeneghaki. Age data in black are from this study and those in red are from (1) Angelier et al. (1977), (2) Fytikas et al. (1976, 1986), (3) Matsuda et al. (1999) and (4) Stewart and McPhie (2006).

We obtained ${}^{\mathrm{40}}\mathrm{Ar}{/}^{\mathrm{39}}\mathrm{Ar}$ ages of 3.41–4.10 and 3.06 ± 0.02 Ma, respectively, from the groundmasses of dacite samples G15M0017 and G15M0015 in the south-west of Milos (Figs. 2 and 13b). Both of these samples are derived from the coherent dacite facies of the rhyolitic Profitis Illias volcano based on Fig. 11 of Stewart and McPhie (2006). Sample G15M0015 yielded much higher radiogenic ⁴⁰Ar (41.77 %) than that of sample G15M0017 (< 10 % of ⁴⁰Ar^∗), and the rhyolite sample M 164 from Fytikas et al. (1986) (23.5 % of ⁴⁰Ar^∗) gave an estimate the eruptive age of 3.08 ± 0.08 Ma to the Profitis Illias volcano which is much younger than that given by our sample G15M0017 (Fig. 12). Therefore, we consider our ${}^{\mathrm{40}}\mathrm{Ar}{/}^{\mathrm{39}}\mathrm{Ar}$ ages of 3.06 ± 0.02 Ma as the best estimate of the emplacement age of the coherent dacite facies of the Profitis Illias volcano.

A basaltic andesite dyke near Kleftiko on the south-western coast of Milos has a K–Ar age of 3.50 ± 0.14 Ma, which only gave 13.9 % of ⁴⁰Ar^∗ (Fytikas et al. 1986). This age is significantly older than the eruptive ages of the Profitis Illias volcano, which the dyke intruded (Stewart, 2003). Although containing relatively low ⁴⁰Ar^∗ (16.87 %), our ${}^{\mathrm{40}}\mathrm{Ar}{/}^{\mathrm{39}}\mathrm{Ar}$ age of 2.66 ± 0.01 Ma with 67.27 % of ⁴⁰Ar^∗ from the groundmass of basaltic andesitic sample G15M0016 of the dyke near Kleftiko is probably an accurate intrusion age.

4.2 The published ages of other volcanic units

Unfortunately, we were not able to date all key volcanic units of the Milos VF. This was due to three factors: (1) we did not collect samples from all units; (2) some of the collected samples were not fresh enough after inspection of thin sections; and (3) some of the ${}^{\mathrm{40}}\mathrm{Ar}{/}^{\mathrm{39}}\mathrm{Ar}$ data indicate that the K–Ar decay system was disturbed. Therefore, we include published age information to establish a complete high-resolution geochronology for the Milos VF.

The published volcanic units that we include are the Profitis Illias volcano (3.08 ± 0.08 Ma with 23.5 ⁴⁰Ar^∗ (%), Fytikas et al., 1986), the Mavro Vouni lava dome (2.50 ± 0.09 Ma with 55.2 ⁴⁰Ar^∗ (%), Angelier et al., 1977) in the south-western part of Milos, the Bombarda volcano (1.71 ± 0.05 Ma with 24.3 ⁴⁰Ar^∗ (%), Fytikas et al., 1986) and the Plakes volcano (0.97 ± 0.06 Ma with 10.2 ⁴⁰Ar^∗ (%), Fytikas et al., 1986, and 0.8–1.2 Ma with 5.4–11.9 ⁴⁰Ar^∗ (%), Matsuda et al. 1999). Scoria deposits that Stewart and McPhie (2006) attributed to an andesitic scoria cone between Milos and Kimolos were produced in a submarine setting and maybe occasionally above sea level. No age data for this deposit have been published so far. However, the stratigraphic position of this scoria deposit is between MIL 365 (2.66 Ma, Stewart and McPhie, 2006) and M103 (1.59 Ma, Fytikas et al., 1986), which is shown in Fig. 10 of Stewart and McPhie (2006). Therefore, this scoria cone was likely active in the north-eastern part of the Milos VF between 2.6 and 1.6 Ma.

Fytikas et al. (1986) also analysed a pumice coming from the Sarakiniko deposits east of Adamas (1.85 ± 0.10 Ma with 13.6 ⁴⁰Ar^∗ (%), Fytikas et al., 1986) (Fig. 2). This unit is reworked pyroclastic sediment of the Adamas lava dome (Rinaldi and Venuti, 2003). Therefore, the K–Ar age from the Sarakiniko unit is not regarded as an eruption age in this study. We did not sample the neighbouring islands of the Milos VF and also did not attempt to date the products of the recent phase of phreatic activity from which Traineau and Dalabakis (1989) obtained ¹⁴C ages of 200 BCE and 200 CE.

4.3 Implications for the stratigraphy of the Milos VF

4.3.1 Start of volcanism in the Milos VF

Figures 13 and 14 summarize our new ${}^{\mathrm{40}}\mathrm{Ar}{/}^{\mathrm{39}}\mathrm{Ar}$ ages in combination with previously published stratigraphic, biostratigraphic, fission track, ¹⁴C, K–Ar and U–Pb age data. We did not consider the Matsuda et al. (1999) data as the fission track ages seem to be offset to other dating technique ages obtained from the same deposits (see Sect. 4.1 above). The exact start of volcanism in the Milos VF is still unclear since these older deposits are strongly hydrothermally altered. Van Hinsbergen et al. (2004) reported five ash layers in the Pliocene sedimentary rocks of southern Milos, ranging between 4.5–3.7 Ma in age, based on biostratigraphy, magnetostratigraphy and astronomical dating. In a slightly wider circle around Milos island, the 6.943 ± 0.005 Ma a1 tephra event recorded in several locations on nearby Crete (Rivera et al., 2011) shows that explosive volcanism along the Aegean arc, possibly on Milos, already occurred during the Messinian. These ash beds cannot be traced to currently exposed centres in the Milos VF and could conceivably be related to volcanic centres further north (Antiparos and Patmos), which were active during this time interval (Vougioukalakis et al., 2019).

Figure 14 Diagram presenting three periods of different long-term volumetric volcanic output rate of the Milos volcanic field based on the new ${}^{\mathrm{40}}\mathrm{Ar}{/}^{\mathrm{39}}\mathrm{Ar}$ data of this study and published data. The location of the different volcanoes is given in Fig. 2 and indicated in the left panel (from left to right: SW, W, NW, N, NE, E, SE and S of Milos). The right panel corresponds to published age data: A – Fytikas et al. (1976); B – Angelier et al. (1977); C – Fytikas et al. (1986); D – Bigazzi and Radi (1981); E – Matsuda et al. (1999); F – Stewart and McPhie (2006); G – Traineau and Dalabakis (1989). Biostratigraphic data of the Neogene sediments (NGs) are from H (Calvo et al., 2012) and I (Van Hinsbergen et al., 2004) calibrated to Raffi et al. (2020) (LCO of Sphenolithus spp. and FO of Discoaster tamalis). The number in the left panel represents the volcanic centres of Milos (see details in Table 5). The start of volcanism (3.08–3.61 Ma) on Milos and the basement of the other islands (Antimilos, Kimolos and Polyegos) are not well constrained and indicated with question marks (see text for discussion). The simplified basement cross section (NS: Neogene sedimentary rock; MB: Metamorphic basement) under Milos volcanic units is based on Fytikas et al. (1989). We used the filled symbols as the best estimate for the eruption ages at the different volcanic centres, and the open symbols are not used as the best estimate due to their relatively large uncertainties.

Biostratigraphy shows that the youngest layer with dateable fossils (bio-event, the last common occurrence of Sphenolithus spp., Van Hinsbergen et al., 2004) in the Neogene sedimentary rocks is 3.61 Myr old (GTS2020, Raffi et al., 2020). The diatomite unit II from Calvo et al. (2012) on top of the oldest volcaniclastic deposit from the north-eastern coast of Milos is constrained within 2.83–3.19 Ma. These data suggest that the oldest products must be older than 2.83 Ma and younger than 3.61 Ma. Our oldest ${}^{\mathrm{40}}\mathrm{Ar}{/}^{\mathrm{39}}\mathrm{Ar}$ ages of this study displayed a wide range of 3.41–4.10 Ma, which is probably not correct due to the alteration of the samples. Alteration might induce Ar loss and that would imply that the age is even older than 3.4–4.1 Ma. The age of 3.50 ± 0.14 Ma given by Fytikas et al. (1986) for an andesitic pillow lava or dyke has been discussed above and probably belongs to a series of basaltic andesite intrusions in the younger dacitic–rhyolitic deposits of Profitis Illias (∼ 3.08 Ma, Fytikas et al., 1986), and therefore the 3.5 Ma age is probably not correct (e.g. Stewart, 2003). Fytikas et al. (1986) measured one sample from Kimolos (Figs. 2 and 3) with an age of 3.34 Ma. Furthermore, Ferrara et al. (1980) reported an age of 3.15 Ma for a lithic clast derived from the Petalia intrusion in the Kastro volcaniclastics of Polyegos. If we assume that this reported age is a cooling age, volcanism in the Milos VF must have started before 3.15 Ma. Although age constraints for the start of volcanism on Milos both from the Neogene sedimentary rocks and the dated volcanic samples are poor, the evidence at this stage would suggest that volcanism in the Milos VF started at ∼ 3.3 Ma.

4.3.2 Periods with different volumetric output

The volume estimates of the Milos VF are hampered by limited exposure of several volcanic units and unknown age relationships. Therefore, not all units can be attributed to a certain volcano. Furthermore, we also do not know how much of the volcanic products was lost through transport by air, sea currents and erosion. Therefore, the discussion here only provides a first-order estimate of the onshore extruded magma volume. Taking into account all these limitations, our age data and the volume estimates by Stewart and McPhie (2006) indicate at least three periods of different long-term volumetric volcanic output rates (Q_e) from ∼ 3.3 to 0.0 Ma. We define a “Period” as a time interval were the Q_e is significantly different from the average output rate (Q_e average = 1.0 × 10⁻⁵ km³ yr⁻¹) of the Milos VF over the last 3.3 Myr. Figure 11 shows that the Q_e can be subdivided into two slow-growth periods (I and III) and one period (II) during which the Q_e was significantly larger.

Table 5 Summary of the eruption ages of the Milos volcanic field.

The lower boundary of Period I is based on our estimate of the oldest volcanic units of Milos at ∼ 3.3 Ma. These oldest units were deposited in the south-west of Milos between ∼ 3.3 and 3.08 Ma and include the BPS of Fytikas et al. (1986) and the felsic pumice-cone/crypto-dome facies of Stewart and McPhie (2006). These deposits have a minimum thickness of 120 m. The estimates of the DRE volume and the Q_e of these earliest volcanic deposits are hampered by the lack of precise age information, the high degree of alteration and structural complexities. Therefore, we only calculated the Q_e of Period I from 3.08 Ma, for which the eruption products are mainly dacitic–rhyolitic in composition (Table 5, Fig. 11), and the first products that can be reliably dated are crypto-domes (3.06 Ma, sample G15M0015) and dykes (2.66 Ma, sample G15M0016) into the BPS of Fytikas et al. (1986) or the units of the Profitis Illias volcano of Stewart and McPhie (2006, 3.08 Ma) in the south-west of Milos. These crypto-domes and dykes were followed by the formation of the submarine Filakopi pumice-cone volcano at 2.66 Ma (Stewart and McPhie, 2006) and the Kalogeros crypto-dome at 2.62 Ma (sample G15M0006) in the north-eastern part of Milos. These two pumice-cone volcanoes contributed 3–11 km³ DRE in volume to the Milos VF. The last two volcanic activities of Period I occurred in the south-west (Mavro Vouni, 2.50 Ma, Angelier et al., 1977) and west of Milos (Mavros Kavos, 2.36 Ma, this study), which produced two high-aspect-ratio andesitic–dacitic lava domes with a total volume of 1–3 km³ DRE (Stewart and McPhie, 2006). During the submarine Period I, which lasted ∼ 1.2 Myr, the estimated Q_e is 0.9 ± 0.5 × 10⁻⁵ km³ yr⁻¹.

The change from periods I to II is based on the sharp increase in the Q_e at 2.13 Ma (Fig. 11). During this period the Q_e (3.0 ± 1.7 × 10⁻⁵ km³ yr⁻¹) increased by a factor of ∼ 3 compared to periods I and III. Period II began with the submarine extrusions of the dacitic–rhyolitic Triades lava dome in the north-west and dacitic Adamas lava dome in the north-east of Milos and was followed by the rhyolitic Dhemeneghaki volcano and the Bombarda volcano in the north-east of Milos. For the Bombarda centre a large age range is reported in the literature (1.71–2.15 Ma, Fig. 13b). We did not successfully date samples from the Bombarda centre, but Rinaldi and Venuti (2003) reported that an age of 1.71 Ma is the best approximation based on other stratigraphic information. For the Dhemeneghaki centre, we obtained an ${}^{\mathrm{40}}\mathrm{Ar}{/}^{\mathrm{39}}\mathrm{Ar}$ age of 1.825 ± 0.002 Ma from obsidian. The Triades, Adamas, Dhemeneghaki and Bombarda centres all developed in submarine settings, as the intercalated sediments from the northern coast of Milos show (Calvo et al., 2012; Fig. 14). The last two volcanic expressions in Period II consist of two submarine-to-subaerial lava dome extrusions: Kontaro (1.59 Ma, Fytikas et al., 1986) and Korakia (1.48 Ma, this study) in the north-west and north-east of Milos, respectively. The products of these two centres are andesitic–dacitic in composition. All volcanic centres of Period II produced 8–30 km³ DRE in volume for the Milos VF.

Period III began with a time interval of 0.4 Ma with no eruptions and has a very low Q_e of 0.25 ± 0.05 × 10⁻⁵ km³ yr⁻¹. The boundary between periods II and III can be placed at the last eruption of Period II, at the start of the first eruption in the low-output interval or halfway in between. The difference between those options is not significant, given the large uncertainties of the volume estimates (Fig. 11), and therefore we have decided to start Period III directly after the last eruption of the high Q_e of Period II. The composition of nearly all Period III volcanic products is rhyolitic, an exception is the dacitic Plakes lava dome (Fig. 11). The Plakes lava dome is probably the last volcano erupting at ∼ 0.97 Ma (Fytikas et al., 1987) in a submarine environment in the north of Milos, whereas the other lava dome in Period III, Halepa, produced rhyolitic lavas in a subaerial setting in the south (Stewart and McPhie, 2006). The Halepa and Plakes domes contributed 1–3 km³ DRE in volume to the Milos VF and were followed by a 0.3 Myr interval with no or limited volcanic eruptions. Two subaerial pumice-cone volcanoes with biotite-bearing rhyolites were constructed during the last 0.6 Myr: the Trachilas and Fyriplaka complexes. The Trachilas complex was active for approximately 300 kyr (0.63–0.32 Ma) in the northern part of Milos. The evolution of this complex began with phreatic eruptions, which became less explosive over time (Fytikas et al., 1986). During the last eruption (0.317 ± 0.004 Ma) of the Trachilas complex rhyolitic pyroclastic deposits filled up the crater area and breached the northern tuff cone walls. The Trachilas complex only added a small volume (1–2 km³ DRE) to the Milos VF. The Kalamos lava dome was also extruded in the south of Milos (Fig. 2) contemporaneously with the Trachilas complex.

The youngest volcanic activity of Milos (0.11 Ma–present) is characterized by subaerial eruptions of biotite–phyric-rhyolite from the Fyriplaka complex in the south of Milos and was studied in detail by Campos Venuti and Rossi (1996). This complex is constructed on a paleosol that developed in a phreatic deposit (green lahar, Fytikas et al., 1986) or lies directly on the metamorphic basement. Campos Venuti and Rossi (1996) indicated that the stratigraphic order is Fyriplaka and Gheraki tuff rings, Fyriplaka lava flow and tuff cone of Tsigrado–Provatas. The total estimated volume of volcanic material is 0.18 km³ DRE. The boundary between the Fyriplaka and Tsigrado tuff cones is characterized by a marked erosive unconformity. The composition of these young volcanic products is very constant (Figs. 10 and 11), as noted by Fytikas et al. (1986) and Campos Venuti and Rossi (1996). The products from Fyriplaka and Tsigrado cones are covered by a paleosol rich in archaeological remains and a phreatic deposit consisting largely of greenschist metamorphic fragments. According to Campos Venuti and Rossi (1996), the Fyriplaka cone was quickly built by phreatic and phreatomagmatic eruptions, as there are no paleosols observed between the different units. However, our data do suggest a large range in ages between 0.11 and 0.06 Ma. Fytikas et al. (1986) also reported a range between 0.14 and 0.09 Ma. These ages are inconsistent with the green lahar age of 27 kyr (Principe et al., 2002), suggesting that the green lahar deposit consists of many different phreatic eruption layers that were formed during a time interval of more than 0.4 Ma, as the Kalamos lava is underlain by a green phreatic eruption breccia (Campos Venuti and Rossi 1996). We, therefore, conclude that phreatic eruptions occurred for more than 400 kyr, predominantly in the eastern part of Milos until historical times (200 BCE–200 CE, Traineau and Dalabakis, 1989).

4.3.3 Temporal evolution of the magma flux and composition

Figure 11 shows temporal major-element variations during the evolution of the Milos VF. The volcanic units of Period III are dominantly rhyolitic in composition, whereas during periods I and II the compositions of volcanic units range between basaltic andesite to rhyolite. However, the ${\mathrm{K}}_{\mathrm{2}}\mathrm{O}/{\mathrm{SiO}}_{\mathrm{2}}$ ratio is constant (0.05 ± 0.02) over the 3.3 Myr evolution of the Milos VF, with one exception: sample G15M0021 collected near Cape Vani, which is altered by hydrothermal processes (e.g. Alfieris et al. 2013). Periods I and III contain large explosive pumice-cone volcanoes, whereas Period II is dominated by effusive dome extrusions. The difference in volcanic structures is not observed in the SiO₂ content and the ${\mathrm{K}}_{\mathrm{2}}\mathrm{O}/{\mathrm{SiO}}_{\mathrm{2}}$ ratio of the volcanic products.

It is noteworthy that the value of the Q_e (0.2–4.7 × 10⁻⁵ km³ yr⁻¹) for the Milos VF is at least 2–3 orders lower than the average for rhyolitic systems (4.0 × 10⁻³ km³ yr⁻¹) and the mean for continental arcs (∼ 70 × 10⁻³ km³ yr⁻¹) (White et al., 2006). Milos overlaps with the lowest Q_e values of the study of White et al. (2006). No data are available for the ratio between intruded magma in the crust below Milos and extruded volcanic units (I:E). White et al. (2006) argued that a ratio of 5:1 (I:E) is probably a realistic estimate for most volcanic centres and that this ratio can be higher in volcanic centres constructed on continental crust. A magma supply rate from the mantle beneath the Milos VF could be estimated in the order of 0.1–3.3 × 10⁻⁴ km³ yr⁻¹. Druitt et al. (2019) reported a long-term average magma supply rate of approximately 1 × 10⁻³ km³ yr⁻¹ beneath the Kameni islands of Santorini, which is comparable to that of Milos. Besides the case of the Santorini VF, no other information on the long-term average magma supply rate of other volcanic centres of the SAVA is available to our knowledge.

Milos is approximately 15 km long (W–E), and a magma production rate of approximately 0.7–22 ${\mathrm{km}}^{\mathrm{3}}{\mathrm{km}}^{-\mathrm{1}}{\mathrm{Myr}}^{-\mathrm{1}}$ can be estimated over the last ∼ 3.3 Myr. Although this magma production rate per kilometre arc length is the onshore estimate for the Milos VF, it is still significantly lower than for oceanic arcs: 157–220 ${\mathrm{km}}^{\mathrm{3}}{\mathrm{Myr}}^{-\mathrm{1}}{\mathrm{km}}^{-\mathrm{1}}$ (Jicha and Jagoutz, 2015). For continental arcs, the long-term magma production rate is more difficult to establish because magmatism is cyclic, and short periods (5–20 Myr) of intense magmatism (“flare-ups”) with 85 ${\mathrm{km}}^{\mathrm{3}}{\mathrm{km}}^{-\mathrm{1}}{\mathrm{Myr}}^{-\mathrm{1}}$ alternate with periods of 25–50 Myr with a low magma production rate of 20 ${\mathrm{km}}^{\mathrm{3}}{\mathrm{km}}^{-\mathrm{1}}{\mathrm{Myr}}^{-\mathrm{1}}$ (e.g. Jicha and Jagoutz, 2015). The periods of low magma production overlap with the magma production rates beneath the Milos VF over the past ∼ 3.3 Myr.

5 Conclusions

This study reports 21 new ${}^{\mathrm{40}}\mathrm{Ar}{/}^{\mathrm{39}}\mathrm{Ar}$ ages and major-element data for 10 volcanic units of the Milos volcanic field.

In combination with previously published age data, geochemistry and facies analysis the following points can be made.

1. The exact age of the start of volcanism in the Milos VF is still unclear due to the high degree of alteration of the oldest deposits. The best estimate based on our new ${}^{\mathrm{40}}\mathrm{Ar}{/}^{\mathrm{39}}\mathrm{Ar}$ ages, published K–Ar data and nannofossil biozones is between 3.5 and 3.15 Ma.

2. Based on the long-term volumetric volcanic output rate, the volcanic history of the Milos VF can be divided into two slow growth periods – periods I (∼ 3.3–2.13 Ma) and III (1.48 Ma–present) – and one relatively fast growth period, Period II (2.13–1.48 Ma).

3. Periods I and II are characterized by andesitic to rhyolitic lavas and pyroclastic units, whereas those of Period III are dominantly rhyolitic. The ${\mathrm{K}}_{\mathrm{2}}\mathrm{O}/{\mathrm{SiO}}_{\mathrm{2}}$ ratio is constant over the 3.3 Myr history of the Milos VF.

4. The long-term volumetric volcanic output rate of Milos is 0.2–4.7 × 10⁻⁵ km³ yr⁻¹, 2–3 orders of magnitude lower than the average for rhyolitic systems and continental arcs.