Project Objective

This module introduces the application of stable water isotopes (δ¹⁸O and δD) to quantify the origin of precipitation and hydrological inputs under varying climatic conditions. Using Bayesian mixing models implemented in MixSIAR (R package), participants will learn to:

Prepare and structure isotopic datasets from precipitation, potential water sources, and discrimination factors.

Apply a Bayesian mixing model to estimate the relative contributions of different hydrological sources (rainfall, snowmelt, groundwater).

Evaluate model outputs through statistical summaries and convergence diagnostics (e.g., MCMC trace plots, Gelman–Rubin statistics) to ensure robustness.

Visualize estimated source proportions and their posterior distributions, enabling clear interpretation in both climatic and hydrological studies.

By the end of the module, participants will be able to trace water fluxes within a watershed, assess the climatic influence on precipitation sources, and communicate their findings effectively using results derived from Bayesian analysis.

In this training context, simulated datasets are employed to provide order-of-magnitude examples and hands-on practice with MixSIAR.

δ¹⁸O (delta-O-18): The ratio of oxygen-18 (¹⁸O) to oxygen-16 (¹⁶O) in water. δD (δ²H, delta-Deuterium): The ratio of hydrogen-2 (²H, deuterium) to hydrogen-1 (¹H) in water.

Step 0: Install and Load Necessary Packages

The following packages are required for the analysis. Install them using the provided instructions or install.packages('package_name'). JAGS must be installed first: JAGS download.

library(rjags)
library(MixSIAR)
require(loo)
require(vegan)
find.package("MixSIAR")
## [1] "C:/R-4.4.1/library/MixSIAR"

Step 1: Generate Synthetic Data

Synthetic data are generated for consumers (observations), sources (rain, snowmelt, groundwater), and isotopic discrimination.

# Consumer data (observations)
mix_data <- data.frame(
  SampleID = paste0("S", 1:4),
  d18O = c(-5.1, -6.2, -4.8, -9.3),
  dD   = c(-40.2, -55.1, -38.7, -65.0)
)
write.csv(mix_data, "precip_mix_data.csv", row.names = FALSE)

# Source data (Rain, Snowmelt, Groundwater)
source_data <- data.frame(
  source = c("Rain", "Snowmelt", "Groundwater"),
  Meand18O = c(-3.0, -10.0, -7.0),
  SDd18O   = c(0.5, 0.8, 0.6),
  MeandD   = c(-20.0, -70.0, -50.0),
  SDdD     = c(2.0, 3.0, 2.5),
  n        = c(5, 5, 5)  # Sample size
)
write.csv(source_data, "precip_source_data.csv", row.names = FALSE)

# Isotopic discrimination (zero if unknown)
discr_data <- data.frame(
  source = c("Rain", "Snowmelt", "Groundwater"),
  Meand18O = c(0, 0, 0),
  SDd18O   = c(0, 0, 0),
  MeandD   = c(0, 0, 0),
  SDdD     = c(0, 0, 0),
  n        = c(5, 5, 5)
)
write.csv(discr_data, "precip_discr_data.csv", row.names = FALSE)

# Display dataset dimensions and preview
dim(mix_data)
## [1] 4 3
head(mix_data)
##   SampleID d18O    dD
## 1       S1 -5.1 -40.2
## 2       S2 -6.2 -55.1
## 3       S3 -4.8 -38.7
## 4       S4 -9.3 -65.0
dim(source_data)
## [1] 3 6
head(source_data)
##        source Meand18O SDd18O MeandD SDdD n
## 1        Rain       -3    0.5    -20  2.0 5
## 2    Snowmelt      -10    0.8    -70  3.0 5
## 3 Groundwater       -7    0.6    -50  2.5 5
dim(discr_data)
## [1] 3 6
head(discr_data)
##        source Meand18O SDd18O MeandD SDdD n
## 1        Rain        0      0      0    0 5
## 2    Snowmelt        0      0      0    0 5
## 3 Groundwater        0      0      0    0 5

Step 2: Load Data into MixSIAR

The synthetic data are loaded into MixSIAR for analysis.

mix <- load_mix_data(
  filename = "precip_mix_data.csv",
  iso_names = c("d18O", "dD"),
  factors = NULL,
  fac_random = NULL,
  fac_nested = NULL,
  cont_effects = NULL
)

source <- load_source_data(
  filename = "precip_source_data.csv",
  source_factors = NULL,
  conc_dep = FALSE,
  data_type = "means",
  mix = mix
)

discr <- load_discr_data(
  filename = "precip_discr_data.csv",
  mix = mix
)

Step 3: Write JAGS Model File

The JAGS model file is created for the Bayesian mixing model.

model_filename <- "MixSIAR_model.txt"
write_JAGS_model(
  filename = model_filename,
  resid_err = TRUE,
  process_err = TRUE,
  mix = mix,
  source = source
)

Step 4: Run Bayesian Model

The Bayesian model is executed using a test run.

jags_output <- run_model(
  run = "test",  # Quick test; replace with "normal" or "long" for production
  mix = mix,
  source = source,
  discr = discr,
  model_filename = model_filename,
  alpha.prior = 1,
  resid_err = TRUE,
  process_err = TRUE
)

Step 5: Produce Results and Plots

Results are generated, including convergence diagnostics and visualizations.

output_options <- list(
  summary_save = TRUE,
  summary_name = "MixSIAR_summary",
  sup_post = TRUE,
  plot_post_save_pdf = TRUE,
  plot_post_name = "posterior_density",
  sup_pairs = TRUE,
  plot_pairs_save_pdf = TRUE,
  plot_pairs_name = "pairs_plot",
  sup_xy = TRUE,
  plot_xy_save_pdf = TRUE,
  plot_xy_name = "xy_plot",
  gelman = TRUE,
  heidel = FALSE,
  geweke = TRUE,
  diag_save = TRUE,
  diag_name = "diagnostics",
  indiv_effect = FALSE,
  plot_post_save_png = FALSE,
  plot_pairs_save_png = FALSE,
  plot_xy_save_png = FALSE,
  diag_save_ggmcmc = FALSE,
  return_obj = TRUE
)

# Generate results
output_JAGS(jags_output, mix, source, output_options)

# Summary of estimated proportions
print(jags_output$BUGSoutput$summary)

# Global posterior proportions
p_global <- jags_output$BUGSoutput$sims.list$p.global

# Boxplot
x11()
boxplot(p_global, main = "Boxplot of Global Proportions", 
        names = colnames(p_global))

# Histogram for the first source
x11()
hist(p_global[,1], breaks = 30, col = "lightblue", 
     main = "Posterior Density: First Source",
     xlab = "Proportion")

# Scatterplot matrix
x11()
pairs(p_global, main = "Matrix Plot of Source Proportions")

Conclusion

The MixSIAR model successfully produced source proportions. Convergence is generally good but could be improved. Results indicate that groundwater is the dominant source in the simulated mixture (~50%). Rain and snowmelt contribute similarly (~23–27%). Wide credible intervals suggest uncertainty in exact proportions. Some convergence diagnostics (Rhat > 1.1) recommend increasing iterations to stabilize chains.

References

Author: “Abdi-Basid ADAN” Email: “”