Kaneki et al.
Progress in Earth and Planetary Science(2024) 11:35
https://doi.org/10.1186/s40645-024-00637-8

Progress in Earth andPlanetary Science

Open Access

METHODOLOGY

An automatic peak deconvolution
code for Raman spectra of carbonaceous
material and a revised geothermometerfor intermediate‑ to moderately high‑grade
metamorphism
Shunya Kaneki1* , Yui Kouketsu2 , Mutsuki Aoya3 , Yoshihiro Nakamura1 , Simon R. Wallis4 ,Yusuke Shimura1   and Ken Yamaoka1   

Abstract
Carbonaceous material (CM) undergoes progressive changes that reflect its thermal history. These changes arein general irreversible and provide valuable information for understanding diagenetic and metamorphic processesof crustal rocks. Among various approaches to quantify these changes, the R2 ratio, area ratio of specific peaks in CM
Raman spectra, is widely used to estimate the maximum temperature of intermediate- to moderately high-grade
metamorphism. The calculation of the R2 ratio requires peak deconvolution of the original spectrum, and the results
depend on the details of how this is carried out. However, a clear protocol for selecting appropriate initial conditions has not been established and obtaining a reliable temperature estimate depends at least in part on the experience and skill of the operator. In this study, we developed a Python code that automatically calculates the R2 ratio
from CM Raman spectra. Our code produces R2 ratios that are generally in good agreement with those of Aoya
et al. (J Metamorph Geol 28:895–914, 2010, https://​doi.​org/​10.​1111/j.​1525-​1314.​2010.​00896.x) for the same Raman
data, with much less time and effort than was the case in the previous studies. We have confirmed that the codeis also applicable to other previous datasets from both contact and regional metamorphic regions. The overall trendof the recalculated data indicates that samples with R2 greater than ~ 0.7 are not sensitive to the changes in CM
maturity and thus should not be used for the calibration of an R2-based geothermometer. We propose a modified geothermometer for contact metamorphism that is strictly applicable to samples with R2 from 0.023 to 0.516,
with the proviso that a laser with a wavelength of 532 nm should be used. A slight extrapolation of the newly
proposed geothermometer up to R2 of 0.57 provides a temperature estimate that is consistent with the geothermometer of Kaneki and Kouketsu (Island Arc 31:e12467, 2022; https://​doi.​org/​10.​1111/​iar.​12467); the boundary
between the two geothermometers corresponds to a temperature of 391 °C.
Keywords Automatic peak deconvolution, Carbonaceous material, Contact metamorphism, Geothermometer,
Raman spectroscopy*Correspondence:
Shunya Kaneki
kaneki.shunya@aist.go.jp
Full list of author information is available at the end of the article© The Author(s) 2024. Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which
permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to theoriginal author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images orother third party material in this article are included in the article’s Creative Commons licence, unless indicated otherwise in a credit line
to the material. If material is not included in the article’s Creative Commons licence and your intended use is not permitted by statutory
regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of thislicence, visit http://creativecommons.org/licenses/by/4.0/.

Kaneki et al. Progress in Earth and Planetary Science(2024) 11:35

Page 2 of 14

1 Introduction
The thermal history of crustal rocks provides important
insights into the diagenetic and metamorphic processesof the Earth’s interior. Carbonaceous material (CM)derived from organic matter is a common constituent of
sediments and metasedimentary rocks. CM undergoes
irreversible changes in its organochemical and crystallographic features—referred to collectively as its maturity—as temperature increases. Quantification of thesefeatures of CM has been widely used as a proxy for thethermal history experienced by the host rock. The maturity of CM has been investigated using numerous different analytical techniques such as X-ray diffraction
analysis, Raman spectrometry, and elemental analysis(e.g., Buseck and Huang 1985; French 1964; Ganz and
Kalkreuth 1987; McKirby and Powell 1974; Seifert 1978;
Sweeney and Burnham 1990; Tuinstra and Koenig 1970;
van Krevelen 1993). The changes in CM maturity are generally thought to depend on time and temperature alone
and be irreversible. However, deformation and in particular brittle deformation may also influence the results(e.g., Nakamura et al. 2015).
Over the past three decades, the use of Raman spectroscopy has dramatically increased in popularity among
researchers in the geosciences (see Fig. 1 in Henry et al.
2019 and Kouketsu 2023). A Raman spectrum of CM
typically exhibits two bands at wavenumbers of 1355 ­cm–1(disordered band) and 1575 ­cm–1 (graphite band), which
are attributed to the disordered and ordered graphitic
structures, respectively (Tuinstra and Koenig 1970). Some
Raman parameters associated with these two peaks, suchas width, position, intensity, and area, have been reported
to show close correlation with the maturity of CM (e.g.,
Beyssac et al. 2002; Kouketsu et al. 2014; Schito et al.
2017; Roberts et al. 1995), as reviewed by Henry et al.
(2019). Beyssac et al. (2002) proposed the first empiricalgeothermometer to estimate the maximum temperature T
during regional metamorphism from the area ratio of specific peaks, i.e., the R2 ratio, as follows:
T [◦ C] = −445(R2) + 641.

(1)R2 was defined as
R2 =

AD1
,
AD1 + AD2 + AG(2)

where AD1, AD2, and AG are the areas of the D1-, D2-,
and G-bands. Following their pioneering work, various
empirical geothermometers based on Raman spectra of
CM have been developed for terrestrial rocks (e.g., Aoya
et al. 2010; Kouketsu et al. 2014; Lahfid et al. 2010; Lünsdorf et al. 2017; Rahl et al. 2005). Mori et al. (2015a, 2017)
performed thermal modeling of contact metamorphism

and concluded that the heating duration required for CM
maturation to attain steady state for a given temperatureis greater than one hundred years. Throughout the present paper, we assume that such steady state conditions
have been achieved for all samples. The idea of estimating T from CM Raman spectra has been extended to fault
rocks (e.g., Hirono et al. 2015; Kaneki et al. 2016; Mukoyoshi et al. 2018) and extraterrestrial rocks (e.g., Busemann et al. 2007; Cody et al. 2008; Homma et al. 2015),
although the direct application of the geothermometers
to these types of rocks may not be appropriate, due tolarge differences in the environments where CM matures(e.g., Homma et al. 2015; Kouketsu et al. 2017).
Of the CM-based geothermometers that have been
proposed (e.g., Aoya et al. 2010; Beyssac et al. 2002;
Kouketsu et al. 2014; Lahfid et al. 2010; Lünsdorf et al.
2017; Rahl et al. 2005), those of Kouketsu et al. (2014)
and Aoya et al. (2010) (hereafter, referred to as K2014
and A2010, respectively) are among the most widelyused. Numerous studies have used one of these geothermometers, or a combination, to estimate T of the samples of interest (e.g., Bonneville et al. 2020; Cavalazzi
et al. 2021; Hickman-Lewis et al. 2020; Mori et al. 2015a,b; Nakamura et al. 2019; Shimura et al. 2021; Yamaoka
et al. 2022). The geothermometers of K2014 and A2010
are based on the widths of specific peaks and R2, respectively, and have distinct but complementary applicable temperature ranges: for K2014 it is 150–400 °C and
for A2010 it is 340–655 °C. Although both approaches
require peak deconvolution of the measured Raman
spectra, the criteria for setting the initial conditions fornonlinear least-squares fitting are not clearly presented
in the original papers, meaning temperature estimates
could vary depending on the skill and experience of theanalysts. Kaneki and Kouketsu (2022) (hereafter, referred
to as KK2022) addressed this issue by developing a
Python code that automatically performs peak deconvolution to calculate the widths of the D1- and D2-bands.
The main objective of the present study is to develop asimilar code to calculate R2.
We first describe the fitting procedures implemented
in our code to reproduce R2 of A2010 for the same dataset. The results obtained are then compared with thosereported by A2010 for the validation of the code. Finally,
we assess the utility of our proposed approach including
discussion of its limitations based on results of applying
it to datasets other than those of A2010.

2 Methods
In the development of our code, we reanalyzed the same
Raman spectral dataset as that measured and analyzed in
A2010. The basic structure of the code is similar to that

Kaneki et al. Progress in Earth and Planetary Science

3

(b)

×103

1.0

Raw
Baseline
Baseline-corrected

2

Intensity

Page 3 of 14

Normalized intensity(a)(2024) 11:35

1

0.8
0.6
0.4
0.2
0

0
1200

1400

1600

1200

Raman shift, cm–1

(c)

0.8
0.6
0.4
0.2

1600

R2 = 0.998
1.0

MeasuredD1-bandD2-bandD3-bandG-bandCalculated

Normalized intensity

1.0

Normalized intensity(d)R2 = 0.991

1400

Raman shift, cm–1

0.8
0.6

MeasuredD1-bandD2-bandD3-bandG-bandCalculated

0.4
0.2
0

0
1200

1400

1600

Raman shift, cm–1

1200

1400

1600

Raman shift, cm–1

Fig. 1 Data processing procedures of our code for calculating R2 from a raw Raman spectrum. One of the spectra of the N29 sample is shown
here as an example. a Raw spectrum, linear baseline, and baseline-corrected spectrum; b normalized baseline-corrected spectrum (referred
to as ‘measured spectrum’ in this paper); c initial conditions of calculated spectrum; and d calculated spectrum after nonlinear least-squares fitting.
Gray areas in a indicate the spectral ranges used to determine a linear baseline. The coefficients of determination R2 in c and d were calculatedusing Eq. (8)of Fitting D in KK2022. Note that we developed our codeby trial and error to reproduce the analytical results ofA2010 as closely as possible for the same dataset. For thispurpose, we tried to follow the analytical procedures ofA2010 as closely as possible (e.g., fitting function, orderof a baseline, number of peaks, etc.).
2.1 Data

We first briefly describe the samples and spectral dataof A2010. Readers can refer to A2010 for more detailed
information, including their experimental setups. Since
the present study focuses on the development, validation, and application of our new code, we will not examine the procedure and methodology used in the original

experiments reported in A2010. A total of ten analyzed
samples was studied, of which two, K08 and N33, were
excluded from the final geothermometer of the present
study (Table 1). The detailed reasons for this selection
are discussed in Sect. 4.3. All samples experienced contact metamorphism, and their independently estimated T
range from 340 ± 25 °C to 655 ± 25 °C. A2010 measured atleast 50 spectra for each sample and a total of 809 spectra. Most of the spectra show two distinct peaks of thedisordered and graphite bands (Fig. 1a), except for highly
matured CM.
A2010 examined the effect of measurement conditions on the analytical results by varying the magnification of the objective lens (50 × or 100 ×), sample type

Kaneki et al. Progress in Earth and Planetary Science(2024) 11:35

Page 4 of 14

Table 1 Summary of temperature, measurement condition, and R2 ratio of each sample used for the geothermometer.
Sample ID

K08*

N33*

N30

N29

N27

N9

N6

N3

N1

K04

T , °C**

340 ± 25387 ± 36410 ± 30443 ± 29475 ± 26550 ± 19584 ± 18604 ± 25625 ± 32655 ± 25

Measurement condition

Aoya et al. (2010)

This study

Mean

SD***

n***

Mean

SD***

SE***

n***

Thin section, ×100, parallel to c-axis

0.678

0.022

21

0.729

0.032

0.007

21

Thin section, ×100, normal to c-axis

0.696

0.012

31

0.743

0.025

0.005

31

Total

0.689

0.019

52

0.737

0.029

0.004

52

Chip, ×100, parallel to c-axis

0.580

0.064

24

0.586

0.063

0.013

24

Thin section, ×50, parallel to c-axis

0.601

0.067

19

0.610

0.072

0.017

19

Thin section, ×100, parallel to c-axis

0.684

0.046

22

0.707

0.043

0.009

22

Thin section, ×100, normal to c-axis

0.615

0.076

29

0.629

0.066

0.013

29

Total

0.616

0.078

96

0.633

0.075

0.008

94

Thin section, ×50, parallel to c-axis

0.483

0.067

23

0.474

0.068

0.015

23

Thin section, ×100, parallel to c-axis

0.470

0.093

25

0.479

0.095

0.019

25

Thin section, ×100, normal to c-axis

0.516

0.056

76

0.533

0.068

0.008

76

Total

0.504

0.064

122

0.516

0.076

0.007

123

Thin section, ×100, parallel to c-axis

0.399

0.035

30

0.389

0.039

0.007

30

Thin section, ×100, normal to c-axis

0.421

0.037

28

0.416

0.040

0.008

28

Total

0.407

0.037

58

0.397

0.040

0.005

57

Thin section, ×100, parallel to c-axis

0.328

0.034

27

0.326

0.033

0.006

27

Thin section, ×100, normal to c-axis

0.348

0.035

27

0.347

0.029

0.006

26

Total

0.336

0.036

55

0.336

0.033

0.005

53

Chip, ×100, parallel to c-axis

0.279

0.047

23

0.226

0.058

0.012

24

Thin section, ×50, parallel to c-axis

0.213

0.058

27

0.173

0.048

0.010

26

Thin section, ×100, parallel to c-axis

0.178

0.068

34

0.151

0.061

0.011

34

Thin section, ×100, normal to c-axis

0.234

0.064

28

0.249

0.050

0.010

27

Total

0.222

0.068

112

0.193

0.068

0.006

112

Chip, ×100, parallel to c-axis

0.156

0.059

20

0.179

0.073

0.016

22

Thin section, ×100, parallel to c-axis

0.124

0.063

29

0.110

0.052

0.010

29

Total

0.137

0.062

49

0.136

0.067

0.010

50

Chip, ×100, parallel to c-axis

0.152

0.083

28

0.087

0.066

0.013

27

Thin section, ×100, parallel to c-axis

0.141

0.094

30

0.089

0.068

0.013

28

Thin section, ×100, normal to c-axis

0.113

0.063

36

0.098

0.069

0.012

36

Total

0.123

0.071

89

0.092

0.067

0.007

91

Chip, ×100, parallel to c-axis

0.073

0.068

28

0.065

0.059

0.011

28

Thin section, ×50, parallel to c-axis

0.073

0.056

25

0.056

0.035

0.007

24

Thin section, ×100, parallel to c-axis

0.132

0.071

30

0.089

0.062

0.012

30

Total

0.094

0.071

83

0.070

0.054

0.006

82

Thin section, ×100, parallel to c-axis

0.010

0.013

27

0.016

0.017

0.003

28

Thin section, ×100, normal to c-axis

0.010

0.023

28

0.030

0.037

0.007

29

Total

0.011

0.020

56

0.023

0.028

0.004

58

*Not used for our regression analysis.
**Temperature and its error are from Aoya et al. (2010), except the K04 sample (Wada and Suzuki 1983).
***SD, SE, and n indicate standard deviation, standard error, and number of data used to calculate them.

(thin section or polished chip), and incident angle of laser
beam (normal or parallel to the c-axis of the crystal structure of graphite). They concluded that these variations in
the experimental conditions had a limited effect on thecalculated R2 (Fig. 8 in A2010).

2.2 Code description
2.2.1 Fitting function

A2010 employed a Voigtian, V , in their fitting procedures, which is defined as the convolution of a Gaussian,G , and a Lorentzian, L, as follows:Kaneki et al. Progress in Earth and Planetary Science(2024) 11:35

Page 5 of 14∞

V (x; ω, σ , γ , A) = ∫ G x′ ; ω, σ , A L x − x′ ; ω, γ , A dx′ ,−∞(3)where ω is the center position, A is the area, and σ and γ
represent the widths of G and L, respectively. G and L areexpressed asA(x − ω)2
exp −,G(x; ω, σ , A) = √(4)2σ 22π σ 2L(x; ω, γ , A) =Aγ.
π (x − ω)2 + γ 2(5)Although we need to determine the values of the arguments of V as an initial condition for nonlinear leastsquares fitting, it is difficult to constrain σ and γ from a
raw Raman spectrum. Some previous studies (e.g., Kouketsu et al. 2014; Lünsdorf and Lünsdorf 2016) employeda pseudo-Voigtian, PV , for peak deconvolution. One candefine a typical PV as follows:

describe the linear baseline using the data in the ranges
1200–1250 ­cm–1 and 1700–1750 ­cm–1 (Fig. S2b). We
defined a ‘mature’ spectrum as one with the maximum
intensity at 1320–1380 ­cm–1 (disordered band) less thanhalf that at 1570–1600 ­cm–1 (graphite band) after thebaseline correction. If a given spectrum satisfies this criterion, the linear baseline is described using the spectral
data in the ranges 1200–1250 ­cm–1 and 1700–1750 ­cm–1,
rather than 1100–1150 ­cm–1 and 1700–1750 ­cm–1.
The baseline-corrected spectrum was then normalized such that its maximum intensity in the range 1100–
1750 ­cm–1 becomes unity (Fig. 1b). This normalizedbaseline-corrected spectrum is referred to as the ‘measured spectrum’ in the present study, and we used thisspectrum for the peak deconvolution.
2.2.3 Initial conditions and peak deconvolution

To perform the nonlinear least-squares fitting of themeasured spectrum, the initial conditions for each peakshould be determined. The approach to setting the initial conditions in our code is similar to that of Fitting D




Ŵ
Ŵ
√
PV (x; ω, Ŵ, A, η) = (1 − η)G x; ω,, A + ηL x; ω, , A ,2
2 2 ln (2)where η is the mixing ratio of G and L, and Ŵ is full widthat half maximum (FWHM). PV in Eq. (6) takes its maximum intensity I at x = ω as follows:ln (2)
12A.
+ηI=(1 − η)(7)Ŵππ
Note that the first approximation to all arguments of
PV can be easily determined from a raw Raman spectrum. Since we are interested in the area ratio and the difference in A between V and PV is less than 1% (Fig. S1),
we concluded that the usage of PV instead of V shouldhave negligible effect on the calculated R2.
2.2.2 Background correction and normalization

Background correction is required for each raw spectrum
prior to peak deconvolution. Although it is obvious that
A2010 employed a linear baseline, the details of the algorithm used to describe it are not clearly presented. In ourcode, we determined the slope and intercept of the linear
baseline by a linear least-squares method using the raw
spectral data in the ranges 1100–1150 ­cm–1 and 1700–
1750 ­cm–1 (Fig. 1a). However, for samples with a highdegree of maturity and thus low effects of fluorescence, asmall bulge sometimes appears around 1100 ­cm–1 (attributable to the overlapping of spectra from quartz or calcite
grains with the spectra from CM), which makes this procedure inapplicable (Fig. S2a). In this case, it is better to(6)

in KK2022. The D1-, D2-, D3-, and G-bands of PV are
assumed to be always present in all the measured spectra. The maximum intensity of the measured spectrum
in the range 1320–1380 ­cm–1 and the associated wavenumber are employed as the initial values for I and ω ofthe D1-band, respectively. The same analysis is applied
to the G-band in the range 1570–1600 ­cm–1. The initial Ŵ of the D1-band is set as the distance between two
wavenumbers where the intensities are half the initial I
in the ranges from 1200 ­cm–1 to initial ω and from initial ω to 1500 ­cm–1. That of the G-band is set to twice thedistance between initial ω and the wavenumber at the
point where the intensity is half the initial I in the range
from 1500 ­cm–1 to initial ω . These approaches to determining appropriate values for the initial Ŵ allow us toanalyze Raman spectra of CM with different maturities
using a single code, eliminating the subjectivity owing
to the selection of the fitting method (a detailed discussion of this point can be found in Sect. 4.3 in KK2022).
The D1- and G-bands are initially assumed to be a pure
Lorentzian, and thus η = 1. To determine the initial conditions of the other two bands, we used the residual spectrum after subtracting the initial D1- and G-bands fromthe measured spectrum (Fig. S3). For the D2-band, 90%of the maximum intensity of the residual spectrum in therange 1610–1630 ­cm–1 and the associated wavenumber
are employed as the initial values for I and ω , respectively. The same approach is employed for the D3-band inKaneki et al. Progress in Earth and Planetary Science(2024) 11:35

the range 1500–1540 ­cm–1. The initial Ŵ of the D2-bandis set to twice the distance between the initial ω and thewavenumber where the intensity is half the initial I in therange from initial ω to 1700 ­cm–1. That of the D3-band is
set to 130 ­cm–1, as in Fitting D of KK2022. The D2- andD3-bands are initially assumed to be half Gaussian andhalf Lorentzian, and thus η = 0.5. We calculated the initial A values of these four bands using Eq. (7).
After determining the initial conditions of the calculated spectrum (Fig. 1c), we performed nonlinear leastsquares fitting on the measured spectrum in the spectral
range 1100–1750 ­cm–1 (Fig. 1d). We used scipy.optimize.
curve_fit function with the Trust Region Reflective algorithm (SciPy v1.10.1), which is applicable even when thebounds on the model parameters are given. We defined
the cost function as half the sum of squares of the residual of the measured spectrum after subtracting the calculated spectrum. The maximum number of iterations was
set to 10,000. All values of tolerances during the iterations were set to ­10–8. Ŵ and A of the four bands can takeany positive values. The value of η can move from zeroto unity. The ω values can vary 1300–1400 ­cm–1, 1600–
1650 ­cm–1, 1450–1550 ­cm–1, and 1550–1600 ­cm–1 forthe D1-, D2-, D3-, and G-bands, respectively. The coefficient of determination, R2, was calculated for each fitting
curve:
2
  meas− IicalcIii
R2 = 1 −(8) ,  calc
calc 2I−I
ii
where Iimeas and Iicalc are the intensities of the measured
and calculated spectra, and I calc is the mean intensity of
the calculated spectrum in the range 1100–1750 ­cm–1.
For each sample, we calculated the values of the mean,
standard deviation (SD), and standard error (SE) of R2 ofthe calculated spectra after optimization. As in A2010,
we then excluded the data with R2 deviating by morethan ± 2SD from the mean value. After this treatment, werecalculated the mean, SD, and SE values.

3 Results
3.1 Performance of the code

Figure 2 shows examples of the initial conditions together
with the measured spectra for each sample. The initial fitting curves show a good match with the measured datawith R2 > 0.97. After optimization, the R2 of the calculated spectra exceed 0.99 regardless of the sample (Fig. 3).
Since the fitting is performed automatically, the same
analytical results are obtained for the same spectral data,irrespective of the analyst. Our code typically analyzes
one spectrum in less than one second using an ordinary

Page 6 of 14

laptop computer and thus greatly reduces the time andeffort required to arrive at a result.
3.2 Comparison with Aoya et al. (2010)

Figure 4 compares R2 of the present study with thosereported by A2010. Error bars indicate the SD. As
described in Sect. 2.1, A2010 investigates the effect ofmeasurement conditions on the calculated R2 by varying the sample type, magnification of an objective lens,and incident angle of laser beam. Figure 4a shows thecomparisons of R2 for each measurement condition. Our
code-generated results are generally in good agreementwith those of A2010 within the range of SD, regardless ofthe measurement conditions. When R2 exceeds ~ 0.6, theD2- and G-bands were not well decomposed for somesamples (e.g., K08 in Figs. 2 and 3), and the results of thepresent study showed significantly larger R2 than those of
A2010 with differences exceeding one SD. These trends
also hold even when the data measured under different
conditions are summed and then analyzed to obtain representative R2 for each sample (Fig. 4b). These results are
consistent with the negligible effects on R2 by varying themeasurement conditions as examined by A2010. All data
are summarized in Table 1.
3.3 Modified geothermometer

Instead of a linear model proposed by Beyssac et al.
(2002) (Eq. (1)), A2010 presents a quadratic equation toexpress the relationship between R2 and T :T [◦ C] = 221.0(R2)2 − 637.1(R2) + 672.3.

(9)

Ideally, R2 calculated using our code should be suitable
to use Eq. (9) for estimating T . However, there are differences in methodology of the present study that makethis simple approach inappropriate. (1) The present study
uses a more limited sample set than A2010 with the data
not only for N33 but also for K08 removed when defining the calibration curve. (2) Our code-generated R2 were
not identical to those of A2010 even for the same spectra(Fig. 4 and Table 1). (3) The uncertainties in both R2 andT should be taken into consideration when formulating
the calibration curve. Due to these differences, we concluded that the previous geothermometer (Eq. 9) needs
to be modified when employing our code, which is themain issue addressed in this subsection.
Since the calibration data have uncertainties in both
R2 and T (Table 1), the coefficients of the quadratic function cannot be calculated as a unique solution to a linear
equation. Therefore, we performed the Deming regression (Text S1), which is applicable to a regression problem that requires iterative treatments to optimize theKaneki et al. Progress in Earth and Planetary Science

Normalized intensity

K08 (T = 340 ± 25 °C)
1.0 R2 = 0.981

N33 (T = 387 ± 36 °C)

N30 (T = 410 ± 30 °C)

N29 (T = 443 ± 29 °C)

R2 = 0.989

R2 = 0.989

R2 = 0.991

N9 (T = 550 ± 19 °C)

N6 (T = 584 ± 18 °C)

N3 (T = 604 ± 25 °C)

R = 0.976

R = 0.996

R2 = 0.992

0.6
0.4
0.2
0

1.0 R = 0.9952

Normalized intensity

Page 7 of 14

0.8

N27 (T = 475 ± 26 °C)

2

2

0.8
0.6
0.4
0.2
0
N1 (T = 625 ± 32 °C)

Normalized intensity(2024) 11:35

K04 (T = 655 ± 25 °C)

1.0 R2 = 0.998

R2 = 0.998

1200

1400

1600

Raman shift, cm–1

1200

1400

1600

Raman shift, cm–1

MeasuredD1-bandD2-bandD3-bandG-bandCalculated

0.8
0.6
0.4
0.2
0
1200

1400

1600

Raman shift, cm–1

1200

1400

1600

Raman shift, cm–1

Fig. 2 Examples of the measured spectra and initial conditions of the calculated spectra for each sample. The coefficients of determination R2 were
calculated using Eq. (8)model parameters. Samples K08 and N33 were excluded
from our regression analysis, owing to their inappropriateness as the calibration samples (see Sect. 4.3 for thedetailed discussion). Following KK2022, we assume thatthe errors in T of the calibration samples can be represented by the SE. We employed SE rather than SD as theuncertainties of each data point, since SE is required toperform the Deming regression (Text S1). Finally, theregression curve obtained is expressed as follows:

T [◦ C] = 393.6(R2)2 − 716.9(R2) + 671.8.

(10)

The geothermometer of Eq. (10) is strictly only applicable to samples with R2 from 0.023 to 0.516 (Table 1).
Figure 5 shows the regression curve (Eq. (10)), its 95%
prediction interval (Text S2), and the data points used
to calculate them (Table 1). Although A2010 defined

the uncertainty of Eq. (9) as the maximum difference
between the estimated and known T of the calibration samples, we employed the prediction interval forthe uncertainty of Eq. (10) owing to its clear statisticaldefinition. Therefore, a direct comparison of the uncertainties between Eqs. (9) and (10) is inappropriate. The
regression curve of Eq. (10) shows a very good fit to thedata with R2 = 0.998, which is similar to that of A2010( R2 = 0.995). The optimization of the quadratic function by the iterative treatment affected the results oftemperature estimation only modestly, with ~ 1 °C difference in T between the regression curves after 0th and10th iterations (Fig. S4). At a given R2, the calculated T
using Eq. (10) differs by a maximum of ~ 10 °C from theresults of using Eq. (9) (Fig. 6), which is similar to the
95% prediction interval of Eq. (10) (~ 15 °C). Note that

Kaneki et al. Progress in Earth and Planetary Science

Normalized intensity

K08 (T = 340 ± 25 °C)
1.0 R2 = 0.996

N33 (T = 387 ± 36 °C)

N30 (T = 410 ± 30 °C)

N29 (T = 443 ± 29 °C)

R2 = 0.997

R2 = 0.998

R2 = 0.998

N9 (T = 550 ± 19 °C)

N6 (T = 584 ± 18 °C)

N3 (T = 604 ± 25 °C)

R = 0.998

R = 0.999

R2 = 0.997

0.6
0.4
0.2
0

1.0 R = 0.998
2

Normalized intensity

Page 8 of 14

0.8

N27 (T = 475 ± 26 °C)

2

2

0.8
0.6
0.4
0.2
0
N1 (T = 625 ± 32 °C)

Normalized intensity(2024) 11:35

K04 (T = 655 ± 25 °C)

1.0 R2 = 0.999

R2 = 0.999

1200

1400

1600

Raman shift, cm–1

1200

1400

1600

Raman shift, cm–1

MeasuredD1-bandD2-bandD3-bandG-bandCalculated

0.8
0.6
0.4
0.2
0
1200

1400

1600

Raman shift, cm–1

1200

1400

1600

Raman shift, cm–1

Fig. 3 Examples of the measured spectra and the calculated spectra after optimization for each sample. The coefficients of determination R2 were
calculated using Eq. (8)Eq. (10) was established using the samples which underwent contact metamorphism, and its direct application to the samples which experienced other processes
including regional metamorphism may not be appropriate, although our code may be useful to establish a new
R2-based geothermometer for these situations.

4 Discussion
4.1 Deviation from Aoya et al. (2010)

R2 of the K08 sample in the present study was larger thanthat of A2010 by an amount that exceeds one SD (Fig. 4
and Table 1). Let us assume that a measured spectrum isperfectly explained by a calculated spectrum, i.e., R2 in
Eq. (8) is strictly unity. It then follows from Eq. (2) that

R2 =

AD1
,
Aall − AD3(11)where Aall represents the entire area of the measured
spectrum in the spectral range 1100–1750 ­cm–1 andthus takes a constant value. The Raman spectra of theK08 sample typically show a clear D4-band at around
1200 ­cm–1 and a strong ‘saddle’ intensity at around1500 ­cm–1 (e.g., Figs. 2 and 3). This indicates that the
optimized AD1 and AD3 values may be strongly influenced
by which algorithm was employed to perform the peakdeconvolution. Therefore, the observed difference in R2
between these two studies may be qualitatively explainedKaneki et al. Progress in Earth and Planetary Science

Page 9 of 14

(b)

0.8
Chip, 100×, parallelThin section, 50×, parallelThin section, 100×, parallelThin section, 100×, normal

0.6

R2 ratio in this study

R2 ratio in this study(a)(2024) 11:35

0.4

0.2

0

0

0.2

0.4

0.6

0.8

R2 ratio in Aoya et al. (2010)

0.8

0.6

0.4

0.2

0

0

0.2

0.4

0.6

0.8

R2 ratio in Aoya et al. (2010)

Fig. 4 Comparison of R2 of the present study with those of Aoya et al. (2010). The results for a each measurement condition and b each sample areshown. All plotted data are summarized in Table 1. When R2 of both studies are identical, the data are plotted on the gray line. Error bars indicate
the standard deviation

700
T [°C] = 393.6(R2)2 – 716.9(R2) + 671.8(R2 = 0.998)

T, °C

600

500
Not used

400
Regression curve
95% prediction interval

300

0

0.2

0.4

0.6

0.8

R2 ratio
Fig. 5 Geothermometer modified from that of Aoya et al. (2010).
R2 calculated using our code are plotted against T experienced
by each of the samples examined by Aoya et al. (2010), together
with the obtained geothermometer (Eq. 10) and its 95% prediction
interval. All plotted data are summarized in Table 1. The black
points indicate the results for the samples K08 and N33, which were
excluded from the regression analysis (see Sect. 4.3 for the detailed
reasons). Error bars for R2 are the standard error, which are almost
hidden by the data points owing to their small values. Error
bars in T were also taken as the standard errors. The coefficientof determination R2 was calculated using Eq. (8)if the optimization algorithm of the present study yields
larger AD1 and AD3 than A2010.
4.2 Applicability of the code to other datasets

Since our code was developed by trial and error to reproduce R2 of A2010 as closely as possible, its applicability todatasets other than those of A2010 should be examined.
We reanalyzed the raw Raman spectral data of the natural
samples measured and analyzed by Kouketsu et al. (2019,2021), Shimura et al. (2021), and Yamaoka et al. (2022).
These studies measured the Raman spectra of natural
CM under the same experimental conditions as A2010,
performed peak deconvolution by manually setting theinitial conditions, and calculated R2. R2 of the samples tobe reanalyzed were reported to range from 0.066 ± 0.041
to 0.655 ± 0.021 (Table S1).
The comparison of R2 obtained in the previous and present studies is shown in Fig. 7. Error bars represent theSD. All the data obtained are summarized in Table S1. R2of the present study generally agree with those of the previous studies within the margin of SD. However, we notethat when R2 is greater than ~ 0.6, our code yielded significantly larger R2 than the previous studies exceeding one
SD, something which was also recognized for the samples used in A2010 (Fig. 4). These results clearly indicatethat our code is applicable to datasets other than those of
A2010, although the code-generated results significantly
differ from the previous analysis when R2 > ∼ 0.6.
4.3 Sensitivity of R2 as a proxy for CM maturity and its
implications for use as a geothermometer

K2014 reported that FWHM of the D1-band, ŴD1, linearly decreases with increasing T while it shows almost
no change after reaching ~ 40 ­cm–1. Therefore, we canconsider ŴD1 > 40 cm−1 as a necessary condition thatŴD1 has its sensitivity to changes in CM maturity. This
subsection focuses on the similar sensitivity test for R2

Kaneki et al. Progress in Earth and Planetary Science(a)(2024) 11:35

Page 10 of 14

(b)

700

10

Aoya et al. (2010) (Eq. (9)):
T [°C] = 221.0(R2)2 – 637.1(R2) + 672.3

Difference in T, °C

T, °C

600

500

This study (Eq. (10)):
T [°C] = 393.6(R2)2 – 716.9(R2) + 671.8

400

300

0

0.2

0.4

0.6

0.8

R2 ratio

5

0

–5

–10

0

0.2

0.4

0.6

0.8

R2 ratio

Fig. 6 Comparison of the present geothermometer with that of Aoya et al. (2010). a Geothermometers proposed by Aoya et al. (2010) (Eq. 9)
and the present study (Eq. 10), and b their difference

and discusses its implications for practical temperature
estimation.
We reexamined the raw spectral data of Aoya et al.
(2010), Kouketsu et al. (2019, 2021), Nakamura et al.
(2019), Shimura et al. (2021), and Yamaoka et al. (2022).
R2 was calculated using our code, while ŴD1 was determined using the code for Fitting E of KK2022. Figure 8
compares R2 and ŴD1 for each sample. Figure 8b showsan enlarged view of the gray area in Fig. 8a. Error bars
indicate the SD. The large SD of ŴD1 of highly matured

R2 ratio in this study

0.8
Kouketsu et al. (2019)
Kouketsu et al. (2021)
Shimura et al. (2021)
Yamaoka et al. (2022)

0.6

0.4

0.2

0

0

0.2

0.4

0.6

0.8

R2 ratio in the previous studies
Fig. 7 Applicability of our code to datasets other than thoseof Aoya et al. (2010). R2 calculated using our code are plottedagainst those reported in the previous studies. All plotted data aresummarized in Table S1. The gray line indicates exact agreement
between the values calculated in the previous and present studies.
Error bars indicate the standard deviation

CM with R2 of ~ 0.1 were attributed to the disappearance of the D1-band (e.g., K04 in Figs. 2 and 3). All thedata obtained are summarized in Table S2. The results
obtained show the same trend regardless of datasets(Fig. 8a). Within the region for low-grade material (upper
right part in Fig. 8a), ŴD1 decreases significantly withincreasing CM maturity while R2 shows almost constantvalues of 0.6–0.8. R2 starts to decrease from ~ 0.7 toward
zero when ŴD1 decreases to ~ 55 ­cm–1 (Fig. 8b). After ŴD1
reaches ~ 40 ­cm–1 at R2 of ~ 0.55, only R2 shows a systematic decrease with CM maturity. Based on these results,
we suggest that R2 < ∼ 0.7 is a necessary condition for
R2 to be sensitive to changes in CM maturity. In otherwords, samples with R2 ≥ ∼ 0.7 should not be employedin the calibration of an R2-based geothermometer. Thisis the reason why our regression curve (Eq. 10) did notinclude the K08 sample (R2 = 0.737 ± 0.029), which enables our geothermometer to yield better results aroundR2 = 0.516 than that of A2010 (Eq. 9) who included theK08 sample (Fig. 6). Since very low-grade CM sometimes
exhibits R2 of ~ 0.7 (e.g., dataset of Nakamura et al. 2019
in Fig. 8a), the condition R2 < ∼ 0.7 alone can be misleading for the assessment of CM maturity. It is important to check the consistency of the results obtained withother observations, including the strong influence of fluorescence (e.g., Schito et al. 2017) and the existence of theD4-band (e.g., Kouketsu et al. 2014), both of which areabsent in the spectra of intermediate-grade CM (Fig. S5).
As described in Sect. 3.3, our modified geothermometer is determined in the range 0.023 ≤ R2 ≤ 0.516 (Fig. 5
and Table 1). However, the geothermometers based onŴD1, including that of KK2022, are applicable only when

Kaneki et al. Progress in Earth and Planetary Science(2024) 11:35(a) 130(b)

FWHM of the D1-band, cm–1

110
100
90

FWHM of the D1-band, cm–1

Aoya et al. (2010)
Kouketsu et al. (2019)
Kouketsu et al. (2021)
Shimura et al. (2021)
Yamaoka et al. (2022)
Nakamura et al. (2019)

120

(b)

80
70
60
50
40
30

Page 11 of 14

0

0.2

0.4

0.6

0.8

R2 ratio

80

Proposed upper limit (R2 = 0.57)
for extrapolation of Eq. (10)

75
70
65
60

Eq. (12)

55
50
45
40
35
0.5

0.6

0.7

0.8

R2 ratio

Fig. 8 Comparison of R2 and ŴD1 for various natural samples. The results for the samples of Aoya et al. (2010), Kouketsu et al. (2019, 2021), Nakamura
et al. (2019), Shimura et al. (2021), and Yamaoka et al. (2022) are shown. All plotted data are summarized in Table S2. R2 were calculated using ourcode, while ŴD1 were determined using the code for Fitting E of Kaneki and Kouketsu (2022). A close-up of the gray area in a is shown in b. The
black-dashed line in b indicates the proposed upper limit R2 = 0.57 for the extrapolation of our geothermometer (Eq. 10). When the calculated(2022) (Eq. (12)) and the present study(R2, ŴD1 ) satisfy Eq. (13) (indicated by the gray line in b), the geothermometers
 of Kaneki and Kouketsu
(Eq. (10)) yield the same T . Star symbol in b represents the point (R2, ŴD1 ) = 0.57, 41.3 cm−1 where the two geothermometers yield T = 391 °C.
Error bars indicate the standard deviation

R2 exceeds ~ 0.55, owing to the insensitivity of ŴD1 tochanges in CM maturity for lower R2 (Fig. 8). Therefore,
in a strict sense, a combination of the code-based geothermometers of KK2022 and the present study is incapable ofestimating T of samples with R2 from 0.516 to ~ 0.55 andŴD1 of ~ 40 ­cm–1. A straightforward solution to this problem would be to refine our modified geothermometer further by considering additional calibration data of samples
whose T are known, ŴD1 are ~ 40 ­cm–1, and R2 are larger
than ~ 0.55 while smaller than ~ 0.7. Although the N33
sample almost satisfies this requirement (Table 1), thebimodal distribution in its R2 histogram (Fig. S6) indicates
its inappropriateness as a calibration sample. An alternative approach is to simply extrapolate the applicable rangeof Eq. (10), as implicitly done in some previous studies (e.g.,
Kaneki and Kouketsu 2022; Kouketsu et al. 2014; Rahl et al.
2005). In this case, we need to consider an appropriate
extent for the extrapolation of our geothermometer. One
potential clue can be obtained when we try to smoothly
connect the geothermometers of KK2022 and the present
study on the data trend shown in Fig. 8. The ŴD1-based geothermometer of KK2022 isT [◦ C] = −2.30 ŴD1 cm−1 + 486.
(12)
By combining Eqs. (10) and (12), we obtain
ŴD1 cm−1 = −171.1(R2)2 + 311.7(R2) − 80.8.

(13)

When the calculated (R2, ŴD1 ) satisfy Eq. (13), thetwo geothermometers yield the same T . For example, if
we assume the connecting T is 391 °C, the corresponding R2 and ŴD1 are 0.57 and 41.3 ­cm–1, respectively(Eqs. (10) and (12)). It is noteworthy that the point(R2, ŴD1 ) = 0.57, 41.3 cm−1 (indicated by star symbol in Fig. 8b) plots well within the data trend, which isdrawn by a number of recalculated data both from contact and regional metamorphisms. Therefore, the slight
extrapolation of our geothermometer up to R2 = 0.57
(within the standard deviation of the lowest-T sample forour calibration curve, N30) provides a consistent switching to the lower-T geothermometer of KK2022, with aboundary corresponding to T = 391 ◦ C. In addition, setting R2 = 0.57 as a boundary between the two geothermometers enables an unambiguous classification of all
the data in Fig. 8 into one of two groups: using our geothermometer (Eq. (10)) if R2 ≤ 0.57 and using KK2022(Eq. (12)) if R2 > 0.57.
4.4 Effect of laser wavelength

Differences in the laser wavelength employed in the
Raman system can significantly affect the spectral intensity and subsequent analytical results (e.g., Wang et al.
1990; Ferrari and Robertson 2001). K2014 and KK2022
demonstrated that for the Raman data of a specific sample measured under different experimental conditions(i.e., spectrometer, laser wavelength, pinhole size, andgrating), some representative Raman parameters (i.e., Ŵ,(2024) 11:35

5 Conclusions
In this study, we have developed a Python code that automatically performs peak deconvolution of the Raman
spectra of carbonaceous material and calculates the R2
ratio. A close comparison with the results of Aoya et al.
(2010) shows very good agreement for the same original
data unless R2 exceeds ~ 0.6. Based on a Deming regression analysis, we propose a modified set of coefficients

R1 and R2 ratios of 633 nm laser

(a)

2.5
2.0
1.5
1.0
0.5
0

R1 ratio
R2 ratio

0

0.5

1.0

1.5

2.0

2.5

R1 and R2 ratios of 532 nm laser(b)
FWHM of 633 nm laser, cm–1

ω , intensity ratio, and area ratio) showed consistent values when the laser wavelength was the same, while theparameters other than Ŵ differed significantly for different laser wavelengths. These results clearly suggest thatthe effect of the laser wavelength on the spectral parameters dominates over the other variable experimental
conditions. Since these two studies analyzed the data of
the sample with T = 165 ◦ C and 330 ◦ C, which is outside
the temperature range of the calibration samples used inour regression analysis, this subsection investigates thedependence of the analytical results using our code on
the laser wavelength for samples with higher T .
We measured the Raman spectra of the N30 sample(T = 410 ± 30 ◦ C) under the same experimental conditions as for A2010, except that we used a laser with awavelength of 633 nm instead of the original 532 nm.
Using our code, we calculated the following six Raman
parameters: intensity ratio of the D1-band to the G-band(i.e., R1 ratio), R2, Ŵ of the D1- and G-bands, and ω ofthe D1- and G-bands. The results obtained are comparedwith those of a 532 nm laser in Fig. 9. Error bars indicate
the SD. All the data obtained are summarized in Table S3.
Of the six parameters examined, only Ŵ of the D1- andG-bands were consistent between two laser wavelengths
within the margin of SD (Fig. 9b), similar to the results
reported in K2014 and KK2022. The other four parameters, including R2, showed significant differences exceeding one SD when the laser wavelength was changed from
532 to 633 nm (Fig. 9a, c). The observed larger R2 for the
longer laser wavelength (Fig. 9a) may be attributed to (i)use of the longer laser wavelength resulting in the larger
R1 ratio (e.g., Ferrari and Robertson 2001), and (ii) the
positive correlation between the R1 and R2 ratios (e.g.,
Aoya et al. 2010). The mean T calculated using Eq. (10)were 407 and 364 °C for a 532 and 633 nm laser, respectively, and their difference exceeded 95% prediction intervals. Considering these facts, geothermometers based on
R2 cannot be reliably applied to spectra that were measured using a different laser wavelength from the one used
to establish it. Therefore, the usage of our geothermometer (Eq. (10)) should be limited to data measured using
a 532 nm laser.

Page 12 of 14

50

40

30
D1-bandG-band

20
20

30

40

50

FWHM of 532 nm laser, cm–1(c) 1600
Center position of 633 nm laser, cm–1

Kaneki et al. Progress in Earth and Planetary Science

1590

1500
1580
1590

1580

1400
D1-bandG-band

1300
1300

1400

1500

1600

Center position of 532 nm laser, cm–1
Fig. 9 Effect of laser wavelength on some Raman parameters.
Comparisons of some Raman parameters of the N30 sample
between the laser wavelengths of 532 and 633 nm are shown. a R1
and R2 ratios, b FWHM of the D1- and G-bands, and c center positionsof the D1- and G-bands are shown. Error bars indicate the standard
deviation. Note that the error bars in c are hidden by the data points
owing to their small values

Kaneki et al. Progress in Earth and Planetary Science(2024) 11:35

for the quadratic calibration curve first proposed by Aoya
et al. (2010) to model the relationship between R2 and the
maximum metamorphic temperature T . This modified
curve is used in our code to calculate the mean T and its
95% prediction interval, although its direct application to
samples experiencing the processes other than the contact metamorphism may not be appropriate. Our code
is applicable to datasets other than those of Aoya et al.
(2010), as long as the Raman spectral data are obtainedusing a 532 nm laser. We document that R2 < ∼ 0.7 is anecessary condition for R2 to be sensitive to changes inthe sample maturity. For this reason, our modified geothermometer was calibrated without the K08 sample(mean R2 of 0.737) and should be preferred over thegeothermometer of Aoya et al. (2010). Since a combination of the geothermometers of Kaneki and Kouketsu(2022) and the present study is incapable of estimating
T of samples with maturity in a particular range, further
refinement or a slight extrapolation of the present geothermometer is required. In the latter case, we propose
R2 = 0.57 as a possible boundary where users shouldswitch from the geothermometer presented here to theone proposed by Kaneki and Kouketsu (2022).
Abbreviations
A2010	Aoya et al. (2010)
CM	Carbonaceous material
FWHM	Full width at half maximum
KK2022	Kaneki and Kouketsu (2022)
K2014	Kouketsu et al. (2014)SD	Standard deviation
SE	Standard error

Supplementary Information
The online version contains supplementary material available at https://​doi.​
org/​10.​1186/​s40645-​024-​00637-8.
Supplementary Material 1.
Supplementary Material 2.
Supplementary Material 3.
Supplementary Material 4.
Acknowledgements
We thank Seira Katagiri, Kazuki Matsuyama, and Shunsuke Ogino of Nagoya
University for their kind help with assessing the utility of the code. SK also
thanks Yukitoshi Fukahata of Kyoto University for helpful comments on the
Akaike Information Criterion and Hiroyuki Noda of Kyoto University for fruitful
discussion on the preliminarily results. The web page created by Toshihide
Ihara and Shin Shigemitsu and maintained by Organic Primary Standards
Group, AIST, greatly helped our regression analysis (https://​unit.​aist.​go.​jp/​mcml/​rg-​orgp/​uncer​tainty_​lectu​re/​index.​html). We thank the associate editor
Madhusoodhan Satish-Kumar for careful handling of the manuscript. Comments from two anonymous reviewers helped us improve the paper.
Author contributions
SK, YK, and MA proposed the initial concept for the project based on work
initiated by SW. YK, MA, YN, YS, and KY provided the raw Raman spectral data
published in the previous papers. YK performed additional measurements
using a 633 nm laser. SK developed, validated, and applied the code. SK wrote
an early version of the first draft, which was subsequently refined by SW. All

Page 13 of 14

authors discussed the results, contributed to the writing of the manuscript,
and approved its final version.
Funding
The present study was supported by JSPS KAKENHI Grant nos. JP20J01284
and JP24K17151 to SK, JSPS KAKENHI Grant no. JP21H05202 to SW, and JSPS
KAKENHI Grant nos. JP20J12701 and JP22H05314 to YS.
Availability of data and materials
Readers can download the code developed in the present study from https://​doi.​org/​10.​5281/​zenodo.​10237​726. In the development, validation, andapplication of the code, we reanalyzed the raw Raman spectral data originallymeasured and analyzed by Aoya et al. (2010), Kouketsu et al. (2019, 2021),
Nakamura et al. (2019), Shimura et al. (2021), and Yamaoka et al. (2022). The
raw Raman spectral data newly measured in the present study are used only
to investigate the dependence of the analytical results on laser wavelength.
Requests for access to the new data should be sent to YK. All analyticalresults obtained are listed in tables either in the main text or in Supporting
Information.

Declarations
Competing interests
The authors declare that they have no competing interests that could haveappeared to influence the work reported in this paper.
Author details
1
Geological Survey of Japan, National Institute of Advanced Industrial Science and Technology, Higashi 1‑1‑1, Tsukuba 305‑8567, Japan. 2 Graduate
School of Environmental Studies, Nagoya University, Furo‑cho, Chikusa‑ku,
Nagoya 464‑8602, Japan. 3 Graduate School of Technology, Industrial
and Social Sciences, Tokushima University, Minami‑Josanjima 2‑1, Tokushima 770‑8506, Japan. 4 Graduate School of Science, The University of Tokyo,
Tokyo 113‑0033, Japan.
Received: 21 December 2023 Accepted: 6 June 2024

References
Aoya M, Kouketsu Y, Endo S et al (2010) Extending the applicability of theRaman carbonaceous-material geothermometer using data from contact
metamorphic rocks. J Metamorph Geol 28:895–914. https://​doi.​org/​10.​1111/j.​1525-​1314.​2010.​00896.x
Beyssac O, Goffé B, Chopin C, Rouzaud JN (2002) Raman spectra of carbonaceous material in metasediments: a new geothermometer. J Metamorph
Geol 20:859–871. https://​doi.​org/​10.​1046/j.​1525-​1314.​2002.​00408.x
Bonneville S, Delpomdor F, Préat A et al (2020) Molecular identification offungi microfossils in a neoproterozoic shale rock. Sci Adv 6:eaax7599.
https://​doi.​org/​10.​1126/​sciadv.​aax75​99
Buseck PR, Huang B-J (1985) Conversion of carbonaceous material to graphite
during metamorphism. Geochim Cosmochim Acta 49:2003–2016.
https://​doi.​org/​10.​1016/​0016-​7037(85)​90059-6
Busemann H, Alexander MO, Nittler LR (2007) Characterization of insolubleorganic matter in primitive meteorites by micro Raman spectroscopy.
Meteorit Planet Sci 42:1387–1416. https://​doi.​org/​10.​1111/j.​1945-​5100.​
2007.​tb005​81.x
Cavalazzi B, Lemelle L, Simionovici A et al (2021) Cellular remains in a ~3.42-billion-year-old subseafloor hydrothermal environment. Sci Adv 7:eabf3963.
https://​doi.​org/​10.​1126/​sciadv.​abf39​63
Cody GD, Alexander CMO, Yabuta H et al (2008) Organic thermometry forchondritic parent bodies. Earth Planet Sci Lett 272:446–455. https://​doi.​
org/​10.​1016/j.​epsl.​2008.​05.​008
Ferrari AC, Robertson J (2001) Resonant Raman spectroscopy of disordered,amorphous, and diamondlike carbon. Phys Rev B 64:075414. https://​doi.​
org/​10.​1103/​PhysR​evB.​64.​075414
French BM (1964) Graphitization of organic material in a progressively metamorphosed Precambrian iron formation. Science 146:917–918. https://​doi.​org/​10.​1126/​scien​ce.​146.​3646.​917

Kaneki et al. Progress in Earth and Planetary Science(2024) 11:35

Ganz H, Kalkreuth W (1987) Application of infrared spectroscopy to theclassification of kerogentypes and the evaluation of source rock and oil
shale potentials. Fuel 66:708–711. https://​doi.​org/​10.​1016/​0016-​2361(87)​90285-7
Henry DG, Jarvis I, Gillmore G, Stephenson M (2019) Raman spectroscopy asa tool to determine the thermal maturity of organic matter: application
to sedimentary, metamorphic and structural geology. Earth Sci Rev
198:102936. https://​doi.​org/​10.​1016/j.​earsc​irev.​2019.​102936
Hickman-Lewis K, Cavalazzi B, Sorieul S et al (2020) Metallomics in deep time
and the influence of ocean chemistry on the metabolic landscapes ofEarth’s earliest ecosystems. Sci Rep 10:4965. https://​doi.​org/​10.​1038/​
s41598-​020-​61774-w
Hirono T, Maekawa Y, Yabuta H (2015) Investigation of the records of earthquake slip in carbonaceous materials from the T aiwan C helungpu fault
by means of infrared and R aman spectroscopies. Geochem Geophys
Geosyst 16:1233–1253. https://​doi.​org/​10.​1002/​2014G​C0056​22
Homma Y, Kouketsu Y, Kagi H et al (2015) Raman spectroscopic thermometryof carbonaceous material in chondrites: four–band fitting analysis andexpansion of lower temperature limit. J Mineral Petrol Sci 110:276–282.
https://​doi.​org/​10.​2465/​jmps.​15071​3a
Kaneki S, Kouketsu Y (2022) An automatic peak deconvolution method for
Raman spectra of terrestrial carbonaceous material for application to
the geothermometers of Kouketsu et al. (2014). Island Arc 31(1):e12467.
https://​doi.​org/​10.​1111/​iar.​12467
Kaneki S, Hirono T, Mukoyoshi H et al (2016) Organochemical characteristics
of carbonaceous materials as indicators of heat recorded on an ancientplate-subduction fault. Geochem Geophys Geosyst 17:2855–2868.
https://​doi.​org/​10.​1002/​2016G​C0063​68
Kouketsu Y (2023) Fusion of spectroscopy and geology: earthʼs interior
environment revealed through light. Jpn Magaz Mineral Petrol Sci
52:gkk.230110a. https://​doi.​org/​10.​2465/​gkk.​23011​0a
Kouketsu Y, Mizukami T, Mori H et al (2014) A new approach to develop theRaman carbonaceous material geothermometer for low-grade metamorphism using peak width. Island Arc 23:33–50. https://​doi.​org/​10.​1111/​iar.​
12057
Kouketsu Y, Shimizu I, Wang Y et al (2017) Raman spectra of carbonaceous
materials in a fault zone in the Longmenshan thrust belt, China; comparisons with those of sedimentary and metamorphic rocks. Tectonophysics
699:129–145. https://​doi.​org/​10.​1016/j.​tecto.​2017.​01.​015
Kouketsu Y, Tsai C-H, Enami M (2019) Discovery of unusual metamorphic
temperatures in the Yuli belt, eastern Taiwan: New interpretation of data
by Raman carbonaceous material geothermometry. Geology 47:522–526.
https://​doi.​org/​10.​1130/​G45934.1
Kouketsu Y, Sadamoto K, Umeda H et al (2021) Thermal structure in subducted
units from continental Moho depths in a palaeo subduction zone, theAsemigawa region of the Sanbagawa metamorphic belt, SW Japan. J
Metamorph Geol 39:727–749. https://​doi.​org/​10.​1111/​jmg.​12584
Lahfid A, Beyssac O, Deville E et al (2010) Evolution of the Raman spectrum ofcarbonaceous material in low-grade metasediments of the Glarus Alps
(Switzerland): RSCM in low-grade metasediments. Terra Nova 22:354–360.
https://​doi.​org/​10.​1111/j.​1365-​3121.​2010.​00956.x
Lünsdorf NK, Lünsdorf JO (2016) Evaluating Raman spectra of carbonaceous
matter by automated, iterative curve-fitting. Int J Coal Geol 160–161:51–
62. https://​doi.​org/​10.​1016/j.​coal.​2016.​04.​008
Lünsdorf NK, Dunkl I, Schmidt BC et al (2017) Towards a higher comparabilityof geothermometric data obtained by raman spectroscopy of carbonaceous material. Part 2: a revised geothermometer. Geostandard Geoanal
Res 41:593–612. https://​doi.​org/​10.​1111/​ggr.​12178
McKirdy DM, Powell TG (1974) Metamorphic alteration of carbon isotopic
composition in ancient sedimentary organic matter: new evidence from
Australia and South Africa. Geol 2:591. https://​doi.​org/​10.​1130/​0091-​
7613(1974)2%​3c591:​MAOCIC%​3e2.0.​CO;2
Mori N, Wallis S, Mori H (2015a) Graphitization of carbonaceous material in
sedimentary rocks on short geologic time-scales: An example from the K
insho-zan area, central J apan. Island Arc 24:119–130. https://​doi.​org/​10.​1111/​iar.​12093
Mori H, Wallis S, Fujimoto K, Shigematsu N (2015b) Recognition of shear
heating on a long-lived major fault using Raman carbonaceous material
thermometry: implications for strength and displacement history of the
MTL, SW Japan. Island Arc 24:425–446. https://​doi.​org/​10.​1111/​iar.​12129

Page 14 of 14

Mori H, Mori N, Wallis S et al (2017) The importance of heating duration for
Raman CM thermometry: evidence from contact metamorphism aroundthe Great Whin Sill intrusion, UK. J Metamorph Geol 35:165–180. https://​doi.​org/​10.​1111/​jmg.​12225
Mukoyoshi H, Kaneki S, Hirono T (2018) Slip parameters on major thrusts at
a convergent plate boundary: regional heterogeneity of potential slip
distance at the shallow portion of the subducting plate. Earth Planets
Space 70:36. https://​doi.​org/​10.​1186/​s40623-​018-​0810-z
Nakamura Y, Oohashi K, Toyoshima T et al (2015) Strain-induced amorphization of graphite in fault zones of the Hidaka metamorphic belt, Hokkaido,
Japan. J Struct Geol 72:142–161. https://​doi.​org/​10.​1016/j.​jsg.​2014.​10.​012
Nakamura Y, Hara H, Kagi H (2019) Natural and experimental structural evolution of dispersed organic matter in mudstones: The Shimanto accretionary complex, southwest Japan. Island Arc 28:e12318. https://​doi.​org/​10.​1111/​iar.​12318
Rahl J, Anderson K, Brandon M, Fassoulas C (2005) Raman spectroscopic carbonaceous material thermometry of low-grade metamorphic rocks: calibration and application to tectonic exhumation in Crete, Greece. Earth
Planet Sci Lett 240:339–354. https://​doi.​org/​10.​1016/j.​epsl.​2005.​09.​055
Roberts S, Tricker PM, Marshall JEA (1995) Raman spectroscopy of chitinozoans
as a maturation indicator. Org Geochem 23:223–228. https://​doi.​org/​10.​1016/​0146-​6380(94)​00126-L
Schito A, Romano C, Corrado S et al (2017) Diagenetic thermal evolutionof organic matter by Raman spectroscopy. Org Geochem 106:57–67.
https://​doi.​org/​10.​1016/j.​orgge​ochem.​2016.​12.​006
Seifert WK (1978) Steranes and terpaues in kerogen pyrolysis for correlation ofoils and source rocks. Geochim Cosmochim Acta 42:473–484. https://​doi.​
org/​10.​1016/​0016-​7037(78)​90197-7
Shimura Y, Tokiwa T, Mori H et al (2021) Deformation characteristics and peaktemperatures of the Sanbagawa Metamorphic and Shimanto Accretionary complexes on the central Kii Peninsula, SW Japan. J Asian Earth Sci
215:104791. https://​doi.​org/​10.​1016/j.​jseaes.​2021.​104791
Sweeney JJ, Burnham AK (1990) Evaluation of a simple model of vitrinitereflectance based on chemical kinetics. Am Assoc Petrol Geol Bull
74:1559–1570
Tuinstra F, Koenig JL (1970) Raman spectrum of graphite. J Chem Phys
53:1126–1130. https://​doi.​org/​10.​1063/1.​16741​08
van Krevelen DW (1993) Coal, 3rd edn. Elsevier
Wada H, Suzuki K (1983) Carbon isotopic thermometry calibrated by dolomitecalcite solvus temperatures. Geochim Cosmochim Acta 47:697–706.
https://​doi.​org/​10.​1016/​0016-​7037(83)​90104-7
Wang Y, Alsmeyer DC, McCreery RL (1990) Raman spectroscopy of carbon
materials: structural basis of observed spectra. Chem Mater 2:557–563.
https://​doi.​org/​10.​1021/​cm000​11a018
Yamaoka K, Wallis SR, Miyake A, Kouketsu Y (2022) Recognition of detrital
carbonaceous material in the Ryoke metamorphic belt by using Raman
thermometry: implications for thermal structure and detrital origin. Lithosphere 2022:3899340. https://​doi.​org/​10.​2113/​2022/​38993​40

Publisher’s NoteSpringer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.