Evaluation of Biasing Factors - Towards a Higher Comparability of Geothermometric Data obtained

Towards a Higher Comparability of Geothermometric Data obtained by Raman spectroscopy of Carbonaceous Material

Part 1: Evaluation of Biasing Factors

Authored by N.K. Lünsdorf, I. Dunkl, B. C. Schmidt, G. Rantitsch, H. v. Eynatten and published as: Towards a Higher Comparability of Geothermometric Data obtained by Raman spectroscopy of Carbonaceous Material. Part 1: Evaluation of Biasing Factors. Geostandards and

Geoanalytical Research, 2014, 38(1), p. 73-94.

2.1. Introduction

Geothermometry by Raman spectroscopy of carbonaceous material (RSCM; Beyssac et al., 2002a) becomes more and more popular due to its ease of applicability and non-destructiveness (e.g. Rantitsch et al., 2004; Guedes et al., 2005; Forer et al., 2009; Huang et al., 2010; Wiederkehr et al., 2011; Endo et al., 2012). The Raman spectrum of carbonaceous material (CM) enclosed in metasediments changes systematically with increasing degree of metamorphism (e.g. Pasteris and Wopenka, 1991; Jehlička and Bény, 1992; Wopenka and Pasteris, 1993). Beyssac et al. (2002a) showed that those changes are mainly controlled by temperature and calibrated a geothermome-ter based on CM Raman spectra. This geothermomegeothermome-ter is based on paramegeothermome-ters calculated by spectral curve-fitting of the Raman bands (Fig. 2.1.1). However, there are multitude of different curve-fitting strategies using variable numbers (2 – 5) of model components (Lorentzian-, Voigt-, Gaussian- or Breit-Wigner-Fano-functions)Voigt-, fitting the acquired spectrum (e.g. Beyssac et al.Voigt-, 2002a; Quirico et al., 2003; Sadezky et al., 2005; Lahfid et al., 2010). Important function pa-rameters (position, height, full width at half maximum FWHM, and area), reflecting the thermal transformation of CM, are derived from the extracted components. This allows the calibration of different parameter-ratios against metamorphic temperature (Beyssac et al., 2002a; Rantitsch et al., 2004; Rahl et al., 2005; Baziotis et al., 2006; Aoya et al., 2010; Lahfid et al., 2010).

As the parameter-ratios are the condensates of many steps (e.g. sample preparation, Raman measurement, and spectrum evaluation) they accumulate several biasing factors. The sources of bias can be grouped into three categories: (1) bias intrinsic to spectral curve-fitting, (2) bias intrinsic to the CM and (3) bias intrinsic to the experimental design and the specific Raman system used.

Factors of the first category are different baseline corrections, different mathematical functions (Gaussian, Lorentzian, Voigt, etc.) used for peak-fitting, and the different number of components used to model the Raman spectra. Examples for the second category are structural anisotropy, sample preparation, and sample-heterogeneity. Factors of the third category include used excitation wavelength, spectral grating and light detection device, among others.

In this paper, many of the above mentioned potential sources of bias are evaluated by a suite of simple experiments in which, at best, only one parameter influences the experimental results.

Figure 2.1.1. A) Representative Raman spectra of the three ’crystallinity levels’. B) An example for the decomposition of a ’crystallinity level 1’ Raman spectrum by five components according to Lahfid et al.

(2010). The two bands of the first-order region are described by D1, D2, D3, D4 and G components. D4 widens the low wavenumber side of the D1 band while D2 appears only as a weak shoulder on the high wavenumber side of the G band.

2.2. Methods and samples

2.2.1. Samples. The sample set (Table 2.2.1) covers a wide structural range from low to high

’crystallinity’ degree of CM. The studied samples derive from the Triassic flysch of the Tethyan Himalayan sequence of SE Tibet (Dunkl et al., 2011), from the Eastern Alps (Rantitsch et al., 2004), and from the Thuringian Forest (Germany) (Kunert, 1999).

In order to enrich CM every sample has been treated chemically. The samples were crushed to particles smaller than 5 to 10 millimeters. The rock-chips were initially placed in a 1:1 solution of 37 % hydrochloric acid to dissolve carbonates and after decantation mixed with 1:1 diluted 48 % hydrofluoric acid to dissolve silicates. After hydrofluoric acid treatment the sample suspensions were decanted and diluted with de-ionized water until a pH-value of 5 to 6 was reached. Remaining fluids were evaporated in a drying oven at 50 °C. About 10 to 20 mg of the dried CM were mixed with 1 - 2 ml de-ionized water in a small glass vial and placed into an ultrasonic bath for about 60 seconds in order to disperse the carbonaceous material. This suspension was deposited on a glass slide.

2.2.2. Raman spectroscopy. All Raman measurements were performed with a Horiba Jobin Yvon HR800-UV spectrometer, with attached Olympus BX41 microscope, if not stated otherwise.

The general measurement configuration used a 488 nm Ar⁺- laser for excitation, a spectral grating with 600 l/mm, a long working distance 100x objective with a numerical aperture of 0.8 and the diameter of the confocal hole was set to 100 µm. If not stated otherwise the laser light was circular polarized, a spectral range of 700 – 2000 cm^-1was recorded in one spectral window in 3-5 accumu-lations of 10 - 30 seconds. The laser power on the sample surface was controlled by density filters to 0.3 – 0.5 mW to exclude thermal alteration of the sample. Per sample 15 measurements were conducted on different sample spots. The Raman system was calibrated against the 520.4 cm^-1 line of a Si-waver.

Table 2.2.1. List of samples. Long: longitude; Lat: latitude; s: standard deviation; n: number of mea-surements; ND: not determined.

Sample Long Lat Region R1 std R2 std RA1 std RA2 std n

L29 87.1829 28.6685 SE Tibet n.d. n.d. n.d. n.d. 0,61 0,01 1,56 0,05 15

DB45 91.1048 29.0483 SE Tibet n.d. n.d. n.d. n.d. 0,58 0,02 1,40 0,10 15

DB21 91.6362 28.9271 SE Tibet n.d. n.d. n.d. n.d. 0,59 0,02 1,47 0,12 15

DB26 92.1574 29.1036 SE Tibet 0,73 0,06 0,49 0,01 n.d. n.d. n.d. n.d. 15

DB28 92.0433 29.1454 SE Tibet 0,65 0,06 0,45 0,02 n.d. n.d. n.d. n.d. 15

DB36 91.6764 28.9881 SE Tibet 0,63 0,03 0,46 0,02 n.d. n.d. n.d. n.d. 15

L1 ca. 87.332 ca. 29.040 SE Tibet 0,40 0,05 0,37 0,02 n.d. n.d. n.d. n.d. 15

DB16 91.1906 28.7015 SE Tibet 0,38 0,03 0,36 0,01 n.d. n.d. n.d. n.d. 15

TU2 92,2606 28,8136 SE Tibet 0,32 0,04 0,33 0,02 n.d. n.d. n.d. n.d. 15

L45 88.0787 28.8651 SE Tibet 0,27 0,03 0,30 0,02 n.d. n.d. n.d. n.d. 15

L57 88.1518 28.8372 SE Tibet 0,15 0,02 0,22 0,02 n.d. n.d. n.d. n.d. 15

Kohl1 47.6617 15.6575 Austria n.d. n.d. n.d. n.d. 0,63 0,01 1,74 0,07 15

MAU 47.0451 13.2452 Austria 0,44 0,03 0,40 0,02 n.d. n.d. n.d. n.d. 15

KL2-2 50.3944 11.4014 Thuringian Forest n.d. n.d. n.d. n.d. 0,61 0,00 1,54 0,02 30

KL2-3 50.4009 11.3529 Thuringian Forest n.d. n.d. n.d. n.d. 0,60 0,00 1,47 0,02 30

KL2-4 50.362 11.4057 Thuringian Forest n.d. n.d. n.d. n.d. 0,60 0,01 1,50 0,04 30

KL2-11 50.2632 11.5218 Thuringian Forest n.d. n.d. n.d. n.d. 0,63 0,00 1,67 0,02 30

KL2-17 50.3549 11.5109 Thuringian Forest n.d. n.d. n.d. n.d. 0,63 0,00 1,73 0,01 30

KL2-18 50.3268 11.3777 Thuringian Forest n.d. n.d. n.d. n.d. 0,57 0,01 1,32 0,03 30

2.2.3. Evolution of the first-order Raman spectrum of CM. After deposition and early diagenesis, the organic content of sedimentary rocks constitutes a heterogeneous mixture of organic compounds. During organic maturation mainly O, H, N and to a lesser degree C are expelled from the organic material, changing the chemical composition and structure of the residual organic material. This process leads to an enrichment of aromatic species (for a review see Vandenbroucke and Largeau, 2007).

The aromatic species form so called ’basic structural units’ (BSU) of polyaromatic (4-10 cycles) layers, isolated or piled up by 2 – 3 units (Oberlin, 1989). The nanometer sized BSU is described by the mean stacking height (L_c) and the mean basal plane diameter (L_a). During the early stages of diagenesis and catagenesis the BSUs are randomly oriented, but start to synchronize their orientation to form molecular orientation domains (Bustin et al., 1995; Vandenbroucke and Largeau, 2007). During graphitization Lc and Laprogressively increase while at the same time the number of defects and the interplaner spacing between the graphene layers is reduced (Buseck and Huang, 1985; Wopenka and Pasteris, 1993). Thus, graphitic material of high ’crystallinity’ has few structural defects, large La and Lc values and a low interplanar spacing.

In the Raman spectra of CM, the above outlined transformation process is reflected by the change in shape (Fig. 2.1.1a) of the most prominent Raman bands in the first order spectrum (ca. 700 – 2000 cm^-1). Overall there are at least five Raman bands in the first order spectrum of CM (Fig. 2.1.1b). Following Sadezky et al. (2005) and Marshall et al. (2010) these bands are denominated as D1 (ca. 1350 cm^-1), D2 (ca. 1620 cm^-1), D3 (ca. 1500 cm^-1), D4 (ca. 1250 cm^-1) and G (ca. 1580 cm^-1). The G-band is assigned to the Raman active E_2g optical phonon in graphite (Tuinstra and Koenig, 1970; Reich and Thomsen, 2004). The D1- and D2-bands are defect-induced (Pimenta et al., 2007) and depend on the excitation energy due to double-resonant Raman scattering (Reich and Thomsen, 2004). For more information see Pócsik et al. (1998), Matthews et al. (1999), Thomsen and Reich (2000), Saito et al. (2001), Reich and Thomsen (2004)

and Pimenta et al. (2007). The D3-band supposedly originates from amorphous carbons and D4-band is attributed to sp²-sp³bonds or C-C and C=C stretching vibrations of polyene-like structures (Sadezky et al. (2005) and references therein).

Generally, the number of Raman bands decrease from low to high metamorphic conditions (Wopenka and Pasteris (1993); Yui et al. (1996); Beyssac et al. (2002b), see Fig. 2.1.1a). As the performed mode of spectral curve-fitting changes with ’crystallinity level’, the recorded spectrum is first evaluated ’by eye’ by a rough qualitative classification (Fig. 2.1.1):

’Crystallinity level 1’: This level describes poorly crystalline CM that exhibits a rather complex spectrum in which two broad, overlapping Raman bands at ca. 1350 cm^-1 (D1) and ca. 1580 to 1600 cm^-1 (G + D2) and a third band at ca. 1250 cm^-1 (D4) as shoulder on the 1350 cm^-1 band are present.

’Crystallinity level 2’: This level describes moderately to well crystalline CM. Here the spectra are less complex, as the band at ca. 1250 cm^-1 (D4) is absent. The intensities of the 1350 cm^-1 band (D1) and the overlapping region between 1350 cm^-1 and 1580 cm^-1 are decreasing while the 1580-1600 cm^-1 band (G) gets more intense and narrow. Moreover, a new band at ca. 1620 cm^-1 (D2) appears as a clear shoulder.

’Crystallinity level 3’: This level describes well crystalline CM and graphite. The spectra are simple with only the 1350 cm^-1 and 1580 cm^-1bands present. The 1350 cm^-1band is broad and of low intensity while the 1580 cm^-1band is intense and sharp (low FWHM). In case of pure graphite, only the G-band appears (Tuinstra and Koenig, 1970).

2.2.4. Spectral Evaluation. Before fitting the first order Raman spectrum of CM a back-ground correction is essential. The backback-ground is usually modeled as a linear, polynomial or spline function. The mode of such baseline is crucial for spectral curve-fitting, as all peak parameters are influenced by the baseline function. As manual baseline correction is very susceptible to subjec-tivity, a linear baseline with two control points is proposed to yield the most reproducible results.

To increase the reproducibility of manual baseline correction, the control points, which define the slope of the linear baseline function, are placed in the spectral region of 800 to 900 cm^-1and 1800 to 1900 cm^-1for all ’crystallinity levels’.

In this study the peak- and curve-fitting software Fityk (Wojdyr (2010); http://fityk.nieto.pl) is used for deconvolution of the Raman spectrum of CM into the different components (D and G bands). However, any other peak fitting software can be used for this purpose. In our approach, the position and shape of the components are detected automatically by the software. If this is not successful, the components are located manually. Because Voigt- and Lorentzian-functions are most commonly used in RSCM-thermometry, all components are modeled here as Voigt- or Lorentzian-functions with unfixed peak parameters (FWHM, height, position, area and shape).

For ’crystallinity level 1’, five components (D1, D2, D3, D4, G) result in a good fit (e.g. Sadezky et al. (2005); Lahfid et al. (2010), see Fig. 2.1.1b). For ’crystallinity level 2’, a good solution is obtained with 3 to 4 components (D1, D2, G, [D3]; see Beyssac et al., 2002a) and for ’crystallinity level 3’ only two components (D1, G) are needed. The components are assigned sequentially to the model, which is also sequentially fitted to the data with by the Levenberg-Marquardt-method (Moré, 1978). If components are displaced during fitting, take unlikely shapes or are in any other way inconsistent, the solution is rejected. Subsequently, the component parameters have to be changed and the model has to be solved again. This procedure is repeated until a satisfying fit is obtained. The complete fitting protocol is available in Appendix S1.

Once the parameters of all components are obtained, different ratios can be calculated which correlate to the maximum metamorphic temperature (Beyssac et al., 2002a; Rantitsch et al., 2004;

Rahl et al., 2005; Baziotis et al., 2006; Aoya et al., 2010; Lahfid et al., 2010). The most common are the R1- and R2-ratio (Beyssac et al. (2002a); Eq. 2.2.1 and Eq. 2.2.2) and the RA1- and RA2-ratio (Lahfid et al. (2010); Eq. 2.2.3 and Eq. 2.2.4). The R1-ratio corresponds to the height of the component divided by the height of the G-component. The integrated area of the D1-component divided by the sum of the integrated areas of the D1-, D2- and G-D1-component gives the R2-ratio. The sum of the integrated areas of the D1- and D4-component divided by the sum of the integrated areas of the D1-, D2-, D3-, D4- and G-components forms the RA1-ratio and the sum of the integrated areas of the D1- and D4-component divided by the sum of the integrated areas of the D2-, D3- and G-component is the RA2-ratio.

(2.2.1) R1 =

2.2.5. Systematic tests of the biasing factors. In order to estimate the impact of different biasing factors, a series of experiments was designed. All experiments focus on selected factors while other factors are kept constant (see Table 2.2.2).

To compare the accuracy of the measurements of the different samples the percental fraction of the standard deviation of the mean, i.e. the relative standard deviations are used as a comparative index. According to the calibration range of the calibration curves of Beyssac et al. (2002a) and Lahfid et al. (2010) it is calculated that an increment of 0.01 in the commonly used parameter-ratios R2, RA1 and RA2 (Beyssac et al., 2002a; Lahfid et al., 2010) is equivalent to 4 °C, 12 °C, and 2 °C, respectively. These values estimate the significance of each experiment with respect to the initial calibration uncertainty which is ±50 °C for Beyssac et al. (2002a).

2.3. Results and discussion 2.3.1. Spectral processing bias.

Test 1 - Influence of curve-fitting strategy on parameter ratios. The factors that influence the fitting are the signal to noise ratio, the position and slope of the baseline, type of function used (Voigt, Lorentzian, Gaussian, etc.) and the start position and shape of the inserted peak. Lahfid et al. (2010) suggested that CM spectra of low-grade metamorphic rocks (’crystallinity level 1’ of this study) should be fitted by Lorentzian functions and not by non-converging Voigt functions.

When ’crystallinity level 2 and 3’ samples are fitted, Voigt functions should be used according to Beyssac et al. (2002a).

In this experiment, the amount of scatter in the parameter-ratios due to the fitting procedure is quantified. From the different ’crystallinity levels’, single spectra with a high signal to noise

Table 2.2.2. Summary of the different tests. The influencing factors are listed in columns and the single experiments are shown in rows. The tested factors are emphasized by gray color and empty fields indicate irrelevant factors in the given test. 1 = single; M = multiple

Test Influenceof... Person Sample(percryst.level) Spot Measurement(perspot) Type(function)ofcomponents Baselinecorrection Software Sampleprep.(chemicallyenrichedCMvs.rockchips) Orientation Spectrometer

1 spectral fitting strategy 1 1 1 1 M 1

2 baseline correction 1 1 1 1 1 M

3 empirical scatter 1 1 1 M 1 1

4 evalutation software 1 M M 1 1 1 M

5 personal fitting strategy M M M 1 1 1 1

6 sample preparation 1 M M 1 1 1 1 M 1 1

7 sample heterogeneity 1 M M M 1 1 1 1 1 1

8 anisotropy 1 1 M 1 1 1 1 1 M 1

9 Raman system M 1 M 1 1 1 1 1 1 M

ratio, are repeatedly (20 times) evaluated, using the above outlined fitting strategy. To exclude a variation due to the baseline correction, baseline-corrected spectra were used for each repetition.

As the system often does not converge with automatically positioned Voigt functions, a strict fitting protocol has been used in which the functions are added manually (see Appendix S1). The same procedure has been performed with Lorentzian functions. Since the signal to noise ratio and slope of baseline is constant, the only variable is given by the initial conditions (i.e. position and shape of the components).

This experiment shows that in the ’crystallinity level 1’, the fitting strategy has a significant in-fluence on the parameter ratios. The repeated evaluation of the same spectrum, results in different RA1 and RA2-ratios when Lorentzian- and Voigt-functions are used (Table 2.3.1). Additionally, the variation in RA1- and RA2-ratios of the 20 evaluations of the same spectrum is greater when only Voigt-functions are used (Table 2.3.1).

In the ’crystallinity levels 2 and 3’, the choice of the function does not influence the results (Table 2.3.1). This is in line with previous results by Lahfid et al. (2010). However, the Lorentzian and Voigt-fits of ’crystallinity level 2 and 3’ samples result in basically the same R1 and R2 ratios.

Thus, there seems to be no justification for changing from Lorentzian- to Voigt-functions for better crystalline samples. Instead, for the sake of linearity, we suggest to use the same function for processing samples of all ’crystallinity levels’. As the fits with Lorentzian functions perform well for low and high crystalline samples, this is the function type to be used for all the fits.

Table 2.3.1. Results of test 1: comparison of R1, R2, RA1 and RA2 of the three ’crystallinity levels’

fitted with Lorentzian and with Voigt functions. r.s.d.: relative standard deviation.

’Crystallinity level 1’ - Lorentzian ’Crystallinity level 1’ - Voigt

Sample R1 R2 RA1 RA2 R1 R2 RA1 RA2

L29 1,701 0,671 0,615 1.599 1.374 0,748 0,643 1,803

r.s.d. [%] 0,1 0,0 0,0 0,0 10,2 3,4 1,7 4,6

DB45 1,281 0,641 0,586 1.414 1.200 0,679 0,618 1,619

r.s.d. [%] 0,0 0,0 0,0 0,0 2,3 0,8 0,2 0,6

DB21 1,423 0,612 0,603 1.521 1.391 0,638 0,633 1,721

r.s.d. [%] 0,0 0,0 0,0 0,0 1,1 0,5 0,2 0,6

’Crystallinity level 2’ - Lorentzian ’Crystallinity level 2’ - Voigt

DB36 0,673 0,477 — — 0,668 0,477 — —