Diagnostic and analytical tools for peak detection and fire-history interpretations using high-resolution sediment-charcoal records
Philip Higuera
Professor, Department of Ecosystem and Conservation Sciences
University of Montana, Missoula, MT, USA
philip.higuera@umontana.edu
https://www.umt.edu/people/phiguera
https://phiguera.github.io/CharAnalysis
Version 2.0 — Updated March 2026
This guide documents CharAnalysis Version 2.0. Where behavior differs from Version 1.1, differences are noted explicitly. The analytical methods are unchanged between versions.
- Part I. Background
- Part II. Using CharAnalysis
- Part III. Understanding CharAnalysis
- Part IV. Acknowledgments
- Part V. Disclaimer
- Part VI. References
CharAnalysis is a set of diagnostic and analytical tools designed for analyzing sediment-charcoal records when the goal is peak detection to reconstruct "local" fire history. The analyses were developed based on widely applied approaches that decompose a charcoal record into low- and high-frequency components (e.g. Clark and Royall 1996; Long et al. 1998; Carcaillet et al. 2001; Gavin et al. 2006), and the program introduced a technique of using a locally-defined threshold to separate signal from noise (Higuera et al. 2009). The program is set up to make explicit the range of choices an analyst must make when implementing this approach. Diagnostic tools help determine if peak detection is warranted and, if so, what parameters are most reasonable. Sensitivity analyses illustrate the impacts of alternative analysis criteria on peak-based fire-history interpretations, and graphical displays and statistical analyses summarize peak-based fire-history metrics. The papers by Higuera et al. (2010) and Kelly et al. (2011) provide the most thorough background for the analyses employed in CharAnalysis, including the assumptions inherent in charcoal peak identification using a local threshold.
CharAnalysis is freely available at https://github.com/phiguera/CharAnalysis. Since its original development in the mid-2000s, the program has been used in dozens of published studies to analyze sediment-charcoal records worldwide. The entire codebase is distributed and well commented — users are encouraged to look under the hood, understand what is going on, and modify the program to suit individual needs.
CharAnalysis is available in both MATLAB (Version 2.0; the reference implementation) and R (Version 2.0; a direct translation). This guide is written primarily for the MATLAB version. R users should consult the R package vignette for R-specific installation, usage, and a description of known numerical differences between the two implementations.
The following two papers provide the most thorough background for the analyses employed in CharAnalysis, including the assumptions inherent in charcoal peak identification using a local threshold:
Higuera, P.E., D.G. Gavin, P.J. Bartlein, and D.J. Hallett. 2010. Peak detection in sediment-charcoal records: impacts of alternative data analysis methods on fire-history interpretations. International Journal of Wildland Fire 19:996–1014.
Kelly, R.F., P.E. Higuera, C.M. Barrett, and F.S. Hu. 2011. A signal-to-noise index to quantify the potential for peak detection in sediment-charcoal records. Quaternary Research 75:11–17.
Selected publications using CharAnalysis (incomplete list, through 2014)
Dunnette, P.V., P.E. Higuera, K.K. McLauchlan, K.M. Derr, C.E. Briles, and M.H. Keefe. 2014. Biogeochemical impacts of wildfires over four millennia in a Rocky Mountain subalpine watershed. New Phytologist 203:900–912.
Brossier, B., F. Oris, W. Finsinger, H. Asselin, Y. Bergeron, and A.A. Ali. 2014. Using tree-ring records to calibrate peak detection in fire reconstructions based on sedimentary charcoal records. The Holocene 24:635–645.
Kelly, R., M.L. Chipman, P.E. Higuera, I. Stefanova, L.B. Brubaker, and F.S. Hu. 2013. Recent burning of boreal forests exceeds fire regime limits of the past 10,000 years. Proceedings of the National Academy of Sciences 110:13055–13060.
Mustaphi, C.J.C., and M.F.J. Pisaric. 2013. Varying influence of climate and aspect as controls of montane forest fire regimes during the late Holocene, south-eastern British Columbia, Canada. Journal of Biogeography 40:1983–1996.
Blarquez, O., M.P. Girardin, B. Leys, A.A. Ali, J.C. Aleman, Y. Bergeron, and C. Carcaillet. 2013. Paleofire reconstruction based on an ensemble-member strategy applied to sedimentary charcoal. Geophysical Research Letters 40:2667–2672.
Barrett, C.M., R.F. Kelly, P.E. Higuera, and F.S. Hu. 2013. Climatic and land cover influences on the spatiotemporal dynamics of Holocene boreal fire regimes. Ecology 92:389–402.
Ali, A.A., O. Blarquez, M.P. Girardin, C. Hely, F. Tinquaut, A. El Guellab, V. Valsecchi, A. Terrier, L. Bremond, A. Genries, S. Gauthier, and Y. Bergeron. 2012. Control of the multimillennial wildfire size in boreal North America by spring climatic conditions. Proceedings of the National Academy of Sciences 109:20966–20970.
Marlon, J.R., P.J. Bartlein, D.G. Gavin, C.J. Long, R.S. Anderson, C.E. Briles, K.J. Brown, D. Colombaroli, D.J. Hallett, M.J. Power, E.A. Scharf, and M.K. Walsh. 2012. Long-term perspective on wildfires in the western USA. Proceedings of the National Academy of Sciences 109:3203–3208.
van Bellen, S., M. Garneau, A.A. Ali, and Y. Bergeron. 2012. Did fires drive Holocene carbon sequestration in boreal ombrotrophic peatlands of eastern Canada? Quaternary Research 78:50–59.
Higuera, P.E., M.L. Chipman, J.L. Barnes, M.A. Urban, and F.S. Hu. 2011. Variability of tundra fire regimes in Arctic Alaska: millennial-scale patterns and ecological implications. Ecological Applications 21:3211–3226.
Briles, C.E., C. Whitlock, C.N. Skinner, and J. Mohr. 2011. Postglacial forest development on different substrates in the Klamath Mountains, northern California, USA. Ecology 92:590–601.
Whitlock, C., C.E. Briles, M.C. Fernandez, and J. Gage. 2011. Holocene vegetation, fire, and climate history of the Sawtooth Range, central Idaho, U.S.A. Quaternary Research 75:114–124.
Higuera, P.E., C. Whitlock, and J. Gage. 2011. Linking tree-ring and sediment-charcoal records to reconstruct fire occurrence and area burned in subalpine forests of Yellowstone National Park, U.S.A. The Holocene 21:327–341.
Hu, F.S., P.E. Higuera, J.E. Walsh, W.L. Chapman, P.A. Duffy, L.B. Brubaker, and M.L. Chipman. 2010. Tundra burning in Alaska: linkages to climatic change and sea-ice retreat. Journal of Geophysical Research — Biogeosciences 115, G04002.
Walsh, M.K., C. Whitlock, and P.J. Bartlein. 2010. 1200 years of fire and vegetation history in the Willamette Valley, Oregon and Washington, reconstructed using high-resolution macroscopic charcoal and pollen analysis. Palaeogeography Palaeoclimatology Palaeoecology 297:273–289.
Chepstow-Lusty, A.J., M.R. Frogley, B.S. Bauer, M.J. Leng, K.P. Boessenkool, C. Carcaillet, A.A. Ali, and A. Gioda. 2009. Putting the rise of the Inca Empire within a climatic and land management context. Climate of the Past 5:375–388.
Marlon, J.R., P.J. Bartlein, M.K. Walsh, S.P. Harrison, K.J. Brown, M.E. Edwards, P.E. Higuera, M.J. Power, R.S. Anderson, C. Briles, A. Brunelle, C. Carcaillet, M. Daniels, F.S. Hu, M. Lavoie, C. Long, T. Minckley, P.J.H. Richard, S.L. Shafer, W. Tinner, C. Umbanhowar, and C. Whitlock. 2009. Wildfire responses to abrupt climate change in North America. PNAS 106:2519–2524.
Walsh, M.K., C. Whitlock, and P.J. Bartlein. 2008. A 14,300-year-long record of fire-vegetation-climate linkages at Battle Ground Lake, southwestern Washington. Quaternary Research 70:251–264.
Huerta, M.A., C. Whitlock, and J. Yale. 2009. Holocene vegetation-fire-climate linkages in northern Yellowstone National Park, USA. Palaeogeography Palaeoclimatology Palaeoecology 271:170–181.
Briles, C.E., C. Whitlock, P.J. Bartlein, and P.E. Higuera. 2008. Regional and local controls on postglacial vegetation and fire in the Siskiyou Mountains, northern California, USA. Palaeogeography Palaeoclimatology Palaeoecology 265:159–169.
Higuera, P.E., L.B. Brubaker, P.M. Anderson, T.A. Brown, A.T. Kennedy, and F.S. Hu. 2008. Frequent fires in ancient shrub tundra: implications of paleorecords for Arctic environmental change. PLoS ONE 3:e0001744.
If you use CharAnalysis in a publication, please cite Higuera et al. (2009), the first study to apply the core analytical tools implemented in CharAnalysis:
Higuera, P.E., L.B. Brubaker, P.M. Anderson, F.S. Hu, and T.A. Brown. 2009. Vegetation mediated the impacts of postglacial climate change on fire regimes in the south-central Brooks Range, Alaska. Ecological Monographs 79:201–219. https://esajournals.onlinelibrary.wiley.com/doi/full/10.1890/07-2019.1
If you used Version 2.0 specifically, please also cite the software:
Higuera, P.E. 2026. CharAnalysis: Diagnostic and analytical tools for peak analysis in sediment-charcoal records (Version 2.0). Zenodo. https://doi.org/10.5281/zenodo.19304064
Updates are documented on the GitHub repository. If you encounter problems, please use the Issues tab at https://github.com/phiguera/CharAnalysis/issues. A GitHub account is required to register a new issue. Before posting, search existing issues to see whether your problem has already been addressed.
Note: Prior to April 2014, CharAnalysis was hosted on Google Code (http://code.google.com/p/charanalysis/). Resolved issues from that period are archived there.
Version 2.0 (March 2026) is the first major update to the codebase. The analytical methods are unchanged from Version 1.1. Key changes include:
- No toolbox dependencies — all calls to the Curve Fitting Toolbox
smooth()function have been replaced with a base-MATLAB implementation (charLowess.m), so CharAnalysis runs on any MATLAB R2019a+ installation without additional licenses. - Deprecated functions replaced (
plotyy,xlsread,xlswrite,hist, and others). - Global variables eliminated, improving reliability when functions are called in non-standard order.
- Computation separated from visualization, making batch processing more straightforward.
- Input validation added, providing clear error messages for common misconfigurations.
- Modular figure interface added — each output figure (Figures 3–9) is implemented as a standalone function callable independently from the MATLAB workspace. A new
'modular'run mode provides interactive and programmatic figure selection. See Section 3.1 for details.
There are three ways to access CharAnalysis, suited to different users and needs.
Install the R package from CRAN. Requires R 4.0 or higher.
install.packages("CharAnalysis")For in-development features ahead of the next CRAN release, install from the
dev branch on GitHub:
# install.packages("devtools")
devtools::install_github("phiguera/CharAnalysis",
subdir = "CharAnalysis_2_0_R",
ref = "dev")For a full worked example, parameter descriptions, output descriptions, and known numerical differences from the MATLAB implementation, consult the package vignette:
vignette("CharAnalysis_intro", package = "CharAnalysis")This is the recommended option for MATLAB users conducting analyses on their own data.
Requirements
- MATLAB R2019a or higher. No additional toolboxes are required.
V1.1 difference: Version 1.1 requires MATLAB 7.0 or higher and the Curve Fitting Toolbox. Version 2.0 has no toolbox dependencies.
Download
Download the repository as a .zip or tar.gz archive, or clone it from the command line:
git clone https://github.com/phiguera/CharAnalysis.git
The archive contains the full repository, including both the MATLAB and R implementations. After unpacking, the only folders needed for MATLAB use are:
CharAnalysis_2_0_MATLAB/— Version 2.0 source code, including the version-comparison scriptz_Compare_CharAnalysis_V1_V2.m.DataTemplates_and_Examples/— template input files and the bundled Code Lake example dataset.
To browse the MATLAB folder online without downloading, see CharAnalysis_2_0_MATLAB on GitHub.
Installation
- Add the
CharAnalysis_2_0_MATLABfolder to your MATLAB search path:addpath '.../CharAnalysis_2_0_MATLAB'
- Save the path:
savepath
Not sure if CharAnalysis is the right tool for your research? You can run the program instantly in your web browser on the bundled Code Lake example dataset — no installation required. Click the badge below or visit the GitHub repository page.
Clicking this link clones the repository to your MATLAB Drive and opens CharAnalysis ready to run. A free MathWorks account is required. University users can log in with their institutional email address for full access.
Note: MATLAB Online is best suited for evaluating CharAnalysis on the bundled example dataset. For analysis of your own data, install locally via Option 1 (R) or Option 2 (MATLAB); managing input files and parameter files is more straightforward in a desktop environment. When running online, set
saveData = 0andsaveFigures = 0in the parameter file to avoid file-write errors.
Open the included template files template_charData.csv and
template_charParams.csv and save them under new names identifying your site
(e.g. CO_charData.csv and CO_charParams.csv). Place these files in the MATLAB
working directory, or specify the full path when calling CharAnalysis.
V1.1 difference: The CSV format is recommended for Version 2.0 and is the only template format distributed with the program. Input files in
.xlsformat remain supported but are no longer documented.
Paste charcoal data into the charData worksheet (Figure 1). The format is the same as for the program Charster (Gavin). For each sample, enter:
| Column | Variable | Description |
|---|---|---|
| 1 | cmTop | Depth at top of sample (cm) |
| 2 | cmBot | Depth at base of sample (cm) |
| 3 | ageTop | Age at top of sample (cal. yr BP) |
| 4 | ageBot | Age at bottom of sample (cal. yr BP) |
| 5 | charVol | Volume of sediment sample (cm³) |
| 6 | charCount | Charcoal count (pieces) |
Enter the site name in cell G1 (column 7, row 1).
Missing values for charVol and charCount (but not for depths or ages) can be indicated by any value < 0 (e.g. −999). CharAnalysis will interpolate across missing data. Note: paste only values into the spreadsheet — formulas are not supported.
Figure 1. The charData worksheet in the template file, where data input occurs.
Parameter choices are entered in column C (column 3) of the CharParams worksheet (.xlsx) or file (.csv) (Figure 2). Parameters are divided into four stages.
Figure 2. The CharParams worksheet in the template file, where analysis parameters are selected.
| Parameter | Description |
|---|---|
zoneDiv |
Years defining the beginning and end of the record, and any zone divisions. Used for plotting and for analyzing fire return intervals. Enter at least two values (beginning and end of record) in ascending order (youngest age first). Additional zone boundaries may be added. |
yrInterp |
Interpolation interval (years). All records are resampled to equal time intervals before analysis. Enter 0 to interpolate to the median sample resolution of the raw record. |
transform |
Data transformation before analysis: 0 = none; 1 = base-10 log; 2 = natural log. For options 1–2, a value of 1 is added to all CHAR values before transformation. |
| Parameter | Description |
|---|---|
method |
Smoothing method for estimating Cbackground: 1 = Lowess; 2 = Robust Lowess (resistant to outliers); 3 = Moving average; 4 = Moving median; 5 = Moving mode. |
yr |
Smoothing window width (years) for estimating Cbackground. |
| Parameter | Description |
|---|---|
cPeak |
Method for calculating Cpeak: 1 = residuals (Cpeak = Cint − Cback); 2 = ratios (Cpeak = Cint / Cback). |
threshType |
Threshold type: 1 = globally defined (based on entire record); 2 = locally defined (sliding window around each sample). |
threshMethod |
Method for determining threshold values: 1 = user-defined values in threshValues; 2 = percentile of a Gaussian noise distribution; 3 = same as 2, but noise distribution determined by a Gaussian mixture model. |
threshValues |
Threshold values to evaluate (up to 4 values). If threshMethod = 1, values are in Cpeak units. If threshMethod = 2–3, values are percentiles of the noise distribution (e.g. 0.95). The last value is used for peak plotting and analysis. |
minCountP |
Cut-off probability for minimum-count analysis (e.g. 0.05). Potential peaks where the preceding minimum count has a probability > minCountP of coming from the same Poisson distribution as the peak count are flagged but excluded from analysis. Set to 0.99 to turn off. |
peakFrequ |
Smoothing window (years) for fire frequency and fire return interval curves. |
sensitivity |
Set to 1 to evaluate sensitivity of results to varying Cbackground window widths; 0 to skip. |
| Parameter | Description |
|---|---|
saveFigures |
Save output figures as .tiff and .pdf files? 0 = no; 1 = yes. Set to 0 when using MATLAB Online. |
saveData |
Append output data to the input file? 0 = no; 1 = yes. Set to 0 when using MATLAB Online. |
allFigures |
Display all diagnostic figures (including Figures 1, 2, and 9)? 0 = no (show only peak analysis results); 1 = yes. |
Open MATLAB and type CharAnalysis into the Command Window. When prompted, enter the name of your parameter file with single quotation marks and the file extension:
>> >> CharAnalysis
>>
>> *****************************************************************
>> CharAnalysis 2.0
>> (c) 2004 - 2026, P.E. Higuera
>> MATLAB(r). (c) 1984 - 2024 The MathWorks, Inc.
>> philip.higuera@umontana.edu
>> Please read documentation:
>> https://phiguera.github.io/CharAnalysis/
>> *****************************************************************
>>
>> CharAnalysis requires an input file in .xls or .csv format.
>> This file includes the selected parameters for peak analysis,
>> and it either includes (.xls) or references (.csv) the input
>> charcoal dataset. The input file must be in the working
>> directory for the program to run.
>>
>> If you choose to save the results (figures and/or data), they
>> will be saved in the directory containing the input file.
>> *NOTE*: you must close the input file before running
>> CharAnalysis for the .xls or .csv file to be updated with new
>> results.
>>
>> *****************************************************************
>>
>> Input the file name OR the full path to the site directory,
>> bounded with single quotations and including the file
>> extension (if file name): 'mysite_charParams.csv'You can also call CharAnalysis directly with the filename as an argument, which is useful for processing multiple records:
>> CharAnalysis('site1_charParams.csv')
>> CharAnalysis('site2_charParams.csv')To return results to the MATLAB workspace:
>> results = CharAnalysis('mysite_charParams.csv')The 'modular' run mode allows individual output figures to be selected interactively or programmatically, without changing the parameter file or the standard workflow. All analytical results and data saving are identical to the standard run — only figure generation differs.
Interactive menu — after the analysis completes, a command-window menu prompts the user to select which figures to generate:
>> CharAnalysis('mysite_charParams.csv', 'modular')Programmatic selection — pass a vector of figure numbers as a third argument to bypass the interactive menu. This is useful for batch processing or scripted workflows:
>> CharAnalysis('mysite_charParams.csv', 'modular', [3 7])Programmatic selection with automatic save — pass true as a fourth argument to save figures without a prompt:
>> CharAnalysis('mysite_charParams.csv', 'modular', [3 7], true)Calling individual figures directly — after any run, individual figure functions can be called from the workspace using the results struct:
>> results = CharAnalysis('mysite_charParams.csv');
>> CharPlotFig7_ContinuousFireHistory(results)
>> CharPlotFig3_CintCbackCpeak(results)The figures available through the modular interface are:
| Figure | Function | Contents |
|---|---|---|
| 3 | CharPlotFig3_CintCbackCpeak |
Cint, Cback, and Cpeak |
| 4 | CharPlotFig4_ThresholdSNI |
Sensitivity to alternative thresholds and SNI |
| 5 | CharPlotFig5_CumulativePeaks |
Cumulative peaks through time |
| 6 | CharPlotFig6_FRIDistributions |
FRI distributions by zone |
| 7 | CharPlotFig7_ContinuousFireHistory |
Continuous fire history |
| 8 | CharPlotFig8_ZoneComparisons |
Between-zone comparisons |
| 9 | CharPlotFig9_ThresholdDetails |
Alternative threshold displays |
Figures 1, 2, and 10 are not part of the modular interface. Figures 1 and 2 are generated during pretreatment and threshold determination and are controlled by the allFigures parameter. Figure 10 is produced only when sensitivity = 1.
| Error | Likely Cause |
|---|---|
| Program quits while interpolating or smoothing | Input file contains formulas rather than plain values; multiple depths have the same age; or smoothing window is longer than the record. |
Unable to find file error when saving results |
The output file does not exist in the current directory. Set saveData = 0 and saveFigures = 0 to skip saving (required for MATLAB Online). |
| Unexpected peak counts vs. Version 1.1 | Verify parameter settings match between versions. Use z_Compare_CharAnalysis_V1_V2.m to compare outputs numerically. |
| Figures appear very large in MATLAB Online | This is a display scaling issue. Drag figure window corners to resize, or use browser zoom (Ctrl/Cmd −) to reduce the interface size. |
Unrecognized field name error when calling a figure function directly |
The results struct was returned by an older version of CharAnalysis.m. Re-run CharAnalysis to generate a fresh results struct that includes the Post, gapIn, and site fields required by the figure functions. |
| Term | Description |
|---|---|
| C | Charcoal accumulation rate (CHAR; pieces cm⁻² yr⁻¹) |
| Craw | CHAR of the raw record |
| Cint | CHAR of the interpolated record |
| Cback | Low-frequency trend in Cint, also termed "background CHAR" or "BCHAR" in the literature |
| Cpeak | High-frequency component of Cint, after Cback is removed |
| Cnoise | One of the two additive components of Cpeak |
| Cfire | The other additive component of Cpeak |
| t | Sample-specific threshold value used to separate Cnoise from Cfire |
| Cthresh | Time series of t, displayed in various ways in output figures |
| SNI | Signal-to-noise index: var(Cfire) / (var(Cfire) + var(Cnoise)) |
The structure of CharAnalysis reflects the main analytical components of most decomposition methods. At each step the user makes one or more parameter decisions (Figure 3). The general steps are:
- Interpolate the raw charcoal series (concentration [pieces cm⁻³], sediment accumulation rate [cm yr⁻¹]) to equal time intervals to define Cinterpolated (pieces cm⁻² yr⁻¹).
- Smooth Cint to model low-frequency trends and define Cback.
- Remove Cback from Cint to create Cpeak, containing only high-frequency variations.
- Define threshold t and apply to Cpeak to separate fire-related samples from non-fire-related samples.
- Screen peaks and remove any that fail to pass the minimum-count criterion.
Figure 3. Decision tree for peak detection in CharAnalysis.
Stratifying Records. Records can be divided into user-defined zones identified by ages. Zone divisions are used for plotting and for statistical comparisons of fire return intervals.
Interpolating. All records must be resampled to equal time intervals to justify the techniques used for peak analysis. Resampling in CharAnalysis determines the relative proportion that each raw sample contributes to each interpolated interval, then weights the raw sample(s) within an interpolated interval based on these proportions. This technique preserves the primary structure of the charcoal data better than binning pseudo-annual values derived by linear interpolation, because it makes fewer assumptions about the pattern of charcoal deposition within sampling intervals. The interpolation interval can be user-defined or set to the median sampling interval of all raw samples.
Data Transformation. Before peak analysis, CHAR data can be transformed via log (base 10) or natural log transformations. In each case, a value of 1 is added to all CHAR values before transformation.
Smoothing refers to the method used to model low-frequency trends in a charcoal record, Cbackground. The analyst must select both a smoothing method and a timescale (window width) that define Cback. Tools to help make this choice are presented in Figure 10 (Cbackground sensitivity analysis). CharAnalysis includes five smoothing methods:
- Lowess — locally weighted scatter plot smooth using least squares linear polynomial fitting.
- Robust Lowess — Lowess smoothing resistant to outliers.
- Moving average filter.
- Moving median — each sample is assigned the median value of Cint within the smoothing window. The resulting series is then smoothed with a Lowess filter.
- Moving mode — each sample is assigned the modal value from Cint within the smoothing window, with values divided into 100 equally-spaced bins. The resulting series is then smoothed with a Lowess filter.
V1.1 difference: Version 1.1 used the
smooth()function from the MATLAB Curve Fitting Toolbox for methods 1–3. Version 2.0 replaces this with a base-MATLAB implementation (charLowess.m), so no toolbox is required.
CharAnalysis includes several methods for peak analysis, reflecting both approaches used in the literature and techniques introduced with the program. The approach involves three steps.
Remove low-frequency trends in Cint to obtain a peak CHAR series, Cpeak. Two options:
- Residuals: Cpeak = Cint − Cback. Assumes an additive relationship between Cback and charcoal introduced by local fires.
- Ratios: Cpeak = Cint / Cback. Assumes a multiplicative relationship.
Determine and apply a threshold value t to each Cpeak sample, flagging the sample as a "peak" if Cpeak > t. Three options:
Option 1 — Define t manually. Threshold values are entered directly in Cpeak units (residual or ratio, depending on the cPeak setting).
Option 2 — Define t based on the Cnoise distribution. The user selects a percentile cut-off of the noise distribution (e.g. 0.95 places the threshold at the 95th percentile of Cnoise). Two approaches for estimating Cnoise:
- Gaussian assumption (
threshMethod = 2): assumes the mean of Cnoise is 0 (residuals) or 1 (ratios) and estimates variance from the lower tail of the Cpeak distribution. - Gaussian mixture model (
threshMethod = 3): uses a two-component Gaussian mixture model to identify the mean and variance of Cnoise more accurately when the noise distribution mean is not exactly 0 or 1.
In both cases, t can be defined globally (based on the entire Cpeak distribution) or locally (based on a sliding window of the same width as the Cback smoothing window). The locally-defined threshold is described in Higuera et al. (2009, 2010). When using a locally-defined threshold, ensure the window width includes at least 30 samples.
Each sample exceeding the threshold is subjected to minimum-count screening before being classified as a charcoal peak. This test calculates the probability that two charcoal counts could arise from the same Poisson distribution (Shiue and Bain 1982). The minimum charcoal count in the 75 years preceding the potential peak is compared to the maximum count within 75 years following it. The test statistic d is:
where X1, V1 are the count and volume of the minimum-count sample, and X2, V2 are those of the maximum-count sample. Peaks with a probability greater than minCountP are flagged but excluded from peak analysis.
Output data are saved to a CharResults worksheet (.xlsx) or file (.csv) when saveData = 1 (Figure 4). The columns are as follows:
| Column | Variable | Description |
|---|---|---|
| A | cmTop_i | Interpolated top depth (cm) |
| B | ageTop_i | Interpolated top age (cal. yr BP) |
| C | charCount_i | Interpolated charcoal counts (pieces) |
| D | charVol_i | Interpolated sample volume (cm³) |
| E | charCon_i | Interpolated charcoal concentration (pieces cm⁻³) |
| F | charAcc_i | Interpolated CHAR (pieces cm⁻² yr⁻¹) |
| G | charBkg | Cbackground (pieces cm⁻² yr⁻¹) |
| H | charPeak | Cpeak (pieces cm⁻² yr⁻¹) |
| I | thresh1 | 1st threshold value |
| J | thresh2 | 2nd threshold value |
| K | thresh3 | 3rd threshold value |
| L | threshFinalPos | Final positive threshold value |
| M | threshFinalNeg | Final negative threshold value |
| N | SNI | Signal-to-noise index |
| O | threshGOF | KS goodness-of-fit p-value for the noise distribution |
| P | peaks1 | Binary: start of identified peak using thresh1 |
| Q | peaks2 | Same as peaks1, but for thresh2 |
| R | peaks3 | Same as peaks1, but for thresh3 |
| S | peaksFinal | Same as peaks1, but for threshFinalPos |
| T | peaksInsig | Peaks flagged by the minimum-count criterion |
| U | peakMag | Peak magnitude (pieces cm⁻² peak⁻¹) |
| V | smPeakFrequ | Smoothed fire frequency (fires per 1000 yr) |
| W | smFRIs | Smoothed fire return intervals (yr fire⁻¹) |
| X | nFRIs | Total number of FRIs per zone |
| Y | mFRI | Mean FRI per zone (yr) |
| Z | mFRI_uCI | Upper 95% CI for mFRI (1000 bootstrapped samples) |
| AA | mFRI_lCI | Lower 95% CI for mFRI |
| BB | WBLb | Weibull scale parameter b (yr) |
| CC | WBLb_uCI | Upper 95% CI for WBLb |
| DD | WBLb_lCI | Lower 95% CI for WBLb |
| EE | WBLc | Weibull shape parameter c |
| FF | WBLc_uCI | Upper 95% CI for WBLc |
| GG | WBLc_lCI | Lower 95% CI for WBLc |
Figure 4. Example of the CharResults worksheet after running CharAnalysis and saving data.
Output figures provide a detailed look at what the program is doing numerically and display publication-quality summaries of charcoal series and peak-analysis results. All time-series plots (except Figure 2) are scaled horizontally based on the amount of time analyzed. Up to 10 figures are produced depending on parameter choices.
Figure 1: Craw and Cinterpolated; smoothing options
Panel (a): raw CHAR displayed as bars with Cinterpolated as a stair-step plot. Panel (b): Cinterpolated with all five smoothing options for a selected window width. Areas with missing values are indicated by grey boxes.
Figure 2: Distribution of Cpeak values
Varies depending on threshold type. For a global threshold, shows the full Cpeak distribution as a histogram with the modeled noise distribution and threshold value. For a local threshold, shows multiple non-overlapping time periods, each with the modeled noise distribution, local threshold values (tyr), SNI, and KS goodness-of-fit (KS p-val).
Figure 3: Cint, Cback, and Cpeak
Panel (a): Cint with Cback overlaid. Panel (b): Cpeak with positive and negative threshold values defining Cnoise, and identified peaks marked as + symbols. Peaks failing the minimum-count criterion are shown as grey dots.
Figure 4: Sensitivity to alternative thresholds and signal quality
Panel (a): CHAR series with peaks from all three threshold values. Panel (b): mean FRI and 95% confidence limits by zone for each threshold. Panel (c): SNI time series. Panel (d): boxplot of all SNI values. Note: y-axes in panels (b) and (c) are log scales.
Figure 5: Cumulative peaks through time
Cumulative sum of identified peaks as a function of time. The slope at any point is the instantaneous fire frequency (fires yr⁻¹). Areas with missing values are indicated by grey boxes.
Figure 6: Fire return interval distributions by zone
Histogram of FRIs within each zone (20-yr bins). If the fitted Weibull model passes the goodness-of-fit test (p > 0.10 if n < 30; p > 0.05 if n ≥ 30), model parameters and 95% confidence estimates are listed along with mFRI, confidence estimates, and n.
Figure 7: Continuous fire history
Top panel: peak magnitude as bars with identified peaks as + symbols. Middle panel: fire return intervals and smoothed FRI curve. Bottom panel: smoothed fire frequency. In all panels, areas with missing values are indicated by grey boxes.
Figure 8: Between-zone comparisons of raw CHAR distributions
Left panel: cumulative distribution functions (CDFs) of raw CHAR values within each zone. Two-sample KS tests compare zones pairwise; a table of p-values is displayed within the plot. Right panel: box plots of raw CHAR values by zone (10th, 25th, 50th, 75th, and 90th percentiles).
Figure 9: Alternative displays of threshold values
Panel (a): Cint with Cback and Cthresh. Panel (b): Cpeak and Cthresh displayed in the ratio domain. Panel (c): Cpeak and Cthresh displayed in the residual domain. This figure illustrates how the selected threshold would appear under both Cpeak definitions.
Figure 10: Sensitivity to Cbackground window width (produced only when sensitivity = 1)
Results from multiple analyses using varying smoothing window widths. For a local threshold: (a) KS goodness-of-fit p-values by window width; (b) SNI distributions by window width; (c) sum of (a) and (b), useful for selecting the optimal window when the two measures show opposing trends. For a global threshold: a three-variable plot of peak count (z) as a function of threshold value (x) and smoothing window (y).
Many features in CharAnalysis are based on analytical techniques from the programs CHAPS (Patrick Bartlein, University of Oregon) and Charster (Daniel Gavin, University of Oregon). The resampling algorithm and minimum-count screening were developed directly from features in Charster. Peak magnitude and smoothed fire frequency displays were developed based on CHAPS. CharAnalysis, like Charster, uses a Gaussian mixture model originally created by Charles Bouman (Purdue University).
Development of the program benefited greatly from discussions with and testing by members of the Whitlock Paleoecology Lab at Montana State University, Dan Gavin, Patrick Bartlein, and Ryan Kelly.
Version 2.0 was developed with the assistance of Claude, an AI assistant by Anthropic. Claude assisted with code modernization, bug fixes, architecture redesign, and documentation. All code was reviewed and validated by the author against Version 1.1 reference outputs.
CharAnalysis was written in MATLAB with resources from the University of Washington, Montana State University, the University of Illinois, the University of Idaho, and the University of Montana. Version 2.0 targets MATLAB R2019a or higher.
THIS SOFTWARE PROGRAM AND DOCUMENTATION ARE PROVIDED "AS IS" AND WITHOUT WARRANTIES AS TO PERFORMANCE. THE PROGRAM CharAnalysis IS PROVIDED WITHOUT ANY EXPRESSED OR IMPLIED WARRANTIES WHATSOEVER. BECAUSE OF THE DIVERSITY OF CONDITIONS AND HARDWARE UNDER WHICH THE PROGRAM MAY BE USED, NO WARRANTY OF FITNESS FOR A PARTICULAR PURPOSE IS OFFERED. THE USER IS ADVISED TO TEST THE PROGRAM THOROUGHLY BEFORE RELYING ON IT. THE USER MUST ASSUME THE ENTIRE RISK AND RESPONSIBILITY OF USING THIS PROGRAM.
THE USE OF THE SOFTWARE DOWNLOADED FROM THE UNIVERSITY OF MONTANA WEBSITE OR GITHUB IS DONE AT YOUR OWN RISK AND WITH AGREEMENT THAT YOU ARE SOLELY RESPONSIBLE FOR ANY DAMAGE TO YOUR COMPUTER SYSTEM OR LOSS OF DATA THAT RESULTS FROM SUCH ACTIVITIES.
Carcaillet, C., Y. Bergeron, P. Richard, B. Frechette, S. Gauthier, and Y. Prairie. 2001. Change of fire frequency in the eastern Canadian boreal forests during the Holocene: does vegetation composition or climate trigger the fire regime? Journal of Ecology 89:930–946.
Clark, J. S., and P. D. Royall. 1996. Local and regional sediment charcoal evidence for fire regimes in presettlement north-eastern North America. Journal of Ecology 84:365–382.
Gavin, D. G., F. S. Hu, K. Lertzman, and P. Corbett. 2006. Weak climatic control of stand-scale fire history during the late Holocene. Ecology 87:1722–1732.
Higuera, P. E., L. B. Brubaker, P. M. Anderson, F. S. Hu, and T. A. Brown. 2009. Vegetation mediated the impacts of postglacial climate change on fire regimes in the south-central Brooks Range, Alaska. Ecological Monographs 79(2):201–219.
Higuera, P. E., D. G. Gavin, P. J. Bartlein and D. J. Hallett. 2010. Peak detection in sediment-charcoal records: impacts of alternative data analysis methods on fire-history interpretations. International Journal of Wildland Fire 19: 996-1014.
Kelly, R.F., P.E. Higuera, C.M. Barrett, and F.S. Hu. 2011. A signal-to-noise index to quantify the potential for peak detection in sediment-charcoal records. Quaternary Research 75: 11-17.
Long, C. J., C. Whitlock, P. J. Bartlein, and S. H. Millspaugh. 1998. A 9000 year fire history from the Oregon Coast Range based on a high-resolution charcoal study. Canadian Journal of Forest Research 28:774–787.
Shiue, W., and L. Bain. 1982. Experiment size and power comparisons for two-sample Poisson tests. Applied Statistics 31:130–134.