Here I have constructed a fan diagram for the Shoal field site from Little Ambergris Cay in the Turks and Caicos Islands.

#first we want to define all variables for our system 
tempC <- 32 #degrees C
#all concentrations here are in mmol/kgw
cCa <- 11
cMg <- 62   
cNa <- 588
cK <-15
cCl <- 611
cSO4 <- 31

#define DIC as a vector of values (units are mmol/kgw) this will be the x axis of the plot 
DIC <- seq(from = 1, to = 3, by = 0.01)

#define Alk as a vector of values (units are mmol/kgw) this will be the y axis of the plot 
Alk <- seq(from = 1, to = 3, by = 0.01)

#pre-allocate space for our output. Note that now we're creating a matrix with number of rows corresponding to the length of our  vector and number of columns corresponding to the length of our [Alk] vector
Omega_arag <- matrix(0,length(DIC),length(Alk))


# here we load the libraries we will need. (metR) is to labeling the contour lines in the plot
library(phreeqc)
library(tidyverse)
## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ dplyr     1.1.4     ✔ readr     2.1.4
## ✔ forcats   1.0.0     ✔ stringr   1.5.1
## ✔ ggplot2   3.5.1     ✔ tibble    3.2.1
## ✔ lubridate 1.9.3     ✔ tidyr     1.3.0
## ✔ purrr     1.0.2     
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag()    masks stats::lag()
## ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors
library(metR)
## 
## Attaching package: 'metR'
## 
## The following object is masked from 'package:purrr':
## 
##     cross
phrLoadDatabaseString(phreeqc.dat)
phrSetOutputStringsOn(TRUE)

# we denote [nn] and [mm] as the variables in used in the forloop

for (nn in 1:length(DIC)) {
  for (mm in 1:length(Alk)){
  input <- c(   
    '  SOLUTION 1 LAC water',
    '  units         mmol/kgw',
    paste('  temp              ',as.character(tempC)),
    paste('  Alkalinity        ',as.character(Alk[mm])),
    paste('  C(4)              ',as.character(DIC[nn])),
    paste('  Ca                ',as.character(cCa)),
    paste('  Mg                ',as.character(cMg)),
    paste('  Na                ',as.character(cNa)),
    paste('  K                 ',as.character(cK)),
    paste('  Cl                ',as.character(cCl)),
    paste('  S(6)              ',as.character(cSO4)),

#here we can use a smaller set of selected outputs since we're only tracking one output variable:
    'SELECTED_OUTPUT       ',
    '  -high precision   TRUE',
    '  -si               aragonite')

  phrRunString(input)
  output <- phrGetSelectedOutput()

#save the Omega value into our output variable
  Omega_arag[nn,mm] <- 10^output$n1$si_aragonite
}
}

#put the data into a dataframe so we can plot it in ggplot. this is trickier than in previous examples because our Omega data is a matrix, not a vector, so we have to do some extra manipulation

#first we use the values in our concentration vectors to name the rows and columns in our matrix of Omega data
rownames(Omega_arag) <- DIC
colnames(Omega_arag) <- Alk

#this function checks to see if Omega_arag is a data frame and makes it one if it isn't already. the "%>%" symbol passes the result of that function to the next line.
as.data.frame(Omega_arag) %>% 
#this changes our row names (which we set to be the values in DIC a few lines above) into a column in the dataframe
  rownames_to_column(var="DIC") %>% 
#this function pivots all the data for which the column name is a variable into two columns, one of that variable (Alk), and one of the data (Omega)
  pivot_longer(
    cols = c(2:length(Alk)+1),
    names_to = 'Alk',
    values_to = "Omega"
    ) %>%
#this function changes the values in the DIC and Alk columns into numeric data so it can be plotted
  mutate(DIC=as.numeric(DIC), 
         Alk=as.numeric(Alk))%>%
#NOW we're ready to plot:
  ggplot() +
  
#these line makes the contour plot with contour labels 
  geom_contour_filled(aes(x = DIC, y = Alk, z = Omega), breaks=seq(from=0, to= 18, by=1)) +
  geom_contour2(aes(x = DIC, y = Alk, z = Omega), color = "black", breaks=seq(from=0, to= 18, by=1)) +
  geom_text_contour(aes(x = DIC, y = Alk, z = Omega), breaks=seq(from=1, to= 18, by=2), min.size = 0.5,
  skip = 0, stroke = 0.2, label.placer = label_placer_random()) +
  
  
#this line adds a point at the ([DIC,[Alk) condition in our base in the system 

  #vector arrow for precipitation 
   geom_segment(
    x = 2.2, #DIC
    y = 2.4,  #Alk
    xend = 2.1, #DIC
    yend = 2.2,  #Alk
    arrow = arrow(length = unit(0.5,"cm"),angle=20),
    color = "white"
    ) +
  annotate(
    "text",
    x=2.2,
    y=2.1,
    label="precipitation",
    size=4,
    color="white"
    ) +
  
    #vector arrow for dissolutuion 
   geom_segment(
    x = 2.2, #DIC
    y = 2.4, #ALk
    xend = 2.3, #DIC
    yend = 2.6, #Alk
    arrow = arrow(length = unit(0.5,"cm"),angle=20),
    color = "white"
    ) +
  annotate(
    "text",
    x=2.1,
    y=2.6,
    label="dissolution",
    size=4,
    color="white"
    ) +
  
  #vector arrow for heterotrophy/respiration   
  
     geom_segment(
    x = 2.2, #DIC
    y = 2.4, #ALk
    xend = 2.3, #DIC
    yend = 2.4, #Alk
    arrow = arrow(length = unit(0.5,"cm"),angle=20),
    color = "white"
    ) +
  annotate(
    "text",
    x=2.5,
    y=2.3,
    label="respiration",
    size=4,
    color="white"
    ) +
  
    #vector arrow for photosynthesis
  
   geom_segment(
    x = 2.2, #DIC
    y = 2.4, #ALk
    xend = 2.1, #DIC
    yend = 2.4, #Alk
    arrow = arrow(length = unit(0.5,"cm"),angle=20),
    color = "white"
    ) +
  annotate(
    "text",
    x=1.8,
    y=2.4,
    label="photosynthesis",
    size=4,
    color="white"
    ) +
  
  #this line adds a point at the ([DIC,[Alk) condition in our base in the system I put this at the end, rater than earlier so the dot appears over the vector arrows 
  geom_point(
    x = 2.2, 
    y = 2.4,
    pch = 21,
    size = 3,
    color = 'black',
    fill = "orange"
    ) +
#customize the axis labels
  labs(y = '[Alk]  (meq/L)',x = '[DIC] (mmol/kgw)',fill = bquote(Omega[arag])) 

#use this to save the plot as an .eps file, or whatever filetype you want (.pdf etc)
  ggsave("shoalfan.eps")
## Saving 7 x 5 in image