My research interests are quite varied, but my publication list should give you a flavour of what I do. Most of my work has been in rock magnetism, but I have also spent time working on problems in climate modelling, time series analysis and statistics. One of the problems I have become interested in is understanding natural mixtures and finding mathematical methods which help us to break up mixed signals into their component parts in order that we can understand them better. I'm currently working as part of the ARC Discovery Project "Australian dust: its response to, and role in, climate change".

Research Project Links

Unmixing magnetic remanence curves without a priori knowledge

Increase in African dust flux at the onset of commercial agriculture in the Sahel region

Characterizing the uncertainty of taxa relative abundances derived from count data

Can oceanic paleothermometers reconstruct the Atlantic Multidecadal Oscillation?

Temperature fingerprinting of the Atlantic meridional overturning circulation response to freshwater forcing

Using non-negative matrix factorization in the unmixing of diffuse reflectance spectra

Free Ebooks

An Introduction to MATLAB for Geoscientists (free Ebook)

Practical Statistics for Geoscientists (free Ebook)

The cooling rate gradient that exists between the glassy rim of a mid-ocean ridge basalt (MORB) pillow and its interior produces a spatial trend in titanomagnetite size and composition, which can be tracked by a number of rock magnetic procedures. A total of 68 remenanece curves were unmixed into a 5 end-member model to study the evolution of the magnetic mineral assemblage through the T787-R1 MORB pillow.

Coercivity spectra and acquisition curves for the five end-members extracted from the MORB data set. The end-members correspond to the different grain-size and compositional components in the basalt. Transitional end-members are required to explain intermediate assemblages, this means that behaviour evolves through the MORB in a non-linear fashion and cannot be modelled as a simple linear combination of SD and MD end-members.

The abundance data follows the expected pattern, with the SD grains that dominate the behaviour at the rapidly cooled margin, gradually being replaced by MD grains towards the interior, via a number of transitional assemblages.

The presented case study demonstrates the ability of the NMF approach to unmix the variations in the magnetic assemblage of the T787-R1 MORB specimen. The description of the spatial variation of the magnetic properties provided by the end-member model corresponds clearly to the intra-pillow processes described by previous authors. In cases such as the T787-R1 where the magnetic mineral inventory is seen to evolve gradually, it is necessary to derive a relatively large number of end-members in order that sufficient transitional assemblages can be defined to produce a good description of the measured data.

Software: IRM-Unmixer is a stand-alone application which will help you to produce a linear mixture model for your own IRM acquisition data. Please download the zip-file, which also contains a manual discussing the ideas behind IRM-Unmixer and a detailed explaination of how to use the software. If you don't have MATLAB on your computer (or you are using a different version of MATLAB to the one I compiled the code with) then you will need to install the MATLAB Runtime package.

Increase in African dust flux at the onset of commercial agriculture in the Sahel region

with Stefan Mulitza

In a proxy based study we reconstructed a record of dust deposition off northwest Africa for the late Holocene. Drought in the Sahel region has been shown to be related to the AMO and thus potentially changes in the strength of the AMOC (Zhang & Delworth, 2006). Via this connection it is feasible that the flux of dust deposited offshore Northwest Africa could contain a record of past AMOC variability. We examined the chemistry and grain-size distribution of terrigenous sediments deposited at a marine site located directly under the West African dust plume. On the basis of the XRF derived chemical composition of the sediments, an end-member unmixing model allowed the proportions of aeolian, river-borne and marine derived material in each sediment sample to be estimated. When combined with sedimentation rates calculated from a high resolution age model, the proportion of aeolian material could be converted into a record of absolute dust flux spanning the last ~3200 years.

A coherent multidecadal variability possibly attributable to variations in the AMOC could not be extracted from the record of dust flux. However it could be shown that dust deposition was related to precipitation in tropical West Africa until the seventeenth century. The most dramatic feature of the record appears at the beginning of the nineteenth century, a sharp increase in dust deposition that appears to parallel the advent of commercial agriculture in the Sahel region. Our findings suggest that human-induced dust emissions from the Sahel region have made a significant contributed to the atmospheric dust load for about 200 years.

Characterizing the uncertainty of taxa relative abundances derived from count data

with Stijn De Schepper and Ulrike Proske

When working with micropalaeontological assemblage data it is necessary to accept that the true population can never be quantified fully. Thus, the investigator must count an assemblage of individuals in the hope of obtaining a statistically representative sample of the population. Once the counting procedure is complete, it is essential to determine the uncertainty associated with the estimated taxa relative abundances. Such uncertainties can be obtained by assuming that the relative abundances originate from a multinomial distribution. With this assumption a confidence region can be easily defined mathematically using existing statistical methods. For assemblages containing more than three taxa, however, the numerical determination and representation of multinomial confidence regions in high-dimensional space becomes an exceedingly cumbersome task. We outlined a simple method that allows diagnostic values of the multinomial confidence region to be calculated for any number of taxa without the need to calculate the confidence region in its entirety. Examples of such values are the most extreme relative abundances, which allow the user to quantify the interval of a given taxon within the confidence region and the maximum/minimum values of various diversity indices. These values can also provide the investigator with important information that can guide decisions concerning the number of individuals that should be counted.

(left) Example of the 95% multinomial condence region (shaded) calculated for the relative abundances A=0.1, B=0.3 and C=0.6 (black point) when N=100, using the approach of Weltje (2002). (right) Determination of extreme values (dashed lines). The span between the lower and upper extremes of each component is labeled using the delta symbol. In addition, the point at which each extreme value is reached (circles) can be determined.

Software: Rcounts is a package which will help you to perform the calculations described above. The package requires the free software R to be installed (versions for Windows, MacOS and Unix can be downloaded from www.r-project.org/). Please download the Rcounts zip-file, which also contains a manual discussing the ideas behind Rcounts and a detailed explaination of how to use the software.

Can oceanic paleothermometers reconstruct the Atlantic Multidecadal Oscillation?

with Andre Paul

The Atlantic Multidecadal Oscillation (AMO) is a leading mode of sea surface temperature (SST) variability with a period of 30-80 years thought to be driven by internal fluctuations of the AMOC. The SST pattern attributed to the AMO plays an important role in modulating key components of the climate system, for example, the African and Indian summer monsoons, North American and European summer climate, and the activity of hurricanes in the Atlantic (Zhang & Delworth, 2006). Therefore understanding the evolution of the AMO and its variability through time is of key importance. A spatial pattern of SST anomalies associated with the AMO, was obtained by calculating the difference in detrended Atlantic SST between 1941-1960 (observed positive AMO phase) and 1971-1990 (observed negative AMO phase) based on instrumental data. The map of AMO related SST anomalies was then converted in anomalous heat fluxes with zero spatial integral. The resulting anomalous heat fluxes form a basin scale dipole about the equator and thus are suitable to impose northward heat transport in the Atlantic, mimicking changes in the AMOC.

A 2000-year integration was employed to bring the UVic ESCM into equilibrium with boundary conditions and insolation forcing appropriate for 1850 CE. The anomalous heat flux was then introduced with a magnitude controlled by a repeating pattern of the AMO index following the time series obtained from the data of Rayner et al. (2003). After the run was completed, the modeled AMO index was calculated based on the SST anomalies with respect to a control run (Figure 1).

The constructed AMO fingerprint was based on the amplitude of temperature anomalies generated in the surface and deeper layers of the ocean as a result of the anomalous heat flux forcing. With this information the question if proxy studies, based on properties such as the composition of foraminiferal tests or the growth-rate / composition of corals, would be able to reconstruct the kind of temperature oscillations expected as part of natural AMO variability? To determine the validity of the fingerprint for any given temperature proxy, factors such as the standard error of estimation and signal filtering had to be taken into consideration.

Observed and modelled AMO index, pre- (lines) and post- (shaded) low-pass filtering. Positive and negative phases of the AMO are shown by red and blue shading, respectively.

The model data and subsequent fingerprinting analysis revealed that, under a best case scenario (annual resolution), existing ocean paleothermometers may be capable of recovering interpretable records of the temperature variations in the North Atlantic associated with the past activity of the AMO. Whilst such a scenario appears to be appropriate for coral-based reconstructions, sediment records have a lower resolution and may have been disrupted by processes such as bioturbation. For the case where sediment records have a decadal resolution, the quoted accuracy of proxies such as Mg/Ca, TEX86 and U37K' suggests that they would yield reconstructions of the temperature anomalies resulting from the AMO with very low signal-to-noise ratios.

Temperature fingerprinting of the Atlantic meridional overturning circulation response to freshwater forcing

with Andre Paul

Hypothesized changes in the strength of the AMOC play a key role in explaining a number of large-scale events recorded in the paleoclimatic record. In particular, cold episodes such as the Younger Dryas and 8.2 ka event have been attributed to reductions in northward heat transport resulting from freshwater induced slowing of the AMOC. In order to understand the underlying processes behind such events it is desirable to design proxy investigations that target regions where reconstructed temperatures potentially carry high fidelity information concerning the strength of the past AMOC. To investigate the relationship between AMOC strength and temperature as freshwater is added to the North Atlantic, a hosing experiment of the type described by Rahmstorf et al. (2005) was performed. A 2000 year integration was employed to bring the model into equilibrium with boundary conditions and insolation forcing appropriate for 1850 CE. Subsequently, an increasing virtual salt flux was applied to the North Atlantic across a uniform band spanning 20-50°N to simulate an uncompensated freshwater input. The resulting reduction in the overturning with increasing freshwater forcing defines the “upper branch” of the major AMOC hysteresis loop. To assess the hysteretic response of the system, a minor hysteresis branch was also determined in a separate integration. This minor branch took the form of a so-called first-order reversal curve (FORC, Mayergoyz (2003)).

The fingerprint of the Atlantic temperature response to a slowdown and recovery of the AMOC was constructed by developing metrics with which to quantify four key factors:

1)    The sensitivity of the temperature response.
2)    The nonlinearity of the temperature response.
3)    The monotonicity of the temperature response.
4)    The level of symmetry of the response to AMOC slowdown verses recovery.

Each of these factors are complementary and necessary in order to identify areas which may appear to be suitable locations for proxy reconstructions based on their overall sensitivity, but suffer from problems such as strong nonlinearity which will limit the ability to relate reconstructed temperature changes to the strength of the AMOC.

The overall structure of the developed fingerprint revealed a number of important insights into the suitability of the temperature record at a given location to provide information concerning the past strength of the AMOC. In particular, studies which reconstruct both sea-surface and seafloor temperatures must consider the existence of (potentially unquantifiable) biases between their records and substantial inconsistencies in how they represent overall strength of the AMOC. Furthermore, the inclusion of a FORC into the experiments demonstrated that large areas of the Atlantic exhibit an asymmetrical temperature response to a declining and then recovering AMOC. This implies that cold events attributed to a decreased AMOC must be considered in detail with an awareness that similar temperatures during the cooling and subsequent warming sections of the event cannot be assumed to correspond to the same AMOC strength.

Software: MATLAB functions are available to perform the analysis, please email me for details.

Using non-negative matrix factorization in the unmixing of diffuse reflectance spectra

Diffuse reflectance spectroscopy (DRS) determines the “colour” of sediment across a given range of wavelengths. A Minolta 2002C instrument is used routinely on research cruises to measure percentage reflectance at 10 nm increments across the visible wavelength band (400 and 700 nm). The visible light wavelength range can provide information on a number of process indicative materials such as iron oxides, clay minerals, carbonates and organic matter. Characteristic spectral peaks in their first derivative reflectance curves help to identify these materials.

Hematite and goethite both have characteristic peaks in the first derivatives of the reflectance spectra. This helps in their identification in sediment samples.

The DRS spectra of natural sediments do however represent a composite signal of the different constituent minerals; therefore it is necessary to “unmix” the data to obtain information on specific materials. Certain conditions must be met if we are to unmix a collection of DRS spectra in a meaningful way and we attempt to characterise two factors:

End-members: represent the reflectance spectra of the constituent materials. The end-member spectra must have the same units as the measured DRS data [% reflectance in the range 0-100]

Fractional abundances: give the proportion of each end-member present in the measured spectra. The abundances must 0 or greater (in other words negative abundances are not allowed) and for each measured spectrum the abundances of the end-members must add to 1 (so-called full additivity).

To obtain such an unmixed representation of a DRS data set we take the matrix of measured data and try to reduce it into the discussed factors (fractional abundances and end-members) that meet the mixing constraints of non-negativity and full additivity. To produce factors that only contain positive values may seem trivial, but in fact it places some tricky non-linear constraints on the unmixing. An analytical solution to the problem is not possible, so instead we must proceed numerically. Starting with a first guess of the solution (in practise just random numbers) the non-negative matrix factorisation algorithm of Lee and Seung (1999) is applied to gradually push the model towards the correct solution.

Unmixed DRS data from ODP Site 967 in the eastern Mediterranean. A four end-member solution provides a good model of the data and an be interpreted in terms of environmental change.

Shown above is an example taken from Mediterranean sediments (ODP Site 967) spanning the period between 900 and 1500 thousand years ago. Four end-members are identified which correspond to organic matter which increases in more humid climates, aeolian dust from Africa transported by the strong winds that form during drier periods, grey ooze that consists of microfossil shells and fluvial material which is predominantly composed of clays transported by the Nile. The unmixing of the DRS data set therefore gives us detailed information concerning climate change and variations in the importance of the different transport mechanisms that fed sediment into the eastern Mediterranean.

Software: DRS-Unmixer is a stand-alone application which will help you to produce a linear mixture model for your own DRS data sets. Please download the zip-file, which also contains a manual discussing the ideas behind DRS-Unmixer and a detailed explaination of how to use the software. If you don't have MATLAB on your computer (or you are using a different version of MATLAB to the one I compiled the code with) then you will need to install the MATLAB Runtime package.

An Introduction to MATLAB for Geoscientists (free Ebook)

I wrote this Ebook a few years ago for a MATLAB course I was teaching. It provides an introduction to the basics of MATLAB and hopefully presents the topics in a digestable way. The pdf can be downloaded from the link below and there is an accompanying zip file which contains the functions and data sets used in the examples.

An Introduction to MATLAB for Geoscientists Ebook (pdf)
An Introduction to MATLAB for Geoscientists, functions and data sets (zip file)

Practical Statistics for Geoscientists (free Ebook)

Statistics are the most powerful tool we have for separating scientific fact from fiction. The aim of this text is to provide an introduction to statistical techniques that are useful in the analysis and characterization of geological data. A focus is placed on conceptual understanding of how specific methods work and the situations in which they can and cannot be applied. A number of practical examples are discussed, providing the opportunity for hands-on learning through the processing of real data sets with statistical software.

An Introduction to MATLAB for Geoscientists Ebook (pdf)
An Introduction to MATLAB for Geoscientists, functions and data sets (zip file)

Random Musings

The following is a small selection of nano-projects that may be of interest. Please note that these ideas are not fully developed and provided as points of interest that might provide you with ideas rather than complete research projects.

A seven line unsupervised unmixer

The unmixer is based on the observation that the vertices of a simplex always lie on the surface on the minimal volume (hyper)ellipsoid which bounds the simplex. Further, the considered simplex is the maximum volume simplex which can be placed inside the bounding (hyper)ellipsoid. Interestingly it also appears that the ratio of the volume of the simplex and the volume of bounding (hyper)ellipsoid is constant for any given number of dimensions.

We can use the above observation to estimate the vertices of a mixing space by finding the minimum volume ellipsoid that bounds the mixture data (this is much simpler than trying to directly find the minimum volume simplex which bounds the data). The important caveat of this approach is that it assumes that a single pure example of each of the end-members is present within your data set and that the measurements are effectively noise free, but hey, how much can you expect from 7 lines of code!

The MATLAB code is setup as follows and assumes that you have the free package YALMIP installed to do the hard work. The routine solves the mixing equation X=M*S for d end-members where X is a matrix of measured spectra (1 row per variable, 1 column per sample), M contains the end-members (1 row per variable, 1 column per end-member) and S contains the mixing proportions (1 row per end-member, 1 column per sample).

Here's what the code actually does (assuming your matrix of data, X, and value of d are defined in advance):

Line 1: performs PCA and represents the data with the first d-1 components.
Line 2: sets up the variables for YALMIP.
Line 3: defines constraints for YALMIP (specifically, each point in X must lie inside the final ellipsoid).
Line 4: find the minimum volume ellipsoid surrounding the data.
Line 5: locate the d points which lie on the surface of the ellipsoid.
Line 6: return the d points to the measurement space (these are the end-members, Mhat).
Line 7: solve the least-squares problem to find the mixing proportions Shat).

So finally, here's the code...

[U,Q,V]=svd(bsxfun(@minus,X,mean(X,2))',0); x=U(:,1:d-1)';
A = sdpvar(d-1,d-1); b = sdpvar(d-1,1); F = set([]);
for i=1:size(x,2), F = F + set(norm(A*x(:,i)+b) <= 1.0); end
solvesdp(F,-logdet(A),sdpsettings('verbose',0));
[Fidx,Fidx]=sort(sqrt(sum((double(A)*x+double(b)*ones(1,size(x,2))).^2)));
Mhat=V(:,1:d-1)*Q(1:d-1,1:d-1)*x(:,Fidx(end-d+1:end))+mean(X,2)*ones(1,d);
Shat=[Mhat;ones(1,size(Mhat,2))]\[X;ones(1,size(X,2))];