Temporal Evolution of Proto-Izu–Bonin–Mariana Arc Volcanism over 10 Myr: Constraints from Statistical Analysis of Melt Inclusion Compositions

International Ocean Discovery Program (IODP) Expedition 351 ‘Izu–Bonin–Mariana (IBM) Arc Origins’ drilled Site U1438, situated in the northwestern region of the Philippine Sea. Here volcaniclastic sediments and the igneous basement of the proto-IBM volcanic arc were recovered. To gain a better understanding of the magmatic processes and evolution of the proto-IBM arc, we studied melt inclusions hosted in fresh igneous minerals and sampled from 30–40Myr old deposits, reﬂect-ing the maturation of arc volcanism following subduction initiation at 52Ma. We performed a novel statistical analysis on the major element composition of 237 representative melt inclusions selected from a previously published dataset, covering the full age range between 30 and 40 Ma. In addition, we analysed volatiles (H 2 O, S, F and Cl) and P 2 O 5 by secondary ion mass spectrometry for a subset of 47 melt inclusions selected from the dataset. Based on statistical analysis of the major element composition of melt inclusions and by considering their trace and volatile element compositions, we distinguished ﬁve main clusters of melt inclusions, which can be further separated into a total of eight subclusters. Among the eight subclusters, we identiﬁed three major magma types: (1) enriched medium-K magmas, which form a tholeiitic trend (30–38Ma); (2) enriched medium-K magmas, which form a calc-alkaline trend (30–39Ma); (3) depleted low-K magmas, which form a calc-alkaline trend (35–40 Ma). We demonstrate the following: (1) the eruption of depleted low-K calc-alkaline magmas occurred prior to 40 Ma and ceased sharply at 35 Ma; (2) the eruption of depleted low-K calc-alkaline magmas, enriched medium-K calc-alkaline magmas and enriched medium-K tholeiitic magmas overlapped between 35 and 38–39Ma; (3) the eruption of enriched medium-K tholeiitic and enriched medium-K calc-alkaline magmas became predominant thereafter at the proto-IBM arc. Identiﬁcation of three major magma types is distinct from the previous work, in which enriched medium-K calc-alkaline magmas and depleted low-K calc-alkaline magmas were not identiﬁed. This indicates the usefulness of our statistical analysis as a powerful tool to partition a mixture of multivariable geochemical datasets, such as the composition of melt inclusions in this case. Our data suggest that a depleted mantle source had been replaced by an enriched mantle source owing to convection beneath the proto-IBM arc from > 40 to 35 Ma. Finally, thermodynamic modelling indicates that the overall geochemical variation of melt inclusions assigned to each cluster can be broadly reproduced either by crystallization differentiation assuming P ¼ 50 MPa ( (cid:2) 2km deep) and (cid:2) 2wt% H 2 O (almost saturated H 2 O content at 50 MPa) or P ¼ 300MPa ( (cid:2) 15km deep) and (cid:2) 6wt% H 2 O (almost saturated H 2 O content at 300MPa). Assuming oxygen fugacity ( f O2 ) of log f O2 equal to þ 1 relative to the nickel–nickel oxide (NNO) buffer best reproduces the overall geochemical variation of melt inclusions, but assuming more oxidizing conditions (log f O2 ¼ þ 1 to þ 2 NNO) probably reproduces the geochemical variation of enriched medium-K and calc-alkaline melt inclusions (30–39 Ma).

The magmatic evolution of the IBM arc-basin system can be reconstructed from the geochemistry of melt inclusions hosted in fresh igneous minerals that are recovered throughout the sequence of coarse (sand to gravel) volcaniclastic sediments, because melt inclusions can be protected from surface processes after entrapment in their host minerals, if these are chemically and physically stable. Brandl et al. (2017) have analysed the major, trace and volatile elements (S and Cl) of 304 melt inclusions hosted in fresh clinopyroxene and plagioclase grains from well-dated volcaniclastic sediments of Unit III of the core ranging from 30 Ma (Rupelian; Lower Oligocene) to 40 Ma (Bartonian; Middle Eocene) (Fig. 2), and discussed the temporal evolution of the proto-IBM arc volcanism. They concluded that (1) volcanism of the proto-IBM arc shifted gradually from calc-alkaline to tholeiitic affinity with time (30-40 Ma) and (2) such a compositional shift is linked to both the volcanic productivity and the maturation of an evolving island arc. We will use the terms 'tholeiitic' and 'calc-alkaline' to refer to rock series or differentiation trends, including andesite and more silicic rocks, using the SiO 2 versus total FeO/MgO diagram of Miyashiro (1974).
In this study, we applied a statistical analysis of the published dataset of melt inclusion compositions by Brandl et al. (2017). Statistical analysis aids in partitioning a mixture of multivariable geochemical datasets and potentially helps us to reconstruct their geological evolution (e.g. Iwamori et al., 2017). The purpose of this study is to differentiate a multivariable geochemical dataset of melt inclusions into several petrogenetically distinct groups and constrain their origins. The geochemical dataset of melt inclusions by Brandl et al. (2017) can represent a 'mixture' of volcaniclastic sediments sourced from volcanic centres of the Kyushu-Palau Ridge (KPR), which is about 100 km east and upslope of Site U1438 (Fig. 1). The KPR is the remnant part of the currently active IBM arc and was active between 25 and 49 Ma (Ishizuka et al., 2011b. In addition, we analysed a subset of 47 carefully selected melt inclusions for their volatiles (H 2 O, S, F and Cl) and P 2 O 5 content using secondary ion mass spectrometry (SIMS) to extend the dataset of Brandl et al. (2017) and better constrain the temporal evolution of the proto-IBM arc volcanism from 30 to 40 Ma.
Gravity flows repeatedly transported material from arc volcanic complexes of the KPR to Site U1438 in the ASB and were probably triggered by large-volume volcanic eruptions and/or flank collapses . Whereas proto-arc basalts (48-52 Ma), boninites (44-48 Ma) and subsequent volcanic rocks representing the initial stages of arc volcanism are exposed in the forearc region of the currently active IBM arc (e.g. Ishizuka et al., 2006Ishizuka et al., , 2011aIshizuka et al., , 2018Reagan et al., 2010Reagan et al., , 2019Straub et al., 2010;Kanayama et al., 2012), few corresponding geological units have been recovered at the rear-arc side to date. Ishizuka et al. (2011b) systematically sampled submarine volcanic rocks along the KPR via dredging; however, their 40 Ar/ 39 Ar ages covered only the youngest stages of the KPR volcanism (25-28 Ma; Ishizuka et al., 2011b; see also Straub et al., 2010), when magmatic activity of the proto-IBM arc ceased through arc rifting and backarc spreading (e.g. Stern, 2004). Studies using drill cores from Site 1201, located at the rear-arc side of the KPR and about 900 km south of Site U1438, found clues to the magmatic history of the proto-IBM arc from the Eocene to the Oligocene (30-35 Ma; Savov et al., 2006). Thus, the drill cores from Site U1438 were expected to complement our knowledge of the magmatic history of the proto-IBM rear-arc between the formation of proto-arc basalts (48)(49)(50)(51)(52) and the cessation of arc volcanism along the KPR (25-28 Ma) (Straub et al., 2010;Ishizuka et al., 2011b).
The recovered 1611 m long core at Site U1438 is composed of a 1461 m thick sedimentary section and  Arculus et al., 2015aArculus et al., , 2015b. Blue is deeper seafloor and red is shallower seafloor. The location of IODP Site U1438 is marked by a star. 150 m of igneous basement rocks. The good recovery (76 %) of the sedimentary section allowed us to reconstruct the magmatic history of the proto-IBM arc in detail (Brandl et al., 2017;Johnson et al., 2017). The sedimentary section is subdivided into four units based on lithology ( Fig. 2; Arculus et al., 2015aArculus et al., , 2015b. The uppermost sedimentary Unit I (160Á3 m thick) is Holocene to late Oligocene in age and thus postdates rifting of the proto-IBM arc. It is composed of hemipelagic fine-grained sediments interbedded with discrete ash layers probably derived from explosive eruptions along the nearby Ryukyu and Kyushu arcs (e.g. Kimura et al., 2015). Unit II (139Á4 m thick) is Oligocene in age and is composed of turbidites (silt-very fine sand). Unit III (1046Á4 m thick) is Oligocene to Eocene in age and is composed of coarser-grained turbidites (medium to very coarse sand to gravel) relative to Unit II. Unit IV (99Á7 m thick) is composed of siliceous pelagic sediments interbedded with tuffaceous sand and volumetrically minor intrusive rocks that possibly reflect the Summary of the lithostratigraphy of Hole U1438 and sampling depths. mbsf, meters below seafloor; cl, clay (<2 À8 mm, dark brown); si, silt (2 À8 -2 À4 mm, brown); vfs-fs, very fine sand-fine sand (2 À4 -2 À2 mm, light brown); ms-vcs, medium sand-very coarse sand (2 À2 -2 1 mm, light gray); gr, gravel (>2 1 mm, dark gray) (after Arculus et al., 2015a;Brandl et al., 2017). Total of 237 melt inclusions from 48 core locations (blue filled circles) were analysed in this study. Ages of sedimentary units are based on the age-depth model of Brandl et al. (2017), which uses shipboard micropaleontological and paleomagnetic studies (Arculus et al., 2015a(Arculus et al., , 2015b. The age range of oceanic igneous crust (46Á8-49Á3 Ma; weighted mean 48Á7 Ma) was determined by Ishizuka et al. (2018).
slow onset of arc volcanism. The age of Unit IV is 42-49 Ma (Fig. 2), overlapping the age of boninitic volcanism at the IBM arc (44-48 Ma). The deposition age of volcaniclastic sediments is based on the age-depth model of Brandl et al. (2017), which is consistent with the U-Pb geochronology (Barth et al., 2017) and the shipboard chronostratigraphy (Arculus et al., 2015a(Arculus et al., , 2015b. Most of the sedimentary units thicken eastwards (Arculus et al., 2015b), indicating that volcaniclastic sediments are sourced from the KPR at the eastern boundary of the ASB (Fig. 1). The drilled interval of the igneous basement rock Unit 1 is 150 m thick with 29 % recovery. The geochemistry of the igneous basement rocks of the ASB (Arculus et al., 2015a;Hickey-Vargas et al., 2018;Yogodzinski et al., 2018) and their age (46Á8-49Á3 Ma; Ishizuka et al., 2018) are similar to those of proto-arc basalts of the forearc region of the IBM arc dated at 50-52 Ma (Reagan et al., 2008(Reagan et al., , 2010Ishizuka et al., 2011aIshizuka et al., , 2014a.

ANALYSIS OF VOLATILES IN MELT INCLUSIONS Samples
The 2-4 cm long half-core sections of volcaniclastic sediments were electrically fragmented using the 'Selfrag Lab' at the Japan Agency for Marine-Earth Science and Technology (JAMSTEC). Fresh silicate minerals containing glassy melt inclusions were handpicked from the fragmented samples under a binocular microscope, mounted on epoxy resin and polished until the melt inclusions were exposed on the surface. Clinopyroxene and plagioclase were the predominant minerals hosting melt inclusions. Orthopyroxene, quartz and amphibole also hosted limited melt inclusions. The shipboard X-ray diffraction analysis (Arculus et al., 2015b) indicated the presence of quartz being restricted to the deeper levels of Unit III [>1120 mbsf (m below sea floor)]. Olivine was not recovered from any samples of Unit III (Brandl et al., 2017). Representative back-scattered electron (BSE) images of melt inclusions are shown in Fig. 3a-c. Plagioclasehosted melt inclusions (<50 lm in diameter) are usually rounded (Fig. 3a). Most of the clinopyroxene-hosted melt inclusions (<100 lm in diameter) are either rounded (Fig. 3b) or slightly angular (Fig. 3c). Shrinkage bubbles are observed in some clinopyroxene-hosted melt inclusions (Fig. 3c). Typically, post-entrapment overgrowth of the host minerals inside the melt inclusions was not identifiable. All analysed melt inclusions are glassy, suggesting that they were rapidly quenched after eruption. We did not analyse altered melt inclusions or melt inclusions with observed daughter minerals. In addition to glass as melt inclusions, tiny minerals (<100 lm long) are sometimes trapped as inclusions. For example, Fig. 3d is a BSE image of apatite and Fe-Ti oxide inclusions hosted by clinopyroxene.
We checked the apparent Fe-Mg partitioning between clinopyroxene host and coexisting melt inclusions to ensure that the composition of clinopyroxene-hosted melt inclusions represents the original melt compositions before entrapment. In this study, we excluded 67 clinopyroxene-hosted melt inclusions for which K D (Fe-Mg) cpx-melt 0Á2 from the dataset of 304 melt inclusions reported by Brandl et al. (2017). We assumed that melt inclusions for which K D (Fe-Mg) cpxmelt 0Á2 are in disequilibrium with host clinopyroxene [K D (Fe-Mg) cpx-melt is 0Á28 6 0Á08; Putirka, 2008] and do not represent original melt compositions, owing to post-entrapment overgrowth of host clinopyroxene and significant compositional modification of melt inclusions. This results in a total of 237 melt inclusions used for this study (Supplementary Data Table S1; supplementary data are available for downloading at http:// www.petrology.oxfordjournals.org).

Analytical methods
We carefully selected a subset of 47 representative melt inclusions (42 clinopyroxene-hosted and five plagioclase-hosted melt inclusions) from the dataset of Brandl et al. (2017), and then analysed them for their volatiles (H 2 O, S, F and Cl) and P 2 O 5 content using SIMS. We used the Cameca IMS-1280HR at the Kochi Institute for Core Sample Research of JAMSTEC. Our samples (wafer containing one-side intersected melt inclusions) were originally mounted in epoxy, polished, and then carbon-coated for major element analyses by electron probe microanalysis (EPMA; Brandl et al., 2017). Prior to the analysis with SIMS, we removed the carbon coatings and mounted the samples in indium metal. The samples were cleaned using acetone and distilled water in an ultrasonic bath, and then were dried in a high-vacuum oven (10 À7 torr) at 90 C for several days (>48 h). After drying, the samples were coated with Au to avoid charge build-up during the SIMS analyses. Samples were then stored in the airlock chamber of the SIMS system at <10 À8 torr for >48 h to improve vacuum conditions before starting the analytical session (Shimizu et al., 2017). The analyses were carried out with an $15 lm defocused Cs þ ion beam and an ion potential of 20 kV (10 kV at the ion source and 10 kV at the sample surface). A normal-incidence electron gun was used for charge compensation of the sample surface. Secondary ions ( 16 OH À , 19 F À , 30 Si À , 31 P À , 32 S À , and 35 Cl À ) were accelerated by À10 kV and were detected by an axial electron multiplier using a magnetic peakswitching method. The mass-resolving power was set to $6000. Further details of the analytical conditions and results for a suite of basaltic reference glasses have been presented by Shimizu et al. (2017). The volatile element content of the basaltic reference glasses (vol-0B, vol-005B, vol-05A, vol-1B, vol-3A, EPR-G3, IND-G1, FJ-G2, MRN-G1, MA42, BCR-2G, BIR-1G and BHVO-1G) is in the range of 0-4Á8 wt% for H 2 O, 8-1018 ppm for F, 12-2833 ppm for Cl, 0-1372 ppm for S, and 0Á027-0Á370 wt% for P 2 O 5 . An internal reference glass EPR-G3 was mounted together with the samples in each indium mount to monitor for potential instrumental drift and check the reproducibility of the analytical results. We excluded CO 2 analyses because of the potential contamination of carbon resulting from the prior carbon coating for analyses with EPMA (e.g. Shimizu et al., 2009).
The matrix effects of SIMS on volatile element content should be considered because we used basaltic reference glasses and analysed melt inclusions ranging widely from basaltic to rhyolitic compositions. The H 2 O content of silicic melt inclusions analysed with SIMS is expected to be lower than the true H 2 O content because of the matrix effects, whereas the S, F and Cl content is not affected (e.g. Hauri et al., 2002;Druitt et al., 2016). Thus, we additionally analysed the H 2 O content of six silicic (dacitic and rhyolitic) melt inclusions using a Fourier transform infrared (FTIR) spectrometer at JAMSTEC. Details of additional analyses and results are described in the Supplementary Appendix.

Results of the volatile element analyses of melt inclusions
The results of the volatile element analyses combined with the results of the major and trace element analyses of Brandl et al. (2017) are summarized in Supplementary Data Table S1. The average values of 2r deviation of the SIMS analyses are 100 ppm for H 2 O, 2 ppm for F, 4 ppm for Cl, 14 ppm for S and 30 ppm for P 2 O 5 based on repeated analyses of EPR-G3 reference glass (Shimizu et al., 2017. Among the elements analysed with SIMS in this study, P 2 O 5 , Cl and S have already been analysed with EPMA and the results have been published by Brandl et al. (2017). We used analytical data of these elements acquired either by SIMS or EPMA in this study. We used P 2 O 5 content analysed with SIMS, because the analytical uncertainty of P 2 O 5 with SIMS (2r ¼ 30 ppm) is much smaller than that with EPMA (2r ¼ 0Á1 wt%). We used Cl content (<6000 ppm) analysed with SIMS, because the analytical uncertainty of Cl with SIMS (2r ¼ 4 ppm) is much smaller than that with EPMA (2r ¼ 200 ppm). However, we consider the analytical results with EPMA to be more reliable for high-Cl melt inclusions (>8000 ppm), because the reference glasses for SIMS analysis have a range of Cl content from 12 to 2833 ppm, thus extrapolating the calibration line of Cl toward >8000 ppm may not be reasonable. We used the S content determined with SIMS if analytical data for both methods were available, because (1) the analytical uncertainty of S with SIMS (2r ¼ 14 ppm) is much smaller than that with EPMA (2r ¼ 200 ppm) and (2) the wavelength of the SKa X-ray changes as a function of the oxidation state of S, which can result in underestimation of the S content when using EPMA. Further details on comparison of P 2 O 5 , Cl and S content analysed by both methods are described in Supplementary Appendix, Fig. S1. Regarding the H 2 O content of six silicic melt inclusions, analytical results with FTIR were higher than those with SIMS because of the matrix effect. The differences in analytical results between the two methods are systematically larger with increasing SiO 2 (Supplementary Appendix, Fig. S2). In this study, we used the H 2 O data analysed with FTIR for these six silicic melt inclusions. Further details are described in Supplementary Appendix. Analytical data of volatile elements of melt inclusions are summarized in Fig. 4, where they are plotted against the deposition age. The analyses cover a wide range of volatile element content. For example, H 2 O ranges from 1Á4 to 6Á9 wt% (Fig. 4a) and S ranges from $10 to 3000 ppm (Fig. 4b). Generally, the upper limit of volatile element content, excluding three high-Cl melt inclusions (>8000 ppm; Fig. 4d), increases with decreasing age, along with that of incompatible element oxides such as P 2 O 5 (Fig. 4e). The H 2 O content decreases from basalt to andesite, but increases again with increasing SiO 2 toward rhyolitic composition (Fig. 5a). Sulphur content monotonously decreases with increasing SiO 2 ( Fig. 4b) and shows a broad correlation with FeO t , which is total iron oxides (FeO þ Fe 2 O 3 , Fig. 5c). Fluorine content correlates with K 2 O ( 1Á5 wt%; Fig. 5d). Chlorine also correlates with K 2 O, except for high-Cl melt inclusions (>8000 ppm; Fig. 5e). Fluorine and phosphorus have a similar incompatibility during melting and crystallization; therefore, the F/P ratio of melt inclusions can reflect that of the mantle source (Saal et al., 2002). The F/P ratio in elemental weight is relatively constant (0Á3-0Á8) at SiO 2 60 wt% and increases to as much as 4Á4 at SiO 2 ! 70 wt% at younger ages (Figs 4f and 5f).

STATISTICAL ANALYSIS OF MELT INCLUSION COMPOSITIONS Methods
Volcaniclastic sediments can be sourced from multiple volcanic centres in the upslope vicinity of Site U1438. Thus, melt inclusion compositions may represent a mixture of volcaniclastic sediments sourced from multiple vent sites, possibly not only from the KPR but, depending on the plate configuration, also from the frontal arc section of the proto-IBM arc. In addition, they could include volcaniclastic material from primary eruptions and reworked sediments.
The statistical analysis of multivariable datasets of such a mixture of materials is thus a useful approach and may help to reconstruct their geological evolution. Statistical approaches commonly used in petrology and geochemistry include principal component analysis (PCA) (e.g. Zindler et al., 1982;Allè gre et al., 1987;Hart et al., 1992;Stracke, 2012;, factor analysis (White & Duncan, 1996), independent component analysis (Iwamori & Albarè de, 2008;Iwamori & Nakamura, 2012Yasukawa et al., 2016), K-means cluster analysis (KCA) (e.g. Temple et al., 2008;Brandmeier & Wö rner, 2016), and implementation of machine learning techniques (e.g. Kuwatani et al., 2014;Petrelli & Perugini, 2016;Ueki et al., 2018). Iwamori et al. (2017)  the structure of a multivariable dataset is best resolved when (1) PCA is applied to the dataset and (2) the eigenvalue normalized-PC scores (which are called 'whitened data') are partitioned by KCA. PCA is commonly used for effectively specifying the uncorrelated base vectors that account for the variance. KCA is also a common statistical method used to partition the multivariate dataset into an assigned number of clusters K, in which the total distance between the mean of a cluster and the individual data points in the cluster is minimized (e.g. MacQueen, 1967, and references therein).
In this study, we applied PCA and KCA to the major element compositions (10 elements) of 237 melt inclusions, following the procedures of Iwamori et al. (2017). We did not include S and Cl content into our dataset to perform statistical analysis because (1) some melt inclusions are devoid of S and Cl content (Table S2 of Brandl et al., 2017) and (2) volatile elements are subject to degassing. First, the raw geochemical data of the 237  (Jochum et al., 2000) with EPMA for in situ monitoring. Analytical uncertainties of volatiles are smaller than symbol size. melt inclusions, on which is imposed a constant-sum (normalized to 100 wt%) constraint, were transformed into logarithmic space by performing centred log-ratio transformation (e.g. Aitchison, 1982Aitchison, , 1984Aitchison, , 1986. Then, the data were processed by primary standardization using the mean and standard deviation. Next, 'whitening' (division of PC scores by the square root of the eigenvalues) was applied to the dataset to decorrelate the variables. Based on the eigenvalues and eigenvectors, a minimum number of variables that account for a sufficiently large proportion (i.e. !90 %) of the sample variance was chosen. The number of eigenvectors that individually account for !5 % of the variance for the dataset of 237 melt inclusions is five (Fig. 6a). In addition, the cumulative contribution of the first five eigenvectors accounts for $90 % of the variance (Fig. 6b), which suggests that the number of variables can be reduced from 10 (number of elements) to five or slightly more eigenvectors . Then, we performed KCA on the pre-processed data from 100 random and different initial conditions by varying K from four to 12 and varying the number of independent components (nic) from three to eight. We finally found that the geochemical dataset of melt inclusions consists of essentially five clusters with five independent components (K ¼ 5 and nic ¼ 5). Results of statistical analysis of the 237 melt inclusions by KCA, including their principal components and independent components, are given in Supplementary Data Table S2.

Results of statistical analyses
The results of KCA are summarized in oxide variation diagrams (Fig. 7). Further details of the results are summarized in Supplementary Appendix, Figs S3 and S4. Downhole distributions of melt inclusions assigned to each cluster are summarized in Fig. 8. Cluster 1 melt inclusions (n ¼ 84) are medium-K mafic melts that form a tholeiitic differentiation trend ( Fig. 7e and g) and occur from 38 to 30 Ma (Fig. 8a). They are characterized by higher TiO 2 , suggesting that they represent melts from a fertile mantle source. Cluster 2 melt inclusions (n ¼ 61) form a calc-alkaline differentiation trend (Fig. 7d) and can be separated further into two subclusters in terms of K 2 O (Fig. 7g): Cluster 2a melt inclusions (n ¼ 2) are high-K melts and Cluster 2b melt inclusions (n ¼ 59) are medium-K melts. Cluster 2a melt inclusions only occur at $30 Ma, whereas Cluster 2b melt inclusions occur at the full range from 30 to 40 Ma (Fig. 8). Cluster 3 melt inclusions (n ¼ 67) are low-K melts that form a calcalkaline differentiation trend ( Fig. 7d and g). They can be further separated into two subclusters: Cluster 3a melt inclusions (n ¼ 2) are characterized by higher Al 2 O 3 (Fig. 7b) and lower MgO (Fig. 7e) at given SiO 2 than Cluster 3b melt inclusions (n ¼ 65). Clusters 3a and 3b can also be distinguished in terms of their ages; the former occur at $30 Ma, whereas the latter occur from 35 to 40 Ma (Fig. 8). Cluster 4 melt inclusions (n ¼ 22) are dacitic and rhyolitic melts that form a calc-alkaline trend. They can be separated further into two subclusters: the Cluster 4a melt inclusion (n ¼ 1) is distinguished in terms of lower Al 2 O 3 (Fig. 7b), higher FeO t (Fig. 7c) and lower K 2 O (Fig. 7g) than Cluster 4b melt inclusions (n ¼ 21) at given SiO 2 . Clusters 4a and 4b can also be distinguished in terms of their ages; the former occur at $37 Ma, whereas the latter occur at $30 Ma (Fig. 8). Cluster 5 melt inclusions (n ¼ 3) are characterized by low P 2 O 5 ($0 wt%; Supplementary Appendix  Table S1).
Twenty-five whole-rock core samples contain melt inclusions assigned to more than one cluster (Fig. 8a), suggesting that the whole-rock core samples are in fact mixtures of volcaniclastic sediments from different volcanic series. Such volcaniclastic sediments might have been mixed during long-distance transport by gravity flow and deposited around Site U1438. The 237 melt inclusions included in this study are hosted in 145 individual host minerals. In the case where two or more melt inclusions are hosted in a single host mineral, these melt inclusions usually fall into an identical cluster. However, we found 11 exceptional clinopyroxene minerals hosting two or more melt inclusions that are assigned to different clusters ( Fig. 8b and Table 1). Most of these cases can be explained by heterogeneity of the host minerals. For example, three host clinopyroxene minerals (D19R3B-min1, D55R3A-min2 and D31R6-min5) are zoned, and the melt inclusions located in the different zones are assigned to different clusters ( Fig. 9a-c). However, the case of the clinopyroxene D55R3B-min3 (Fig. 9d) (Jochum et al., 2000) with EPMA. The tholeiitic-calc-alkaline dividing line in (d) is from Miyashiro (1974), where FeO t is total iron oxides (FeO þ Fe 2 O 3 ). The definition of low-K, medium-K and high-K series in (g) is after Gill (1981). discussed below. In another case, two melt inclusions, which are assigned to different clusters, are hosted in a homogeneous clinopyroxene mineral (D27R3B-min4; Fig. 9e). This shows that heterogeneity of the host minerals is not the sole explanation for melt heterogeneity.
The temporal evolution of the melt as reflected by SiO 2 , TiO 2 , Al 2 O 3 , FeO t , MgO, CaO and K 2 O content is shown in Fig. 10 and interpreted in Fig. 11. It is important to note that the ages of the first appearance of Clusters 1 (38 Ma), 2b (39 Ma) and 3b (40 Ma) are significantly different, although the data are sparse for this time period (38-40 Ma).

Characteristics of volatile elements in each cluster
The results of the volatile element analyses of melt inclusions, combined with the results of KCA, are presented in Fig. 12. Analytical data of volatile elements are available for Cluster 1, 2b, 3b and 4b melt inclusions. We distinguished two coeval subgroups in Cluster 3b melt inclusions: a low-H 2 O subgroup and a high-H 2 O subgroup (Figs 11 and 12). Both subgroups represent magmas that erupted simultaneously in the period from 40 to 35 Ma. Some Cluster 1 and 3b (high-H 2 O subgroup) melt inclusions (such as D20R4B-min2-mi1 and D27R3B-min4-mi2) contain 6-7 wt% H 2 O, which is the maximum range of H 2 O content of melt inclusions reported from island arcs (e.g. Plank et al., 2013;2-6 wt%). Generally, the H 2 O content of melt inclusions decreases from 50 to 60 wt% SiO 2 but increases again from 60 to 70 wt% SiO 2 (Fig. 12a). In terms of sulphur content, systematic differences are observed among the four clusters. Some Cluster 1 melt inclusions enriched in H 2 O are also enriched in sulphur, and the maximum sulphur content of Cluster 1, 2b, 3b and 4 melt inclusions is $3000 ppm (D27R3B-min4-mi2), $760 ppm, $1200 ppm and $260 ppm, respectively ( Fig. 12b and c). The sulphur content decreases monotonically with increasing SiO 2 (Fig. 12b).
The F content linearly correlates with K 2 O ( 1Á5 wt%), except for the Cluster 4b melt inclusions when K 2 O ! 1Á5 wt% (Fig. 12d). The Cl content also positively correlates with K 2 O, except for three high-Cl (up to $10 000 ppm) melt inclusions (D55R3B-min3-mi2, D55R3B-min3-mi3 and D55R3B-min3-mi5), all of which are hosted in a single clinopyroxene D55R3b-min3 (Figs  9d and 12e). The F and Cl content of Cluster 3b melt inclusions is lower than that of the other clusters ( Fig. 12d and e). The deviation of the F/P ratio is relatively small (0Á3-0Á8) at SiO 2 60 wt% (Clusters 1, 2b and  Table 1 and Supplementary Data Tables  S1 and S2. 3b), which reflects that of the mantle source (Saal et al., 2002). The differences in the F/P ratio (0Á3-0Á8) are due to the differences in P 2 O 5 content of the melt inclusions; higher F/P ratios ($0Á8) correspond to low P 2 O 5 (around 0Á1 wt%) and lower F/P ratios ($0Á3) correspond to high P 2 O 5 (0Á2-0Á3 wt%) (Fig. 12f). The F/P ratio significantly increases from one to five at SiO 2 ! 70 wt% (Cluster 4b; Fig. 12f).  Table S1). Trace element patterns of 34 selected melt inclusions, assigned to Clusters 1, 2b, 3b, 4a and 4b, are normalized to the depleted mid-ocean ridge basalt (MORB) mantle (DMM;

Characteristics of trace elements in each cluster
Workman & Hart, 2005) and presented in Fig. 13a et al., 2015) from the Ogasawara (Bonin) islands, an uplifted segment of the proto-IBM forearc, are also shown for comparison ( Fig. 13e-g). It is important to note, however, that the latter suites from the frontal proto-IBM arc volcanoes (40-48 Ma; Kanayama et al., 2014;Umino et al., 2015) may not be directly comparable with our rear-arc melt inclusions (30-40 Ma). The trace element patterns of Cluster 1 and 2b melt inclusions ( Fig. 13a and b) overlap as a whole, with a  Table 1 and Supplementary Data Tables S1 and  S2. (See main text for a detailed description.) Host mineral: 52 53 Melt inclusion: Host mineral: 56 60 Melt inclusion:

High-K melts
Calc-alkaline trend

Cluster 3a
Low-K and high-Al melts Calc-alkaline trend n = 65

Cluster 3b
Low-K melts Calc-alkaline trend

Cluster 4b
Low-to high-K silicic melts

Cluster 5
Low-P melts n = 3 low-H2O subgroup high-H2O subgroup Fig. 11. Schematic diagram showing an overview of the occurrence of each cluster with age (in Ma).

Cluster 4b melt inclusions
Low   Brandl et al. (2017). (e-f) DMM-normalized trace earth element patterns for volcanic rocks from the proto-IBM arc. Data of trace element compositions are from Kanayama et al. (2014) and Umino et al. (2015). Normalizing values for DMM are from Workman & Hart (2005). respectively ( Fig. 13d and e). The trace element patterns of Cluster 3b and 4a melt inclusions are almost identical (Fig. 13c). Combined with major element composition, Cluster 4a lies on the trend of Cluster 3b in SiO 2 -K 2 O space (Fig. 7g). This observation indicates that melt inclusions assigned to Clusters 3b and 4a can be derived from the same mantle source. Cluster 3b melt inclusions (35-40 Ma; Fig. 13c) are not directly comparable with boninites (44-48 Ma; Fig. 13g). Silicic melt inclusions assigned to Cluster 4b show higher trace element abundances with TiO 2 depletion (Fig. 13d).

Characteristics of clusters from a volatiles perspective
The primary melts of Cluster 1 melt inclusions can be enriched in H 2 O (!6-7 wt%) and S (!2000-3000 ppm) inferred from high content of H 2 O (6-7 wt%) and S (2000-3000 ppm) in some Cluster 1 basaltic melt inclusions ( Fig. 12a and b). Cluster 2b melt inclusions contain lower H 2 O (2-3 wt%) and lower S (<1000 ppm) than Cluster 1 melt inclusions ( Fig. 12a and b), suggesting that the H 2 O and S contents of the primary melts corresponding to Cluster 2 may be lower than those of Cluster 1. Experimental studies have shown that H 2 Orich and H 2 O-poor primitive melts control magmatic differentiation along a calc-alkaline or tholeiitic differentiation trend (e.g. Hamada & Fujii, 2008;Zimmer et al., 2010). However, these experimental constraints are in contrast to our observations from Cluster 1 (tholeiitic affinity and possibly H 2 Orich primary melt) and Cluster 2b (calc-alkaline affinity and possibly H 2 O-poor primary melt) melt inclusions. We infer that Cluster 1 melt inclusions form a tholeiitic differentiation trend because H 2 O is effectively lost from the melt during differentiation (Fig. 12a). Cluster 2b melt inclusions form a calc-alkaline differentiation trend, although they have a lower H 2 O content than Cluster 1 melt inclusions. Thus, other factors than the differentiation of a H 2 O-rich primitive melts may explain the origin of calc-alkaline Cluster 2 melt inclusions, which may include slab melting and/or crustal assimilation (e.g. Francis et al., 1980;Yogodzinski et al., 1995;Benito et al., 1999). The F and Cl contents of Cluster 3b melt inclusions (both high-H 2 O and low-H 2 O subgroups) are lower than those in Cluster 1 and 2b melt inclusions (Fig. 12d and  e). The lower F and Cl contents are consistent with our inference that the mantle source of Cluster 3b melt inclusions is more depleted in incompatible elements. Partial melting of such a depleted and possibly refractory mantle source may be possible through the addition of slab-derived fluids (e.g. Pearce et al., 1992b;Morishita et al., 2011).
Cluster 4b melt inclusions are characterized by the highest F content reported in our study (Fig. 12d). The F content linearly correlates with K 2 O ( 1Á5 wt%) but deviates from this correlation when K 2 O ! 1Á5 wt%. In addition, the F/P ratio of Cluster 4b melt inclusions increases at SiO 2 ! 70 wt% (Fig. 12f). The deviation of F and subsequent increase in F/P may result from the crystallization of F-bearing apatite (e.g. Green & Watson, 1982) and is consistent with our petrographic observations. Indeed, apatite occurs as mineral inclusions in clinopyroxene in silicic core samples (e.g. U1438E-50R, $40 Ma; Fig. 3d).

Characteristics of clusters from an igneous petrology perspective
In addition to geochemical constraints, the conditions under which each cluster was formed can also be evaluated from the perspective of igneous petrology. Pearce element ratio plots can be used to interpret basaltic suites that experienced various degrees of magmatic differentiation (e.g. Pearce, 1968;Ernst et al., 1988;Russell & Nicholls, 1988). Among the proposed Pearce element ratio plots, an Al/K versus (2Ca þ Na)/K plot was applied to the mafic melt inclusions assigned to Clusters 1, 2b and 3b (Fig. 15a). Cluster 1, 2b and 3b (low-H 2 O subgroup) melt inclusions show an identical trend, implying that the geochemical variation of these melt inclusions can be explained by the crystallization of clinopyroxene and plagioclase-a result not surprising given that our melt inclusions are hosted in these minerals. Melt inclusions assigned to Cluster 3b (high-H 2 O subgroup) are offset from the main trend (Fig. 15a). Clinopyroxene and plagioclase crystallize synchronously in the low-H 2 O melts (Fig. 15b), resulting in approximately constant Al 2 O 3 with increasing SiO 2 . In the case of the high-H 2 O melts (Fig. 15b), clinopyroxene crystallizes first (increasing Al 2 O 3 with increasing SiO 2 ), to be followed later (at $50 wt% SiO 2 ) by the crystallization of plagioclase (decreasing Al 2 O 3 with increasing SiO 2 ), because the onset of plagioclase crystallization is suppressed under H 2 O-rich conditions (Fig. 15a and c). Brandl et al. (2017) analysed melt inclusions hosted in clinopyroxene and plagioclase and did not find olivine in any of the core samples from Unit III. However, this does not necessarily mean that olivine did not crystallize from the corresponding parental melts at all, because olivine is easily altered by seawater and/or hydrothermal fluids (e.g. Pokrovsky & Schott, 2000;Ueda et al., 2017) and therefore may not be preserved in volcaniclastic sediments deposited at the seafloor. We assume that olivine crystallized from mafic melts, in addition to clinopyroxene and plagioclase. This assumption is necessary to apply projection of melts saturated with olivine þ clinopyroxene þ plagioclase in the basalt tetrahedron as a geobarometer (e.g. Walker et al., 1979;Grove & Bryan, 1983;Grove & Baker, 1984;Baker & Eggler, 1987). The normative composition of melt inclusions was calculated following the procedures of Tormey et al. (1987) and Grove (1993) and projected from the plagioclase apex, with olivine þ clinopyroxene þ plagioclase cotectics at 0Á1 MPa (Walker et al., 1979) and 200 MPa (Berndt et al., 2005) and with orthopyroxene þ clinopyroxene þ plagioclase cotectics at 400 MPa (Hamada & Fujii, 2008) as pressure references (Fig. 16). A caveat is that a post-entrapment overgrowth of host clinopyroxene could slightly shift the projected position of melt inclusions off the clinopyroxene (CPX) apex, resulting in apparently higher pressure. Cluster 1 and 2b melt inclusions are projected between pressure references of 0Á1 and 200 MPa ( Fig. 16a and b). Two levels are identified regarding Cluster 3b melt inclusions: lower pressures (0Á1-200 MPa) for the low-H 2 O subgroup and higher pressures (200-400 MPa) for the high-H 2 O subgroup (Fig. 16c). Considering the experimental constraint that clinopyroxene crystallizes earlier than plagioclase from hydrous basaltic melt at pressure conditions corresponding to middle to lower crust (!300 MPa; Hamada & Fujii, 2008), it is reasonable to attribute the distinction between the low-and high-H 2 O subgroups of Cluster 3b to a crystallization differentiation under lower (0Á1-200 MPa) and higher ($400 MPa) pressures, respectively. We argue that (1) a new magma chamber was formed at shallower levels as the magma-plumbing system evolved (e.g. Ushioda et al., 2018) and therefore several magma chambers at different depths in the crust may have been present, (2) the melt was saturated with H 2 O as most arc melts are (e.g. Plank et al., 2013;2-6 wt%), and (3) the saturated-H 2 O content of Cluster 3b melts decreased with decreasing pressures (e.g. Hamada et al., 2011Hamada et al., , 2014. There seems to be no systematic difference in the pressure conditions under which crystallization differentiation proceeded (0Á1-200 MPa) among Cluster 1, 2b and 3b (low-H 2 O subgroup) melt inclusions (Fig. 16).

Origin of each cluster constrained by thermodynamic modelling
Clusters 1, 2b, 3b (low-H 2 O subgroup) and 3b (high-H 2 O subgroup) are composed of larger numbers of melt inclusions when compared with other clusters, and range in composition from basalt to andesite as a result of crystallization differentiation. Here, using thermodynamic modelling, we test whether geochemical variation of these four clusters can be explained solely by crystallization differentiation. Among the thermodynamic models designated to simulate crystallization differentiation, we use COMAGMAT 3.72 (e.g. Ariskin et al., 1993;Ariskin 1999;Ariskin & Barmina, 2004), because it simulates crystallization differentiation of hydrous arc magmas more reliably than MELTS (e.g. Almeev et al., 2004;Hamada, 2006;Kimura & Ariskin, 2014). COMAGMAT 3.72 is based on the two-lattice melt model, a kind of sub-ideal solution model after Bottinga & Weill (1972), Drake (1976) and Nielsen & Drake (1979), combined with a series of experimentally determined, mineral-melt geothermometers with empirical terms to compensate for the non-ideality of silicate melts. When simulating crystallization differentiation using COMAGMAT 3.72, the starting conditions must be given a priori, and we performed forward simulation of 'fractional crystallization' starting from the most undifferentiated melt in each cluster. After performing multiple simulations, we found that geochemical variation of Cluster 1 melt inclusions can be reproduced from their most undifferentiated melt with 2 wt% H 2 O (U1438E-27R5W56-I11) at a pressure of 50 MPa (Fig. 17). This pressure condition is in agreement with the estimated pressure range (0Á1-200 MPa; Fig. 16a), and the H 2 O content agrees with the analytical results of H 2 O measurements in most of our melt inclusions (2-3 wt%) (Fig. 12a). Variation in Al 2 O 3 (Fig. 17b), which is not reproduced perfectly, can be explained by changes in initial H 2 O content (0-2 wt%). Oxygen fugacity (f O2 ) is the main control on the stability field of magnetite, the crystallization of which affects the TiO 2 and FeO t content of the evolving melts. Thus, we performed our simulations at variable oxygen fugacities of log f O2 equal to 0, þ1 and þ2 relative to the nickel-nickel oxide (NNO) buffer and found that log f O2 ¼ þ1 NNO best reproduces the geochemical variation of Cluster 1 melt inclusions (Fig. 17).
In a similar manner, the geochemical variation of Cluster 2b melt inclusions can be explained by crystallization differentiation at P ¼ 50 MPa from their most undifferentiated melt with 2Á0 wt% H 2 O (U1438D-60R4-I1; Fig. 18). These pressure and H 2 O conditions are in agreement with the estimated pressure range (0Á1-200 MPa; Fig. 16b) and analytical results of H 2 O in melt inclusions (2-3 wt%; Fig. 12a). Variation in Al 2 O 3 (Fig. 18b) can be reproduced by changing the initial H 2 O content (0-2 wt%). Trends in K 2 O are not perfectly reproduced (Fig. 18g), probably because the initial K 2 O content was set too high. A more oxidizing condition (log f O2 ¼ þ1 to þ2 NNO) than that of Cluster 1 probably reproduces geochemical variation of Cluster 2b melt inclusions (Fig. 18). We previously suggested that the calc-alkaline differentiation trend of Cluster 2b may be explained by slab melting and/or crustal assimilation. Here, however, we demonstrate that the tholeiitic differentiation trend of Cluster 1 and the calc-alkaline differentiation trend of Cluster 2b may result from relatively lower and higher f O2 conditions, respectively (e.g. Hamada & Fujii, 2008).
The origin of Cluster 3b should be discussed with respect to the H 2 O content of the individual subgroups. The low-H 2 O subgroup can be explained by crystallization differentiation at P ¼ 50 MPa from their most undifferentiated melt with 2Á0 wt% H 2 O (U1438D-63R1W-I6; Fig. 19) at log f O2 ¼ þ1 NNO. The Cluster 4a melt inclusion, which extends from the geochemical range of Cluster 3b melt inclusions towards higher SiO 2 content, can be reproduced by crystallization differentiation from the undifferentiated Cluster 3b melt inclusions (low-H 2 O subgroup) (Fig. 19). This argument is consistent with similar trace element patterns observed in Cluster 3b and 4a melt inclusions (Fig. 13c). Trends in TiO 2 (Fig. 19a) and K 2 O (Fig. 19g)     (2) such a compositional shift is linked to both the volcanic productivity, as expressed by deposition rates of volcanic sediments, and the maturation of an evolving island arc. As an extended study of Brandl et al. (2017), this study separated the geochemical dataset of melt inclusions into five statistically robust clusters, which can be further separated into a total of eight subclusters, and their origins were discussed. We separated the 'calc-alkaline affinity' of Brandl et al. (2017) into two distinct clusters: a cluster of medium-K calc-alkaline melt inclusions (Cluster 2b) and another cluster of low-K calc-alkaline melt inclusions (Cluster 3b) representing geochemically enriched and depleted mantle sources, respectively. We have demonstrated the following: (1) the eruption of depleted low-K calc-alkaline magmas (Cluster 3b) occurred prior to 40 Ma and ceased sharply at  (2) eruption of depleted low-K calc-alkaline magmas (Cluster 3b), enriched medium-K calc-alkaline magmas (Cluster 2b) and enriched medium-K tholeiitic magmas (Cluster 1) overlapped temporally between 35 and 38-39 Ma; (3) the eruption of enriched medium-K tholeiitic and enriched depleted medium-K calc-alkaline magmas became predominant thereafter (Figs 10 and 11). Our findings thus present a different, more detailed overview of the temporal evolution of the proto-IBM arc volcanism when compared with those of Brandl et al. (2017).

New insights based on Site U1438 melt inclusion study
Of particular interest is the question regarding the nature of the mantle source that generated depleted low-K calc-alkaline magmas (Cluster 3b) in the time period of 35-40 Ma. Previous modelling of rare earth elements has shown that both tholeiitic and calc-alkaline arc basalts of the Ogasawara (Bonin) Islands (38-45 Ma) were generated by partial melting of fertile lherzolitic mantle more enriched than DMM ( Fig. 13e-g; Kanayama et al., 2014;Umino et al., 2015). This argument suggests that mantle convection was  established soon after the beginning of the Pacific slab subduction at 52 Ma (e.g. Arculus et al., 2015a;Reagan et al., 2017) and replaced the depleted residual mantle after extraction of proto-arc basalts by more fertile mantle through convection by 45 Ma. Such mantle convection was induced by the drag force of the subducting slab lithosphere and was accompanied by mantle upwelling as counter-flow from the deeper rear-arc regions (e.g. Iwamori, 1998). Mantle convection has been an essential driving force to activate the arc volcanism of the IBM arc and the geochemical evolution of the source mantle (e.g. Straub et al., 2010). The volcanism of the proto-IBM arc shifted from the eruption of boninites on Chichijima-Mukojima-Guam (44-48 Ma) to the eruption of arc magmas on Hahajima Island and the western scarp of the Bonin Ridge (38-45 Ma; Ishizuka et al., 2006Ishizuka et al., , 2011aKanayama et al., 2014;Umino et al., 2015). Arc volcanism along the KPR followed (25-49 Ma; Ishizuka et al., 2011bIshizuka et al., , 2018, on which we have focused in this study. Cluster 3b melt inclusions (35-40 Ma; Fig. 13c) are less depleted than boninites (Fig. 13g), and we thus infer that the mantle source of Cluster 3b melt inclusions may reflect the residue of proto-arc basalts at the rear-arc side (Reagan et al., 2017). High H 2 O content (6-7 wt%) in some Cluster 3b melt inclusions (high-H 2 O subgroup) is consistent with the idea of partial melting of a depleted and possibly refractory mantle source through the addition of slab-derived fluids to generate primary Cluster 3b melts. Numerical simulation suggests that fluids derived from a subducted slab hydrate a limited area of overriding lithosphere at the shallower level when subduction initiates, and that dehydration reaches a steady state where the hydration of serpentine becomes predominant at deeper levels thereafter (Arcay et al., 2005). The high H 2 O content in some Cluster 3b melt inclusions may be explained by a focused release of slab-derived fluids following subduction initiation at 52 Ma. We infer that the rocks with enriched medium-K tholeiitic (Cluster 1) and enriched medium-K calc-alkaline affinities (Cluster 2b) represent differentiated rock series from primary melts generated by partial melting of replenished, more enriched mantle sources. We also infer that slight differences in oxygen fugacities of such primary melts resulted in the formation of coeval tholeiitic and calc-alkaline affinities, as discussed in this study based on thermodynamic modelling. We have demonstrated the following: (1) the eruption of depleted low-K calc-alkaline magmas (Cluster 3b) occurred from >40 Ma and ceased sharply at 35 Ma; (2) the eruption of depleted low-K calc-alkaline magmas (Cluster 3b), enriched medium-K calc-alkaline magmas (Cluster 2b) and enriched medium-K tholeiitic magmas (Cluster 1) overlapped temporally between 35 and 38-39 Ma; (3) the eruption of enriched medium-K tholeiitic magmas (Cluster 1) and enriched medium-K calc-alkaline magmas (Cluster 2b) became predominant thereafter. Such temporal evolution of the proto-IBM arc volcanism reflects a replenishment of enriched mantle into depleted mantle through convection. Identification of such distinct magma types becomes possible by (1) using drilled core samples from the deep sea of the rear-arc side, because no corresponding geological record of the proto-IBM rear-arc volcanism has been recovered to date, and (2) applying statistical analysis on multivariable major element composition of melt inclusions and interpreting the findings combined with their trace and volatile element compositions. Thermodynamic modelling indicates that the overall geochemical variation of melt inclusions assigned to each cluster can be broadly reproduced by crystallization differentiation assuming P ¼ Assuming an oxygen fugacity (f O2 ) of log f O2 equal to þ1 relative to the nickel-nickel oxide (NNO) buffer best reproduces the geochemical variation of melt inclusions assigned to Clusters 1 and 2b, but assuming a more oxidizing condition (log f O2 ¼ þ1 to þ2 NNO) probably reproduces geochemical variation of Cluster 3b melt inclusions (low-H 2 O subgroup). We infer that slight differences in oxygen fugacities resulted in the formation of coeval tholeiitic and calc-alkaline affinities.