Quantifying the Tetrad Effect, Shape Components, and Ce–Eu–Gd Anomalies in Rare Earth Element Patterns

Plots of chondrite-normalised rare earth element (REE) patterns often appear as smooth curves. These curves can be decomposed into orthogonal polynomial functions (shape components), each of which captures a feature of the total pattern. The coefficients of these components (known as the lambda coefficients—λ\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\lambda $$\end{document}) can be derived using least-squares fitting, allowing quantitative description of REE patterns and dimension reduction of parameters required for this. The tetrad effect is similarly quantified using least-squares fitting of shape components to data, resulting in the tetrad coefficients (τ\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\tau $$\end{document}). Our method allows fitting of all four tetrad coefficients together with tetrad-independent λ\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\lambda $$\end{document} curvature. We describe the mathematical derivation of the method and two tools to apply the method: the online interactive application BLambdaR, and the Python package pyrolite. We show several case studies that explore aspects of the method, its treatment of redox-anomalous REE, and possible pitfalls and considerations in its use.


Introduction
Rare earth element (REE) patterns are plots of measured REE abundances in a sample, in which the concentration of each element is divided by its corresponding concentration in a reference material, a process termed normalisation (Coryell et al. 1963). The most common reference material is an average of ordinary chondrite meteorites (CI), often taken to represent an unmodified primordial composition of the solar system and the Earth. Once chondrite-normalised (CN), any deviations from a horizontal array record some geochemical fractionation process. As the geochemical properties of the REE generally vary in a smooth manner from La to Lu (Anenburg 2020 and references therein), the pattern itself is a smooth line whose shape is informative of the processes that led to the sample formation.
Quantification of various shape aspects of patterns is usually achieved using normalised element ratios. For example, La N /Lu N (where subscript N indicates normalised values) measures the overall slope of a pattern, and La N /Sm N might provide information on the degree of light REE (LREE) fractionation from the middle REE (MREE). However, the choice of elements is often arbitrary and chosen to support a certain hypothesis, or designed for use in simpler environments (e.g., mid-ocean ridge basalts), and can never provide information on the shape on a shorter wavelength than the one used for the ratio (Anenburg 2020). For example, the three different REE patterns shown in Fig. 1 (A, B, C) are clearly different and represent different petrogenetic processes, yet their La N /Lu N are nearly identical. Additional measures, such as Dy/Dy* (the deviation of Dy from a straight line connecting La and Yb) that attempts to quantify pattern concavity, are only reliable for simple patterns (Davidson et al. 2013). For example, Dy/Dy*= 1 for pattern D in Fig. 1, suggesting a simple straight pattern. Yet, a quick glance at the pattern reveals it to be an artefact of its curvature, in which Dy is located precisely at the La-Yb interpolating line, simply by chance.
In cases of such complexity, descriptions of REE pattern are often supplemented by qualitative terms such as "sinusoidal", "U-shaped", or "spoon-shaped". As these are loosely defined qualitative terms, they are somewhat ambiguous (i.e., one person might understand "spoon-shaped" differently from another) and unquantifiable (to what degree is a pattern spoon-shaped?). The problem is compounded when two or more qualitative shape descriptions are combined. Using normalised element ratios to describe an overall sloping REE pattern which is simultaneously spoon-shaped and sinusoidal becomes a futile endeavour. In such cases a simple visual inspection of the pattern is more informative, as "a picture is worth a thousand words". However, visual inspections become challenging when many analyses are present and the problem of overplotting becomes serious as many individual lines are entangled into an incomprehensible jumble, from which it is not straightforward to discern trends. For instance, trends in a data set representing a rock suite in which pattern sinusoidalities  (4.6, −15, 225, 7800, 24000). The linear line connecting La and Yb in pattern D is plotted in light green increase whereas spoon-shapednesses decrease are likely to be obscured by having "a thousand pictures".
The above-described difficulties are greatly reduced if REE patterns are taken as mathematical functions constructed from simpler functions, each describing a shape, which we term "shape components". The method was originally developed by O'Neill (2016) and employs least-squares fitting of coefficients to orthogonal polynomials of increasing powers, termed the "lambda shape coefficients". Here, we explicate the method of using λ shape components and coefficients to describe REE patterns, and expand it to include τ coefficients that describe the tetrad effect-subtle variation in REE behaviours that affects individual groups of four consecutive elements.

Lambda Shape Coefficients
The function used by O'Neill (2016) to describe REE pattern shapes is where ln [REE i ] N is the natural logarithm of the normalised element concentration. REE plots are plotted using the common logarithm (i.e., log 10 [REE i ] N ), which makes the choice of natural logarithm seem odd; the reasoning for the use of the natural logarithm is presented in O'Neill (2016). In any case, transformations between the two logarithm bases are trivial (see appendix in Anenburg 2020). In Eq. 1, f orth n are orthogonal polynomial functions of eightfold ionic REE radii r i (Shannon 1976). The functions of increasing order are as follows The parameters β, γ n , δ n , and n are pre-calculated constants, such that each f orth n describes a fixed, predetermined shape (see Appendix). The use of ionic radii is preferred to atomic number, as it widens the distance between each LREE relative to HREE. First, each LREE often varies in its properties from its neighbour more than the HREE: in most natural samples, the fractionation degree of Ce differs from Pr more than Ho differs from Tm. This behaviour follows that of ionic radii, and by using them, the greater fractionation degree is accommodated by lower-order shape components. Second, O'Neill (2016) empirically demonstrated that using ionic radii results in a superior fit compared to using atomic numbers. The drawback is that the mathematical treatment becomes slightly more complex, but this is not an issue in nearly all applications if the available software described below is made use of. We follow the convention used by O'Neill (2016) of using the eightfold coordinated "ionic radius" from Shannon (1976), itself based on an earlier structure determination of simple lanthanide fluorides (Greis and Petzel 1974). Whether these values are sufficiently accurate and admissible as representative values for the entire field of geochemistry is still debatable, while being mindful that the "ionic radius" of an element varies based on the local environment and that ions are not rigid balls (Liu et al. 2013;Gibbs et al. 2013). There is ongoing work attempting to refine REE ionic radius values (Gagné 2018), and it is possible that slightly different radii will be used in future versions of the tools presented here.
A complete REE pattern is thus described by the linear combination of the shape coefficients (λ n ) acting as scaling factors for their respective shape components ( f orth n ). Specifically, λ 0 represents the average REE abundance, λ 1 represents the linear slope, λ 2 represents quadratic curvature, λ 3 represents cubic curvature (i.e., sinusoidality), and λ 4 represents the quartic curvature (i.e., W-shape), essentially providing finer control that, when combined with the previous coefficients, allows for various spoon shapes and the like. Figure 2 shows an example of a REE pattern constructed by the gradual addition of shape coefficients. An interactive online app, ALambdaR (https:// lambdar.rses.anu.edu.au), allows easy visualisation of the effect of shape coefficients on REE patterns, and readers are encouraged to experiment with it in order to gain an intuitive understanding of the method.

Background and Previous Studies
In addition to the smooth variation in REE geochemical properties across the lanthanide row, there are cases where the REE are divided into four segments ("tetrads"), each consisting of four elements: (i) La-Ce-Pr-Nd, (ii) Pm-Sm-Eu-Gd, (iii) Gd-Tb-Dy-Ho, and (iv) Er-Tm-Yb-Lu (Kawabe 1992). Within each tetrad, the middle two elements are fractionated relative to the elements at the edges (e.g., Tb and Dy are fractionated relative to Gd and Ho). The fractionation can be positive or negative, which is manifested as a series of ridges or troughs on a REE pattern, occasionally called M-type and W-type, respectively. Various combinations of these types can also be observed . The tetrad effect stems from quantum mechanical effects on the electronic structure and energies of the 4f subshell (Kawabe 1992). Note that these M-type and W-type tetrads are distinct from the quartic f orth 4 of a smooth pattern described by λ 4 , often referred to as M-shaped or W-shaped.
Previous attempts to quantify the tetrad effect included the ratio of the middle elements from their values on an interpolated line connecting the edge elements (Irber 1999), or calculation of the standard deviation of this ratio (Monecke et al. 2002;Abedini et al. 2020). These methods suffer from three major problems. Firstly, the magnitude of the tetrad effect can only be calculated if all four elements are measured and do not present anomalies. This eliminates tetrad 2, as it includes Pm, which does not occur naturally, and Eu, which is commonly anomalous because of its two oxidation states (Eu 2+ and Eu 3+ ). This can also lead to elimination of tetrad 1, as Ce is also occasionally anomalous because of its two oxidation states, Ce 3+ and Ce 4+ , with the latter being common in low-temperature environments; this is especially problematic because the tetrad effect is most likely to occur in surface environments and highly evolved magmatic systems (Bau 1997) up to the low hundreds of degrees Celsius. Secondly, these methods are not suitable for REE patterns which exhibit curvature, namely the vast majority. These interpolation methods will report non-zero tetrads even if the pattern had not been influenced by tetrad effects, simply because curved lines deviate from straight lines (Bau 1997), as these methods essentially calculate the degree of precisely that deviation. Thirdly, tetrads are not symmetric. They arise from the electronic structure of the elements and the stability of empty (La), half-full (Gd), and full (Lu) f shells (Kawabe 1992;Peppard et al. 1969). Additionally, quarter-full f shells should also contribute to the tetrad effect, but the number of f orbitals (14) is not divisible by four, and true quarter-filled shells will only be realised in non-existing elements that occur between Nd and Pm, and Ho and Er (Nugent 1970). This mismatch leads to asymmetry which is not captured by simple linear interpolation on a normalised REE pattern plot. Finally, the two methods based on standard deviations (Monecke et al. 2002;Abedini et al. 2020) may erroneously quantify analytical noise as a tetrad effect if the middle elements are shifted to the same direction.
A superior method for tetrad effect quantification was developed by Minami and Masuda (1997). It involves least-squares fitting of four quadratic curves to the four tetrads, mathematically constraining the positions of the curves such that they intersect correctly at the quarter-filled non-existent element positions using the method of Lagrange undetermined multipliers. The method also solves the problem of missing elements-they are simply excluded from the fitting procedure, using the remaining elements instead. Finally, the method partially solves the issue of curvature by dividing the entire REE pattern into four linear slopes unto which parabolas are superimposed. Unfortunately, the method has largely remained unadopted by the geochemical community, likely due to its mathematical complexity and the lack of accessible computational tools to apply it.

Tetrad Functions Compatible with λ Shape Coefficients
A major limitation for tetrad fitting is that the method of Minami and Masuda (1997) cannot be easily applied in conjunction with the REE shape coefficients described above, as the least-squares matrix construction is different. Minami and Masuda (1997) only use a subset of the REE for the fitting of each tetrad, whereas λ fitting considers all elements simultaneously. Furthermore, as the tetrad fitting already includes some degree of curvature (albeit in low resolution), there is a conflict between that and some of the λ coefficients. Here we adapt their method so that it is compatible with the O'Neill (2016) λ shape coefficients, and allows compatible fitting of both shape and tetrad coefficients.
Each tetrad is initially represented by a negative parabola that rises above the x-axis at its respective location and is zero elsewhere, essentially the function Four such parabolas have to be generated, each shifted laterally so that their x-axis intersections occur at the following positions on an atomic number axis: (i) La and Nd-Pm midpoint: 57 and 60.5, (ii) Nd-Pm midpoint and Gd: 60.5 and 64, (iii) Gd and Ho-Er midpoint: 64 and 67.5, and (iv) Ho-Er midpoint and Lu: 67.5 and 71. Their maxima are located between the intersections at 58.75, 62.25, 65.75, and 69.25, respectively. Consequently, 1−x 2 is transformed to form all four tetrads via f tet (g(z)) where g(z) = (z−z 0 ) 1.75 , z 0 are the maxima positions, and 1.75 is half of the parabola width (3.5/2) at f tet = 0. Finally, as the λ method expects ionic radii rather than atomic number, a third function that transforms r to z is required. Previous attempts to parameterise r as a function of z have found agreement with Slater's model for atomic shielding (Raymond et al. 2010;Seitz et al. 2007), but here we use a simple polynomial function h(r ) = z that maps the Shannon (1976) radii to atomic number. The polynomial order can be arbitrarily large, and we find that a seventh-order polynomial is sufficient to accurately generate parabolas at the correct locations when the parabolas are calculated using f tet (g (h (r ))).
Each tetrad is then described by τ n f tet n , where τ n are each of the four tetrad coefficients, with positive values indicating a negative parabola (i.e., a ridge, or M-type tetrad), and negative values indicating a positive parabola (i.e., a trough, or W-type tetrad). Dark blue lines are the smooth patterns derived from λ n terms, light blue lines are the tetrads derived from τ n terms, light green lines are the quadratic curves that underlie the tetrad model, and purple lines are the complete REE patterns (i.e., + T)

Combining λ and τ Shape Coefficients
The sum of all functions described so far leads to the complete function that describes REE patterns Two examples of complete REE patterns that include both λ and τ terms are given in Fig. 3.
Fitting data to this form is a straightforward procedure in most modern data analysis software packages. There are currently two tools available for fitting a function of this form (Eq. 4) to REE patterns. The first is BLambdaR, an interactive online tool written in the Shiny package for R, using R's native linear algebra routines (Chang et al. 2020; R Core Team 2020; Anderson et al. 1999;Wilkinson and Rogers 1973;Chambers 1992). It is accessible at https://lambdar.rses.anu.edu.au/blambdar/. The second is pyrolite, a Python package which includes functions for fitting lambdas and tetrads amongst a suite of tools tailored to transformation and visualisation of geochemical data (Williams et al. 2020).

Treatment of Uncertainties
Naturally, the fitting of often noisy real-world data using λ and τ coefficients is never perfect. Several numerical indicators can be used to assess the fit quality. The first indicator is the coefficient of determination (r 2 ), a number ranging from 0 to 1 which is commonly interpreted as the proportion of variability in the data explained by the model, with higher r 2 values indicating lower unexplained variance. Whilst useful for linear fits, r 2 has the problem that it increases for every fitted parameter, even if said parameter does not actually improve the fit. As our model may include up to nine parameters (five λ and four τ ), this may be an issue of concern. Therefore, we use the adjusted r 2 where n is the number of fitted REE and p is the number of fitted coefficients. Note that the denominator often includes a −1 term in the statistical literature to account for an intercept parameter, but in our case this parameter (λ 0 ) is included in p, rendering decrementation by one unnecessary. r 2 adj attempts to account for the potentially spurious increase of r 2 when many parameters are fitted, and provide a less biased estimate of unexplained variance (Freund et al. 2010).
The second indicator, which measures the goodness of fit, is the reduced chi-squared where ν indicates degrees of freedom (i.e., fitted variables minus fitted parameters), the numerator is the sum of squared residuals, the hat symbol indicates expected value from fit, and s indicates the assumed uncertainty on REE measurements as a fraction (i.e., 5% uncertainty should be taken as s = 0.05). See Bevington and Robinson (2003) and the Appendix in O'Neill (2016) for more details. In general, lower χ 2 ν is better. However, given that we know that there are errors in the data because analytical methods are not foolproof and unavoidable uncertainties exist-s(ln[REE]) * 2 > 0it is very unlikely for χ 2 ν to be too low, as that would imply much lower uncertainties in the REE data than is possible according to instrumental and other sources of error. χ 2 ν is then an indicator of this issue: χ 2 ν = 1 indicates that the fitted residuals of the model are in accordance with the assumed data uncertainty, χ 2 ν > 1 indicates that the model does not do a good job of fitting or that we have underestimated data uncertainties, and χ 2 ν < 1 indicates that we are overfitting the model (i.e., it fits a curve with no actual significance) or that we have overestimated data uncertainties. Geochronologists and other geoscientists might be more familiar with an alternative term for χ 2 ν , the mean squared weighted deviation, or MSWD (Powell et al. 2002). Individual χ 2 ν values may not be particularly informative themselves, as a single measurement may be exceptionally accurate or inaccurate just by chance. However, when integrated along a large database, all individual χ 2 ν values are themselves χ 2 -distributed with a theoretical maximum that corresponds to the actual s, allowing one to estimate an overall representative uncertainty on REE data across the database (O'Neill 2016).
We also include statistics on individual parameters: the standard error (SE) and p-value. The SE is calculated by taking the square root of the diagonal of the variancecovariance matrix derived from the measured data (Freund et al. 2010;Bevington and Robinson 2003). The ratio of the coefficients and their respective SE is taken, from which a double-sided probability of that value or larger under a t distribution with ν degrees of freedom is computed as the p-value. Better fits result in lower pvalues. In hypothesis testing, p-values smaller than 0.05 are often taken as a rubber stamp for statistical significance, or whether there is a real influence by a certain fitted parameter. We would caution against such a simplistic interpretation in this case, as no hypothesis testing is conducted per se, and the physical interpretation of the shape coefficients is debatable (see below). There is an ongoing discussion on the use and abuse of p-values (Amrhein et al. 2017;McShane et al. 2019;Greenland 2019), and we recommend that they are taken as a semi-quantitative indicator of whether a certain pattern exhibits a certain shape. Values as low as 0.0001 can safely be interpreted as "the shape definitely exists", and values as high as 0.3 may require more scrutiny. This is particularly important for tetrads: if they cannot be consistently and significantly fitted, it is best to leave them out of the model, as their inclusion breaks λ orthogonality, as discussed in the next section.

Orthogonality and Lack of Thereof During Tetrad Fitting
One principal strength of the λ shape parameters and their respective shape components is predicated upon their orthogonality. All components are independent of each other. For instance, bending of a linearly sloping pattern in a fashion consistent with f orth 2 would increase the magnitude of λ 2 with no influence on λ 1 . Another striking example, which we encourage readers to attempt on their own using the online tools, is to generate REE data with = (0, 0, 0, 0, −20000) in ALambdaR so that the pattern consists essentially of an M-shaped f orth 4 . Next, feed the REE data to BLamb-daR, which correctly extracts the initial vector (permitting minor deviations due to rounding errors, etc.). Unticking the option to fit λ 4 leads BLambaR to counterintuitively report = (0, 0, 0, 0), as if the fitted pattern is a simple horizontal line. This behaviour provides an excellent demonstration of orthogonality: any f orth n polynomial is essentially invisible to all other orthogonal polynomials, and allows the method to be used consistently across the range of all naturally occurring REE patterns. A certain shape (e.g., quadratic curvature: λ 2 ) is intercomparable across all other patterns, regardless of the degree of their slope and sinusoidality, and other shapes described by λ =2 . The above-described property of orthogonality breaks down when tetrad components (τ n f tet n ) are added to the least-squares fitting. The quadratic f tet n functions are only defined in their respective four-element tetrads and are consequently independent of each other. That is, excluding or including a certain tetrad will not directly affect the other tetrads. However, they are not orthogonal with respect to the f orth n functions, and cannot be made so. The consequence is that when both λ and τ coefficients are sought in a single least-squares fit, the fitted coefficients do not follow the expected mutual orthogonal independence. This can be easily seen for the above pattern (λ 4 = −20000) by ticking all τ n while λ 4 remains unticked. The pattern contains no tetrads and no other λ-shapes-because that is how we prepared it, yet the fitting leads to the following spurious parameters: = (−0.115, −1.59, −12.1, 666) and T = (0.57, 0.043, −0.01, 0.15). The problem leads to a concern: how can the fitted parameters be trusted and interpreted when attempting to quantify the tetrad effect? An attempt to address these issues is given in the case studies below.

Physical Significance of λ Shape Coefficients
Accurately measured REE patterns represent actual fractionation processes experienced by the analysed sample. An underlying assumption for least-squares fitting is that real data are well represented by the fitted model. This leads to the question: are the predetermined f orth n shape functions meaningful, and do they directly represent actual natural processes? The physical significance of λ 0 is trivial, as it represents the overall average chondrite-normalised abundance of REE in a pattern and does not affect the shape. REE patterns essentially reflect REE partitioning and segregation between two or more phases. This partitioning is invariably related to crystal-chemical constraints, and is often modelled using the lattice strain model (Blundy and Wood 2003). Quadratic fits to measured REE partition coefficients in D − r i space are usually accurate (where D is the concentration ratio between two phases), particularly for minerals that contain a single crystallographic site that can host REE. High-quality partition data often reveal minor asymmetry, which is successfully modelled by an additional cubic r 3 i term. The physical significance of λ 1 and λ 2 is therefore a proxy for simple lattice strain partitioning, as they similarly describe quadratic curves in ratio −r i space. The parabola defined by λ 2 f orth 2 , on its own, has an extremum close to Eu (pattern A in Fig. 4), and it can be shifted horizontally by varying λ 0 and λ 1 (as any parabola defined by y(x) = ax 2 +bx +c can be transformed in x − y space by varying b and c), demonstrated by patterns B and C in Fig. 4. However, this is a simplification because the hypothetical patterns in Fig. 4 would have formed by fractionation of a single-REE-site mineral from an unfractionated CI chondrite pattern. Real REE patterns are unlikely to be that simple, because they are formed by fractionation of numerous minerals, each having more than one REE-hosting crystallographical site, from an alreadyfractionated source. In essence, a currently measured REE pattern reflects the entire lattice strain partitioning pedigree since the formation of the solar system. Modelling this history is theoretically possible, but arduous and most likely non-unique, leading to accurate reconstruction being unachievable. Nevertheless, fractionated natural patterns will inevitably be smoothly curved, leading to the use of λ 3 and λ 4 . There is arguably no direct physical significance to the cubic f orth 3 and quartic f orth 4 , but these higher-order polynomials represent fine-tuning and approximate small variations that stem from the cumulative parabolas, of varying heights, widths, and positions, encapsulated in a single REE pattern. The lack of physical significance does not, however, detract from the usefulness of the method, as it does an excellent job of describing various shape features of REE patterns. The method also has practical applications in petrogenetic modelling by parameterising REE crystal/melt partition coefficients as D = δ 0 + δ 1 f orth 1 + · · · and using these coefficients in petrogenetic process vectors ( ) and coefficients (ψ n ). See O'Neill (2016) for more details and examples.

Quantification of Cerium and Europium Anomalies
A common feature of REE patterns is the presence of anomalies: the deviation of the redox-sensitive elements Ce and Eu from the smooth pattern exhibited by all other elements. The presence, sign, and magnitude of Ce and Eu anomalies are useful indicators of various petrogenetic processes. For example, Ce and Eu anomalies in zircon are indicative of the magma oxidation state from which the zircon crystallised (Loucks et al. 2020 and references therein). Anomalies are typically quantified by taking the ratio of the observed normalised concentration of Ce and Eu to a calculated normalised concentration derived by the geometric mean of the two neighbouring elements, often termed Ce* and Eu* (geometrically being the linear interpolation between the two elements on a logarithmic plot): Ce/Ce * LI = √ La N × Pr N and Eu/Eu * LI = √ Sm N × Gd N , where subscript LI indicates linear interpolation. A

Fig. 5
A zircon-like REE pattern generated using the parameters = (4.6, −63, −545, −1400, 10500) in light blue, and the same pattern with Ce and Eu anomalies in dark blue. Lightest blue are linear lines connecting La and Pr, and Sm and Gd. The pattern is truncated at Tb for clarity, and the full pattern can be obtained by inputting the vector into ALambdaR major disadvantage of this method is the assumption of linearity, namely the question of whether the geometric mean indeed represents the true concentration of Ce and Eu had there been no redox effects. Whilst this method may give acceptable results for relatively straight patterns, it rapidly becomes inaccurate when patterns are strongly curved. For example, the REE pattern shown in Fig. 5 is a typical zircon-like pattern, with positive Ce and negative Eu anomalies. Ce anomalies in zircon are a useful petrogenetic indicator (Burnham and Berry 2014;Trail et al. 2015), and many approaches to zircon Ce anomaly quantification have been published in recent decades (Zhong et al. 2019). Visually, the difference between the interpolated line and the non-anomalous pattern looks negligible-the difference is smaller than the symbol size and the line plots on top of the circles. However, this is misleading, as the y-axis covers six orders of magnitude. We indicate the expected value of Ce from the smooth curve constructed from shape components as Ce * λτ . The pattern has been generated with a Ce/Ce * λτ = 40, but calculating the anomaly using the interpolated line results in Ce/Ce * LI = 47.5, an overestimation of about 19%. Likewise, Eu/Eu * λτ = 0.6 (the true Eu anomaly), whereas the calculated anomaly via interpolation is Eu/Eu * LI = 0.667 , an underestimation of about 11%.
Following the above example, more accurate anomaly quantifications are conducted by taking the ratio of measured Ce and Eu to that value expected from the leastsquares fit, instead of linear interpolation. Clearly, this requires that anomalous Ce and Eu are not used for the actual fitting. It should be noted that for this method to succeed, accurate values for all REE are crucial, particularly for elements at the extreme edges (i.e., La, Yb, and Lu) when calculating a Ce anomaly. In some materials, this proves difficult, with Lu having the lowest chondritic abundance of all the REE. Additionally, La is notoriously difficult to measure in minerals that reject LREE, most notably zircon (Zou et al. 2019;Burnham 2020). In such cases, we recommend a visual inspection of the measured pattern, fitted line, and statistical indicators as given by pyrolite and interactively by BLambdaR, to assess whether an anomaly determination is reliable.

Gadolinium Anomaly and Tetrad Effect in Seawater
REE compositions and patterns of seawater are a useful proxy for studying various ocean processes such as mixing of reservoirs, sediment load, hydrothermal vent input, etc. Our understanding of these processes might be improved if shape coefficients are used to describe seawater. As an example, we took the representative seawater data "BATS" of van de Flierdt et al. (2012) and fitted λ coefficients to it, treating Ce and Eu as anomalous. The match between the fitted line and the data is close, but not perfect (adjusted r 2 is 0.95 and 0.90 for shallow and deep waters, respectively, Fig 6a, Online Resource A1). Adding τ coefficients markedly improved the fit, with adjusted r 2 increasing to 0.997 and 0.985, respectively (Fig 6b, Online Resource A2).
It may be argued that the improved fit is least-squares trickery, because increasing the amount of parameters invariably leads to better fits and greater r 2 . However, there are several lines of evidence to suggest the effect is real. First, we are using the adjusted r 2 which accounts for spurious increases in determination caused by additional parameters. Second, all τ coefficients are consistently negative (i.e., reflecting a W -type tetrad effect) and showing roughly the same magnitude, except τ 1 which is known to be stronger in geological environments (Fig. 7). Had the additional τ coefficients described noise and spuriously improved the least-squares fit, this consistency would have been unlikely. Third, the uncertainty range of the τ coefficients is well away from zero (Fig. 7), and the τ significance p-values for deep water-consisting of more abundant REE, hence more likely to be accurately and precisely measured-are all under 0.012. For less accurate shallow water, the τ 1 p-value is 0.068, whereas the others are under 0.04.
Whether seawater has a tetrad effect or not has been previously debated (McLennan 1994;Kawabe et al. 1998), but our modelling shows that, when reliable REE data are normalised using reliable CI values, a significant tetrad effect is indeed observed, in agreement with recent observations (Ernst and Bau 2021). Seawater is commonly normalised to PAAS (Post-Archean Australian shale) rather than CI. Fitting τ coefficients to accurately measured PAAS data of Pourmand et al. (2012) . 7 Tetrad coefficients for seawater (BATS) and PAAS. Uncertainty shown is 1 standard error not contain any statistically significant tetrad effect, with p-values ranging from 0.24 to 0.96 (Fig. 7, Online Resource A3). Therefore, the choice of either CI or PAAS for normalisation is unlikely to generate, obscure, or erase tetrad effects. The tetrad fitting also shows that the minor positive Gd anomaly observed in Fig. 6 is the result of the tetrad effect and represents a natural phenomenon. A positive Gd anomaly in seawater, lakes, or the rivers that feed them is interpreted as anthropogenic contamination, and the anomaly is usually quantified using deviations of Gd from linear interpolation of other REE on a PAAS-normalised plot (Kulaksız and Bau 2011a), or by using nonorthogonal linear or polynomial least squares (Kulaksız and Bau 2011b;Müller et al. 2002). As these methods do not take overall pattern curvature or tetrad effects into consideration, subtle anomalies which derive from the tetrad effect might be erroneously ascribed to anthropogenic contamination. We recommend that least-squares fits to REE patterns using λ and τ parameters be examined to check for consistency with the linear interpolation or fitting methods.

Tetrad Effects in the Toongi Rare Metal Deposit
The Toongi Deposit in NSW, Australia, hosts substantial resources of Zr, Hf, Nb, Ta, REE, and Y. The mineralisation is hosted in eudialyte and various later alteration products such as bastnäsite, vlasovite, milarite, catapleiite, and gaidonnayite (Spandler and Morris 2016). These minerals often exhibit unusual and complex REE patterns, and Spandler and Morris (2016) qualitatively noted the presence of the tetrad effect in the patterns they obtained from the Toongi minerals. Figure 8 shows τ coefficients derived from the Spandler and Morris (2016) data (see Online Resource A4 and A5). As discussed above, there is a clash between τ and λ shape components, particularly when the patterns are highly curved at the edges, in which case the same features are fitted by both τ 1 , τ 4 , and λ 4 . To investigate the issue, we fitted the data with λ 4 (light green) and without λ 4 (dark blue). It is immediately evident that eudialyte has small positive tetrads, which are negligible compared to the later minerals. This is consistent with the inference of Spandler and Morris (2016) that it is a primary mineral, unaffected by lower-temperature hydrothermal processes. For minerals other than eudialyte, there is a clear difference between the λ 4 inclusive and exclusive values of τ 1 and τ 4 . Which values are correct? The two middle tetrads-τ 2 and τ 3 -are informative as they remain mostly unchanged irrespective of whether λ 4 is fitted or not. Therefore, we expect the external tetrads to follow the same behaviour, which is only achieved when λ 4 is not fitted. In contrast, τ 1 and τ 4 are diminished in magnitude and reversed in sign when λ 4 is included, a result inconsistent with the middle tetrads. Therefore, in this case of strongly curved LREE and HREE, we proffer that λ 4 is spuriously preventing the correct identification of the tetrads.
The LREE tetrads, τ 1 and τ 2 , are strongly positive in all secondary minerals. One possibility is that LREE in the Toongi secondary assemblage are dominated by bastnäsite and other REE carbonates. Although Spandler and Morris (2016) do not provide full REE data for these minerals, partial patterns based on La, Ce, Pr, Nd, and Sm suggest the presence of a negative τ 1 for REE carbonates. The secondary minerals shown in Fig. 8  The above discussion is not meant as a comprehensive study of REE partitioning and the tetrad effect in Toongi, but serves as an example of a range of possibilities made available by the use of quantitative λ and τ determinations. It also shows an example in which it is best to exclude λ 4 from the least-squares fitting.

λ 4 and Ce Anomalies in the Troodos Ophiolite
In our final case study, we examine REE patterns from oceanic plagiogranites collected from the Troodos ophiolite in Cyprus by Anenburg et al. (2015). A representative pattern is shown in Fig. 9. All τ coefficients are highly insignificant, regardless of whether λ 4 is included in the fit, and therefore there is patently no tetrad effect exhibited by the pattern (Online Resources A6 and A7), although by simple visual inspection of the pattern, it might appear as if there is a tetrad effect. This pattern is remarkable by how well it is described by λ 4 (Online Resource A8). Its magnitude is rather large (>10,000) compared to some of the lower-order coefficients (λ 2 ≈ 4 and λ 3 ≈ 200). λ 4 is often less precisely determined than the lower-order coefficients (O'Neill 2016), but in this case SE(λ 4 ) < 20%, whereas SE(λ 2 ) > 40% and SE(λ 3 ) > 20%.
The peculiar dominance of λ 4 can be explained as a combination of the three most REE-rich minerals that occur in the plagiogranites: hornblende, zircon, and epidoteallanite. REE patterns for the Troodos hornblende are characterised by two positively straight segments: a steep segment from La to Sm, and an almost horizontal segment from Gd to Lu (Gillis 2002), such that it would be well described by λ 1 < 0 and λ 2 < 0. Epidotes and allanites are LREE-rich (Anenburg et al. 2015), and zircon is HREErich (Chen et al. 2020), such that small amounts of these accessory minerals allow the lightest LREE and heaviest HREE to poke out of the hornblende-dominated pattern. This serves as an example of the lack of distinct physical significance for λ 4 . The central wide ridge records the maximum of hornblende curvature, and the two steep peripheral segments record allanite and zircon (Fig. 9). However, there is no reason to assume that λ 4 is particularly efficient in describing hornblende, allanite, or zircon patterns. It is just coincidental that the combinations of three lower-order f orth n and their lateral translations were accurately described by λ 4 .
The presumably coincidental yet excellent λ 4 fit raises an interesting question. What is the nature of the Ce anomaly in the pattern shown in Fig. 9? Is it a negative anomaly, as Ce plots below the La-Pr interpolation line (Ce/Ce * LI < 1), or is it positive, as Ce plots above the fitted line (Ce/Ce * λτ > 1)? The REE patterns of a rock sample reflect the REE patterns of its constituent minerals. The plagiogranites described by Anenburg et al. (2015) are LREE-poor, with Ce concentrations of around 5 ppm. This low Ce content is very sensitive to Ce-rich materials, such as zircon. Plagiogranite-hosted zircon from Troodos is exceptionally Ce-rich, with Ce/Ce* often reaching values greater than 100 (Chen et al. 2020;Morag et al. 2020), and as such is expected to contribute substantial Ce to the whole rock budget (along with trace contributions from epidote and allanite). This anomalous contribution of Ce would not be part of a smooth curve, and propagate the anomaly to the whole rock, under the assumption that Ce-rich zircon did not fractionate earlier and impart a measurable anomaly to the melt. This case study shows-in addition to the examples shown earlier-that Ce anomalies obtained using the deviation from the polynomial fit are more accurate where REE patterns exhibit curvature and better reflect geological and mineralogical observations relative to the deviation from the interpolation line. In cases of high curvature such as the one presented here, the inferred anomaly even has the opposite sign.

Concluding Remarks
We present a method for parameterising REE patterns, combining the orthogonal polynomials of O'Neill (2016) with adapted tetrad polynomials based on Minami and Masuda (1997). Our method provides uncertainty estimates and is implemented in two forms: an interactive online app written in R and Shiny (BLambdaR), and a Python package that includes our λ and τ fitting algorithm with other geochemical data analysis tools (pyrolite). The case studies show that there is no one single fitting recipe that applies universally, and we recommend that different fitting options be explored. The case studies also show some of the considerations that apply when treating real-world data and how they might be approached.