Resolving the Interpretation of Magnetic Coercivity Components From Backfield Isothermal Remanence Curves Using Unmixing of Non‐Linear Preisach Maps: Application to Loess‐Paleosol Sequences

Unmixing of remanent magnetization curves, either isothermal remanent magnetization (IRM) or backfield IRM, is widely used in rock magnetic and environmental magnetic studies to discriminate between magnetic coercivity components of different origins. However, the wide range of physical properties of natural magnetic particles gives rise to an ambiguous interpretation of these components. To reduce this ambiguity and provide a straightforward interpretation of coercivity components in terms of domain state, interactions, and constituent magnetic phases, we combined backfield IRM unmixing with unmixing of nonlinear Preisach maps for two typical mid‐latitude northern hemisphere loess‐paleosol sequences. Both backfield IRM and nonlinear Preisach maps unmixing are based on the same non‐parametric algorithm, and provide similar endmembers (EMs) in the two sections studied. The first EM (EM1) has a low median coercivity (∼21 mT) and is a non‐interacting single domain (SD) magnetite/maghemite of pedogenic origin. The second EM (EM2) has a moderate median coercivity (∼60 mT) and is a mixture of pseudo‐single domain/multidomain, SD magnetite/maghemite and non‐interacting SD hematite, all of eolian origin. The same EM1 found in both sections suggests that this component's grain size and coercivity are independent of pedogenesis intensity. The same EM2 indicates that a similar magnetic population is being transported and deposited, irrespective of the dust source area and loess granulometry. The approach outlined here provides strong evidence that non‐parametric backfield IRM unmixing isolates physically realistic EMs. Unmixing nonlinear Preisach maps elucidates these EMs in terms of domain states and their constituent magnetic phases.

Two types of approaches have been used to differentiate between various types of magnetic components.The first approach is based on a parametric method in which the remanent curve (either isothermal remanent magnetization (IRM) or backfield IRM), or its first derivative (known as the coercivity distribution), is fitted by a linear combination of individual component distributions modeled by known mathematical functions, such as lognormal (Kruiver et al., 2001;Robertson & France, 1994;Stockhausen, 1998), skewed generalized Gaussian (Egli, 2003(Egli, , 2004a(Egli, , 2004b)), skewed normal (SN) (Maxbauer et al., 2016b), or Burr type XII (Zhao et al., 2018) distributions.The parametric method has several drawbacks.First, the unmixing results depend strongly on the shape of basis functions chosen to model the data (e.g., Egli, 2021;Li et al., 2020;Zhao et al., 2018).Second, the quality of measurements and the initial guess for the individual components highly influence the results of the IRM unmixing (e.g., He et al., 2020;Zhao et al., 2018).Third, the overlap of coercivity spectra adds to the ambiguity of determining individual components and their number in the fitting process (Egli, 2021;He et al., 2020).Moreover, micromagnetic simulations have shown that the IRM curves for an assemblage of noninteracting magnetic particles with a simple grain size distribution need to be fitted by multiple Gaussian functions, due to the distribution of grain elongation and the different types of switching of particles with different shapes (Bai et al., 2021).
The second approach was proposed by Heslop and Dillon (2007) and uses a dimensionality reduction nonparametric method, namely non-negative matrix factorization (NMF).This method determines the shape of the coercivity distributions and their contributions by using only the measured data.The main advantage of the non-parametric method is that it overcomes the problems associated with selecting the number of components and types of mathematical functions to fit.Several studies have used this method to identify the main processes responsible for magnetic properties of loess-paleosol deposits (Necula et al., 2013(Necula et al., , 2015;;Nie et al., 2014), soils (Hu et al., 2020), lacustrine sediments (Fabian et al., 2016;Nie et al., 2017), or marine sediments (Just et al., 2012), as well as to diagnose remagnetization in natural samples (Dekkers, 2012).The disadvantages of this method are that the endmembers (EMs) provided may not represent single components (He et al., 2020;Heslop, 2015;Just et al., 2012), and that if the overlap of the coercivity distributions is significant, the recovery of the pure EMs and their corresponding abundances may fail (He et al., 2020;Heslop, 2015).
Several studies (e.g., He et al., 2020;Hu et al., 2021;Maxbauer et al., 2017;Roberts et al., 2018) have used a complementary approach, combining results from IRM unmixing with another non-parametric unmixing method, principal component analysis of first-order reversal curves (FORCs) (Harrison et al., 2018;Lascu et al., 2015).FORC-PCA has proven to be a powerful method that can discriminate between signatures of domain states, magnetostatic interactions, or mineral phases, providing an improved understanding of magnetic particle mixtures.However, FORCs are in-field measurements, and, in addition to remanence contributions, also record transient, induced and thermal relaxation phenomena (e.g., Harrison et al., 2019;Zhao et al., 2017).Therefore, a direct link between components resulting from IRM unmixing and domain states revealed by FORC-PCA unmixing is difficult to establish (e.g., Egli, 2021;Zhao et al., 2017).
Here we explore whether IRM coercivity components reflect one or multiple domain states by using a nonparametric unmixing method of individual backfield curves combined with a non-parametric unmixing of nonlinear Preisach maps (Church et al., 2016).Nonlinear Preisach maps are built from a sequence of backfield IRM curves, each starting from a different conditioning field, covering the entire remanence space.Therefore, integration of the Preisach function over this space is equivalent to an IRM curve measured from an initially AF demagnetized state (Church et al., 2016).In addition, no smoothing is performed on the Preisach maps, thus avoiding the artifacts associated to this type of post-measurement processing.Furthermore, nonlinear Preisach maps perform better than FORC diagrams for complex compositions where the remanent signal is masked by non-remanent contributions (Church et al., 2016).
We test this combined approach on two Romanian loess-paleosol sequences chosen to reflect different climatic conditions and different source areas, one in the Southern Wallachian Plain, and the other on the Black Sea shore (Necula et al., 2015).These sites are typical for mid-latitude northern hemisphere loess deposits, known as "Chinese loess type" (Marković et al., 2015).Magnetic properties of loess-paleosol sequences from Southeastern Europe record glacial-interglacial cyclic climatic changes on a global scale, allowing the reconstruction of past climatic conditions for the last 1 million years (Marković et al., 2015;Obreht et al., 2019;Song et al., 2018;Zeeden et al., 2017).Rock magnetic studies have shown that Southeastern European loess-paleosol deposits, like their Chinese counterparts, are characterized by so-called "magnetic enhancement" in the paleosols, because of pedogenesis during interglacials (Buggle et al., 2014;Fitzsimmons et al., 2012;Jordanova & Petersen, 1999;Marković et al., 2015;Necula et al., 2013Necula et al., , 2015;;Panaiotu et al., 2001;Zeeden et al., 2016).Loess deposits from Serbia, Bulgaria and Romania have demonstrated that this magnetic enhancement is mainly due to the presence of fine-grained, superparamagnetic and single domain (SD) ferrimagnetic minerals (magnetite and/or maghemite) produced in situ during wet and warm climate conditions.During glacials, when it is generally dry and cold, the accumulation of aeolian dust prevails, resulting in a significant input of coarse detrital magnetite in the pseudosingle domain (PSD) to multidomain (MD) range (e.g., Jordanova & Jordanova, 2020;Laag et al., 2021;Necula et al., 2015).These domain states are inferred either from bulk magnetic parameters measurements, or from IRM unmixing coercivity distributions, and typically suffer from ambiguity in interpretation (Egli, 2021).The dual IRM-nonlinear Preisach maps approach proposed here provides a straightforward interpretation of IRM unmixing EMs in terms of domain state, magnetostatic interactions, and composition, casting a new light on the populations of grains controlling the magnetic properties characteristic of the Chinese loess type.

Study Sites
The first loess-paleosol sequence is in the Southern Wallachian Plain near the village of Lunca, on the left bank of the Sâi River, approximately 15 km north of its confluence with the Danube River (43˚50′56″N, 24˚45′56″E) (Figure 1).The section is about 34 m-thick and contains six loess layers (L1-L6), the recent soil (S0) and six intercalated paleosols (S1-S6).The section is continuous, without visible erosion or hiatuses, and is the longest loess-paleosol sequence from Romania accessible in outcrop.A chronostratigraphy of the Lunca deposit was established using optically-stimulated luminescence and a magnetic age-depth model (Constantin et al., 2015).The luminescence ages obtained on fine and coarse quartz have shown that loess L1 developed during marine isotopic stages (MIS) 2-4.The magnetic age-depth model based on the correlation of a preliminary magnetic susceptibility record with the stack of 57 globally distributed benthic δ 18 O records (Lisiecki & Raymo, 2005) is in agreement with the luminescence ages, indicating that palaeosols S1-S4, and S6 can be correlated with the interglacials in the interval MIS 5-11, and MIS 17, respectively, while paleosol S5 spans MIS 13-15 (Constantin et al., 2015).According to this age-depth model, loess accumulation started at Lunca around 700 ka, making this sequence one of the oldest loess-paleosol sections from Romania.These findings agree with the Danube loess stratigraphy (Marković et al., 2015) and the recent luminescence chronology of the northeastern Bulgarian loesspaleosol sequences (Balescu et al., 2019).
The second loess-paleosol sequence is located about 300 km east of Lunca, on the Black Sea shore near the village of Costineşti (Figure 1), and contains five loess layers (L1 to L5), the recent soil (S0) and five interbedded paleosols (S1 to S5).This section is about 12 m thick and covers the last ∼580 ka (Constantin et al., 2013), exhibiting a sedimentation rate approximately three times lower compared to Lunca.Necula et al. (2015) showed that during the last ∼580 ka the palaeosols, formed during interglacials, experienced pedogenic alteration, resulting in high amounts of superparamagnetic, SD and pseudosingle domain magnetite/maghemite grains and hematite.The loess layers, formed during glacial periods, are mainly dominated by MD and/or pseudosingle domain oxidized magnetite and some hematite, all of aeolian origin.

Methods
For this study we measured backfield IRM, FORCs, and grain size on the Lunca samples, and performed nonlinear Preisach measurements on both Lunca and Costineşti samples.The median grain size for Costineşti was determined from the measurements performed by Necula et al. (2015).Backfield IRM curves for Costineşti and their NMF unmixing are from Necula et al. (2015).The magnetic susceptibility data for Lunca section are from Constantin et al. (2015).

Lunca Backfield IRM Measurements and NMF Unmixing
The Lunca section was cleaned, prepared for sampling by avoiding bioturbated areas, and sampled at every 10 cm for magnetic susceptibility (Constantin et al., 2015) and granulometry measurements (a total of 336 samples).Based on magnetic susceptibility measurements we selected samples for remanence measurements at every ∼20 cm (a total of 168 samples).Both backfield remanent moment and direct moment curves were measured as part of the same protocol using a PMC MicroMag 3900 vibrating sample magnetometer.Samples were first magnetized at 1 T to obtain an IRM and then the field was reversed and nonlinearly incremented from 0 to 1 T in 30 steps.To align IRM unmixing with Preisach map analysis we measured backfield IRM instead of IRM acquisition curves: the NLP protocol devised by Church et al. (2016) is built on a backfield protocol even if the curves are measured in forward fields.Secondly, using backfield IRM measurements allows for reproducibility of the experiment.
For backfield IRM unmixing we used the non-parametric algorithm of Heslop and Dillon (2007), which requires IRM acquisition curves as input data.To simulate IRM acquisition, we rescaled (division by a factor of 2), reversed and inverted our backfield curves (Heslop & Dillon, 2007).The algorithm also requires all IRM curves to be measured using identical measurement fields.Therefore, we used linear interpolation to anchor all the curves to the same acquisition fields.The unmixing algorithm assumes that all IRM curves can be explained by a linear mixture of a small number of EMs.The estimation of the number of EMs included in the unmixing model was based on the calculation of the coefficient of determination between the data and the model, versus the number of EMs (Heslop & Dillon, 2007).We also performed standard parametric remanence decomposition on representative soil and loess samples from Lunca (S0 and L1) using MaxUnmix (Maxbauer et al., 2016b) to facilitate the interpretation of the EMs derived from NMF analysis (He et al., 2020).

Nonlinear Preisach Maps Measurements and Unmixing
Nonlinear Preisach maps were generated from remanence curve data sets collected using the Lake Shore Cryotronics 8604 vibrating sample magnetometer at the Biogeomagnetism Lab, National Museum of Natural History, Smithsonian Institution.We prepared 24 samples, one from each loess and paleosol layer at both sites.For each sample, 50 Preisach map curves were measured at 50 steps per curve, on a nonlinear scale from 0.1 mT to 1 T, with a 1 T saturation field and an averaging time of 1 s.Nonlinear Preisach maps were calculated following the protocol of Church et al. (2016).
Unmixing of the nonlinear Preisach maps was performed using independent component analysis (ICA), as it is implemented in the HyperSpy Python library (de la Peña et al., 2011Peña et al., , 2021)).ICA belongs to the "blind source separation" class of unsupervised analysis methods, which aim to recover EM source signals from their linear mixture without any knowledge about the mixing process or the source signals themself (Hyvärinen et al., 2001).
ICA assumes that the EMs are statistically independent, and that they have non-gaussian distributions (Hyvärinen & Oja, 2000;Hyvärinen et al., 2001).In matrix notation the mixing model can be written as where x is the observations matrix, A in the mixing matrix and s are the independent components.ICA model estimates the unmixing matrix such that the recovered sources are statistically independent.Thus a recovered independent component, y, consists of a linear combination of the observations contained in the matrix x.Mathematically can be written as where the vector w should be one of the rows of the A 1 matrix (Hyvärinen & Oja, 2000).According to the central limit theorem, a sum of two independent random variables is more Gaussian than the original variables.Thus y is least Gaussian when it is identical to one independent component s, and as such w can be chosen so as to maximize the nongaussianity of y (Hyvärinen & Oja, 2000).Details of the ICA unmixing algorithm and the selection of EMs can be found in Text S1 and Figure S1 in Supporting Information S1.As a pre-processing step, the data are "whitened" by transforming them into a new vector, in which the measurements are uncorrelated and their variances equal unity.The most popular method to perform the whitening is to apply singular value decomposition (SVD) directly on the data.Singular value decomposition plays two roles.First, SVD is used to determine the number of EMs, an operation ICA is not designed to perform (de la Peña et al., 2011;Hyvärinen & Oja, 2000;Hyvärinen et al., 2001).Second, once the number of EMs has been established, the first SVD factors, corresponding to the number of EMs, are used as initialization for the ICA.The latter step is done to avoid the problem of overfitting (de la Peña et al., 2011;Hyvärinen & Oja, 2000;Hyvärinen et al., 2001).To determine the EMs for the nonlinear Preisach data set, and their contributions to each sample, the paired SVD + ICA approach was applied on the Preisach curves.Prior to the application of SVD + ICA, the remanence curves were normalized to the saturation remanent magnetization.The entire procedure was performed directly on the curves rather than on the Preisach maps images (cf., Harrison et al., 2018).The SVD + ICA approach was also applied on the backfield IRM data set, to directly compare with the results obtained from the nonlinear Preisach maps unmixing.SVD + ICA was applied to the coercivity distributions with each coercivity spectrum being normalized to sum to one.Details about the unmixing procedure and selection of EMs can be found in Text S2 and Figure S2 in Supporting Information S1.

FORC Measurements and Unmixing
Thirteen sister samples from Lunca were used to generate FORC diagrams for the loess and paleosol layers.Firstorder reversal curves were collected using the PMC MicroMag model 3900 vibrating sample magnetometer at the Paleomagnetic Laboratory, University of Bucharest.For each sample, 150 FORCs were measured using a saturating field of 1 T, a field increment of 0.875 mT, and an averaging time of 500 ms.FORC data were processed with FORCinel 3.07 (Harrison & Feinberg, 2008), using the VARIFORC variable smoothing algorithm (Egli, 2013).FORC-PCA was conducted according to the method of Lascu et al. (2015) using the algorithms of Harrison et al. (2018).

Grain Size Measurements
Grain-size analyses of 336 samples from Lunca collected at 10 cm sampling interval were performed using a Horiba laser instrument model LA950 dispersed with sodium hexametaphosphate.Samples were additionally dispersed using the ultrasonic device of the Horiba instrument.The grain-size distribution was calculated using the Mie theory, using refractive indices of 1.33 (fluid) and 1.54 (sample) (Mie, 1908).The loess-paleosol backfield IRM data from Lunca can be explained using a two-EM mixing model, with a R 2 coefficient of determination of 0.99924 (Figure 2).A three-EM mixing model, with a R 2 coefficient of determination of 0.99936, provides only a minor improvement in terms of the model fit.The remanence gradient curves of the two EMs have median coercivities of 21 mT (EM1) and 40 mT (EM2) (Figure 2a).Typical examples of individual representative samples fits are presented in Figures 2c and 2e.For comparison, the NMF unmixing end-members and individual representative samples fits for Costineşti section are displayed in Figures 2d and 2f (Necula et al., 2015).Coefficients of determination between 0.998 and 0.999 for all the Lunca and Costineşti samples prove that the backfield IRM data are well described by the NMF model.

Results and Discussion
The variation of the relative contributions of the two EMs with depth is plotted in Figure 3a.Similar relative contributions resulted from NMF unmixing of backfield IRM curves from Costineşti section (Figure 3b) (Necula et al., 2015).The absolute contributions of EM1 and EM2, along with the granulometry and magnetic susceptibility, χ lf (Constantin et al., 2015) for Lunca section are plotted in Figure 3c.EM1 and the pedogenic grain size fraction (<5 μm) correlate positively with χ lf , whereas EM2 and the airborne fraction (>16 μm) correlate negatively with χ lf .

Nonlinear Preisach Maps
The nonlinear Preisach maps of representative loess and paleosol samples from both loess sections are displayed in Figure 4.The amplitude of the maps of all paleosols are dominated by a non-interacting SD ridge along the main diagonal (upper left to bottom right) of the diagram (Church et al., 2016).The SD ridge extends to 1 T, with the signal intensity attaining maximum between 20 and 30 mT.In addition to the SD ridge, there is a weaker amplitude contribution to the background of the Preisach map, showing spreading about the diagonal, with an asymmetrical, downward shift indicating the presence of PSD/MD particles (Church et al., 2016).
The nonlinear Preisach maps of the loess samples contain the same features, but showing a higher intensity of the background signature, relative to the SD ridge.The weak signal extension along the diagonal, with coercivities up to 1T in both loess and paleosol samples, suggests an antiferromagnetic contribution.

Unmixing of Nonlinear Preisach Maps
The variability in the nonlinear Preisach maps for the Lunca and Costineşti data sets is accounted for by two principal components (PCs), which explain more than 99.99% of the variance in both sections (Figure S3 in Supporting Information S1).The EMs and their relative contributions are depicted in Figure 5 (Lunca) and Figure 6 (Costineşti).Both sections present virtually identical EMs.EM1 (Figures 5a and 6a) is dominated by the SD diagonal ridge with low coercivities (maximum between 20 and 30 mT) (Church et al., 2016).EM2 (Figures 5b and 6b) is dominated by the PSD/MD signature (Church et al., 2016).The range of coercivities spanned by EM2 is relatively large (from a few mT to 1 T).In addition to the coarse grain signature, there is a strong signal along the diagonal, at coercivities of 50-100 mT, indicating a contribution from SD particles that are different from the ones dominating EM1.There is an additional, weaker, high coercivity signal along the diagonal, at coercivities up to 1 T, suggesting an antiferromagnetic contribution of non-interacting SD particles.
The effects of measurements noise on the calculated nonlinear Preisach maps are clearly visible in Figure 4. Calculating statistical significance levels would be necessary to distinguish between signal and noise (see e.g., Heslop & Roberts, 2012).However, the use of SVD as a first step in the unmixing procedure presents the advantage of performing denoising on the EMs, compared to individual Preisach maps (Harrison et al., 2018).This smoothing effect can be seen for EM2 (Figures 5b and 6b), which shows a negative region below the main Journal of Geophysical Research: Solid Earth positive signal, between 30 and 100 mT, which is a characteristic feature of PSD/MD magnetite (Church et al., 2016;Egli, 2021;Lascu et al., 2018) that is hard to distinguish otherwise in individual Preisach maps.
The distribution of coercivities from a nonlinear Preisach map can be calculated directly through summation over the magnetization of the remanent Preisach map (Church et al., 2016).The resulting distributions of coercivities for the two EMs show peaks at ∼22 mT for EM1 and ∼60 mT for EM2 (Figures 5c and 6c) for both sections.The EM1 peak value is similar to that of EM1 from the NMF IRM unmixing, but the distribution is restricted to coercivities <100 mT.The EM2 peak is at a value that is higher than that of EM2 from NMF IRM unmixing.The relative contributions of the two EMs are plotted in Figures 5d and 6d.Lunca (this study) and Costineşti (Necula et al., 2015) sections respectively; (c) and (e) Individual representative samples fits of loess and paleosol using the two end members model for Lunca section; (d) and (f) the same as (c) and (e) but for Costineşti section (Necula et al., 2015).

Origin of the Magnetic Minerals
The NMF unmixing model for the backfield IRM curves for the Lunca section (Figure 2a) shows the presence of two components, with median coercivities of 21 mT (EM1) and 40 mT (EM2).EM1 is interpreted as pedogenic because it covaries with the magnetic susceptibility and has an important contribution in all paleosols, and only a minor contribution in the loess layers (Figures 3a and 3c).Using the same unmixing method, the pedogenic component was also found dominant in the paleosol units from the Dobrogea Plateau at Costineşti and Mircea Vodă (see Figure 2b; Necula et al., 2013Necula et al., , 2015)), as well as in the Chaona loess section (Chinese Loess Plateau, Nie et al., 2014).Other IRM deconvolution methods showed similar median coercivities for pedogenic magnetite/  (Necula et al., 2015); (c) Absolute remanence contributions of the two endmembers, together with magnetic susceptibility, χ lf (Constantin et al., 2015), and the relative abundances of the clay (<5 μm) and airborne (>16 μm) granulometric fractions for Lunca section.maghemite worldwide (Egli, 2004a(Egli, , 2004b;;Geiss et al., 2008;Hu et al., 2013;Jordanova & Jordanova, 2020;Jordanova et al., 2011;Spassov et al., 2003).The currently accepted interpretation of this component is that it is pedogenic magnetite, with a grain size spanning the SP/SD boundary (20-25 nm) and a wide enough distribution to carry a remanent magnetization.Our IRM unmixing results corroborate previous findings that coercivity analysis of modern soils and paleosols highlights a pedogenic magnetite component which is relatively consistent across different continents and environmental conditions (e.g., Maxbauer et al., 2016a and references therein).EM2 is dominant in loess units, with only minor contributions in palaeosols, and is interpreted as a detrital component.This component has also been found in the Dobrogea Plateau, but with a higher median coercivity: ∼60 mT at Mircea Vodă (Necula et al., 2013) and ∼50 mT at Costineşti (Figure 2b; Necula et al., 2015).A comparable detrital component has also been reported in Lower Danube Basin loess, as well as in the Chinese Loess Plateau (e.g., Hu et al., 2013;Jordanova et al., 2011;Nie et al., 2014;Spassov et al., 2003).All these studies have interpreted this component as being a detrital magnetite or partially oxidized detrital magnetite (magnetite core with maghemite rim resulting from weathering).

Journal of Geophysical
As Church et al. (2016) pointed out, nonlinear Preisach maps are more robust than IRM curve unmixing because they allow both the separation of components with distinct coercivity distributions, and the identification of magnetization processes associated with different magnetic phases found in a bulk sample.The Preisach EM1 coercivity distribution is identical both at Lunca and Costineşti, with a median coercivity of ∼22 mT (Figures 5c  and 6c).However the EM1 component from backfield IRM NMF unmixing shows a high-coercivity tail which is not present in the corresponding Preisach EM1.The Preisach EM2 distribution has a median coercivity of ∼60 mT (Figures 5c and 6c), a higher value than that of EM2 from NMF IRM unmixing, especially for Lunca.These discrepancies may come from using different experimental protocols, number of data points (30 for IRM acquisition vs. 2,500 for Preisach), conditioning field steps, measurement noise, unmixing algorithms, etc (e.g., Heslop, 2015;Heslop & Dillon, 2007).To resolve these discrepancies, we applied the SVD + ICA unmixing method to the IRM backfield curve data sets from both Lunca and Costineşti.The analysis yielded two EMs in both sections, which explain about 99.96% of the variance at Costineşti and 99.98% of the variance at Lunca (Figure S4 in Supporting Information S1).The two EMs and their relative contributions are illustrated in Figure 7.At both locations the median coercivity for EM1 is 21 mT and for EM2 is 60 mT, matching closely the median coercivities of the Preisach EMs (Figures 5c and 6c).This analysis provides strong evidence that both sets of EMs describe the same magnetic components.The backfield IRM EM1 is controlled by the same population of grains that controls the Preisach EM1 component.Based on the Preisach maps of EM1 (Figures 5a and 6a), we can determine unambiguously that this population of grains is dominated by low-coercivity non-interacting SD magnetite.Thus, using Preisach unmixing, we can clearly state that EM1 contains a single magnetic population, exclusively of pedogenic origin.On the other hand, Preisach unmixing demonstrates that EM2 does not represent a single magnetic phase, but is comprised of a mixture of PSD/MD and SD magnetite, as well as non-interacting SD hematite, all with the same origin, likely eolian.These multiple components cannot be further extracted because their coercivity spectra overlap substantially (e.g., Egli, 2021;He et al., 2020).
The relative abundances of EM1 and EM2 differ according to the unmixing method used.Thus, NMF IRM unmixing provides higher EM1 abundances (Figures 3a and 3b) than SVD + ICA IRM unmixing (Figures 7b and  7d), and Preisach maps unmixing (Figures 5d and 6d).In case of nonparametric unmixing methods the shape of EMs and the resulting relative abundance are closely related (He et al., 2020;Heslop, 2015).Thus, the main cause of discrepancy between the results of NMF and SVD + ICA backfield IRM unmixing is related to the specific protocols each method uses: different shapes of the EM coercivity distributions provided by each unmixing protocol will give rise simultaneously to different abundances.For example, NMF considers only non-negative, sum-to-1 solutions for the mixing matrix, while SVD + ICA does not have such constraints.In case of a two EM model, considering only sum-to-1 solutions for the mixing matrix forces the relative contributions of the EMs to take 0 and 1 values (e.g., Figure 3a).This implies that the EMs themselves are contained in the data set, assumption that is not always realistic.More commonly, all samples are mixtures between the two EMs in different proportions, with none of the samples representing a pure EM.On the other hand, SVD + ICA takes as EMs the most non-Gaussian solutions without imposing sum to 1 condition.Therefore, it is worth remind readers to carefully take into consideration the sensitivity of the unmixing methods in determining the EMs reflecting some inherent uncertainty associated to each specific analysis.
On the other hand the coercivity distributions from Preisach maps unmixing (Figures 5c and 6c) and SVD + ICA backfield IRM unmixing match closely (Figures 7a and 7c).Furthermore, the EM proportions from both methods show a strong linear dependence (R 2 = 0.99 for Costineşti and R 2 = 0.96 for Lunca) and slopes close to 1 for both sections (Figure 8).This reinforces the fact that the backfield IRM coercivity distributions and Preisach maps capture the same information.Small differences between the results from IRM curves and Preisach maps could be due to the fact that the measurements were done on sister samples, that the backfield IRMs were measured for fields from 1 mT to 1 T, while Preisach maps contain a wider range of fields (from 0.1 mT to 1 T), and that there was a significant difference between the number of backfield IRM samples and Preisach maps samples (168 and 131 vs. 13 and 11 respectively).
However, results from both unmixing methods provide the same message: loess layers are dominated by detrital magnetic particles of eolian origin, while paleosols are dominated by the pedogenic fraction consisting mainly in non-interacting SD magnetite/maghemite.Some exceptions are paleosol S1 from Lunca and loess units L4 and L5 from Costineşti, which exhibit a more even proportion of EMs.The Lunca S1 paleosol has probably not been developed enough to lead to a predominance of the pedogenic fraction over the detrital one.This is supported by the bulk magnetic susceptibility and pedogenic fraction (<5 μm) variations (Figure 3c) which contain the lowest values in this unit at Lunca.The Costineşti L4 and L5 loess layers, on the other hand, have a significant EM1 fraction suggesting that these units are affected by pedogenesis.This is also supported by all the concentrationdependent magnetic parameters from Costineşti, which, in these units, show values as high as in the S0, S1 and S2 paleosols (Necula et al., 2015).Therefore, Preisach maps can be used as a straightforward tool for explaining the IRM coercivity distributions in terms of domain state for the main ferromagnetic minerals.

Single Sample Parametric Backfield IRM Modeling
The unmixing results for parametric remanence decomposition using MaxUnmix (Maxbauer et al., 2016b) are depicted in Figure 9.For S0, the remanence gradient has a single peak and it is almost symmetric and very narrow suggesting overlap between coercivity components (Figure 9).We tested the parametric unmixing in two steps.Firstly, we checked if decomposition depends on the type of model distributions used to fit the data (Egli, 2003;Zhao et al., 2018).MaxUnmix application includes two of these namely SN and simple Gaussian (G) distributions.At least three components are needed to fit the coercivity distribution for both types of curves (Figures 9a  and 9b), keeping the same initial parameters for both.For the first two components the median coercivity remains the same (4 and 22 mT using SN, and 5 and 21 mT using G respectively).The third component, however, has different coercivities in the two cases (88 mT for SN and 70 mT for G).Moreover, the amplitude of the components varies depending on the type of curve used, the most noticeable being the second component (Figures 8a  and 8b).Secondly, we tested the influence of initial conditions on the unmixing results.Therefore, we changed the coercivity of the third component to ∼100 mT and only slightly for the first two components.The new fitting solution seems to be as good as the previous one showing that a third component with a median coercivity of 102 mT is also valid (Figures 9c and 9d).However, the first component's median coercivity is modified substantially when using G, showing the dependence of curve type.Overall, only the second component (21-24 mT median coercivity) is more or less stable, regardless of initial conditions and types of functions used for fitting the data.
The interpretation of the three components can be ambiguous.The first component may be interpreted as having several origins: it could represent MD or SP particles signatures (Egli, 2004a(Egli, , 2004b)), could be the result of magnetostatic interactions between magnetic particles (Egli, 2004a(Egli, , 2004b;;Heslop et al., 2002Heslop et al., , 2004)), could be simply a filler because the model distributions cannot mimic the real data (Egli, 2003;Zhao et al., 2018), or any combination of the above.The third component can be interpreted as partially oxidized magnetite if coercivity is in the 70-88 mT range, but can also be interpreted as maghemite and/or hematite if its coercivity is ∼100 mT (e.g., Spassov et al., 2003).For the second component (21-24 mT) the interpretation is more or less straightforward: low coercivity pedogenic magnetite most probably in SD state (e.g., Maxbauer et al., 2016a).
For L1 we also obtained three components without constraining the coercivity values.For the SN model the median coercivities of the components are 7, 24, and 78 mT (Figure 9e), while for the G model they are 8, 25, and 78 mT (Figure 9f), matching the components obtained for the unconstrained S0 models.However, the third component is consistent regardless of the type of theoretical curves used to fit the data, and can be interpreted as partially oxidized magnetite.The amplitude of this component is larger compared to the paleosol sample, where the second component has the largest amplitude, indicating a dominance of particles of eolian origin.
Given the ambiguity of the first component, there are probably only two major components that contribute to the loess and paleosol representative samples from Lunca: a pedogenic component and one of eolian origin.The pedogenic component is similar to that provided by the NMF and SVD + ICA methods, whereas the eolian one has notably higher coercivity.If additional information about the sample (e.g., domain state, magnetic mineralogy, magnetostatic interactions) is lacking it is challenging to interpret these solutions.Therefore, using Preisach maps we can straightforwardly interpret the coercivity components given by any of these remanence curves unmixing methods (parametric and/or non parametric).

FORC-PCA Unmixing
FORC-PCA yielded two EMs along the first PC (eigenvector 0), which explains nearly 98% of the variability in the data set (Figures 10a and 10b).EM1 (Figure 10c) has a signature dominated by a low coercivity central ridge (maximum amplitude at about 10 mT), suggesting the presence of non-interacting SD particles.The background EM1 signature could indicate strong magnetostatic interactions or the presence of a PSD/MD component.EM2 is dominated by a PSD/MD signal (Figure 10d), with a weak central ridge signature.The FORC diagrams of the EMs suggest they are likely mixtures of pedogenic and eolian components, indicating that the FORC-PCA separation is not optimal.Moreover, there is no high coercivity component captured.For this type of highresolution measurement the results are limited to a narrow range of coercivities and interaction fields due to  also not smoothed, thus avoiding artifacts associated to this post processing operation.FORC unmixing results cannot be directly compared to the remanence-based unmixing results because, FORCs capture a combination of induced, transient, remanent, and thermal relaxation magnetizations.An effect of this is that FORC EM1 maximum coercivity has a lower value (10 mT) than EM1 maximum coercivity from Preisach maps unmixing (∼21 mT).Moreover the Preisach unmixing provides better separation of the true pedogenic and eolian components.Preisach EM1 is easily identifiable as noninteracting SD particles (dominance of ridge signal) of pedogenic origin, while Preisach EM2 contains a mix of particles of eolian origin (PSD/MD + higher coercivity SD ferrimagnetic, as well as antiferromagnetic contributions).Finally the entire FORC-PCA procedure (collecting data, processing FORCs, and setting up the FORC-PCA) is at least as time consuming as the Preisach unmixing.All these indicate Preisach maps as the most suitable to help the interpretation of IRM unmixing in term of domain states and interactions, at least for loess-paleosol deposits.

Pedogenesis and Provenance in Lower Danube Loess-Paleosol Deposits
The two loess sections were developed under different climatic conditions.At Lunca, pedogenesis intensity decreases gradually from S5 to S1, suggesting progressively drier interglacial conditions with time (Figure 3c).Obreht et al. (2016) found a similar aridification pattern in the Stalač loess section in central Serbia, documenting a transition from Mediterranean climate to more continental conditions in the Balkans from MIS 9 to present.The same aridification trend was described in the Zeemun loess sequence in the Middle Danube Basin starting with MIS 11 (Laag et al., 2021).A similar trend of increased continentality over time is also observed around the Azov Sea, suggesting a supra-regional climatic trend (Liang et al., 2016;Obreht et al., 2016).The Lunca section documents the highest intensity of pedogenesis the S5 paleosol suggesting that this aridification trend could have started in the area as early as MIS13-15, 500-600 ka ago (Figure 3c).
Preisach maps unmixing revealed identical EM1 in both sections irrespective of the degree of pedogenesis and the evolution of climate during the last 700 ka.This implies that the pedogenesis produce the same type of magnetic mineral, namely non interacting SD magnetite/maghemite, with the same coercivity distributions.The variations in EM1 abundances at the two sites must then reflect the intensity of pedogenesis, which is in turn controlled by climate, with more magnetite being produced during wetter conditions.
The source areas for the loess deposits in the Middle and Lower Danube Basin, as well as for the Black Sea shore, are predominantly local (Lehmkuhl et al., 2021).In the Lower Danube Basin, the main source of dust is thought to be the very alluvium of the Danube River (Fenn et al., 2022;Lehmkuhl et al., 2021).In the Dobrogea Plateau and Black Sea shore, the main source of dust is the continental shelf of the Black Sea, which is exposed during glacials as a result of sea level fall (Hoyle et al., 2021;Jipa, 2014).Loess granulometry shows a systematic finer grain size at Costineşti compared to Lunca (e.g., more than 4 times finer in unit L5), supporting the model of different source areas for the two locations (Figure 11).
Finer dust could also imply transport from greater distances.Indeed, it has been documented that additional sources of dust for this area are Ukrainian glaciofluvial sediments (Buggle et al., 2008), and the Aral and Caspian arid lands (Stephens et al., 2003).Preisach EM2, which is of eolian origin, can be used to characterize the magnetic population of the source area, as well as transport processes.EM2 shows identical magnetic mineral content both at Costineşti and Lunca, suggesting that the same magnetic components are transported and deposited in the region irrespective of the dust provenance and transport processes.Moreover SD and/or PSD magnetite is probably oxidized, which gives its higher coercivity compared to paleosols (e.g., Ge et al., 2014;Özdemir & Dunlop, 2010;van Velzen & Dekkers, 1999).
Journal of Geophysical Research: Solid Earth 10.1029/2024JB029004

Conclusions
Our rock magnetic investigation of two loess-paleosol sections from Southeastern Europe showed that the combination of non-parametric backfield IRM unmixing and Preisach maps unmixing can be successfully used for identifying magnetic components and processes in loess-paleosol deposits, providing a considerable advantage in their interpretation.Our findings can be summarized as follows.
• The nonparametric unmixing method can provide more robust insight into the difference in magnetic properties of constituent magnetic components than parametric unmixing method.However, uncertainty still remains as the number of EMs, their shapes and thus abundances may still depend on the unmixing method used.
Preisach maps provide useful information on the internal structure of the EMs that is critical for interpreting the mineral phases and domain states of the components.• Our results demonstrate that EM1 is dominated by pedogenic, non-interacting SD magnetite/maghemite with low coercivity.EM2 is a mixture of magnetic phases comprising PSD/MD and SD magnetite/maghemite with higher coercivity than EM1, as well as non-interacting SD hematite, all of eolian origin.• Pedogenesis produces the same type of magnetic mineral, namely non-interacting SD magnetite/maghemite, with the same coercivity distribution regardless of location.Since this coercivity component has been consistently found worldwide in loessic soils, nonlinear Preisach measurements can provide a way to unambiguously diagnose its presence.• The eolian EM contains both coarse (PSD/MD) and fine (SD of higher coercivity than EM1 most probably due to low temperature oxidation) ferrimagnetic particles, as well as an antiferromagnetic component.Thus, nonlinear Preisach maps allow the identification of the domain states of the mineral phases making up this complex composition.The same EM2 in the two loess sequences indicates that the same magnetic dust fraction is transported and deposited irrespective of provenance and transport processes.• There are drawbacks to using non-linear Preisach maps.Because they are based on remanence curves, Preisach maps will not capture the superparamagnetic component of pedogenic enhancement, which may be significant.The measurements are much more time consuming than IRM curves but not than FORCs and, like FORCs, they are generally impractical for a large number of samples.Non-linear Preisach measurements are thus ideal for characterization of representative samples, at a resolution that captures all the magnetic features of interest.

Figure 1 .
Figure 1.Location of the loess-paleosol sections investigated in this study: Lunca and Costineşti.

Figure 9 .
Figure 9. Parametric unmixing of backfield isothermal remanent magnetization curves for representative samples of paleosol (S0) and loess (L1) from Lunca section: (a) S0 modeling using skewed normal curves; (b) S0 modeling using Gaussian curves; (c) and (d) the same as (a) and (b) with modified initial parameters; (e) and (f) the same as (a) and (b) for L1.Gray dots represent the data and yellow line represents the modeled coercivity distribution.Shaded areas represent 95% confidence intervals associated with each component.