Integrating plant wax abundance and isotopes for paleo-vegetation and paleoclimate reconstructions: a multi-source mixing model using a Bayesian framework

Plant wax n-alkane chain length distribution and isotopes have been studied in modern ecosystems as proxies to reconstruct vegetation and climate of the past. However, most paleo-proxies focus on either concentrations or isotopes, whereas both carry complementary information on the mixing sources. We propose a multi-source mixing model in a Bayesian framework that evaluates both chain length distributions and isotopes simultaneously. The model consists of priors that include user-defined source groups and their associated parametric distributions of n-alkane concentration and δ¹³C. The mixing process involves newly defined mixing fractions such as fractional leaf mass contribution (FLMC) that can be used in vegetation reconstruction. Markov Chain Monte Carlo is used to generate samples from the posterior distribution of these parameters conditioned on both data types. We present three case studies from distinct settings. The first involves n-C₂₇, n-C₂₉, and n-C₃₁ alkanes in lake surface sediments of Lake Qinghai, China. The model provides more specific interpretations on the n-alkane input from aquatic sources than the conventional P_aq proxy. The second involves n-C₂₉, n-C₃₁, and n-C₃₃ alkanes in lake surface sediments in Cameroon, western Africa. The model produces mixing fractions of forest C₃, savanna C₃, and C₄ plants, offering additional information on the dominant biomes compared to the traditional two-end-member mixing regime. The third couples the vegetation source model to a hydrogen isotope model component, using biome-specific apparent fractionation factors (ε_a) to estimate the δ²H of mean annual precipitation. By leveraging chain length distribution, δ¹³C, and δ²H data of four n-alkane chains, the model produces estimated precipitation δ²H with relatively small uncertainty limits. The new framework shows promise for interpretation of paleo-data but could be further improved by including processes associated with n-alkane turnover in plants, transport, and integration into sedimentary archives. Future studies on modern plants and catchment systems will be critical to develop calibration datasets that advance the strength and utility of the framework.

1 Introduction

Plant wax n-alkyl compounds, including n-alkanes, n-alkanoic acids, n-alkanols, and n-esters are important biomarkers in paleoenvironmental reconstructions. Among them, mid- and long-chain n-alkanes (23–35 carbons) are the most frequently studied lipid compound class due to their great abundance in higher plants and excellent preservation in sedimentary archives (e.g., Sachse et al., 2012; Bush and McInerney, 2013; Freeman and Pancost, 2014; Diefendorf and Freimuth, 2017; Liu and An, 2020). The use of n-alkanes as paleoenvironmental proxies relies on our understanding of their biosynthesis and extensive surveys of their distributions in modern plants. These n-alkanes of various chain lengths are synthesized via the same acetogenic biochemical pathway in all higher plants, which results in the characteristic of odd-over-even carbon chain length predominance (e.g., Eglinton and Hamilton, 1967; Kolattukudy et al., 1976; Cheesbrough and Kolattukudy, 1984; Hayes, 1993; Kunst and Samuels, 2003; Chikaraishi et al., 2004). The persistence of such chain length predominance in sedimentary archives has been used to inform the state of preservation in plant-derived n-alkanes (e.g., Bray and Evans, 1961; Eglinton and Hamilton, 1967; Buggle et al., 2010; Zech et al., 2013; Brittingham et al., 2017).

Despite the shared odd-over-even chain length predominance, the absolute amounts and relative abundances of n-alkanes differ between major plant groups, and the distribution of chain lengths in sediments has been used as a chemical fingerprint to reconstruct vegetation composition. In terms of absolute amounts, angiosperms generally produce greater quantities of long-chain n-alkanes than gymnosperms (Diefendorf et al., 2011, 2015; Diefendorf and Freimuth, 2017), while terrestrial plants produce greater quantities of long-chain n-alkanes than submerged or emergent aquatic macrophytes (Ficken et al., 2000; Mead et al., 2005; Aichner et al., 2010; Liu et al., 2015; Liu and Liu, 2016). In terms of relative abundances, n-alkane chain length distributions differ between plant functional types such as trees, shrubs, forbs, grasses, and succulents and between deciduous and evergreen plants (Rommerskirchen et al., 2006; Vogts et al., 2009; Diefendorf et al., 2011; Bush and McInerney, 2013; Carr et al., 2014; Garcin et al., 2014; Badewien et al., 2015; Diefendorf et al., 2015; Magill et al., 2019). However, both absolute amounts and relative abundances of n-alkanes have been found to be highly variable within groups, as they are influenced by physiological and environmental factors such as photosynthetic pathways, leaf age, temperature, elevation, light intensity and water stress (e.g., Koch et al., 2006; Diefendorf et al., 2011; Macková et al., 2013; Tipple et al., 2013; Bush and McInerney, 2015; Feakins et al., 2016; Liu et al., 2017; Suh and Diefendorf, 2018; Griepentrog et al., 2019). The large uncertainties associated with characteristic n-alkane chain length distributions have limited their application in quantitative vegetation reconstruction (Bush and McInerney, 2013). Several models have been developed to reconstruct vegetation composition from n-alkyl lipid chain length distribution (Jansen et al., 2010; Gao et al., 2011; Peaple et al., 2021), but they all require a wide spectrum of chains to be analyzed, which may not be feasible depending on the type of sediment analyzed.

Stable carbon isotopes of chain-length-specific n-alkanes have also been used extensively to indicate changes in vegetation of the past. This is based on the well-established empirical evidence that the distribution of carbon isotope ratios (using the δ notation) of plant n-alkanes largely mirrors that of bulk plant tissue, primarily reflecting differences in the photosynthetic pathways among terrestrial plants (e.g., C₃ vs. C₄) and the source of carbon in aquatic plants (e.g., Collister et al., 1994; Mead et al., 2005; Aichner et al., 2010; Diefendorf and Freimuth, 2017; Liu and An, 2020). For example, terrestrial C₄ plants produce more ¹³C-enriched n-alkanes than C₃ plants (e.g., Bi et al., 2005; Rommerskirchen et al., 2006; Vogts et al., 2009; Kristen et al., 2010; Badewien et al., 2015). Such a distinction has been used to identify large-scale changes between C₃ and C₄ vegetation based on δ¹³C values of n-alkanes extracted from sedimentary archives (e.g., Bird et al., 1995; Huang et al., 2000; Freeman and Colarusso, 2001; Schefuß et al., 2003; Hughen et al., 2004; Castañeda et al., 2009; Niedermeyer et al., 2010; Tipple and Pagani, 2010; Schefuß et al., 2011; Vogts et al., 2012; Uno et al., 2016; Zhou et al., 2017; Andrae et al., 2018; Polissar et al., 2019). At a finer scale, carbon isotope compositions of n-alkanes from C₃ plants are influenced by both biological and environmental factors, such as taxonomy, plant functional types, water availability, temperature, and altitude (e.g., Diefendorf et al., 2010; Cernusak et al., 2013; Diefendorf et al., 2015; Wu et al., 2017; J. Wang et al., 2018; Liu and An, 2020). Compared to the uncertainties associated with n-alkane chain length distribution, the uncertainties of δ¹³C values in C₃ and C₄ plants are relatively constrained and have been incorporated into reconstructions of vegetation cover (e.g., Garcin et al., 2014; Uno et al., 2016; Andrae et al., 2018; Polissar et al., 2019). However, when the reconstructions are based on a linear mixing relationship with a careful selection of one specific n-alkane compound (often n-C₂₉ or n-C₃₁), they are associated with the challenge that n-alkane production of the selected chain could vary between the end-members. If the uncertainties associated with n-alkane production are not accounted for, it is possible that the uncertainty associated with paleo-vegetation reconstructions could be mischaracterized and that the interpretations could be biased.

Traditional interpretations of n-alkane proxies often rely on one line of evidence: either chain length distribution or δ¹³C values of one chain. While many studies have reported δ¹³C values and relative abundance of multiple n-alkane chains with their respective uncertainties, they are often interpreted independently, using the interpretation of one to qualitatively support that of another. A more rigorous approach to interpreting these data could be developed through the use of a proxy system model (Evans et al., 2013; Dee et al., 2018; Konecky et al., 2019), in which the common and unique assumptions underlying the interpretation of both chain-length and isotope data are explicitly represented. Such a model would provide a quantitative framework within which one proxy type can be used to constrain the uncertainties involved in the interpretation of the other, and vice versa. In the case of plant n-alkanes, both chain length distribution and δ¹³C values of multiple chains carry complementary information related to the mixing of compounds from different vegetation sources. The combination of both proxies, therefore, may help to better constrain uncertainties, reduce biases, and refine our interpretations. This study describes a generic multi-source n-alkane mixing model that incorporates both lines of evidence (chain lengths, carbon isotopes) and their associated uncertainties in the interpretation. Using case studies based on published sedimentary n-alkane records, the primary goals of this study are (1) to demonstrate model utility by proposing and evaluating multi-source n-alkane mixing regimes, (2) to demonstrate model potential by providing informed interpretations on the mixing regimes with well characterized uncertainties, and (3) to explore the use of the model in paleoclimate reconstruction by linking the mixing model to a simple model describing the ecosystem-specific H-isotope fractionation between n-alkanes and local precipitation.

2 Methods

2.1 Model structure

The proposed approach is achieved in a Bayesian hierarchical modeling framework (Fig. 1), which can leverage information from multiple proxies to provide a robust statistical basis for proxy integration. The hierarchical model is then inverted using Markov Chain Monte Carlo (MCMC) methods (Geman and Geman, 1984) to obtain posterior parameter estimates that are conditioned simultaneously on all proxy data. Similar modeling approaches have been applied to meta-analyses of paleoclimatic and vegetation proxies (e.g., Garreta et al., 2010; Li et al., 2010; Tingley et al., 2012; Bowen et al., 2020) but have not been specifically proposed for proxy interpretation of n-alkanes.

Figure 1 Proposed model structure of the Bayesian hierarchical framework for interpretation of n-alkane records, where model components are represented by boxes. Links between components are numbered and referenced in Sect. 2. RA stands for relative abundance of n-alkane chains. FLMC_i is the fractional leaf mass contribution from source i. FSC_n,i is the fractional contribution from source i to chain n.

2.1.1 Prior distributions

Our understanding of species-level n-alkane chain length distribution and chain-specific δ¹³C values is based on empirical evidence from extant plants (Sachse et al., 2012; Bush and McInerney, 2013; Diefendorf and Freimuth, 2017). The extensive record of literature and the published empirical data form the basis for prescribing the prior chain length distributions used in this work.

The n-alkane concentrations (µg g⁻¹ of dried leaves) are influenced by the taxonomy, growth forms, and growing conditions of the plant (Han et al., 1968; Ficken et al., 2000; Rommerskirchen et al., 2006; Vogts et al., 2009; Diefendorf et al., 2011; Bush and McInerney, 2013, 2015; Andrae et al., 2019; Magill et al., 2019). Within a source group that is associated with a specific combination of these factors, we expect n-alkane concentration to follow a group-specific chain length distribution pattern. The distribution of concentration values for a given chain within a well-defined source group typically follows a log-normal distribution (Garcin et al., 2014). Therefore, we assume that for each sample drawn from source group i (Fig. 1a), n-alkane concentrations are derived from a multivariate log-normal distribution:

$$\begin{array}{}\text{(1)}& {\mathrm{Conc}}_{\mathrm{1}:n,i}\sim \mathrm{multi}-\ln -N\left({\mathit{\mu }}_{\mathrm{1}:n,i}{\mathbf{\Omega }}_{\mathrm{1}:n,i}\right),\end{array}$$

where ${\mathit{\mu }}_{\mathrm{1}:n,i}$ are the means of the natural log-transformed n-alkane concentration values for chains 1 to n, while ${\mathbf{\Omega }}_{\mathrm{1}:n,i}$ is a n × n variance–covariance matrix of the natural log-transformed n-alkane concentration values calculated from empirical data for each source group i.

Similarly, n-alkane δ¹³C values are primarily influenced by taxonomy, photosynthetic pathways, and growing conditions of the plant (e.g., Bi et al., 2005; Rommerskirchen et al., 2006; Vogts et al., 2009; Diefendorf et al., 2010; Kristen et al., 2010; Cernusak et al., 2013; Badewien et al., 2015; Diefendorf et al., 2015; Wu et al., 2017; J. Wang et al., 2018; Liu and An, 2020). Within a source group i that is associated with a specific combination of these factors, we expect n-alkane δ¹³C to follow a group-specific distribution pattern. Because n-alkanes of different chain lengths are synthesized via the same biochemical pathway, ${\mathit{\delta }}^{\mathrm{13}}{\mathrm{C}}_{\mathrm{1}:n,i}$ of all chains should be correlated (Eglinton and Hamilton, 1967; Cheesbrough and Kolattukudy, 1984; Hayes, 1993; Kunst and Samuels, 2003; Chikaraishi et al., 2004). Therefore, we assume that for each sample drawn from source group i (Fig. 1a), ${\mathit{\delta }}^{\mathrm{13}}{\mathrm{C}}_{\mathrm{1}:n,i}$ are derived from a multivariate normal distribution:

$$\begin{array}{}\text{(2)}& {\mathit{\delta }}^{\mathrm{13}}{\mathrm{C}}_{\mathrm{1}:n,i}\sim \mathrm{multi}-N\left({\mathit{\mu }}_{\mathrm{1}:n,i}^{\mathrm{C}}{\mathbf{\Omega }}_{\mathrm{1}:n,i}^{\mathrm{C}}\right),\end{array}$$

where ${\mathit{\mu }}_{\mathrm{1}:n,i}^{\mathrm{C}}$ are the means of δ¹³C for chains 1 to n, while ${\mathbf{\Omega }}_{\mathrm{1}:n,i}^{\mathrm{C}}$ is a n × n variance–covariance matrix calculated from empirical data for each source group i.

Calculation of the variance–covariance matrices does not allow for missing values in the empirical dataset. Therefore, data entries with any missing value are removed before calculation of the variance–covariance matrices. The prior distribution parameters allow random samples to be drawn from the prior distributions and used in the process model calculation (Step 1 in Fig. 1). By using prior distributions derived entirely from modern plant n-alkane data to inform the interpretation of sedimentary alkane data, the version of the model implemented here implicitly assumes no change in the δ¹³C value of atmospheric CO₂ and minimal incorporation of reworked and potentially pre-aged n-alkanes into sedimentary archives.

2.1.2 Process model

The model consists of a generic mixing process with multiple sources, following the principle of isotope mass balance. First, a fixed number of random draws (k) from the prior distributions are used to generate samples of Conc_n,i and δ¹³C_n,i to characterize chain n for each source group i (similar to Fig. 1a). Concentration values are in units of micrograms of n-alkane per gram of dried leaves, and the random draws thus give a distribution of alkane yields per unit leaf mass sampled from the prior. Both concentration and isotope data used to specify the prior distributions typically represent measurements from individual plants, and the process of drawing these k random samples is intended to represent the production of an integrated sample from k individuals. Although in many cases we would expect more than 50 individuals of a particular source group to contribute to a sedimentary sample, k = 50 is chosen to balance sample representation of prior distributions, process stochasticity, and algorithmic complexity. In general, larger values of k are more prescriptive in that the random sample will conform more closely to the prior distribution, and we consider k = 50 to be appropriately conservative in terms of allowing some variation in the samples drawn from the priors.

Second, we calculate the mixture of alkanes from all groups as a function of the fractional leaf mass contribution (FLMC), where FLMC_i gives the relative leaf biomass from group i contributing alkanes to the sedimentary sample. Because $\sum _i{\mathrm{FLMC}}_i=\mathrm{1}$; ${\mathrm{FLMC}}_{\mathrm{1}},\mathrm{\dots },{\mathrm{FLMC}}_i$ follows a Dirichlet distribution:

$$\begin{array}{}\text{(3)}& {\mathrm{FLMC}}_{\mathrm{1}},\mathrm{\dots },{\mathrm{FLMC}}_i\sim \mathrm{Dir}\left({\mathit{\alpha }}_{\mathrm{1}},\mathrm{\dots },{\mathit{\alpha }}_i\right),\end{array}$$

where ${\mathit{\alpha }}_{\mathrm{1}},\mathrm{\dots },{\mathit{\alpha }}_i$ are the Dirichlet concentration parameters. Using an uninformative prior assuming equal mixture of the sources, the Dirichlet concentration parameters are set to be the same. The values ${\mathit{\alpha }}_{\mathrm{1}}=\mathrm{\dots }={\mathit{\alpha }}_i=\mathrm{1}$ are chosen assuming a uniform probability of FLMC_i within its range (0, 1). It is important to note that potential differences in leaf mass or lipid turnover are not incorporated into this formulation: the values of FLMC_i may not be directly related to the standing biomass of source i if the alkane and/or leaf turnover rates for the groups differ.

Third, the weighted average relative abundance (RA) of each n-alkane n in the mixture can be calculated from the concentrations and mixing fractions:

$$\begin{array}{}\text{(4)}& {\mathrm{RA}}_{n,\mathrm{mix}}={\displaystyle \frac{\sum _i\left({\mathrm{FLMC}}_i\sum {\mathrm{Conc}}_{n,i}\right)}{\sum _i^n\left({\mathrm{FLMC}}_i\sum {\mathrm{Conc}}_{n,i}\right)}},\end{array}$$

where in the numerator, $\sum {\mathrm{Conc}}_{n,i}$ is the sum of concentrations from the random draws (k) for chain n and source group i, while $\sum _i\left({\mathrm{FLMC}}_i\sum {\mathrm{Conc}}_{n,i}\right)$ is the fraction-weighted sum of concentrations for chain n of all source groups. In the denominator, ${\sum }_i^n\left({\mathrm{FLMC}}_i\sum {\mathrm{Conc}}_{n,i}\right)$ is the total n-alkane concentration for all chains and all source groups. Equation (4) dictates that $\sum {\mathrm{RA}}_{n,\mathrm{mix}}\equiv \mathrm{1}$.

Lastly, for chain n and its δ¹³C_n,mix, according to isotope mass balance

$$\begin{array}{}\text{(5)}& {\mathit{\delta }}^{\mathrm{13}}{\mathrm{C}}_{n,\mathrm{mix}}={\displaystyle \frac{\sum _i\left[{\mathrm{FLMC}}_i\sum \left({\mathrm{Conc}}_{n,i}\times {\mathit{\delta }}^{\mathrm{13}}{\mathrm{C}}_{n,i}\right)\right]}{\sum _i\left({\mathrm{FLMC}}_i\sum {\mathrm{Conc}}_{n,i}\right)}},\end{array}$$

where in the numerator $\sum \left({\mathrm{Conc}}_{n,i}\times {\mathit{\delta }}^{\mathrm{13}}{\mathrm{C}}_{n,i}\right)$ is the sum of the products of n-alkane concentration and δ¹³C in each round of random draws (k) for source i, while $\sum _i\left[{\mathrm{FLMC}}_i\sum \left({\mathrm{Conc}}_{n,i}\times {\mathit{\delta }}^{\mathrm{13}}{\mathrm{C}}_{n,i}\right)\right]$ is the fraction-weighted sum of the sums for chain n of all source groups. In the denominator, $\sum _i\left({\mathrm{FLMC}}_i\sum {\mathrm{Conc}}_{n,i}\right)$ is the same as in the numerator in Eq. (4). Equation (5) is essentially a concentration-weighted isotope mass balance equation for δ¹³C.

The process model specifies the numerical relationships between random samples from the prior distributions and the simulated metrics of interest (RA_n,mix, δ¹³C_n,mix). The simulated metrics are then used in the data model for further evaluations (Step 2 in Fig. 1).

2.1.3 Data model and model inversion

All proxy data are subject to errors that are associated with the proxy observations themselves (Evans et al., 2013). Therefore, RA_n,mix and δ¹³C_n,mix are modeled with their respective measurement error (Fig. 1c). Because each of the measured RA values are ratios by nature, they are assumed to be associated with a Cauchy error term while centering around the true value RA_n,mix:

$$\begin{array}{}\text{(6)}& {\displaystyle }{\mathrm{RA}}_{n,\mathrm{mea}}\sim \mathrm{Cauchy}\left({\mathrm{RA}}_{n,\mathrm{mix}},{\mathit{\tau }}_{n,\mathrm{mea}}\right),\text{(7)}& {\displaystyle }{\left({\mathit{\tau }}_{n,\mathrm{mea}}\right)}^{\mathrm{2}}\sim \mathrm{Gamma}(\mathrm{shape}=\mathrm{5},\mathrm{rate}=\mathrm{0.05}),\end{array}$$

where τ_n,mea is the measurement error of RA_n,mea, and the hyper parameter shape and rate are set to 5 and 0.05, respectively, corresponding to the small error margins of n-alkane concentration measurements and the associated RA values, which are often reported to the third decimal digit.

Because the uncertainty of the measured δ¹³C values typically follows a normal distribution, they are assumed to be associated with a Gaussian error term while centering around the true value δ¹³C_n,mix:

$$\begin{array}{}\text{(8)}& {\mathit{\delta }}^{\mathrm{13}}{\mathrm{C}}_{n,\mathrm{mea}}\sim N\left({\mathit{\delta }}^{\mathrm{13}}{\mathrm{C}}_{n,\mathrm{mix}},{\mathit{\sigma }}_{n,\mathrm{mea}}\right),\end{array}$$

where σ_n,mea is the measurement error of δ¹³C_n,mea, which is often reported as 1 standard deviation (SD) in the literature.

In the model inversion step (Steps 3 to 5 in Fig. 1), MCMC is used to propose samples of all model parameters conditioned on the measured δ¹³C and RA results.

2.1.4 Model implementation

The model structure described above is coded in the BUGS (Bayesian inference Using Gibbs Sampling) language (Lunn et al., 2012), and implemented in R version 4.0.5 (R Core Team, 2021), using the “rjags” package with the standalone JAGS (Just Another Gibbs Sampler) encoder installed separately (Plummer, 2021). Three chains are run in parallel, and the number of iterations is set at 800 000 to ensure model convergence, with the first 200 000 interactions as burn-ins. Chain thinning is set at once per 240 iterations. Convergence is assessed visually via trace plots (R package “mcmcplots”, Curtis, 2018) and with reference to the convergence factor “rhat” (Gelman and Rubin, 1992) and effective sample sizes reported by the “rjags” package. The iteration parameters are chosen to ensure complete convergence with rhat values smaller than 1.01. The average run time for a three-chain δ¹³C evaluation is approximately 2 h.

2.1.5 Model output

Once the model inversion is completed for each model parameter, its posterior density is summarized via the kernel density estimation function “density” in the R package “stats” with default settings. For selected parameters, posterior density summaries such as the maximum a posteriori estimation (MAPE), median density estimation, and the 89 % highest density interval (HDI) are reported using functions in the R package “bayestestR” (Makowski et al., 2019). Out of all the parameters, the posterior densities of two parameters are of special interest in the mixing process. One of these parameters is FLMC_i (as in Eq. 3, Fig. 1e). The other parameter is derived from the posterior samples generated by the MCMC process, which we term as the fractional contribution of source i to each specific chain n (FSC_n,i, Fig. 1e):

$$\begin{array}{}\text{(9)}& {\mathrm{FSC}}_{n,i}={\displaystyle \frac{{\mathrm{FLMC}}_i\sum {\mathrm{Conc}}_{n,i}}{\sum _i\left({\mathrm{FLMC}}_i\sum {\mathrm{Conc}}_{n,i}\right)}},\end{array}$$

where the numerator is the fraction-weighted sum of randomly drawn (k) concentrations for chain n of source i. The denominator is the same as the numerator of Eq. (4). Bivariate density plots and contour lines are used to explore covariation patterns of FLMC_i among the sources.

2.2 Case studies

Here we provide three case studies that demonstrate model characteristics and offer alternative interpretations of previously published n-alkane proxy data. We considered that n-alkanes extracted from lake surface sediments fulfill most of the assumptions associated with the model structure, and we therefore chose lake surface sediments as the archive of interest for the first two case studies. The third case study is based on n-alkanes extracted from marine core sediments to demonstrate model utility in paleorecords. The case studies are not intended to formally test or validate the modeling approach but to illustrate how the model framework can offer an accessible and statistically informative approach for interpretations of n-alkanes and similar proxy data.

2.2.1 Evaluating aquatic and terrestrial n-alkane input in sediment samples of Lake Qinghai

Background

Lacustrine sediments can incorporate n-alkanes that are sourced from both terrestrial and aquatic plants (Ficken et al., 2000; Mead et al., 2005; Aichner et al., 2010). Quantifying contributions from the different sources to specific n-alkanes is important because it guides paleoenvironmental interpretation of n-alkane proxies. For example, n-alkane chain length distribution and stable isotopes (δ¹³C and δ²H) have been used to investigate the composition of sedimentary organic matter and regional hydrological conditions through time (e.g., Ficken et al., 1998; Schwark et al., 2002; Huang et al., 2004; Sachse et al., 2004; Mead et al., 2005; Aichner et al., 2010). Many studies use the proxy $P_{\mathrm{aq}}=({\mathrm{C}}_{\mathrm{23}}+{\mathrm{C}}_{\mathrm{25}})/({\mathrm{C}}_{\mathrm{23}}+{\mathrm{C}}_{\mathrm{25}}+{\mathrm{C}}_{\mathrm{29}}+{\mathrm{C}}_{\mathrm{31}})$ (Ficken et al., 2000) as an indicator of terrestrial or aquatic source dominance to interpret measured isotope ratios of sedimentary n-alkanes. It quantifies the relative proportion of mid-chain (n-C₂₃ and n-C₂₅) vs. long-chain homologs, based on the evidence that mid-chain alkanes are largely produced by aquatic plants, while long-chain alkanes (n-C₂₉ and n-C₃₁) are largely produced by terrestrial plants (e.g., Ficken et al., 2000; Aichner et al., 2010; Gao et al., 2011; Duan and Xu, 2012; Wang and Liu, 2012; Liu and Liu, 2016). However, long-chain alkanes (e.g., n-C₂₇ and n-C₂₉) have been found in aquatic macrophytes, while mid-chain alkanes (e.g., n-C₂₅) have been found in some terrestrial plants (e.g., Aichner et al., 2010; Gao et al., 2011; Duan and Xu, 2012; Wang and Liu, 2012; Liu et al., 2015; Dion-Kirschner et al., 2020). These findings add some uncertainty to the use of P_aq and similar indices as indicators of n-alkane contributions from aquatic plants.

The δ¹³C values of n-alkanes have also been used as an indicator of terrestrial or aquatic source dominance (e.g., Mead et al., 2005; Aichner et al., 2010; Liu et al., 2015). This is based on the observation that some aquatic macrophytes that grow in shallow waters display δ¹³C values that are distinctive from those of terrestrial C₃ plants and algae due to their ability to utilize dissolved bicarbonate ions as their carbon source (e.g., Allen and Spence, 1981; Prins and Elzenga, 1989; Keeley and Sandquist, 1992). However, the interpretation of δ¹³C values for a specific n-alkane homologue using a linear two-end mixing relationship will be sensitive to differences in n-alkane concentrations between the end-member vegetation types, leading to contradictory interpretations from different compounds. The presence of additional (beyond two) potential sources further complicates the interpretation of such data.

To demonstrate how our model can provide quantitative estimates of n-alkane mixing fractions of multiple aquatic and terrestrial plant sources, simultaneously conditioned on isotopic and concentration data for all measured chain lengths, we apply the model using published lake surface sediment data from three locations in Lake Qinghai.

Data compilation

The region of interest is the Qinghai–Tibetan Plateau in western China, where both freshwater and saline lakes are abundant. The land cover of the region is dominated by alpine meadow, steppe, and shrubland, which consist of almost exclusively C₃ plants (Yu et al., 2001; Duan and Xu, 2012; Liu et al., 2015). Most of the lakes in the region are shallow and densely populated with submerged aquatic macrophytes and green algae during the short growing season (Aichner et al., 2010; Liu et al., 2015; Liu and Liu, 2016). Based on these potential sources of n-alkanes in the region, we define three source groups: (1) terrestrial plants, (2) aquatic macrophytes, and (3) algae.

Table 1 Lake surface sediment samples used in case study 1, with measured chain-specific n-alkane concentrations and δ¹³C values; these data were originally published by Liu et al. (2015).

Figure 2 Empirical data of per sample n-alkane concentrations (µg g⁻¹ of dried sample, blue histograms) and δ¹³C (‰, VPDB, red histograms) of n-C₂₇, n-C₂₉, and n-C₃₁ alkanes and their corresponding prior distributions (blue Gaussian curve overlays) of published plants in the terrestrial, aquatic macrophyte, and algae sources; raw data are compiled in supplementary data EA-2 (Yang, 2022). Estimated prior parameters of each source, including the means and variance–covariance matrices, are reported in Tables S1 and S2.

Per-sample n-alkane concentrations (µg g⁻¹ of dried sample) and δ¹³C (‰, VPDB) of n-C₂₇, n-C₂₉, and n-C₃₁ alkanes are compiled from the literature for terrestrial plants in the region and aquatic macrophytes and algae in the lakes (Aichner et al., 2010; Liu et al., 2015; Liu and Liu, 2016). The compiled dataset (supplementary data EA-2, Yang, 2022) is used to provide parameter estimates for the prior distributions (as explained in Sect. 2.1.1) for each source. Parameters of the multivariate normal distributions are estimated using the R packages “stats” (R Core Team, 2021). To evaluate the assumption that n-alkane concentrations follow a log-normal distribution for each chain in each source, the distribution parameters are estimated using the R package “EnvStats” (Millard, 2013). Goodness of fit is visually assessed using quantile–quantile plots (Fig. S1 in the Supplement). The n-alkane concentrations and δ¹³C values of selected lake surface sediments from Lake Qinghai are compiled from Liu et al. (2015). Other data that characterize the composition and preservation of n-alkanes in the geological archives, such as the P_aq proxy, are also referenced in the data compilation (Table 1). Because the data are modeled to follow multivariate normal and log-normal distributions (Sect. 2.1.1), the prior distributions should be theoretically assessed in a multivariate space, but such an assessment would be difficult to visualize. Instead, for each variable in the multivariate space, the multivariate density is “collapsed” onto a single dimension by integrating the densities in other dimensions (as in Fig. 2) using the function “omxMnor” in the R package “OpenMx” (Neale et al., 2016).

2.2.2 Evaluating n-alkane input along a vegetation gradient in lake sediments in western Africa

Background

Biome composition and its change over time have profound implications on climatic shifts and the evolution of the biosphere. The tropical grassland biome has been of particular interest because of its unique δ¹³C signature and its important association with climate of the past and mammalian evolution (e.g., Cerling, 1992; Quade and Cerling, 1995; Latorre et al., 1997; Janis et al., 2002; Bobe and Behrensmeyer, 2004; Strömberg, 2011; Kaya et al., 2018; Polissar et al., 2019). Many studies used δ¹³C values of long-chain n-alkanes to provide estimates of C₃ $/$ C₄ vegetation cover, based on linear mixing calculations using data from one or more chain lengths (Garcin et al., 2014). In many cases, n-C₃₁ alkane is used, as its concentration is often similar in C₃ and C₄ plants (e.g., Uno et al., 2016; Polissar et al., 2019). However, this approach does not take advantage of the distinct chain length dominance patterns associated with different biomes. For example, C₄ grasses are associated with much higher production of the n-C₃₃ alkane than most C₃ plants (Krull et al., 2006; Rommerskirchen et al., 2006; Bush and McInerney, 2013; Garcin et al., 2014), while tropical rainforest plants are associated with the very high production of the n-C₂₉ alkane (Vogts et al., 2009). Using the information from chain length distribution in addition to δ¹³C values has the potential to produce more detailed reconstructions of biome composition beyond the binary C₃ $/$ C₄ ratio.

To demonstrate how the model can provide quantitative estimates of n-alkane contributions from multiple broadly defined biomes, we apply the model using published lake surface sediment data from three locations in a vegetation gradient in Cameroon, western Africa.

Data compilation

The region of interest is Sub-Saharan Africa where rainfall amount and seasonality are the primary determinants of biome types (Sankaran et al., 2005; Aleman et al., 2020). In western Africa, vegetation cover is dominated by rainforest close to the Equator, wooded grassland to the north, and a transitional zone in between (Huang et al., 2000; Rommerskirchen et al., 2003; Garcin et al., 2012, 2014; Schwab et al., 2015). In particular, tropical forest and savanna woody plant species have been shown to co-occur in this region (Aleman et al., 2020), making it an ideal place to investigate the potential of using n-alkane proxies to reconstruct the dominant C₃ biome. Based on the above information, we define three source groups: (1) tropical C₄ plants, (2) savanna C₃ plants, and (3) rainforest C₃ plants. Succulent plants are excluded in the analysis due to their low presence in the region.

Table 2 Lake surface sediment samples used in case study 2, with measured chain-specific n-alkane concentrations and δ¹³C values; these data were originally published by Garcin et al. (2014).

Figure 3 Empirical data of per-sample n-alkane concentrations (µg g⁻¹ of dried sample, blue histograms) and δ¹³C (‰, VPDB, red histograms) of n-C₂₉, n-C₃₁, and n-C₃₃ alkanes and their corresponding prior distributions (blue Gaussian curve overlays) of published plants in the C₄ plants, savanna C₃, and rainforest C₃ sources. Raw data are compiled in supplementary data EA-3 (Yang, 2022). Estimated prior parameters of each source, including the means and variance-covariance matrices, are reported in Tables S3 and S4.

Per sample n-alkane concentrations (µg g⁻¹ of dried sample) and δ¹³C (‰, VPDB) of n-C₂₉, n-C₃₁, and n-C₃₃ alkanes are compiled from the literature for terrestrial plants in western, southwestern, and eastern Africa (Ali et al., 2005; Rommerskirchen et al., 2006; Vogts et al., 2009; Kristen et al., 2010; Badewien et al., 2015; Diefendorf and Freimuth, 2017; Magill et al., 2019). The compiled dataset (supplementary data EA-3, Yang, 2022) is used to provide parameter estimates for the prior distributions for each source (as explained in Sect. 2.1.1). Prior parameters are estimated using the same methods as in Sect. 2.2.1. Goodness of fit is visually assessed using quantile–quantile plots (Fig. S2). The n-alkane concentrations and δ¹³C values of selected lake surface sediments are compiled from Garcin et al. (2014), with estimated fractional C₄ vegetation cover ($f_{{\mathrm{C}}_{\mathrm{4}}}$) referenced from the same publication (Table 2). Following the same methods as in Sect. 2.2.1 (“Data compilation”), univariate prior distributions are plotted in Fig. 3 as density curves.

2.2.3 Estimating precipitation δ²H from n-alkane records in marine sediments off the Zambezi River mouth

Background

Hydrogen isotope ratios (δ²H, ‰, VSMOW) of terrestrial plant wax lipids have been extensively used to interpret changes in hydroclimate and rainfall regimes in the African tropics (e.g., Tierney et al., 2008; Konecky et al., 2011; Schefuß et al., 2011; Berke et al., 2012; Feakins, 2013; Kuechler et al., 2013; Collins et al., 2014; Costa et al., 2014; Shanahan et al., 2015; Castañeda et al., 2016; Niedermeyer et al., 2016; Caley et al., 2018; Norström et al., 2018). This is based on the common assumption that variation in precipitation δ²H is the primary driver of δ²H variation in terrestrial plant wax lipids (Huang et al., 2004; Sachse et al., 2006; Hou et al., 2008). Investigations into the apparent fractionation factor (ε_a) between plant wax lipids and source water revealed that ε_a of grasses (both C₃ and C₄) is much more negative than that of other C₃ plants (Krull et al., 2006; Liu et al., 2006; Smith and Freeman, 2006; Sachse et al., 2012; Kahmen et al., 2013a; Gamarra et al., 2016; Z. Wang et al., 2018). A strong positive correlation is also found between aridity and ε_a in C₃ plants, while a much weaker correlation is present in C₄ grasses (Sachse et al., 2004; Hou et al., 2008; Feakins and Sessions, 2010; Polissar and Freeman, 2010; McInerney et al., 2011; Sachse et al., 2012; Kahmen et al., 2013a, b; Schwab et al., 2015; Gamarra et al., 2016). These findings suggest that both vegetation type (e.g., C₃ vs. C₄ plants) and aridity (e.g., in rainforest vs. savanna C₃ biomes) can impact the interpretation of sedimentary plant wax δ²H records (Douglas et al., 2012; Collins et al., 2013; Feakins, 2013; Nelson et al., 2013; Wang et al., 2013a; Z. Wang et al., 2018). In addition, as mentioned in Sect. 2.2.2, C₄ grasses produce more n-C₃₁ and n-C₃₃ alkanes than C₃ plants in general. As a result, when the n-alkanes are sourced from a C₃ $/$ C₄ mixed-vegetation region, the magnitude of change through time in δ²H of n-C₃₁ and n-C₃₃ alkanes can be higher than that of n-C₂₉ alkane (Schefuß et al., 2011; Wang et al., 2013a). While different methods have been proposed to account for the vegetation effect on lipid δ²H (Collins et al., 2013; Feakins, 2013; Magill et al., 2013; Konecky et al., 2014), most interpretations are based on data of one selected chain (e.g., either n-C₂₉ or n-C₃₁ for n-alkanes) while data of other chains are not directly interpreted. With the model structure in Sect. 2.1 and prior distributions of broadly defined biomes in Sect. 2.2.2, we assume that the mixing fractions among the vegetation types can be resolved using δ¹³C and chain length distribution of multiple n-alkane chains. One possible solution to account for the effects of aridity to n-alkane δ²H is to use biome-specific ε_a between n-alkanes and mean annual precipitation (MAP). However, challenges remain in the high parametric uncertainty of observed ε_a (e.g., in n-C₂₉, 1σ ranges from 22 ‰ to 32 ‰ in different plant functional types; Konecky et al., 2019), which can introduce a substantial amount of uncertainty in the reconstructed precipitation δ²H values from a single n-alkane chain (Polissar and D'Andrea, 2014). While the parametric uncertainties of individual n-alkane chains are likely to be similar, simultaneous evaluation of δ²H in multiple chains can theoretically reduce the overall uncertainty of estimated δ²H of MAP by leveraging their common sensitivities to precipitation.

To demonstrate how the reconstructed biome compositions in Sect. 2.2.2 can be used to further assist interpretation of n-alkane δ²H and quantitative estimation of precipitation δ²H, we expand the existing model structure outlined in Sect. 2.1 to permit evaluation of associated n-alkane δ²H, δ¹³C, and relative abundance of multiple chains simultaneously. We apply the expanded model to published marine core data (GIK 16160-3, lat. −18.2412^∘, long. 37.8685^∘) off the Zambezi River mouth, Indian Ocean (Wang et al., 2013a). This dataset also allows for a qualitative comparison (results reported in Supplement) to the pollen record from another marine core nearby (GeoB9311-1, lat. −21.55^∘, long. 36.4166^∘; Dupont and Kuhlmann, 2017).

Prior distributions of ε_a

To estimate δ²H in MAP, we first use the definition of the apparent fractionation factor (ε_a) between δ²H of n-alkane and MAP:

$$\begin{array}{}\text{(10)}& {\mathit{\epsilon }}_{\mathrm{a}n,i}=\left({\displaystyle \frac{{\mathit{\delta }}^{\mathrm{2}}{\mathrm{H}}_{n,i}+\mathrm{1000}}{{\mathit{\delta }}^{\mathrm{2}}{\mathrm{H}}_{\mathrm{MAP}}+\mathrm{1000}}}-\mathrm{1}\right)\times \mathrm{1000},\end{array}$$

where the δ²H terms are numerically expressed with a factor of 1000.

Solving Eq. (10), δ²H_n,i can be expressed as follows:

$$\begin{array}{}\text{(11)}& {\mathit{\delta }}^{\mathrm{2}}{\mathrm{H}}_{n,i}=\left({\displaystyle \frac{{\mathit{\epsilon }}_{\mathrm{a}n,i}}{\mathrm{1000}}}+\mathrm{1}\right)\times \left({\mathit{\delta }}^{\mathrm{2}}{\mathrm{H}}_{\mathrm{MAP}}+\mathrm{1000}\right)-\mathrm{1000},\end{array}$$

where δ²H of alkane chain n and source group i is determined by ε_a of chain n and source group i, as well as δ²H_MAP.

Similar to Eq. (2), we expect ε_a among the chains of source group i to follow a group-specific distribution pattern and assume that for each sample drawn from source group i, ${\mathit{\epsilon }}_{\mathrm{a}\mathrm{1}:n,i}$ are derived from a multivariate normal distribution:

$$\begin{array}{}\text{(12)}& {\mathit{\epsilon }}_{\mathrm{a}\mathrm{1}:n,i}\sim \mathrm{multi}-N\left({\mathit{\mu }}_{\mathrm{1}:n,i}^{\mathit{\epsilon }}{\mathbf{\Omega }}_{\mathrm{1}:n,i}^{\mathit{\epsilon }}\right),\end{array}$$

where ${\mathit{\mu }}_{\mathrm{1}:n,i}^{\mathit{\epsilon }}$ are the means of ε_a for chains 1 to n, while ${\mathbf{\Omega }}_{\mathrm{1}:n,i}^{\mathit{\epsilon }}$ is a n × n variance–covariance matrix calculated from empirical data for each source group i.

δ²H of MAP is modeled to follow a normal distribution with the following mean and precision parameters:

$$\begin{array}{}\text{(13)}& {\mathit{\delta }}^{\mathrm{2}}{\mathrm{H}}_{\mathrm{MAP}}\sim N\left({\mathit{\mu }}_{\mathrm{MAP}}^{{\mathit{\delta }}^{\mathrm{2}}\mathrm{H}},{\mathit{\tau }}_{\mathrm{MAP}}^{{\mathit{\delta }}^{\mathrm{2}}\mathrm{H}}\right),\end{array}$$

where ${\mathit{\mu }}_{\mathrm{MAP}}^{{\mathit{\delta }}^{\mathrm{2}}\mathrm{H}}$ is modeled with an uninformative prior using a uniform distribution with the following bounding values:

$$\begin{array}{}\text{(14)}& {\mathit{\mu }}_{\mathrm{MAP}}^{{\mathit{\delta }}^{\mathrm{2}}\mathrm{H}}\sim U(-\mathrm{200},\mathrm{20}),\end{array}$$

and ${\mathit{\tau }}_{\mathrm{MAP}}^{{\mathit{\delta }}^{\mathrm{2}}\mathrm{H}}$ is modeled with a target error of ±10 ‰, using the following shape and rate parameters:

$$\begin{array}{}\text{(15)}& {\left({\mathit{\tau }}_{\mathrm{MAP}}^{{\mathit{\delta }}^{\mathrm{2}}\mathrm{H}}\right)}^{\mathrm{2}}\sim \mathrm{Gamma}(\mathrm{shape}=\mathrm{50},\mathrm{rate}=\mathrm{5000}).\end{array}$$

Process model for n-alkane δ²H

The process model for n-alkane δ²H follows the same model structure outlined in Sect. 2.1.2 and is based on the same assumptions associated with n-alkane concentration distribution among the chains. This provides a concentration-weighted isotope mass balance equation similar to Eq. (5), for chain n and its δ²H_n,mix, as follows:

$$\begin{array}{}\text{(16)}& {\mathit{\delta }}^{\mathrm{2}}{\mathrm{H}}_{n,\mathrm{mix}}={\displaystyle \frac{\sum _i\left[{\mathrm{FLMC}}_i\sum \left({\mathrm{Conc}}_{n,i}\times {\mathit{\delta }}^{\mathrm{2}}{\mathrm{H}}_{n,i}\right)\right]}{\sum _i\left({\mathrm{FLMC}}_i\sum {\mathrm{Conc}}_{n,i}\right)}}.\end{array}$$

Data model and model inversion for n-alkane δ²H

The measured δ²H data are assumed to relate to the modeled δ²H_n,mix with a Gaussian error term:

$$\begin{array}{}\text{(17)}& {\mathit{\delta }}^{\mathrm{2}}{\mathrm{H}}_{n,\mathrm{mea}}\sim N\left({\mathit{\delta }}^{\mathrm{2}}{\mathrm{H}}_{n,\mathrm{mix}},{\mathit{\sigma }}_{n,\mathrm{mea}}^{\ast }\right),\end{array}$$

where ${\mathit{\sigma }}_{n,\mathrm{mea}}^{\ast }$ is the measurement error of δ²H_n,mea, which is often reported as 1 standard deviation (SD) in the literature.

In the model inversion step, MCMC is used to propose samples of all model parameters conditioned on the measured δ²H, δ¹³C, and RA data for multiple n-alkane chains. Five chains are run in parallel at 1 million iterations with 200 000 burn-ins. Average run time for the dual-isotope evaluation is approximately 8 h.

Data compilation

The sedimentary record of interest (GIK16160-3) is located off the Zambezi River mouth, which has a temporal range of ca. 0.1–36.7 ka BP (Wang et al., 2013a, b). The dominant biomes of the catchment area today are the Zambezian woodland savanna and grassland savanna, with coastal forest and Afromontane biomes also present (White, 1983; Dupont et al., 2011; Wang et al., 2013a; Dupont and Kuhlmann, 2017). Previous studies have interpreted the sedimentary record to reflect hydrological and vegetation changes in the catchment area of the Zambezi River, as well as patterns of sediment transport from the catchment vs. regions north of the river mouth (Wang et al., 2013a; Khon et al., 2014; van der Lubbe et al., 2014; Kasper et al., 2015; van der Lubbe et al., 2016; Lattaud et al., 2017). The published n-alkane record of this core comprises n-alkane concentrations; δ²H; and δ¹³C values of n-C₂₇, n-C₂₉, n-C₃₁, and n-C₃₃ alkanes (Wang et al., 2013a). Due to the high computational demand of the method, only 10 sets of data in the sedimentary sequence are used in this case study. The selected data span across late Marine Isotope Stage 3 (MIS 3), the Last Glacial Maximum (LGM), Heinrich Stadial event 1 (HS1), the Younger Dryas (YD), and the Holocene. Ice-volume-corrected δ²H values (Wang et al., 2013a) are used to account for baseline changes in δ²H of source moisture.

For the prior distributions of n-alkane concentration and δ¹³C, we use the same empirical dataset as in Sect. 2.2.2, while adding a distribution for n-C₂₇ alkane. To account for the post-industrial change in the δ¹³C of atmospheric CO₂, corrections are made to the δ¹³C means of the prior distributions, using the δ¹³C of atmospheric CO₂ of the year 2010 (−8.3  ‰ VPDB, Graven et al., 2017) and the pre-industrial values from an Antarctic ice core (Schmitt et al., 2012). For the prior distribution of ²H fractionation between plant wax and MAP (ε_a), location-specific plant wax δ²H data are compiled from the literature. We use a subset of the global compilation by Sachse et al. (2012) by setting the latitudinal limits to −30 and 30^∘ to reflect the tropical and subtropical latitudes of the catchment area of the Zambezi River. We also include more recently published plant n-alkane δ²H data (Kahmen et al., 2013b; Feakins et al., 2016; Griepentrog et al., 2019), with the same latitudinal limits. The δ²H of MAP is obtained from the Online Isotopes in Precipitation Calculator (OIPC version 3.1, Bowen and Revenaugh, 2003; IAEA/WMO, 2015; Bowen, 2022). The OIPC data product is chosen here because studies on modern plant wax lipid δ²H are often not associated with data of plant source water (soil or xylem water). Moreover, OIPC data have been widely applied to global models of plant ecology and paleoclimate and specifically to plant wax n-alkanes (e.g., Sachse et al., 2006; West et al., 2008; Sachse et al., 2012; Nusbaumer et al., 2017; Konecky et al., 2019). Since n-alkanes in marine sedimentary records often integrate climatic signals across large spatial areas and time intervals at the millennial scale, we consider the interpolated data of OIPC a reasonable set of prior information in the estimation of MAP δ²H. ε_a for each n-alkane chain in each data entry is calculated using Eq. (10). The same source groups (rainforest C₃ plants, savanna C₃ plants, C₄ plants, as in Sect. 2.2.2) are assigned to the data compilation based on either taxonomic identification or description of sample location from the original publications (supplementary data EA-5, Yang, 2022). Prior parameters of group-specific ε_a are estimated using the same methods as in Sect. 2.2.1 (“Data compilation”; Table S5). Univariate prior distributions of ε_a are plotted in Fig. S3.

2.3 Sensitivity tests

The model results are to some degree sensitive to the prior parameter estimates, which are derived from empirical data but imperfectly known in our case studies. Here we use data associated with case study 2 to explore the influence of prior parameter estimates on model output. To produce a different set of prior parameter estimates, plant samples are selected from western Africa only, as a subset of the sub-Saharan empirical dataset (123 entries out of the 301 entries from the original dataset, supplementary data EA-4, Yang, 2022). Prior parameters are estimated from this dataset using the same methods as described in Sect. 2.2.1 (“Data compilation”)​​​​​​​. The resultant prior distributions differ from those of the sub-Saharan dataset primarily in the C₄ source, in that they display lower n-alkane concentrations and less negative δ¹³C values (Fig. 4). Model sensitivity is assessed via a comparison between the model outputs based on the different priors.

Figure 4 Comparisons of prior distributions based on wester African plant samples (thick red curves) vs. sub-Saharan African plant samples (thin blue curves). The three left-most columns illustrate n-alkane concentrations. The three right-most columns illustrate n-alkane δ¹³C values.

3 Results

3.1 Case study 1

For each of the three data points of the Lake Qinghai case study (Fig. 5a–c), the posterior densities of FLMCs of terrestrial plants and aquatic macrophytes vary substantially between the samples (Fig. 5d–f), while the distributions of the algal FLMC are almost the same between the samples. FLMC from aquatic macrophytes is the highest in the high-P_aq sample (QHS13-5S, P_aq = 0.34, Table 1), with the MAPE at 0.53 (summary statistics in Table 3). In comparison, FLMC from aquatic macrophytes is the lowest (MAPE close to 0, Table 3) in the mid-P_aq sample (QHS13-7S, P_aq = 0.22, Table 1). FLMC from aquatic macrophytes is intermediate in the low-P_aq sample (QHS13-9S, P_aq = 0.19, Table 1) at 0.09 (MAPE, Table 3).

Figure 5 (a–c) Sedimentary long-chain n-alkanes (n-C₂₇, n-C₂₉, and n-C₃₁), their δ¹³C (‰, VPDB) and relative abundance values (excluding other chains) in published lake surface sediment samples from Lake Qinghai (Liu et al., 2015). (d–f) Posterior densities of fractional leaf mass contribution (FLMC) of terrestrial plants, aquatic macrophytes, and algae as model output conditioned on the δ¹³C and relative abundance values.

Table 3 The maximum a posteriori probability estimates (MAPE), the medians, and the 89 % highest density intervals (HDI) of posterior densities of the mixing fractions as model output using data in published lake surface sediment samples from Lake Qinghai (Fig. 5).

Figure 6 Bivariate density plots of the posterior densities of fractional leaf mass contribution (FLMC) of terrestrial plants, aquatic macrophytes, and algae in published lake surface sediment samples from Lake Qinghai (Liu et al., 2015).

Although the algal FLMC is not well constrained (Fig. 5), a strong “trade-off” correlation (the increase of one correlates with the decrease of another) between algal FLMC and those of the other two sources is apparent in bivariate density plots (Fig. 6). The high-P_aq (QHS13-5S) and low-P_aq (QHS13-9S) samples display multiple local modes for FLMCs, represented by multiple bivariate density peaks (as indicated by the contour lines and color gradient, Fig. 6).

Figure 7 Posterior densities of fractional source contribution (FSC_n) of terrestrial plants, aquatic macrophytes, and algae to each specific alkane chain (n-C₂₇, n-C₂₉ and n-C₃₁) in published lake surface sediment samples from Lake Qinghai (Liu et al., 2015).

The posterior densities of the algal FSC_n to the sedimentary long-chain n-alkane pool are consistently low (Fig. 7). In comparison, the FSC₂₇ of aquatic macrophytes can be substantial in both the high-P_aq (QHS13-5S) and low-P_aq (QHS13-9S) samples, while FSC₂₉ of aquatic macrophytes is substantial only in the high-P_aq sample (QHS13-5S). The terrestrial source shows consistently high FSC₃₁s across all samples.

3.2 Case study 2

For the western African transect case study (Fig. 8a–c), the posterior densities for fractional leaf mass contribution of C₄ plants, savanna C₃ plants, and rainforest C₃ plants vary substantially between samples (Fig. 8d–f). The C₄ plant contribution is the highest in the high-C₄ sample (Rhum, $f_{{\mathrm{C}}_{\mathrm{4}}}$ = 0.72, Table 2), at 0.59 (MAPE, Table 4). The FLMC of C₄ plants is lowest in the low-C₄ sample (Baro, $f_{{\mathrm{C}}_{\mathrm{4}}}$ = 0.05, Table 2) at 0.10 (MAPE, Table 4). FLMCs from the savanna C₃ and rainforest C₃ plants are generally less strongly constrained than those of C₄ plants, except for in the mid-C₄ sample (Asso), which shows a low FLMC from the savanna C₃ plants (MAPE close to 0, Fig. 8e, Table 4). The MAPEs of FLMCs of C₄ plants of the high-C₄ (Rhum) and the low-C₄ (Baro) samples fall within the reported mean ±1σ ranges for $f_{{\mathrm{C}}_{\mathrm{4}}}$ at these sites based on satellite imaging (Garcin et al., 2014, Table 2). The MAPE of FLMC of C₄ plants of the mid-C₄ sample (Asso) falls outside the reported $f_{{\mathrm{C}}_{\mathrm{4}}}$ mean ±1σ range.

Figure 8 (a–c) Sedimentary long-chain n-alkanes (n-C₂₇, n-C₂₉ and n-C₃₁), their δ¹³C (‰, VPDB), and relative abundance values (excluding other chains) in published lake surface sediment samples from Cameroon (Garcin et al., 2014). (d–f) Posterior densities of fractional leaf mass contribution (FLMC) of tropical C₄ plants, savanna C₃ plants, and rainforest C₃ plants as model output conditioned on the δ¹³C and relative abundance values.

Table 4 The maximum a posteriori probability estimates (MAPE), the medians, and the 89 % highest density intervals (HDI) of posterior densities of the fractional leaf mass contributions (FLMCs) as model output using data in published lake surface sediment samples from Cameroon (Fig. 8).

Bivariate density plots show relatively weak “trade-off” correlation patterns between the FLMCs of savanna and rainforest C₃ sources in the High C₄ sample (Rhum) and Low C₄ sample (Baro, Fig. 9). All three samples display a single mode of most likely vegetation contributions, represented by a single peak in the bivariate density plots (as indicated by the contour lines and color gradient, Fig. 9), although a secondary mode with higher rainforest and lower savanna C₃ FLMCs exists at Rhum.

Figure 9 Bivariate density plots of the posterior densities of fractional leaf mass contribution (FLMC) of tropical C₄ plants, savanna C₃ plants, and rainforest C₃ plants in published lake surface sediment samples from Cameroon.

Figure 10 Posterior densities of fractional source contribution (FSC_n) of terrestrial plants, aquatic macrophytes, and algae to each specific alkane chain (n-C₂₉, n-C₃₁ and n-C₃₃) in published lake surface sediment samples from Cameroon.

The distributions of FSC_n of all three sources to each specific chain (n-C₂₉, n-C₃₁ and n-C₃₃) vary substantially between samples. Within each sample, the rainforest C₃ source consistently shows a higher FSC₂₉ value over FSC₃₁ and FSC₃₃ (Fig. 10). In comparison, the C₄ source shows a lower FSC₂₉ value over FSC₃₁ and FSC₃₃ with the exception of the high-C₄ sample (Rhum), in which the C₄ FSC₂₉ and FSC₃₁ are almost identical (Fig. 10). Within each sample, the uncertainties of FSC_n vary among the chains, with the high-C₄ sample (Rhum) exhibiting the highest overall uncertainty in FSC₂₉ and FSC₃₁, and the low-C₄ sample (Baro) exhibiting the highest overall uncertainty in FSC₃₃. The high-C₄ sample (Rhum) also shows bimodal distributions in FSC₂₉ and FSC₃₁ for the rainforest C₃ and savanna C₃ sources. Among the samples, the FSC_n uncertainty of the C₄ source is the lowest.

3.3 Case study 3

The MAPEs of rainforest C₃ FLMC are consistently low, ranging from 0.01 to 0.25 in the GIK 16160-3 sequence (Fig. 11). MAPEs of savanna C₃ FLMC are consistently high, ranging from 0.51 to 0.80. MAPEs of C₄ FLMC range from 0.17 to 0.39. MAPEs of C₄ FLMC were relatively stable at ∼ 0.4 from the end of MIS3 to the LGM and gradually decreased to ∼ 0.2 through HS1 and the YD entering the early Holocene. By contrast, MAPEs of rainforest C₃ FLMC gradually increased within the same time span to its peak value at ∼ 0.25 but experienced an abrupt drop by the end of HS1 to a consistently low value through the YD and the Holocene. MAPEs of savanna C₃ FLMC were relatively steady at ∼ 0.55 from MIS 3 to HS1 but experienced an abrupt increase to ∼ 0.8 by the end of HS1 and remained relatively stable through the YD and the Holocene.

Figure 11 (a) Model estimates of fractional leaf mass contribution (FLMC) of three broadly defined biomes through time, based on n-alkane data of the GIK 16160-3 core (Wang et al., 2013a). Solid lines represent maximum a posteriori estimates (MAPE) of the three sources. Dashed lines represent the 89 % highest density interval (HDI) of the posterior distributions. (b) The n-alkane relative abundance and (c) n-alkane δ¹³C of the same core (Wang et al., 2013a). (d) Global sea level change (Lambeck et al., 2014). (e) BIT index indicating terrestrial sedimentary input (Kasper et al., 2015). (f) The εNd of clay indicating relative influence of the Zambezi River (van der Lubbe et al., 2016). YD represents the Younger Dryas (∼ 13–11 ka). HS1 represents Heinrich Stadial 1 (∼ 18–15 ka). LGM represents the Last Glacial Maximum (∼ 26–19 ka).

Figure 12 (a) Model estimates of mean annual precipitation δ²H based on n-alkane δ²H (ice-volume-corrected data, IVC), δ¹³C, and relative abundance of four n-alkane chains from the GIK 16160-3 core (Wang et al., 2013a). Solid lines represent maximum a posteriori estimates (MAPE) of the three sources. Horizontal dashed lines represent the 89 % highest density interval (HDI) of the posterior distributions. (b) The δ²H (IVC) of n-C₂₉ and n-C₃₁ chains reported by Wang et al. (2013a) for reference. YD represents the Younger Dryas (∼ 13–11 ka). HS1 represents Heinrich Stadial 1 (∼ 18–15 ka). LGM represents the Last Glacial Maximum (∼ 26–19 ka).

The MAPEs of MAP δ²H range from −42 ‰ to −30 ‰, broadly consistent with the range of the modern MAP δ²H record (−37 ± 15 ‰) at the GNIP Ndola station (lat. −13.0000^∘, long. 28.6500^∘), Zambia (IAEA/WMO, 2022). From the end of MIS 3 to the LGM, MAP δ²H decreased from a local peak value of −35 ‰ at 30 ka to the lowest value of −42 ‰ at 18.6 ka (Fig. 12a). During the glacial–interglacial transgression, MAP δ²H briefly spiked to −32 ‰ during HS1 and immediately dropped to a local low value of −40 ‰ between HS1 and the YD, and it then went up again to −34 ‰ during the YD. The overall range of estimated MAP δ²H (∼ 12 ‰) is also similar to that of the n-C₂₉ alkane but is smaller than that of the n-C₃₁ alkane (Fig. 12b). The estimated MAP δ²H sequence is positively correlated with δ²H of both n-C₂₉ and n-C₃₁ alkanes (p values = 0.0206 and 0.0387, respectively, Fig. 12b), although n-C₂₉ shows a higher correlation coefficient (r = 0.713) than n-C₃₁ (r = 0.658). The overall uncertainty (range of the 89 % HDI) of estimated δ²H of MAP is ∼ 13 ‰ (Fig. 12a), which is much smaller than the parametric uncertainty (−1σ to 1σ, equivalent to a range of 68 % HDI assuming a normal distribution) of ε_a for, e.g., the n-C₂₉ alkane of trees at > 40 ‰ (Sachse et al., 2012; Konecky et al., 2019).

3.4 Model sensitivity

3.4.1 Sensitivity to parameters of the prior distributions

Using the western Africa prior dataset (Fig. 4), with lower n-alkane concentration and less negative δ¹³C values in the C₄ source, resulted in consistent shifts to higher values of FLMC of the C₄ source relative to the results obtained with the sub-Saharan prior dataset (Fig. 13). With the published data on fractional C₄ vegetation cover as a reference (Garcin et al., 2014), this set of priors produces FLMC of the C₄ source similar to the reported fractional C₄ vegetation cover ($f_{{\mathrm{C}}_{\mathrm{4}}}$) at the high-C₄ site (Rhum) but produces much higher FLMCs than published $f_{{\mathrm{C}}_{\mathrm{4}}}$ values at the mid- and low-C₄ sites (Asso and Baro).

Figure 13 Comparisons of posterior densities of fractional leaf mass contribution of the C₄ source (FLMC${}_{{\mathrm{C}}_{\mathrm{4}}}$) based on two sets of prior distributions: Prior 1 (solid lines) using sub-Saharan Africa dataset and Prior 2 (dashed lines) using the western Africa dataset. Detailed prior distributions are illustrated in Fig. 4, and $f_{{\mathrm{C}}_{\mathrm{4}}}$ values are from Garcin et al. (2014).

3.4.2 Sensitivity to proxy type

The model shows different sensitivity to chain length distribution and carbon stable isotopes. It is relatively insensitive to chain length distribution: only minor changes in the posterior densities are observed when the likelihood evaluations of RA are removed (left column, Fig. 14). By contrast, the model is much more sensitive to the n-alkane δ¹³C data: the central tendencies of the posterior densities shifted substantially when the likelihood evaluations are removed (right column, Fig. 14). Substantial increases in the dispersion of the posterior densities are also observed, indicating more weakly constrained FLMCs (right column, Fig. 14).

Figure 14 Model sensitivity to proxy type in model inversion (Eqs. 6 and 8), using case study 2 as an example. The right column (No RA) shows the model output with the likelihood evaluations of RA completely removed from model evaluation. The left column (No δ¹³C) shows the model output with the likelihood evaluations of δ¹³C completely removed from model evaluation.

4 Discussion

4.1 Interpretations of results

4.1.1 Case study 1

Long chain n-alkanes (especially 27–31 carbons) in lake sediments are possibly a mixture of both terrestrial and aquatic sources (e.g., Aichner et al., 2010; Duan and Xu, 2012; Wang and Liu, 2012; Liu and Liu, 2016; Dion-Kirschner et al., 2020). Based on model output, both FLMC and FSC₂₇ of aquatic macrophytes are not negligible in the low-P_aq sample (QHS13-9S, Figs. 5 and 7). By contrast, FLMC and FSC₂₇ of aquatic macrophytes are minimal in the mid-P_aq sample (QHS13-7S, Figs. 5 and 7). These results suggest that P_aq alone may not be a reliable indicator of the relative aquatic macrophyte input into the sedimentary archive. In addition, for the high-P_aq sample (QHS13-5S), FSC₂₉ of aquatic macrophytes is not negligible, while FSC₃₁ of the terrestrial source is consistently close to 1 among the samples (Fig. 7). This suggests caution should be taken when interpretations of terrestrial paleoenvironment are based on the isotope ratios of n-C₂₉ alkane alone, while previous interpretations based on the n-C₃₁ alkane are likely more reliable.

The relatively unconstrained FLMC of the algae source (Fig. 5) and its “trade-off” with that of the terrestrial source (Fig. 6) suggest that the possibility of algae being an important biomass source of the lake surface sediment cannot be eliminated, and would impact the interpreted contributions from other sources. This possibility is consistent with the observations that the deep-water lake bottom is covered mainly by green algae in Lake Qinghai (Liu et al., 2015). The model approach successfully identified such a possibility, despite the consistently low FSC₂₇ of the algal source (Fig. 7), most likely due to its limited production of the n-C₂₇ alkane (Liu et al., 2015). Algae produce greater amounts of short-chain n-alkanes than mid-chain and long-chain homologues (Han and Calvin, 1969; Gelpi et al., 1970; Cranwell et al., 1987). Theoretically, the inclusion of n-C₂₅ or even shorter chains would provide better constraints on both the aquatic macrophyte and the algae sources. However, this option proved impractical due to the many missing values of n-C₂₅ in the empirical records of terrestrial plants. The other concern would be that mid-chain n-alkanes are also produced by microorganisms (Grimalt et al., 1988; Park, 2005; Ladygina et al., 2006; Brittingham et al., 2017), the influence of which is difficult to assess when only odd-chain alkanes are analyzed.

The Tibetan Plateau has been argued to be an ideal region to investigate the input from aquatic macrophytes to the organic matter and n-alkane pool in lake sediments, due to the minimal presence of terrestrial C₄ plants in the region (Wang, 2003; Aichner et al., 2010; Duan and Xu, 2012; Liu et al., 2015). While the carbon isotope values of these two n-alkane sources largely overlap (Aichner et al., 2010; Liu et al., 2015), they display different chain length dominance patterns: aquatic macrophytes produce relatively high proportions of the n-C₂₇ alkane, while terrestrial C₄ plants produce relatively high proportions of the n-C₃₁ alkane (Rommerskirchen et al., 2006; Bush and McInerney, 2013; Badewien et al., 2015; Liu et al., 2015; Magill et al., 2019). The consistently low δ¹³C values of the n-C₃₁ alkane among the samples confirm the minimal contribution from terrestrial C₄ plants. In a hypothetical situation where both aquatic macrophytes and terrestrial C₄ plants are present, increased aquatic macrophyte contribution would increase the δ¹³C value of the n-C₂₇ alkane more than the n-C₂₉ alkane, while increased terrestrial C₄ contribution would increase the δ¹³C value of the n-C₃₁ alkane more than the n-C₂₉ alkane. The proposed modeling approach is theoretically capable of distinguishing the n-alkane input from aquatic macrophytes and terrestrial C₄ plants by leveraging both chain length distribution and δ¹³C values of these n-alkane chains. While we do not have a case study to clearly demonstrate such a possibility, future studies are encouraged to explore this further.

4.1.2 Case study 2

This case study demonstrates that the proposed framework can leverage chain length distribution and δ¹³C of multiple n-alkane chains and offer interpretations on the mixing ratios of multiple sources. The mid-C₄ sample (Asso) shows a high relative abundance of n-C₂₉ alkane and contrasting δ¹³C values among the n-alkane chains (Fig. 8b, which suggests that the most likely mixing regime of the sample is a much higher FLMC of rainforest C₃ plants than savanna C₃ plants. For the high-C₄ sample (Rhum) the model also recovered bimodal mixing possibilities between the rainforest and savanna C₃ sources (Fig. 10b1). Such information on the possible mixing patterns of the C₃ biomes cannot be achieved with the traditional two-end-member mixing model using n-alkane δ¹³C values alone.

FLMC as a metric for vegetation reconstruction is not directly comparable with $f_{{\mathrm{C}}_{\mathrm{4}}}$ estimated from remote sensing imagery (e.g., MODIS Vegetation Continuous Fields). A carefully selected calibration dataset could support a rigorous comparison of the two metrics. On the other hand, assuming that $f_{{\mathrm{C}}_{\mathrm{4}}}$ estimated from satellite imagery is reflecting the true vegetation distribution contributing to the study samples, the differences between FLMC and $f_{{\mathrm{C}}_{\mathrm{4}}}$ could inform possible biases associated with FLMC. As presented, the central tendencies of the C₄ FLMC at the low- and mid-C₄ sites are higher than $f_{{\mathrm{C}}_{\mathrm{4}}}$, while C₄ FLMC at the high-C₄ site is lower than $f_{{\mathrm{C}}_{\mathrm{4}}}$ (Fig. 8). Such a pattern of biases cannot be explained by a higher leaf mass to vegetation cover ratio of the C₃ biomes. One possible explanation is with a spatial pattern of n-alkane sourcing and integration. Tropical African lake shores often host a gradient of mixed vegetation from C₄ grasses and sedges to C₃ trees and shrubs (Vesey-Fitzgerald, 1963; Greenway and Vesey-Fitzgerald, 1969; Howard-Williams and Walker, 1974). If n-alkanes in lake sediments are sourced consistently more from the immediate lakeshore region than from further away, the sediments would more likely integrate a mixed vegetation signal than the surrounding regions, which is consistent with the observed differences between FLMC and $f_{{\mathrm{C}}_{\mathrm{4}}}$.

4.1.3 Case study 3

The model recovered a prominent shift in the vegetation source of the n-alkane record from the LGM to the YD: a gradual increase in rainforest FLMC to its peak during HS1 and an abrupt decrease by the end of HS1 (Fig. 11a). This pattern is driven by both n-alkane chain length distribution and δ¹³C of the sequence, with an n-C₂₉ dominance and low δ¹³C values for the n-C₂₉ and the n-C₃₁ alkanes within this time interval (Fig. 11b, c; Wang et al., 2013a). This time interval coincided with the lowest sea level (Fig. 11d; Lambeck et al., 2014) and a high terrestrial organic matter input as indicated by the BIT index (Fig. 11e; Kasper et al., 2015). Neodymium isotopes (εNd, ‰), an indicator of sedimentary clay source, suggest that the Zambezi outflow was the dominant source at the core location during the LGM and HS1, but the outflow from small coastal rivers draining the highlands north of the Zambezi River mouth was relatively high during the early Holocene (Fig. 11f; van der Lubbe et al., 2016). The correlations between n-alkane chain length distribution, δ¹³C, sea level, and other sedimentary source proxies (see also Just et al., 2014; van der Lubbe et al., 2014; Lattaud et al., 2017) suggest that the source area of n-alkanes at the core location may have changed substantially. Variations in plant wax δ²H in the tropics have been interpreted to primarily reflect variations in local rainfall amount (e.g., Tierney et al., 2008; Schefuß et al., 2011; Berke et al., 2012; Collins et al., 2014; Costa et al., 2014; Shanahan et al., 2015; Castañeda et al., 2016; Niedermeyer et al., 2016; Norström et al., 2018). However, with evidence that supports a substantial change in plant wax sourcing area, alternative explanations for the n-alkane δ²H variation should be considered.

Prior to 15 ka BP, the contribution of Zambezi River sediments to the core location was relatively stable (Fig. 11d, e). The transient rise in rainforest FLMC during the deglaciation is likely associated with alkanes sourcing from near-coastal forests. At the peak of rainforest FLMC during HS1, the estimated MAP δ²H reached its highest point, which is consistent with expectations that coastal rainfall, proximal to ocean moisture sources, would have high δ²H values (Dansgaard, 1964; Rozanski et al., 1993). Immediately after 15 ka BP, coincident with abrupt sea level rise (Fig. 11d), rainforest FLMC dropped to near zero and savanna FLMC became dominant. At the same time, a transient drop in estimated MAP δ²H is consistent with the lower values expected for inland moisture contributing to n-alkanes mobilized from the savanna biome. This pattern aligns well with interpretations of the changing properties of terrigenous sediments in the nearby GeoB9307-3 core (Just et al., 2014). Going into the early and mid-Holocene with rising sea level and warmer conditions, there is an overall rise in MAP δ²H leading to the mid-Holocene peak, which can be explained by a gradual increase in moisture availability in the Zambezi catchment. This interpretation is also supported by patterns of sedimentary transport (Just et al., 2014; van der Lubbe et al., 2016). As such, when interpreting variations in plant wax δ²H record from a marine core, it is important to consider moisture source and transport history as alternative mechanisms to temporal change in rainfall isotope values at fixed locations.

4.1.4 Sensitivity tests

The results of sensitivity test 1 (Fig. 13) show that using a different set of prior distributions (Fig. 4) can produce somewhat different central tendencies in FLMCs with the same sedimentary n-alkane data. This demonstrates that the outputs are model dependent: they provide a basis for interpreting proxy data in the context of a specific model with its associated priors and assumptions. Using the published $f_{{\mathrm{C}}_{\mathrm{4}}}$ values as a reference, Prior 2 (the western Africa prior) is associated with overall higher estimated C₄ FLMCs than those with Prior 1 (the sub-Saharan Africa prior, Fig. 13). This unidirectional shift in FLMCs is consistent with the expectation of a systematic difference between the two prior distributions. Based on the comparison with satellite-based vegetation cover, this further suggests that Prior 1 is perhaps a better representation of the vegetation sources that produced the n-alkane mixtures in the sedimentary samples than Prior 2 (given the caveats of this comparison, discussed above).

The results of sensitivity test 2 (Fig. 14) show that δ¹³C places stronger constraints on the results than chain length distribution. This demonstrates the importance of stable isotopes as tracers in compound-specific mixing models. In contrast to other carbon-isotope based methods, the model proposed here might still provide reasonable resolution of sources when the chain-length distribution of the sources is substantially different.

4.2 Model achievements

The proposed model framework offers a more integrative approach to interpretation of n-alkane records that could better leverage these rich data to support more nuanced and specific reconstructions of paleo-vegetation and paleoclimate. It has potential advantages over traditional approaches to the interpretation of sedimentary n-alkanes.

First, the model offers a numerical solution (Sect. 2.1.2) that accounts for the uncertainty associated with n-alkane production among the sources (Figs. 2 and 3), which is coupled with the mass balance equation of n-alkane δ¹³C (Eq. 5). This allows δ¹³C to be used to constrain the uncertainty associated with n-alkane chain length distributions, which have previously been used in isolation to reconstruct vegetation (Jansen et al., 2010; Gao et al., 2011; Peaple et al., 2021). The model also addresses a common assumption in the interpretation of lipid δ¹³C data via linear mixing relationships, in which n-alkane production for a selected chain is assumed to be the same among all sources (Bush and McInerney, 2013; Garcin et al., 2014).

Second, the case studies demonstrate that the model can be used to explore mixing regimes of multiple sources, which offers alternative interpretations compared to the traditional two-end-member mixing regime using δ¹³C of one n-alkyl lipid. Case studies 1 and 2 showcase the flexibility of the model framework since the vegetation sources and their prior distributions are user-defined. They also highlight the strengths of the Bayesian framework in disentangling multivariate mixing scenarios in plant wax lipids. The multi-source framework can also be used to assist the interpretations of associated n-alkane δ²H as demonstrated in case study 3.

Moreover, the new metric FLMC has the potential to evaluate leaf mass integration patterns in sedimentary archives. The molecular distribution and δ¹³C of n-alkanes have been used extensively in the assessment of sedimentary organic matter input (e.g., Wiesenberg et al., 2004; Mead et al., 2005; Aichner et al., 2010; Kristen et al., 2010; Seki et al., 2010; Hockun et al., 2016; Ankit et al., 2017; Liu et al., 2020). Because n-alkane concentrations (µg g⁻¹) involved in the mixing model are ultimately derived from dried leaf mass (Sect. 2.1.2), the model can potentially make an explicit connection between FLMC (as in Figs. 5 and 8) and patterns of organic matter sourcing. Although the model presented here does not yet support these interpretations, addition of components to enable e.g., parameterizing of biome-specific primary productivity would be relatively straightforward in future applications.

Lastly, by leveraging compound-specific stable isotope data from multiple n-alkane chains, the model framework is capable of filtering out noise from true environmental signals. Many studies use plant wax δ²H data of one selected chain to estimate source water δ²H. However, the reported uncertainty of the source water δ²H estimation can be miscalculated (Polissar and D'Andrea, 2014). The advantage of the Bayesian framework is in its straightforward propagation of uncertainties (Tingley et al., 2012) in the mixing process of multiple plant wax compounds. For example, in case study 3, the relatively narrow uncertainty bounds of the estimated MAP δ²H values (Fig. 12a) demonstrates that the model framework is less sensitive to the parametric uncertainty associated with any one n-alkane chain than conventional approaches. This can reduce the risks of potentially interpreting data anomalies as environmental signals.

4.3 Future directions

The proposed mixing model is the bare core of a potentially more comprehensive proxy system model, an approach that has gained traction in recent developments of paleoenvironmental reconstruction (e.g., Garreta et al., 2010; Li et al., 2010; Tingley et al., 2012; Evans et al., 2013; Dee et al., 2018; Konecky et al., 2019; Bowen et al., 2020). A proxy system model is a representation of the complete proxy system that ideally includes four components: the environment, the sensor, the archive, and the observation (Evans et al., 2013). The mixing model here primarily describes the sensor and observation components of the complete proxy system, while the environment and the archive components have not been incorporated (see Dee et al., 2015, 2018; Konecky et al., 2019 as examples for other model components). Future efforts that elaborate on the model structure will provide updated model assumptions, which should be based on systematic investigations of the specific proxy system components. Here are three categories of model improvements to consider.

4.3.1 Improvements on prior distributions

A better characterization of n-alkane chain length distribution and δ¹³C in modern plants can potentially expand model application, due to the prior distributions' reliance on empirical data. At the same time, substantial gaps in our knowledge of n-alkanes still exist, especially from certain underrepresented regions of the world, such as western, central, and southeastern Asia, as well as certain biomes, such as Afromontane and coastal forests (see global compilations by, e.g., Bush and McInerney, 2013; Diefendorf and Freimuth, 2017). Future studies should aim at bridging the gaps with systematic surveys of the concentrations, δ¹³C, and ε_a of n-alkanes and potentially other plant wax compounds in modern plants. The need for a centralized database for plant wax lipids is also worth mentioning, as it will facilitate easy access to published records of both lipid concentrations and stable isotope values. From a different perspective, the development of mechanistic ways to represent the prior covariance patterns in concentrations, δ¹³C, and ε_a of n-alkanes can reduce the model's reliance on complete data records and improve its tolerance of missing values. These efforts can improve the flexibility of model application and the robustness of model results.

4.3.2 Better understanding of the processes and the incorporation of additional tracers

Processes such as n-alkane turnover and transportation have not been incorporated into the model framework primarily because they are still poorly understood (Sachse et al., 2012; Diefendorf and Freimuth, 2017). For example, factors associated with plant growth, such as primary productivity, leaf lifespan, and leaf turnover patterns (e.g., Burnham, 1989; Greenwood, 1991; Hauke and Schreiber, 1998; Jetter and Schäffer, 2001; Ellis and Johnson, 2013; Tipple et al., 2013; Suh and Diefendorf, 2018; Suh et al., 2019), can influence the overall n-alkane turnover rate of a specific plant source. The n-alkane transportation can be another source of systematic biases in the interpretation of n-alkane records (e.g., Schefuß et al., 2004; Yamamoto et al., 2013; Nelson et al., 2018; Freimuth et al., 2019; McFarlin et al., 2019; Suh et al., 2019). Catchment-specific integration effects such as spatial and temporal lipid integration (see case study 3) can also affect the interpretation of n-alkane records (e.g., Schefuß et al., 2004; Seki et al., 2010; Vogts et al., 2012; Douglas et al., 2014; Garcin et al., 2014; Herrmann et al., 2016; French et al., 2018; Freimuth et al., 2019, 2021). Understanding these processes would ideally involve additional tracers or proxies. Future studies that investigate n-alkane turnover patterns may consider tracking the total amount of n-alkanes produced and its relationship with the primary productivity of the plants in a year. The incorporation of, e.g., pollen, sedimentary tracers, and multiple isotopes (e.g., Seki et al., 2010; Feakins, 2013; Just et al., 2014; Kasper et al., 2015; van der Lubbe et al., 2016; Feakins et al., 2018; Freimuth et al., 2019), as process indicators can help to inform the relative importance of different transportation processes.

4.3.3 Incorporation of environmental factors into the proxy system

Environmental factors, such as temperature and precipitation, are primary determinants of vegetation composition, but they also influence n-alkane concentration, δ¹³C and δ²H (sensor values) through plant physiological response. The inclusion of environmental factors into the model framework would naturally involve n-alkane δ²H in a more mechanistic way (e.g., Liu and Yang, 2008; Sachse et al., 2012; Sessions, 2016; Liu and An, 2018). For example, ε_a of n-alkanes can be partitioned into a biosynthetic fractionation term (ε_bio) and a evaporative enrichment term (ε_eva). The multi-level influence from environmental factors to, e.g., ε_a of n-alkanes, include a relatively static element, such as a taxon-specific biosynthetic fractionation term (ε_bio) (although, see Newberry et al., 2015, for counterexamples; Cormier et al., 2018), and dynamic elements such as environment-dependent evaporative enrichment terms (ε_eva) for leaf water and soil water (e.g., Douglas et al., 2012; Sachse et al., 2012; Kahmen et al., 2013a; Tipple et al., 2015; Cernusak et al., 2016; Sprenger et al., 2016; Liu and An, 2018; Konecky et al., 2019). Future studies should consider developing such model components, which will ultimately permit proxy system interpretation at the environmental level.

5 Conclusion

Traditional interpretations of n-alkane proxies often rely on either chain length distribution or δ¹³C values of one chain to reconstruct vegetation composition. An alternative approach is to combine the information from both lines of evidence to refine our interpretations. We presented a Bayesian modeling approach that simultaneously evaluates both chain length distribution and stable carbon isotope measurements of sedimentary n-alkanes. The model incorporates the uncertainties associated with n-alkane concentration and carbon isotope ratios of multiple chains, using a generic mixing process with isotope mass balance. We presented three case studies to illustrate how the model can be applied to the interpretation of sedimentary n-alkane records. The first one involves published long-chain n-alkane records from lake surface sediments of Lake Qinghai, China. The second one involves published long-chain n-alkane records from lake surface sediments along a vegetation gradient in Cameroon, western Africa. The third one involves published long-chain n-alkane records with associated δ²H data from a marine core off the Zambezi River mouth. The first two case studies demonstrate that compared to traditional two-end mixing models, our approach can resolve mixing fractions (FLMC) of multiple n-alkane sources, providing alternative interpretations of vegetation composition to the same n-alkane record. With biome-specific ε_a component added to the framework, the third case study demonstrates that the model is capable of estimating precipitation δ²H with a reasonable range of uncertainty by leveraging δ²H of multiple n-alkane chains. Despite these achievements, several processes associated with n-alkane integration in sedimentary archives are still not accounted for. Future studies on n-alkane turnover and transportation will improve our understanding of the biases and constraints associated with n-alkane records. The Bayesian model framework could be further improved by adding more mechanistic modeling components such as n-alkane integration processes and additional environmental processes involving δ²H. The modeling approach represents a continuously evolving framework that can incorporate new understandings and leverage additional proxies in the future.