Investigating the Lid Effect on the Generation of Ocean Island Basalts: 1. Geochemical Trends

Ocean island basalts (OIBs) are generated by mantle plumes, with their geochemistry controlled by a combination of source composition, temperature, and thickness of overlying lithosphere. For example, OIBs erupting onto thicker, older oceanic lithosphere are expected to exhibit signatures indicative of higher average melting pressures. Here, we quantitatively investigate this relationship using a global data set of Neogene and younger OIB compositions. Local lithospheric thicknesses are estimated using theoretical plate‐cooling models and Bayes factors are applied to identify trends. Our findings provide compelling evidence for a correlation between OIB geochemistry and lithospheric thickness, with some variables (SiO2, Al2O3, FeO, Lu) showing linear trends that can be attributed to increasing average melting pressure, whereas others (CaO, La, λ0, and λ1) require a bi‐linear fit with a change in gradient at ∼55 km. Observed variations in highly incompatible elements are consistent with degrees of melting that decrease with increasing lithospheric thickness, as expected. Nevertheless, at thicknesses beyond ∼55 km, the implied degree of melting does not decrease as rapidly as is suggested by theoretical expectations. This observation is robust across different lithospheric thickness estimates, including those derived from seismic constraints. We infer that at thicknesses exceeding ∼55 km, weak plumes fail to effectively thin overlying lithosphere and/or produce insufficient melt to erupt. This is supported by independent estimates of plume buoyancy flux, indicating that OIB magmatism on older lithosphere may be biased toward hotter plumes. In addition, we find evidence for a “memory effect” of incomplete homogenization of melts during their ascent.

The lithospheric mantle is cool and refractory.Accordingly, it is unlikely to melt and generate magmas (e.g., Katz et al., 2003).In addition, the lithosphere is highly viscous and is therefore difficult to mechanically deform (e.g., Burov & Gerya, 2014;Burov et al., 2007;Campbell, 2007;Duvernay et al., 2021Duvernay et al., , 2022;;Jones et al., 2017).As a consequence, it is expected to act as a lid that limits plume ascent and thereby dictate the lowest melting pressure for plume-derived melts (Figure 1).This behavior is the so-called "lid effect," first proposed by Watson and McKenzie (1991) and subsequently examined in several studies at both global (e.g., Dasgupta et al., 2010;Ellam, 1992;Haase, 1996;Humphreys & Niu, 2009;Niu, 2021;Niu et al., 2011) and regional scales (e.g., D. R. Davies, Rawlinson, et al., 2015;Gibson & Geist, 2010;Hole & Millett, 2016;Klöcking et al., 2018;Liu et al., 2016).Despite this extensive body of work, a complete and statistically rigorous assessment of the relationship between lithospheric thickness and the geochemistry of plume-derived magmas has not yet been established: previous studies have either described this relationship qualitatively or only made use of simple linear statistics (e.g., D. R. Davies, Rawlinson, et al., 2015;Ellam, 1992;Humphreys & Niu, 2009;Niu, 2021;Niu et al., 2011).Several important questions remain, including: 1. Do available geochemical data statistically support existence of a lid effect?2. Are observed trends consistent with theoretical expectations for partial melting at different pressures?3. What other processes might be affecting observed trends?
The last few years have seen progress in several areas that are pertinent to answering these questions.There has been a steady increase in the quantity and accessibility of high-quality data available on melt geochemistry, improvements in the accuracy and resolution of lithospheric thickness data sets, and the advent of comprehensive statistical techniques to examine any potential relationship between the two.There is, therefore, an opportunity to undertake a critical revaluation of evidence for the lid effect.
Our study exploits an extensive and carefully curated data set of geochemical analyses for OIBs, extracted from the ever-growing open-source GeoRoc database (https://georoc.eu).The data set is filtered to eliminate those samples whose geochemistry has been strongly altered after initial magma generation.Alongside the geochemical Figure 1.Schematic cartoon illustrating how oceanic lithosphere acts as a lid, hindering the ascent of mantle plumes.The dashed line represents the spinel-garnet transition.When a plume rises beneath thin lithosphere, large melt volumes will be produced with more melts generated within the spinel stability field, thus exhibiting a low-pressure signature.Conversely, when a plume rises beneath thick lithosphere, melt volumes are smaller and melting will principally occur within the garnet stability field, displaying a high-pressure signature.
parameters examined by previous studies, we analyze geochemical diagnostics on Rare Earth Elements (REEs) that have been recently proposed by O'Neill (2016) and are expected to show a clear pressure signal owing to their sensitivity to the degree of melting and the spinel-garnet phase transition.The latter, a pressure-sensitive aluminum-rich phase change, induces a substantial change to the peridotite mineral assemblage, with different REEs exhibiting varying compatibility between the two phases (e.g., Sun & Liang, 2013;Wood et al., 2013).Furthermore, we exploit new estimates of lithospheric thickness, based upon both theoretical models of oceanic spreading and observational constraints from seismic data (Hoggard, Czarnota, et al., 2020;Richards, Hoggard, Crosby, et al., 2020).Using a probabilistic Bayesian approach that is capable of detecting sharp changes in geochemical trends, we investigate the role of lithospheric thickness in controlling OIB geochemistry and explore the mechanisms that underpin the trends that we observe.
The remainder of our paper is structured as follows.In Section 2.1, we introduce our OIB database, our approach to filtering this data, and the geochemical diagnostics examined.In Section 2.2, we describe the lithospheric thickness estimates at each individual island, derived using both plate-cooling models and local constraints from surface-wave tomography models.In Section 2.3 we present a probabilistic Bayesian approach developed and utilized to analyze relationships between geochemistry and lithospheric thickness.Our results are presented in Section 3, with their sensitivities, implications for our understanding of the lid-effect, the role of the lithosphere in modulating plume melting, and other processes that may affect OIB chemistry and contribute to the trends, discussed in Section 4.

Geochemical Data Set
In compiling our geochemical database of the products of plume melting in oceanic settings, we have chosen to focus solely on OIB data and neglect data associated with LIPs.This omission is because LIPs are the meltproducts of plume heads and are also often associated with continental break-up.They regularly occur in the vicinity of the continent-ocean boundary and consequently often display a strong crustal signature (e.g., Chung & Jahn, 1995;J. H. F. L. Davies et al., 2021;Owen-Smith et al., 2017).It is also therefore difficult to estimate lithospheric thickness at the time of eruption (e.g., Courtillot et al., 1999;Hill, 1991).

Source of Analyses
Geochemical data for major and trace element concentrations are compiled from OIB data in the GeoRoc database.As the number of high-quality glass samples is limited, the data are derived principally from analyses of bulk rocks (with some additional glass analyses where available).The GeoRoc database contains geochemical information from over 20,000 OIB samples from the Atlantic, Indian and Pacific Oceans, with their locations mapped in Figure 2 and listed in Tables S1 and S2 of the Supporting Information S1.Our database incorporates concentrations of major (SiO 2 , Al 2 O 3 , MgO, FeO, TiO 2 , Na 2 O, K 2 O, CaO, P 2 O 5 ) and trace elements (REEs, U, Nb, Ba, Th), as well as derivative parameters describing REE patterns (λ 0 , λ 1 , and λ 2 , from O'Neill, 2016).Major elements with high concentrations are likely influenced by the stabilities of minerals under varying pressure and their compatibilities in mantle peridotite.We expect that major elements with lower concentrations (usually <5 wt.%) and trace elements are sensitive to phase changes and the degree of melting, which, in turn, are sensitive to pressure.The combined use of both major and trace element parameters can therefore offer a more complete picture of the impact of the lithospheric lid on mantle melting processes.

Database Filtering
Melts generated from peridotite melting are subject to various physiochemical processes during their ascent and whilst residing in magma chambers, such as fractional crystallization and crustal assimilation (e.g., Class & Goldstein, 1997;Sisson & Grove, 1993;Straub et al., 2013;Ubide et al., 2022).Additionally, post-eruptive hydrothermal alteration can substantially alter the original chemical signature of basalts (e.g., Khogenkumar et al., 2016;Saito et al., 2015).Some previous studies of the lid effect have chosen to use all available OIB geochemical data without attempting to screen samples that are heavily impacted by these additional processes (e.g., Humphreys & Niu, 2009).In our analyses, however, we have filtered OIB samples to isolate those that exhibit a composition most similar to that of the primitive magma.We therefore restrict our data set to samples  2020) with locations of selected ocean island basalt samples (black dots).Only the name of the archipelago for each island group is displayed, but each individual island's lithospheric age and thickness are considered separately during the analysis.(b) Present-day oceanic lithospheric thickness based on a global plate-cooling model (Richards, Hoggard, Crosby, et al., 2020).(c) Present-day oceanic lithospheric thickness constrained by surface-wave tomography (Hoggard, Czarnota, et al., 2020).
that have not undergone excessive alteration or fractional crystallization after initial generation.We do so by applying the following filters to the data: 1.Only those samples with SiO 2 43-54 wt.% are accepted in order to exclude melts that fall outside of the basalt field (Figure 3a); 2. Only samples with MgO 7-16 wt.% are accepted.Values with MgO < 7 wt.% are likely to have been subjected to extensive fractional crystallization (e.g., Sisson & Grove, 1993) and may contain clinopyroxene and/or plagioclase phenocrysts or have experienced clinopyroxene and plagioclase crystallization, complicating interpretation of major element trends.Samples with MgO > 16 wt.%are rejected because they are likely to contain olivine phenocrysts (Figures 3a and 3b, e.g., Albarède et al., 1997); 3. Samples with a loss on ignition > 3 wt.%are rejected to eliminate basalts subjected to excessive levels of posteruptive hydrothermal alteration (e.g., Greenberger et al., 2012); 4. Samples with Nb/U < 30, La/Nb > 1.2, or La/Ba and Nb/U values outside of the ellipse of Fitton et al. (1991) are rejected because they are likely to have been contaminated by continental crust (e.g., Condie, 1999;Hofmann, 2003;Rudnick, 1995;Figures 3c and 3d).
Applying these filters to the global OIB data set results in a subset of 1,737 samples, each consisting of concentrations of major elements, trace elements and REEs.

Correction for Fractional Crystallization
When magma travels through the lithosphere or remains in a magma chamber, any fractional crystallization that occurs alters the concentration of major and trace elements in the remaining melt (e.g., Jackson et al., 2012; Ubide  et al., 2022).Provided that the mineral phases that have crystallized are not complex, we can "revert" this process to estimate concentrations of both major and trace elements in the primary magma.To do so, we use the Pet-rolog3 software (Danyushevsky & Plechov, 2011) to reintroduce olivine into evolved OIBs until MgO concentrations reach 16 wt.%,which is the assumed MgO content of magma that is in equilibrium with the mantle (e.g., Norman & Garcia, 1999).Despite some studies showing that minerals fractionate throughout magma ascent (e.g., Liu et al., 2016;Lundstrom et al., 2003), we make the simplifying assumption that this olivine did not crystallize until melts reached a magma chamber at ∼0.3 GPa (∼9 km depth).This choice of depth roughly coincides with the Moho, where the drop in density from mantle to crustal rocks results in melts becoming neutrally buoyant, allowing magma to remain in the chamber for a more extended period of time (Ryan, 1988(Ryan, , 1994)).In the continuous, pure fractional crystallization process, we assume that partition coefficients for trace elements in olivine remain constant.For each individual OIB sample, we use the major element calculations of Petrolog3 to determine how much olivine to "add back in" to obtain the composition of the primitive magma.
Accordingly, the concentration of each trace element in the primitive magma (c p ) is calculated via where c l is the measured concentration of each element in the sample, X is the fraction of olivine crystallized, and D is the associated partition coefficient (Shaw, 1970).

Shape of REE Patterns
Due to their high charge and large ionic radii, REEs behave as incompatible elements in most mantle minerals.Moreover, the consistency of REE chemical valence makes their ionic radius systematically decrease with increasing atomic number (so-called lanthanide contraction; Ahrens, 1952).Since REEs occupy identical crystal lattice positions, their partition coefficients therefore exhibit a systematic dependence on atomic number, with lower atomic number REEs (Light Rare Earth Elements; LREEs) possessing larger radii and being more incompatible.Accordingly, during partial melting, REEs with a smaller atomic number more preferentially enter the melt than their heavier counterparts, an imbalance that is particularly pronounced at low degrees of melting.A caveat to this basic behavior is that heavy rare earth elements (HREEs) readily substitute for Al 3+ in garnet and, hence, can be compatible in garnet.As such, low-fraction melts generated within the garnet stability field will have lower HREE concentrations than equivalent melts generated in the spinel stability field.Many laboratory experiments have been conducted to constrain the partition coefficients of REEs, with results consistent with these aforementioned theoretical predictions (e.g., Fujimaki et al., 1984;Johnson, 1994Johnson, , 1998;;McKenzie & O'Nions, 1991).It is also worth noting that, due to the general incompatibility of REEs in all low-pressure mineral phases, their relative proportions are generally unaffected by fractional crystallization at low pressure., where positive and negative λ 2 indicate concave or convex REE patterns, respectively.λ 2 exhibits a non-linear, non-monotonic relationship as a function of changing source composition and degree of melting, so we therefore limit ourselves to only a brief discussion of its results.In contrast to the simple ratios between two REEs, such as Ce/Y and La/Sm, that have been extensively used in previous studies (e.g., Ellam, 1992;Humphreys & Niu, 2009;Niu, 2021), λ i considers all REEs except Eu and is more robust to the idiosyncrasies of individual element behavior.

Lithospheric Thickness and Eruptive Age
The thickness of oceanic lithosphere as a function of ocean floor age is commonly approximated through one of two theoretical cooling models: (a) the half-space model, in which lithospheric thickness is proportional to the square root of lithospheric age (Turcotte & Oxburgh, 1967); and (b) the plate model, where lithospheric thickness increases with plate age, but asymptotes toward a constant value beyond a certain age due to heat resupply from below (McKenzie, 1967).Plate-model predictions have been shown to provide an improved match to heat flow and bathymetry observations in older ocean floor and also inferences of lithosphere-asthenosphere boundary depth obtained from seismology (McKenzie, 1967;Parsons & Sclater, 1977;Richards, Hoggard, Crosby, et al., 2020;Richards et al., 2018).Accordingly, the plate model is our preferred reference and we test two different versions of it: one derived from globally averaged subsidence and heatflow data and the second providing optimal fits to subsets of these data from each individual oceanic basin (Atlantic, Indian and Pacific Oceans; Richards, Hoggard, Crosby, et al., 2020).We present results for the latter in the main text and also conduct assessments of the sensitivity of our results to this choice, with a summary presented in Supporting Information S1.In all cases, the potential temperature in the model is fixed to 1333°C and the base of the lithosphere is assumed to follow the 1175 ± 50°C isotherm (Richards, Hoggard, Crosby, et al., 2020).
A limitation of theoretical cooling models is that they assume oceanic lithospheric thickness varies solely as a function of ocean-floor age and, hence, cannot capture local deviations away from this average behavior.Seismological observations, particularly from surface-wave tomography, provide a way of mapping these local variations in lithospheric thickness, including those potentially induced by the impingement of mantle plumes (Ballmer et al., 2011;Duvernay et al., 2022;Richards, Hoggard, White, & Ghelichkhan, 2020;Schaeffer & Lebedev, 2013).Accordingly, to complement our plate-model derived estimates of lithospheric thickness and explore the sensitivity of our results to regional lithospheric thickness variations, we also make use of a seismologically derived model of lithospheric thickness from Hoggard, Czarnota, et al. (2020).
We separate ocean islands into two categories: products from off-axis and on-axis plumes.For off-axis islands, we estimate lithospheric thickness using the aforementioned plate-cooling and seismologically derived models.Unfortunately, neither theoretical cooling models nor global-scale seismic estimates are good at constraining lithospheric thickness above on-axis plumes.The former do not capture the consequences of increased melt generation and hence thicker crust above plumes, while the latter suffer from the limited resolution of surface waves at depths shallower than ∼75 km (Priestley & McKenzie, 2006;White & McKenzie, 1989).For on-axis islands, we therefore obtain lithospheric thickness from local estimates of crustal thickness, assuming that melting extended to the top of the underlying mantle as is observed in ophiolites (e.g., Pallister & Hopson, 1981).Seismic estimates for Moho depths are as follows: Iceland ∼20-30 km (White et al., 1996); Ninetyeast Ridge, ∼15-25 km (Grevemeyer et al., 2001); Walvis Ridge ∼10-25 km (for lithosphere that is now aged between 60 and 100 Ma; Goslin & Sibuet, 1975;Graça et al., 2019).At each of these sites, we calculate average lithospheric thickness according to 1 2 (h max + h min ) , where h max and h min are the maximum and minimum estimates of Moho depth, respectively.
Estimating lithospheric thickness at the time of eruption requires knowledge of lithospheric age at that time, which can be obtained by subtracting the OIB age from the present-day lithospheric age (Figure 2a).Present-day lithospheric age for each island is obtained from the global grid of Seton et al. (2020), with the age range of OIBs on each island constrained, where possible, by the onset and termination of the shield stage of volcanism or, in cases where geological constraints on the shield period are unavailable or unclear, the maximum and minimum age of OIB samples (Tables S1 and S2 in Supporting Information S1).
To estimate lithospheric thickness at the time of eruption for off-axis plumes, we assume that both the present-day lithospheric age (t crust ) and the OIB age (t OIB ) on each island follow a Gaussian distribution as where μ 1 is the oceanic crustal age, σ 1 is half of the age misfit, μ 2 is the mean of maximum and minimum OIB eruption ages, and σ 2 is a quarter of the length of the OIB major eruption period.t crust and t OIB can be considered as independent random variables, thus the age of oceanic lithosphere at the time of OIB volcanism (t erupt ) should also follow a Gaussian distribution given by Lithospheric thickness is estimated from the theoretical cooling models by assuming that it lies between the 1125 and 1225°C isotherms.We assume that lithospheric thickness (z) at a given time follows a Gaussian distribution according to in which μ 3 is the mean of the lithospheric thickness obtained from the 1125 and 1225°C isotherms and σ 3 is a quarter of the difference in depth between them.For each island, we randomly choose a t erupt based on Equation 4and calculate the corresponding lithospheric thickness using Equation 5. Iteratively repeating this process until reaching a stable distribution of thickness estimates yields the plate-model-derived mean value of lithospheric thickness beneath each ocean island.
For the seismically constrained estimates of lithospheric thickness, we test two end-member scenarios: (a) lithospheric thickness at the time of eruption is equivalent to that of the present day; and (b) following eruption and movement away from the location of the plume tail, the lithosphere has re-thickened to its present-day value by conductive cooling following a half-space model.The true scenario likely falls between these two assumptions.Both assumptions yield similar results, likely because the majority of OIBs in our data set are young (<10 Ma) and the lithosphere cannot substantially rethicken over such a short time frame.Correcting for this process makes no appreciable difference to our results (<5 km thickness change; see Tables S4 and S5 in Supporting Information S1) and the size of this correction is generally smaller than the depth range covered by the 1175 ± 50°C isotherms.When using seismically derived estimates of lithospheric thickness, we therefore adopt the first option above.
Estimated lithospheric thickness at the time of eruption, based on either the basin-specific plate models (Richards, Hoggard, Crosby, et al., 2020) or seismological constraints (Hoggard, Czarnota, et al., 2020), are tabulated in Tables S4 and S5 of the Supporting Information S1.Plate model thicknesses for the Atlantic basin are slightly greater than those derived from the global-average model, whereas in the Indian and Pacific basins, lithospheric thickness estimates from basin-based models are similar to, or thinner than, those of the global model (Figure S1 in Supporting Information S1).

Bayesian Model Selection
To investigate the variation of each geochemical parameter with lithospheric thickness, we have plotted and parameterized OIB geochemical data against lithospheric thickness at the time of eruption.To understand whether a particular data set suggests a trend, or a change in gradient, we make use of Bayes factors: the ratio of the evidence or marginal likelihood between two competing statistical models (Jeffreys, 1935;Kass & Raftery, 1995).The evidence represents the integral of the likelihood over the prior for a given model choice.In our case, it quantitatively evaluates how likely it is to generate the observed geochemical data set, based on a specified model (i.e., a function that describes the trend of the geochemical parameters against lithospheric thickness).Therefore, given two or more competing models, the model with the larger evidence is preferred.Computing the evidence is difficult, particularly for large dimension models, but for this problem we use Dynamic Nested Sampling (Skilling, 2006;Speagle, 2020), which gives both posterior and evidence estimates in a single analysis.
The geochemical data include the raw and fractional crystallization-corrected concentrations of major elements, trace elements and λ values calculated from REE concentrations.To determine whether a given geochemical parameter is sensitive to lithospheric thickness or influenced by any potential sudden changes in mantle composition, such as the phase change from spinel to garnet peridotite, three models were compared: (a) a constant value model (which would imply no sensitivity to lithospheric thickness); (b) a linear model (which suggests a lid effect); and (c) a bi-linear model that permits a change in gradient at some depth determined by the data (Figure 4).We choose not to examine exponential models since they are monotonic, so incapable of describing a reverse in a trend or detecting the depth of a potential trend change.High-order functions are also not considered, to avoid over-fitting and a loss of generalization.To estimate posterior probability densities of the model parameters for each candidate model, we choose an independent Gaussian likelihood which is written as where M is the number of islands, N i the number of samples for each island i, p i is the model prediction of the geochemical concentration for island i, d j is the observed data value for each sample from that island, and σ is the overall standard error.This formulation assumes that the data follows the standard normal distribution at each lithospheric thickness.
We fix values for the minimum and maximum lithospheric thicknesses (z min and z max in Figure 4), resulting in one, two and four unknown parameters for the constant, linear and bi-linear models, respectively.The use of Bayes factors to test the relative support of competing models is subtly affected by the choice of priors.Regarding priors for the y values (i.e., the geochemical data), we adopt an empirical Bayes approach and set the prior to be Gaussian with mean and standard deviation equal to that of the overall data.The mean and standard deviation of z (i.e., the lithospheric thickness of a putative transition in the trend for the bi-linear model) are assumed to be 60 and 5 km respectively, close to the average of all thickness data.To test the sensitivity of evidence calculations, we ran repeated tests that changed the standard deviation of the prior by ±10%, which resulted in an average change of log 10 evidence values of ±0.14.Similarly, changing the standard deviation by ±50% (a comparatively large change in the prior) resulted in an average change in the log 10 evidence of ±0.18.We are therefore confident that the choice of priors for our Bayesian evidence calculations are reasonable and that sensitivity to the choice of priors is minor.Nonetheless, in the Bayesian evidence results herein, a reasonable error bound on the numerical log 10 evidence values would be ±0.2.
The evidences for constant, linear and bi-linear models are denoted as E 0 , E 1 , and E 2 , respectively.Since evidence values are typically vanishingly small numbers, they are usually represented by their logarithms.A candidate model with a larger evidence value is to be preferred, for example, model "A" with a log 10 evidence of 1,000 is a hundred times more likely than a competing model "B" with a log 10 evidence of 1,002.Generally, a difference in the log 10 evidence greater than 2 is taken to be statistically significant (Jeffreys, 1935;Kass & Raftery, 1995).We note that model evidence values can only be compared for the same geochemical parameter, as they are influenced by the range of values in the data and sample sizes.As a reminder, the constant model implies that a geochemical parameter is insensitive to changes in lithospheric thickness.The linear model can detect an overall trend but is incapable of describing a change or reversal in trend.The bi-linear model can be useful for identifying a change point in a trend and even detecting a reversal of the trend, but is more sensitive to outliers.For a given geochemical parameter, if log 10 E 1 log 10 E 0 > 2, we are confident in saying that it varies with lithospheric thickness.Furthermore, if log 10 E 2 log 10 E 1 > 2, we can say that a change point or kink can be found in the data trend.In these cases, we provide histograms of the depth of the likely kink in the model and calculate its mean and standard deviation.

Sensitivity to Sites With Large Numbers of Samples
Due to the form of our likelihood function in Equation 6, clusters of large numbers of measurements from a single site could potentially bias the results.Two notable examples of this are the large OIB sample sizes of Iceland and Hawaii.To test the robustness of our results to potential biasing from these two localities, we repeat the calculation of posterior probability densities and evidence values for each geochemical parameter using: (a) all data (i.e., our reference case); and (b) the data set with samples from both Iceland and Hawaii excluded.Removal of Hawaiian samples is of particular relevance because they represent the only OIBs located on thick lithosphere that are dominated by tholeiites (e.g., MacDonald & Katsura, 1964).All other OIBs at and beyond these lithospheric thicknesses consist predominantly or exclusively of alkali basalts (e.g., Fisk et al., 1988;Gautier et al., 1990;Schmincke, 1982).

Results
To provide a relatively simple overview that gets at the essence of our results, we have chosen to focus in the main text on a preferred reference case.This case includes the initial correction of geochemical concentrations for the effects of fractional crystallization, uses data from all OIB localities within our database, and adopts lithospheric thicknesses from basin-specific plate-cooling models.While we discuss any important differences that arise from changes to this reference setup in the main text, the full suite of associated figures and results are presented in Supporting Information S1.

Geochemical Histograms
Raw histograms of major element concentrations for all OIB data, before application of sample filters, display slightly skewed Gaussian-like distributions with peaks at approximately 7 wt.% for MgO, 48 wt.% for SiO 2 , 3 wt.% for TiO 2 , and 14 wt.% for Al 2 O 3 (blue bars in Figures 5a-5d).The MgO peak at 7 wt.%broadly coincides with the minimum in magma density at 7-8 wt.% MgO calculated using Petrolog3 at 0.1 GPa, which is consistent with expectations that the lightest magmas are the most likely to erupt at the surface (see Figure S2 in Supporting Information S1; Danyushevsky & Plechov, 2011).The continuous distribution of major element concentrations is consistent with expectations for mixing of distinct, end-member reservoirs to varying extents, which is also supported by isotopic evidence (e.g., Hart et al., 1992).Filtering the raw data according to the criteria outlined in Section 2.1.2has limited impact on distributions for SiO 2 , Al 2 O 3 , and TiO 2 , but the filtered MgO histogram retains only the right-hand side of the distribution due to the sharp cut-off of samples with MgO < 7 wt.%(green bars in Figures 5a-5d).Histograms of the REE shape parameters for filtered OIB samples exhibit more scatter and less clean unimodal behavior (Figures 5e-5g).Nevertheless, λ 0 has a clear peak at ∼3.3.λ 1 is left skewed, with a peak around 10 and more than 80% of samples have λ 1 > 5. λ 2 is somewhat bimodal, with a central peak at approximately 15 and a subsidiary peak at 40.

Example of Statistical Results
To illustrate our procedure for quantifying the relationship between lithospheric thickness and various geochemical parameters, we present two examples for Al 2 O 3 and λ 1 in Figure 6.Both use our reference setup, in which the global OIB data set is filtered and corrected for fractional crystallization, with lithospheric thickness evaluated via the basin-specific plate model.Black crosses represent individual OIB samples and log 10 of the evidence is provided for each of the three types of model.
For both Al 2 O 3 and λ 1 , we find that the evidence increases by ∼160 when moving from constant to linear models (compare Figure 6a with Figure 6b, and Figure 6e with Figure 6f), supporting the existence of a lid effect for both Al 2 O 3 and λ 1 .However, we see contrasting results when the bi-linear model is introduced.For Al 2 O 3 , evidence values for linear and bi-linear models are similar (Figures 6b and 6c), implying the absence of any obvious transition in the trend as a function of lithospheric thickness.The resulting probability distribution of potential kink depths is therefore broad and poorly constrained in Figure 6d, and we infer that Al 2 O 3 in OIBs decreases linearly with increasing lithospheric thickness, with no definitive kink.On the other hand, λ 1 shows a clear preference for a bi-linear model, with an increase in the log 10 evidence value of ∼87 over a linear model (Figures 6f and 6g).The associated probability distribution for the kink is tightly constrained in the depth range of 49-56 km, with an average of ∼52 km (Figure 6h).Based on this preferred bi-linear model, the most likely trend for λ 1 is that it increases with lithospheric thickness until a depth of ∼52 km, before subsequently remaining approximately constant.

Summary of Evidence Evaluation Results
Values of log 10 E 1 log 10 E 0 and log 10 E 2 log 10 E 1 have been determined for each geochemical parameter, under our reference setup.As a reminder, when greater than a key threshold value of two (i.e., more than hundred-fold increase in the likelihood), the former indicates statistical preference for a linear model over a constant one, while the latter indicates a bi-linear rather than linear relationship.
The preferred model for each geochemical parameter is shown in Figure 7, with further details in Figures 8-11, and can be summarized as follows: 1.All geochemical parameters prefer either a linear or bi-linear model over a constant model, indicating universal sensitivity to lithospheric thickness.2. For major elements Al 2 O 3 , FeO, and SiO 2 , data are optimally fitted by linear models (Figure 8).Conversely, TiO 2 , Na 2 O, K 2 O, CaO, and P 2 O 5 data are optimally fitted by bi-linear models (Figures 9c-9g); 3.For trace elements, the highly incompatible elements La and Th are best fitted by bi-linear models (Figures 9a  and 9b), whereas the less incompatible Yb and Lu are best fitted by linear models (Figure 10); 4. For parameters describing REE patterns, λ 0 and λ 1 are optimally fitted by bi-linear models (Figures 11a and  11b), whereas λ 2 prefers a linear model (Figure 11c). 5.For geochemical parameters that prefer a bi-linear model, kink depths generally occur at lithospheric thicknesses of 50-60 km.

Existence of a Lid Effect
Lithospheric thickness dictates the minimum pressure of plume melting through the so-called "lid effect."It affects OIB chemistry in two ways (e.g., Humphreys & Niu, 2009;Niu, 2021;Watson & McKenzie, 1991).First, by inhibiting upwelling beyond a certain depth, lithospheric thickness limits the maximum degree of melting (expressed as fraction, F).We therefore expect F to be inversely proportional to lithospheric thickness, which will have a substantial impact on the concentrations of highly incompatible trace elements.Second, the pressure at which melting occurs has strong implications for the mineral phases present in the residue following partial melting.In particular, over the depth range of interest here, the stable aluminum-rich phase converts from garnet (Mg 3 Al 2 Si 3 O 12 ) to spinel (MgAl 2 O 4 ) with decreasing pressure, subsequently becoming plagioclase (CaAl 2 SiO 8 ) at shallow depths beneath mid-oceanic spreading centers (e.g., Masaaki, 1980).Despite our analyses being subject to uncertainty, particularly in relation to estimates of lithospheric thickness and assumptions on uniform source composition, the data support a linear or bi-linear trend between all geochemical parameters and lithospheric thickness, providing universal evidence for the lid effect and corroborating the conclusions of, for example, Humphreys and Niu (2009), Dasgupta et al. (2010), andNiu (2021).
Nonetheless, it is clear from our results that different geochemical parameters exhibit distinct responses to the lid effect.Some trends (e.g., Al 2 O 3 ) show a linear relationship with lithospheric thickness, whereas others show a bilinear relationship with an abrupt change at a certain depth (e.g., λ 0 , λ 1 ).In the following sections, we discuss potential explanations for these behaviors.
We start with major element trends that are best fitted by linear models, with an emphasis on the relationship to pressure-dependent mineral assemblages.
We then discuss the remaining major and trace elements, relating observed trends to the influence of variations in F and the spinel-garnet phase transition.Although Yb and Lu are best fitted by linear models, we include them in this section because their behavior is associated with an interplay between F  and the spinel-garnet phase transition.We finish by discussing REE trends, described by λ i , drawing on the lessons learned from the interpretation of trace elements trends.

Major Elements With Linear Trends
Concentrations of the major elements SiO 2 , FeOT, and Al 2 O 3 in OIBs show a linear dependence on lithospheric thickness (Figure 8).In mantle melts, these components are known to be buffered by the mineral assemblage of the mantle residue and the observed trends are consistent both with experimental studies on the Calcium, Magnesium, Aluminum, Silicon system (e.g., Walter & Presnall, 1994) and with the results of previous observational studies (e.g., Humphreys & Niu, 2009;Niu, 2021;Niu et al., 2011).SiO 2 exhibits a moderate decrease with increasing lithospheric thickness (Figure 8a).Its concentration in mantle melts is buffered by the two most abundant minerals in the upper mantle, olivine and orthopyroxene, according to the reaction Increasing pressure drives this reaction to the right, expanding the stability field of orthopyroxene at the expense of olivine (e.g., Campbell & Nolan, 1974;Walter & Presnall, 1994).As a consequence, as the average melting pressure increases beneath thicker lithosphere, the residue contains more SiO 2 -rich orthopyroxene and the corresponding melts produced are increasingly SiO 2 -poor (e.g., Bohlen & Boettcher, 1981;Bohlen et al., 1980).We note that Herzberg (1992) further proposed that the decrease in SiO 2 with increasing melt pressure stops at ∼45 wt.% SiO 2 due to low degrees of melting in the presence of garnet, but this cut-off behavior is not observed in either our analyses or in previous studies (e.g., Dasgupta et al., 2010;Scarrow & Cox, 1995).We therefore suggest that the spinel-garnet transition has limited influence on the SiO 2 content of OIBs, with Reaction 7 and associated buffering of the silica content by olivine and orthopyroxene being the key control.
We also attribute the linear increase in FeOT with increasing lithospheric thickness (Figure 8b) to the relative stabilities of olivine and orthopyroxene as a function of pressure.Olivine contains more Fe than orthopyroxene and increasing the pressure stabilizes orthopyroxene at the expense of olivine.As a consequence, for similar degrees of melting, high-pressure melts contain more Fe than low-pressure melts.The relative abundance of olivine and orthopyroxene in the residue was also used by Niu (2016) to explain the increase in FeOT in mid-ocean ridge basalts (MORB) with increasing ridge axial depth.Analysis of our OIB data set suggests that this trend can be extended over a greater depth range than is possible with the MORB data alone.
Al 2 O 3 linearly decreases with increasing lithospheric thickness (Figure 8c), which we believe can be attributed to an increase in the Al content of clinopyroxene and, to a lesser extent, orthopyroxene, with increasing pressure.Al 3+ can occupy either the tetrahedral or octahedral sites within the pyroxene crystal lattice.The two tetrahedral sites are characterized by a central cation (usually Si 4+ ) surrounded by four oxygen atoms, whereas the two larger octahedral sites are positions in which the central cation is surrounded by six oxygen atoms and are usually occupied with cations that have greater ionic radii, such as Ca 2+ .Increasing pressure shrinks the octahedral M1 and M2 sites in both pyroxenes and allows more Al 3+ to enter the M1 site, which is the smaller of the two octahedral sites (Colson & Gust, 1989).The octahedral Al 3+ can either be charge balanced by additional Al 3+ replacing Si 4+ in an adjacent tetrahedral site according to the reaction or by Na + or K + replacing a divalent ion on the larger M2 site (Campbell & Borley, 1974;Safonov et al., 2011), according to reaction As a consequence, the Al 2 O 3 content of the residual pyroxenes increases with increasing pressure and the Al 2 O 3 concentration in the melt correspondingly  decreases.This simple interpretation may be complicated, however, by Al 3+ being buffered by reactions between spinel and pyroxene in the spinel stability field and between garnet and pyroxene in the garnet stability field.We might therefore have also expected a dependence on the spinel-garnet transition.The fact that our Al 2 O 3 trends do not require a bi-linear model in our reference setup suggests that this is not the case, perhaps because both spinel and garnet contain two Al 3+ ions and the increasing Al content of pyroxenes with pressure is therefore not affected by the spinel-garnet transition.
In summary, we infer that variations in the concentration of major elements SiO 2 , FeOT, and Al 2 O 3 in OIBs are dominated by gradual changes in mineral assemblage as a function of pressure rather than variations in F or effects arising from the spinel-garnet phase transition.

Major and Trace Elements With Bi-Linear Trends
The behavior of trace elements, which do not form stoichiometric components in minerals, can be understood using distribution coefficient, D, for the partitioning of the element between a mineral and the melt.During partial melting of mantle peridotite, the concentration of a given element in the aggregate melt (C l ) during batch melting is given by and, in the case of fractional melting, is where D′ is the bulk partition coefficient, F is the degree of melting and C s is the concentration of the element in the source before melting (Shaw, 1970(Shaw, , 1979)).C l is therefore controlled by the combined effect of D′ and F. Nevertheless, for incompatible trace elements where D′ is low (usually <0.01), these equations can be simplified to m indicating that C l is proportional to 1 F , regardless of the melting mechanism.The lower the value of D′, the more reliable this approximation becomes.Similarly for moderately low values of D′ (i.e., <0.2), given that degrees of melting for OIBs are never higher than ∼0.2, we can simplify Equations 10 and 11 to In cases where D′ > 0.2 and differences in the partition coefficients for different minerals are large, which can occur, for example, across the spinel-garnet phase transition, C l is influenced by both pressure and the phase change.The case of spinel-garnet can be represented by a reaction between spinel and pyroxene to give garnet and olivine according to The transition is abrupt and temperature dependent (e.g., Klemme & O'Neill, 2000).For temperatures appropriate for mantle plumes, the transition occurs over a ∼5 km depth range somewhere between 70 and 85 km, depending on the mantle composition (e.g., Klemme & O'Neill, 2000;Robinson & Wood, 1998;Tomlinson & Holland, 2021;Wood et al., 2013).Note that the shallower plagioclase-spinel transition in peridotite is not relevant to this study because the plagioclase stability field extends only to pressures of 0.8 GPa in fertile lherzolite and to 0.6 GPa in depleted lherzolite, corresponding to depths of 24 and 18 km, respectively (Borghini et al., 2010).
For geochemical parameters that are best fitted by bi-linear models, we divide them into two groups: (a) elements exhibiting low partition coefficients (D′ < 0.2), including Th, La, Ti, P, K, and Na, which we propose can be interpreted primarily in terms of F (although Na and K may be further influenced by D values for pyroxenes at high pressure); and (b) Ca, which requires consideration of both F and the spinel-garnet phase change.These bilinear trends were not identified in previous studies (e.g., Dasgupta et al., 2010;Ellam, 1992;Humphreys & Niu, 2009;Niu, 2021;Niu et al., 2011) and we discuss their likely origin.

Incompatible Elements
D′ values for the incompatible elements investigated in this study decrease in the following order: Ti ≈ P ≈ Na > La > K > Th.All are optimally fitted by bi-linear models, in which their concentrations initially increase rapidly with increasing lithospheric thickness, before remaining flat or increasing at a significantly reduced rate in the cases of Th, La, Ti, P, or slightly decreasing in the cases of K and Na.The kinks in slopes all occur at 50-60 km depth (Figure S3 in Supporting Information S1).
La, K, and Th are highly incompatible in peridotites, regardless of whether the major aluminum-rich phase is spinel or garnet (D′ ≤ 0.01; Table 1).Following on from our interpretation of Equations 12 and 13 for such elements, at constant potential temperature, we expect variations in their concentration to be proportional to 1 F as a function of lithospheric thickness and insensitive to the spinel-garnet phase transition (dashed line in Figure 12a).This behavior should impart an increase in incompatible trace element concentrations at larger thicknesses, with a steeper rate of increase at greater thicknesses.This prediction is consistent with the observed increase in incompatible element concentrations with increasing lithospheric thickness beneath thinner lithosphere.When lithospheric thickness exceeds 50-60 km, however, it is not consistent with the slightly increasing, flat or decreasing concentrations observed.We can further demonstrate this aspect by converting our observed concentrations of La into estimates of F as a function of lithospheric thickness (solid purple line in Figure 12a) and comparing it to the predicted F curves (note that the resulting F curve is insensitive to the choice of La or Th).There is an agreement between the shapes of the two curves for lithospheric thickness <55 km but they become inconsistent at larger thicknesses, implying that another process modulates  concentrations of incompatible elements beyond thicknesses of ∼55 km.The most important conclusion that can be drawn from the analyses of highly incompatible Th and La is that F remains nearly constant for lithosphere thickness >55 km.
The moderately incompatible elements Na, P and Ti have D′ ∼ 0.06-0.08 in both the spinel and garnet stability fields and follow similar trends.For these elements, D′ cannot be neglected and Equation 13should be used to interpret changes in their concentrations.Since D′ varies little with mineralogy for these elements, it can be regarded as a constant and F becomes the dominant variable.As a consequence, experimental and theoretical constraints imply that these moderately incompatible element concentrations should again increase with increasing lithospheric thickness, with steeper rates of increase at greater thicknesses.The contribution from D′ should reduce the effect of F, diluting the concentration ratio of these elements between the melts and residue at higher pressures without altering the underlying trend.This prediction is consistent with observations for lithospheric thicknesses less than 55 km, but it is inconsistent with thicker lithosphere trends, which again suggest minimal changes in F at larger thicknesses.This aspect is important to keep in mind for the following interpretations.
Our analyses demonstrate that the concentrations of Na and K differ from other incompatible elements in that they show a slight decrease with increasing lithospheric thickness beyond the kink (Figures 9e and 9f).At these pressures, F is expected to be small and to remain nearly constant with increasing lithospheric thickness.The observed trends in Na and K may therefore indicate that variations in D′ are playing a role.As discussed in Section 4.2, increasing pressure allows entry of Al 3+ into the clinopyroxene M1 site.Via Reaction 9, this substitution can be charge balanced by replacing a X 2+ cation in the larger M2 site with Na + and/or K + , resulting in an increase of D′ for K + and Na + with increasing pressure.With minimal changes in F, this effect may explain their observed decrease in concentration with increasing lithospheric thickness, relative to highly incompatible La and Th, although we note that Na + is expected to have a stronger affinity for the M2 site than K + because its ionic size is closer to the size of the site (Safonov et al., 2011).

CaO
Calcium is the element most likely to be affected by the spinel-garnet transition because garnet contains stoichiometric Ca, whereas spinel does not.The principal repositories for Ca 2+ in garnet peridotites are, in order of decreasing affinity, clinopyroxene > garnet > orthopyroxene > olivine.Beneath shallow lithosphere, there is a steady decrease in CaO concentration with increasing lithospheric thickness up to ∼55 km (Figures 9m and 9n).
We attribute this behavior to the continuous decrease in F, previously deduced from analyses of incompatible element trends: as F decreases, less clinopyroxene melts and the resulting melt has a lower Ca concentration.Our interpretation of the incompatible element trends suggests that, beyond the kink, F should remain approximately constant or continue to decrease, but at a reduced rate.Therefore, we would expect a further decrease in CaO (albeit at a lower rate), rather than the increase that is observed.The cause of this increase is unclear, but assuming that F is not changing (as suggested by the most incompatible elements), it requires that, with increasing pressure, a Ca-rich phase (presumably Ca-rich pyroxene) melts in preference to moderately Ca-poor garnet and orthopyroxene.We note that previous studies, which applied linear regression to the data, found no discernible trend between CaO and lithospheric thickness (e.g., Humphreys & Niu, 2009), probably because the reversal in trends from decreasing to increasing CaO counteract each other.

Yb and Lu
Concentrations of Yb and Lu change little with increasing lithospheric thickness, showing only a slight, linear decrease (Figure 10).This behavior occurs even though these elements exhibit an order of magnitude difference in compatibility between spinel (incompatible) and garnet (compatible; Table 1), from which we might expect to see a kink in their trends.
Within the spinel stability field, the decrease in F with increasing pressure, required by the incompatible trace element trends, is offset by increasing D′ due to an increasing amount of clinopyroxene in the residue at higher pressures (e.g., Green & Ringwood, 1967).Within the garnet stability field, the approximately constant F with increasing lid thickness is initially offset by an increase in D′ as garnet replaces spinel and, subsequently, as pressure continues to increase, by garnet partially replacing pyroxene according to Reaction 14. Increasing D′ may contribute to reducing the rate of increase for Yb and Lu within the spinel stability field, however it cannot produce the observed decreasing trend since these elements are incompatible in the spinel zone.
Superimposed on these changes is the migration of low Yb-Lu melt from the garnet zone into the spinel zone, where it partially offsets the potential increased concentrations of Yb and Lu due to their lower D′ in the spinel zone.This behavior is termed the "memory effect," whereby erupted melts preserve high-pressure signatures.Such "memory" could be the result of high and low-pressure melts mixing in a continuous melting column (e.g., Elliott et al., 1991), or direct sampling of melts across a range of depths without mixing.In the former scenario, relative proportions of melt from the spinel and garnet zone are critical to determining erupted melt composition: increasing pressure can increase Yb-Lu concentrations in melts within the spinel zone while simultaneously reduce the melt volume within that zone, thus limiting low-pressure melts' impact on the final average Yb-Lu concentration, and resulting in a decrease of Yb-Lu with increasing lithospheric thickness.The latter scenario is also reasonable because of the large scatter of Yb and Lu recorded in our data set.The true scenario is likely in between these two end-members.Overall, the combined impact of the competing influences of increasing F and the systematic decrease in the proportion of melt coming from the spinel zone with increasing pressure, sampling melts from various depths, and changes in D′ as garnet is replaced by spinel, can plausibly explain the slight decrease in Yb and Lu concentrations with increasing lithosphere thickness.Taken together, this result implies that the absence of a kink in Yb-Lu trends, which might be expected from the spinel-garnet phase transition, is obscured by the memory effect.

Shape of REEs
We next discuss the shape parameters, λ i , for REE concentration patterns (O'Neill, 2016).λ 0 is the average of the logarithmic concentration of all REEs except Eu, normalized by each element's concentration in chondrites.
Increasing lid thickness reduces the degree of melting, thereby elevating the concentration of highly incompatible LREEs (Figure 9) while having a limited effect on the concentration of HREEs (Figure 10).It is therefore not surprising that λ 0 follows trends defined by the highly incompatible elements, with a kink at ∼55 km (Figure 11a).
λ 1 measures the enrichment of LREEs relative to HREEs and has a bi-linear trend, similar to λ 0 (Figure 11b).Previous studies described variations in REE trends using ratios, such as La/Yb and Sm/Yb, and related observed changes to an increasing abundance of garnet in the residue as the melting pressure increases (e.g., Ellam, 1992;Humphreys & Niu, 2009).However, these studies did not recognize either a kink or the influence of changes in F on LREE and, hence, the slope of the REE pattern.As noted in connection with λ 0 , the LREEs initially increase with increasing lithospheric thickness to ∼55 km (Section 4.3.1),driven by changes in F, then remain nearly constant, while HREE concentrations change little throughout (Section 4.3.3).Therefore, the combined behavior of LREEs and HREEs accounts for the observed variation in λ 1 .
λ 2 measures the curvature of the REE pattern, and its value depends on the relative volume of melts from the spinel and garnet zones, as well as the source composition (Table S6 in Supporting Information S1).Even when considering melts at a constant pressure, increasing degrees of melting make λ 2 first decrease and then increase.Therefore, due to its non-monotonic relationship with multiple factors, its complexity offers limited information, and we avoid detailed interpretation of λ 2 trends.

Trend Robustness
To analyze the robustness of our results to initial OIB processing steps, potential bias toward heavily sampled localities, and/or choice of lithospheric thickness model, we evaluate Bayes factors for a suite of additional scenarios including: (a) filtered OIB data prior to-and post-corrections for fractional crystallization; (b) a subsample of the OIB data set, where a percentage of samples are randomly removed; (c) data sets where Iceland and Hawaii samples are included/excluded; and (d) lithospheric thickness obtained from the basin-specific platecooling models versus a model derived from seismic tomography.A summary of these results are presented in Figure 13, with further plots presented in Figures S4 and S5 of the Supporting Information S1.

Correction for Fractional Crystallization
Our reverse-fractionation calculations in Petrolog3 suggest that primitive OIB melts commonly undergo 5%-25% fractional crystallization in the magma chamber (Figure S6 in Supporting Information S1).Correcting for fractional crystallization does alter absolute geochemical concentrations, but the preference for a linear over constant dependence on lithospheric thickness remains unchanged for all major and minor elements, and REE shape parameters (Figure 13b).In most cases, the preference for either linear or bi-linear models also remains unchanged, the exception being SiO 2 , for which the non-corrected data shows a preference for a bi-linear model.This behavior is a direct consequence of imposing the SiO 2 >43 wt.% filtering criteria: for the non-fractionation corrected data, this results in a hard cut-off of all data below this value and preferentially increases preferred model values in thicker lithosphere, whereas following fractionation correction, the kink is smoothed out and SiO 2 linearly decreases with increasing lithospheric thickness (Figures S7 and S8 in Supporting Information S1).

Data Subsets
The possibility of bias from specific samples was tested by removing some data to see if the observed trends remain-so-called bootstrapping.We randomly removed 20% and 40% of all OIB samples while keeping all other variables consistent.In most of these tests, the preferred model and the trend for each of the geochemical parameters did not change (Figure 13c).Exceptions occurred when the preferred model was close to the edge of a threshold before removing some samples.Nonetheless, our overall interpretation remains valid, and the observed geochemical trends are not strongly affected by sampling bias.

Heavily Sampled Localities
Hawaii and Iceland are heavily sampled in comparison to other localities in our OIB database, yielding highdensity data clusters that may potentially introduce bias into the results.Nevertheless, we find that excluding these sites does not generally alter evidence in favor of the lid effect.For all geochemical parameters, there is still strong preference for either a linear or bi-linear model over a constant one.The two exceptions occur for Ca and Na, where the evidence with respect a constant model is still greater than 2 but less than 20 (Figure 13d).
With regards to preference for a bi-linear versus linear fit, removal of these localities has more of an influence on results.For geochemical parameters where all sites are optimally fitted by bi-linear models, excluding Iceland and Hawaii can either increase or reduce values of log 10 E 2 log 10 E 1 , depending on the parameter (Figures S4a and S4e in Supporting Information S1).In general, more parameters transition to preferring a bi-linear model (e.g., Al 2 O 3 , FeOT, Yb, and Lu), with Si maintaining preference for a linear model and only Th switching from a bilinear to linear model.This indicates that evidence in favor of bi-linear trends is not attributable to potential sampling bias from Hawaii and Iceland.
Nevertheless, it is interesting to note that excluding Icelandic and Hawaiian samples has an impact on the slope of the trend in thick lithosphere (beyond the kink depth of ∼55 km).For incompatible elements Th, La, TiO 2 , and P 2 O 5 , the slope after the kink increases (Figure S9 in Supporting Information S1), which is also the case for λ 0 and λ 1 (Figure S10 in Supporting Information S1), albeit still at a lower rate than would be expected from theoretical arguments for melting at constant potential temperature and composition.For Na 2 O and K 2 O, it reduces the rate of concentration decrease with increasing lithospheric thickness (Figure S9 in Supporting Information S1).All of these differences can likely be attributed to the concentrations of incompatible elements in Hawaiian basalts being lower than those of other OIBs on lithosphere of similar thickness.As noted in Section 2.4, Hawaiian basalts are dominated by tholeiites, whereas all other OIBs on thick lithosphere have alkali basalt affinities.Tholeiites are produced by higher degrees of partial melting (e.g., Yoder & Tilley, 1962) and concentrations of incompatible elements are therefore expected to be diluted, which is consistent with the Hawaiian plume being hotter and stronger than other plumes beneath thick lithosphere (e.g., Hoggard, Parnell-Turner, & White, 2020).Hawaiian OIBs may also originate from a source that is more depleted in highly incompatible elements if, as suggested by Hofmann and Jochum (1996) and Pietruszka et al. (2013), it contains a significant amount of recycled oceanic gabbro.Regardless of the exact nature of the Hawaiian plume, its distinctive characteristics, coupled with the large number of samples available, can influence the slope of geochemical trends in thick lithosphere but does not refute evidence for the lid effect.

Alternative Estimates of Lithospheric Thickness
As introduced in Section 2.2, a limitation of theoretical cooling models is that they cannot capture local deviations in lithospheric thickness away from the average value for ocean floor of a given age (e.g., D. R. Davies et al., 2019).By comparing expected values with local estimates obtained from seismological constraints, we demonstrate that there are systematic differences between the two at the sites of OIBs in our database (Figure 14).For ocean islands on lithosphere younger than ∼30 Ma, such as at Easter Island or the Azores, seismically inferred estimates of local lithospheric thickness systematically exceed expectations from plate-cooling models.This offset is likely artificial, being a consequence of surface wave tomography having limited resolution at depths shallower than ∼75 km and therefore smearing shallow velocity structure into greater depths in regions of thin lithosphere (see Section 2.2).For older lithosphere on the other hand, seismically inferred estimates of presentday lithospheric thickness beneath each ocean island are consistently thinner than expectations from plate models.
In contrast to the artifacts in regions of thin lithosphere, this observation is likely real and probably reflects destabilization and thinning of the lithosphere by the underlying mantle plume (e.g., G. F. Davies, 1994;Dumoulin et al., 2001).As a consequence of these two effects, the majority of lithospheric thickness estimates beneath OIBs from the seismic model fall in the 40-100 km range, which is slightly narrower than the associated range of 30-120 km from plate-cooling models.
When switching to estimates of lithospheric thickness inferred from the seismic model, we find that all geochemical trends are best fitted by bi-linear models, including those that display linear trends when using lithospheric thickness from basin-specific plate models (Figure S4f in Supporting Information S1).In particular, there is a moderate preference (5 < log 10 E 2 log 10 E 1 < 10) for bi-linear models in the cases of Al 2 O 3 , FeO, and Yb, as well as a slight preference (2 < log 10 E 2 log 10 E 1 < 5) in the cases of SiO 2 and Lu.Nevertheless, constant models still perform poorly (Figure 13e), and our data set displays robust evidence for existence of the lid effect.

Processes Contributing to Observed Geochemical Trends
In our reference setup, bi-linear trends generally fall into two categories: (a) the highly incompatible elements (Th, K, La, P, Ti, Na), in which concentrations are entirely controlled by the degree of melting and therefore indicate that F remains approximately constant at lithospheric thicknesses greater than that of the kink; and (b) for CaO and parameters describing REE patterns (λ 0 and λ 1 ), where the depth of the kink is being determined by the combined effects of changing F interacting with the spinel-garnet phase transition.For the latter, the role of this phase change is to maintain HREE concentrations with increasing lithospheric thickness, while allowing the concentration of the LREE to continue to increase as F reduces.
It is important to note that the kinks in our bi-linear trends may not, in reality, reflect sharp change points, but rather a gradual transition in the trend as a function of increasing lithospheric thickness.While our findings show that the kink is generally identified at a depth of 50-60 km (Figure S3 in Supporting Information S1), since these trends are sensitive to both variations in F and the spinel-garnet phase transition, it is incorrect to infer the phase transition depth directly from the kink depth.This point is further emphasized by the fact that we observe garnet signatures in some OIBs that are generated beneath thin lithosphere (e.g., through trends of Yb and Lu), which can be attributed to the memory effect of high-pressure melts from the garnet stability field incompletely mixing with lower pressure melts from the spinel stability field.
Our resulting inferences on the degree of melting as a function of lithospheric thickness (Figure 12a) suggest that, beyond a certain lithospheric thickness, F becomes approximately constant.This behavior is unexpected since, based on the lid effect and theoretical models of plate cooling, we would expect lithosphere to continue to thicken and, all other aspects being equal, cause a continuous reduction in F. One potential explanation for this behavior could be progressive thinning of overlying lithosphere by upwelling plume material.Small-scale convection above mantle plumes is known to be more prevalent beneath thicker lithosphere (e.g., Ballmer et al., 2011;D. R. Davies et al., 2016;Dumoulin et al., 2001;Duvernay et al., 2021;Le Voci et al., 2014;van Hunen et al., 2003), making it more likely that the base of older lithosphere would become unstable upon plume impingement.This argument is supported by our observation in Section 4.5.4 and Figure 14 that seismically inferred estimates of lithospheric thickness are consistently thinner than those predicted by plate-cooling models in older lithosphere.Accordingly, beyond the ∼55 km kink depth, lithospheric thickness above mantle plumes is unlikely to increase at a rate consistent with cooling model expectations, thereby reducing the rate of expected reduction in the degree of melting.A further contributing factor is that, since the solidus temperature increases with pressure (Figure 12b), weaker plumes with lower excess temperatures may fail to cross the solidus and generate melt beneath thick lithosphere.This effect would be compounded by the fact that weaker plumes generate smaller melt volumes that are more likely to get trapped at depth and not erupt onto the seafloor.Variations in lithospheric thicknesses are insufficient to fully account for observed trends, thus we believe that temperature is an essential additional factor in shaping OIB composition (as also advocated by e.g., Lustrino et al., 2022).We refer to this behavior as the "temperature effect" and have investigated two lines of independent evidence that might support it.
First, we have explored potential relationships between lithospheric thickness and the potential temperature of OIB sources as estimated from geochemical or geophysical arguments (e.g., Ball, White, et al., 2021;Bao et al., 2022;Putirka, 2008, Figure S11 in Supporting Information S1).No clear patterns have emerged (although such estimates are known to be uncertain; e.g., Bao et al., 2022;Herzberg et al., 2007).Second, we have compared lithospheric thickness to recent analyses of plume buoyancy flux from Hoggard, Parnell-Turner, and White (2020).Here, we find that magmatic plumes beneath thicker lithosphere generally have higher buoyancy fluxes, potentially indicative of higher excess temperatures (Figure 15).This observation is consistent with our suggestion that, beyond the kink depth, degrees of melting are approximately constant due to preferential sampling of progressively hotter plumes from regions of thicker lithosphere (i.e., the rate of decrease in F is at least partially offset by the increase in plume temperature).
Taken together, local variations in lithospheric thickness away from average expectations from theoretical cooling models, sampling biases associated with progressively hotter plumes in regions of thicker lithosphere, and source region heterogeneities, are all plausible contributors to observed incompatible element trends.

Limited Evidence for Melt Re-Equilibration at Base of Lithosphere
Both Iceland and Hawaii have a large number of samples and exhibit a wide spread of compositions (e.g., Figure 16).Previous studies have attributed these ranges to variations in the fertility of the mantle source (e.g., Humphreys & Niu, 2009;Jones et al., 2017;Niu et al., 2011).Furthermore, Niu (2021) have also suggested that none of this spread can be attributed to differences in the initial melting pressure (i.e., there is no memory effect), since OIB melts re-equilibrate with the surrounding mantle during their ascent to the surface.Melting in an ascending mantle plume is expected to occur over a depth range of several tens of kilometers.If such re-equilibration reactions  do occur, however, we would expect major element concentrations buffered by olivine and pyroxene to be strongly homogenized, while highly incompatible trace elements (e.g., Th, K, and La) that have D′ < 0.01 in both the garnet and spinel stability fields should retain their original spread.
Our analyses have found no evidence to support such a process for re-equilibration of plume-derived melts: in other words, we find robust evidence for preservation of geochemical signatures across a range of depths in erupted melt products (i.e., the memory effect).OIBs from Hawaii and Iceland, for example, show a negative correlation between SiO 2 and FeOT (Figures 16a and 16b), which can be attributed to melts generated at a range of different pressures and has previously been suggested to occur in many OIBs (e.g., Scarrow & Cox, 1995).In the case of Iceland, most of the data cluster at ∼46.5 wt.% SiO 2 , but some samples extend toward ∼42 wt.% SiO 2 .The high SiO 2 samples could relate to melts generated at low pressure, while samples with lower SiO 2 are generated by smaller degrees of melting at higher pressure.There is also a cluster of samples at the high SiO 2 end of the Hawaiian array, albeit at 48.5 wt.% SiO 2 , with the data more evenly spread across the array.As with Iceland, we suggest that the low-SiO 2 , high-FeOT basalts at Hawaii are produced by melts separating from the mantle at depth (i.e., far below the lithospheric lid) and that they have subsequently erupted without reequilibrating during their ascent.
Importantly, whilst variations in FeOT at a given SiO 2 can potentially be explained by mantle source heterogeneity, the correlation between SiO 2 and λ 1 provides convincing support for melting across a range of pressures without complete homogenization and mixing (Figures 16c and 16d).Melt re-equilibration would be expected to bound major element concentrations within a limited range, but have little impact on incompatible trace element concentrations.Accordingly, following re-equilibration, limited correlation would be expected between major and trace elements, which is not borne out by our observations.

Conclusions
Our study yields insights into the role of lithospheric thickness variations in influencing the geochemical characteristics of OIBs.Our results support existence of the lid effect, in which lithospheric thickness limits the lowest melting pressure of upwelling mantle plumes and has an important influence on OIB geochemistry.Our statistical analyses suggest that REE patterns, and major and trace element concentrations, are influenced by lithospheric thickness, with some geochemical parameters best fitted by linear trends and others by bi-linear trends with a kink at thicknesses of 50-60 km.Although other factors such as source heterogeneity, melts separating from the mantle at various depths below the lid, and bias from heavily sampled localities are accepted to influence OIB geochemical trends, the observed trends remain overall consistent with expectations from the lid effect.Such trends can be explained by a combination of pressure-driven changes in the degree of melting and mineral assemblage, especially the spinel-garnet transition.The behavior of highly incompatible elements suggests that the degree of melting decreases rapidly with increasing lithospheric thickness until thicknesses reach ∼55 km, but subsequently decreases at a significantly lower rate with increasing thicknesses.This behavior is inconsistent with theoretical expectations based solely on the lid effect and suggests that: (a) plumes impinging beneath thicker lithosphere may be more effective at thinning the overlying lid, thereby modulating changes in the degree of melting; and (b) only melts from plumes with higher potential temperatures can penetrate thick lithosphere and reach the seafloor, consistent with solidus temperatures increasing with pressure and evidence that magmatic plumes under thicker lithosphere have higher buoyancy fluxes.Although our explanation may not be the sole contributor to the observed geochemical trends, with other factors likely to play a role, our results demonstrate an association between OIB geochemistry and lithospheric thickness.Quantitatively assessing the relative roles of other contributions requires further research.
The depth of the spinel-garnet transition zone cannot be directly identified from the trends observed herein.Nonetheless, the signature of this phase transition is evident in observed trends for Yb and Lu.These trends require that a signature of melt produced within the garnet zone is carried by many OIBs originating beneath thin lithosphere, indicative of a memory effect within plume-derived melts.This interpretation is further supported by geochemical trends from different samples generated in the same plume: it is therefore likely that OIB melts, generated at varying pressures, can ascend to the surface without mixing or re-equilibrating at the base of the lithosphere.
Taken together, our results have implications for magma generation, migration and mixing beneath OIBs, which will be vital for connecting these intricate processes to the larger-scale dynamics of upwelling mantle plumes.Our study provides quantitative constraints on the relationship between modern OIB geochemistry and lithospheric thickness, which will underpin future efforts to invert the geochemical composition of volcanic lavas for the temperature and pressure of their mantle source (e.g., Ball, Czarnota, et al., 2021;Klöcking et al., 2020), including those preserved from earlier periods of Earth's history, revealing changes in lithospheric thickness through space and time.Such point-wise constraints are also required by a new class of data-driven geodynamical simulation that aim to recover the spatial and temporal evolution of the mantle and its impact at the surface (e.g., Bunge et al., 2023;Ghelichkhan et al., 2023).This research was undertaken with the assistance of resources from the National Computational Infrastructure (NCI Australia), an NCRIS enabled capability supported by the Australian Government, and was partially supported by the Australian Government through the Australian Research Council's Discovery Projects funding scheme (project DP200100053) and Discovery Early Career Researcher Awards (DE220101519).The authors are also grateful for funding provided by the Chinese Scholarship Council (CSC) scholarship and Shen-su Sun scholarship.We thank Anthony Burnham, Thomas Duvernay, Hugh O'Neill and Fred Richards for help and feedback.We thank Boda Liu and Antonio Manjon-Cabeza Cordoba for constructive and valuable reviews, alongside the editor, Paul Asimow.Open access publishing facilitated by Australian National University, as part of the Wiley -Australian National University agreement via the Council of Australian University Librarians.

Figure 2 .
Figure 2. (a) Present-day oceanic lithospheric age from Seton et al. (2020) with locations of selected ocean island basalt samples (black dots).Only the name of the archipelago for each island group is displayed, but each individual island's lithospheric age and thickness are considered separately during the analysis.(b) Present-day oceanic lithospheric thickness based on a global plate-cooling model(Richards, Hoggard, Crosby, et al., 2020).(c) Present-day oceanic lithospheric thickness constrained by surface-wave tomography(Hoggard, Czarnota, et al., 2020).

Figure 3 .
Figure 3. Ocean island basalt database and sample filtering criteria.(a) SiO 2 versus MgO; colored dots = original samples colored by Gaussian kernel density estimation, normalized from 0 to 1; dashed lines = filtering criteria corresponding to SiO 2 43-54 wt.% and MgO 7-16 wt.%; white circles = subset of data that pass all filtering criteria.(b) Same for TiO 2 versus MgO.(c) Same for Nb/U versus La/Nb, where criteria of >30 and <1.2, respectively, are applied.(d) Same for La/Ba versus La/Nb, where only samples inside the ellipse of Fitton et al. (1991) are accepted.

Figure 4 .
Figure 4. Schematic cartoon showing (a) a constant model with one unknown parameter; (b) a linear model with two unknown parameters; and (c) a bi-linear model with four unknown parameters.Associated model variables are labeled.

Figure 5 .
Figure 5. Concentration histograms of (a) MgO; (b) SiO 2 ; (c) TiO 2 ; and (d) Al 2 O 3 in our ocean island basalt (OIB) data set.Original, unfiltered data are colored blue, while data in green represent the subset of data remaining following application of screening filters outlined in Section 2.1.2.For simplicity, histograms of (e) λ 0 , (f) λ 1 , and (g) λ 2 values are shown only for filtered OIB samples.

Figure 6 .
Figure 6.Statistical evidence evaluation results for Al 2 O 3 and λ 1 under our reference setup.(a) Al 2 O 3 as a function of lithospheric thickness fitted using a constant model; black crosses = individual samples; blue shading = probability density; yellow line = mean model; red dotted lines = expected spinel-garnet transition depths for typical mantle potential temperatures expected in plumes (e.g., Klemme & O'Neill, 2000; Robinson & Wood, 1998; Tomlinson & Holland, 2021); inset gives log 10 evidence value.(b) Same for a linear model.(c) Same for a bi-linear model; gray band = prior distribution for kink depth with one standard deviation width.(d) Probability distribution of kink depths, shown as a green histogram; gray band = prior; yellow line = mean value.(e-h) Same as (a-d), albeit for λ 1 .

Figure 7 .
Figure7.Optimal model type for each geochemical parameter under our reference setup.Ticks denote optimal model; strength of color fill indicates level of preference for that model type (i.e., when a simpler model has an evidence value that is within 20 but less than 2 of the optimal model, it is filled with color that linearly increases in intensity).

Figure 9 .
Figure 9. Statistical evidence evaluation results for geochemical parameters optimally fitted by bi-linear models.Data and panel contents same as for Figure 8 but for (a) Th, (b) La, (c) TiO 2 , (d) P 2 O 5 , (e) K 2 O, (f) Na 2 O, and (g) CaO.The horizontal bar below panel (d) shows the probability distribution of the likely kink depth, with more opaque colors indicating that a kink is more likely at that depth.

Figure 10 .
Figure 10.Statistical evidence evaluation results for heavy rare earth elements, which are optimally fitted by linear models.Data and panel contents same as for Figure 9 but for (a) Yb and (b) Lu.

Figure 11 .
Figure 11.Statistical evidence evaluation results for Rare Earth Element shape parameters.Data and panel contents same as for Figure 9, but for (a) λ 0 , (b) λ 1 and (c) λ 2 .Note that λ 0 and λ 1 are optimally fitted by bi-linear models, whereas λ 2 is optimally fitted by a linear model.

Figure 12 .
Figure 12.(a) Solid line = degree of melting (F) as a function of lithospheric thickness inferred from Equation10and the bilinear trends for La concentrations, with basin-specific plate model-derived estimates of lithospheric thickness; dashed lines = theoretical degree of melting for decompression melting of dry, primitive peridotite at different potential temperatures from the parameterization ofKatz et al. (2003), as modified byBall et al. (2022); dash-dotted line = same for a wet 1400°C source with H 2 O = 500 ppm, which is thought to be an upper bound for water content in plume source regions (e.g.,Asimow & Langmuir, 2003;Wallace, 1998).(b) Solidi for peridotite with 0-500 ppm water contents.

Figure 13 .
Figure 13.(a) Optimal model type for each geochemical parameter under our reference setup (filtered ocean island basalt samples from all localities corrected for fractional crystallization, with lithospheric thicknesses from basin-specific platecooling models; as in Figure 7).(b) Same as the reference setup except data are not corrected for fractional crystallization.(c) Same as the reference setup except 40% of the data are randomly removed.(d) Same as the reference setup except Icelandic and Hawaiian samples were excluded.(e) Same as the reference setup albeit with lithospheric thickness taken from the model based on seismic tomography.

Figure 14 .
Figure 14.Difference between local lithospheric thickness beneath each island from the seismic model and that predicted by the basin-specific plate models, as a function of lithospheric age at time of ocean island basalt eruption.Localities in the Atlantic, Pacific and Indian Oceans are represented by green, blue and red circles, respectively.

Figure 15 .
Figure 15.Relationship between buoyancy flux of magmatic plumes from Hoggard, Parnell-Turner, and White (2020) and lithospheric thickness estimated from the basin-specific plate models (for values greater than the kink depth of ∼55 km).

Figure 16 .
Figure 16.Co-variation of pressure-sensitive geochemical parameters for ocean island basalt samples located at Iceland and Hawaii.(a) FeOT as a function of SiO 2 for Icelandic samples, following application of filtering and corrections for fractional crystallization.(b) Same for Hawaiian samples.(c-d) Same for λ 1 values as a function of SiO 2 .

Table 1
Partition Coefficients for Incompatible Elements in the Main Peridotite Minerals