Paquete hydrochem de R: tratamiento de datos de hidroquímica

Guía práctica, conceptos químicos y manual comentado con gráficos

1 Sobre esta página

Actualmente el paquete hydrochem está disponible en CRAN, en https://cran.r-project.org/web/packages/hydrochem/index.html. Ahí se encuentra el manual en html y pdf.

Tiene el formato estándar de los paquetes del CRAN, sin gráficos. Como una funcionalidad del paquete es la generación de diagramas, para ver su potencial en este HTML se reproducen los ejemplos con las gráficas resultantes. El script de Quarto base ha sido generado con Claude Opus 5.5.

2 Qué es hydrochem

hydrochem (versión 0.1.0, publicada en CRAN el 24-06-2026; autores Guillaume Cinkus y Raphaël Bondu) es un paquete de R para procesar, analizar y representar datos hidroquímicos dentro de un flujo de trabajo programático y reproducible. Su objetivo declarado es sustituir el clásico ciclo hoja de cálculo + software gráfico propietario (AquaChem, Diagrammes, etc.) por funciones que se pueden encadenar en un script.

Cubre cuatro bloques:

Bloque	Funciones
Preparación y control de calidad	replace_bdl(), convert_unit(), ionic_balance(), stat_summary()
Caracterización del agua	water_type(), tds(), total_hardness(), total_alkalinity()
Índices hidrogeoquímicos y de riego	sar(), percent_sodium(), kelly_ratio(), permeability_index(), magnesium_hazard(), rsc(), cai1(), cai2()
Diagramas	plot_ternary(), plot_piper(), plot_durov(), plot_schoeller(), plot_stiff(), plot_collins(), plot_stabler(), plot_gibbs(), plot_biplot(), save_plot()

Algunas ideas de diseño que conviene tener claras desde el principio:

• Todas las funciones reciben un data.frame (o data.table) en formato ancho: una fila por muestra y una columna por ion.
• Las funciones de cálculo devuelven un vector (un valor por fila), pensado para añadirlo como nueva columna. Encaja de forma natural con dplyr::mutate().
• Las funciones de conversión y depuración (convert_unit(), replace_bdl(), stat_summary()) devuelven un data.table.
• Las funciones de diagramas devuelven objetos ggplot (o un cowplot en el caso de los biplots), así que se pueden modificar o guardar como cualquier gráfico de ggplot2.
• Internamente el paquete trabaja con data.table y ggplot2; no depende del tidyverse, pero convive bien con él.
• Muchas funciones emiten mensajes informativos ("TDS computed with: ...") que indican qué iones han entrado en el cálculo. Conviene leerlos: son la primera auditoría del resultado.

2.1 Instalación

install.packages("hydrochem")

Requiere R ≥ 4.1 y versiones recientes de sus dependencias: data.table ≥ 1.16.2, ggplot2 ≥ 3.5.1, cowplot ≥ 1.2.0 y scales ≥ 1.4.0. Si tenéis una instalación de R antigua, actualizad antes esas librerías o la instalación fallará al cargar el namespace.

library(hydrochem)

Warning: package 'hydrochem' was built under R version 4.6.1

library(tidyverse)

packageVersion("hydrochem")

[1] '0.1.0'

AdvertenciaColores por grupo en R antiguos (R 4.3 o anterior)

Los diagramas colorean los grupos con la paleta Okabe-Ito de grDevices::palette.colors(). En las versiones actuales de R esa función devuelve un vector sin nombres (argumento names = FALSE), pero en R 4.3 y anteriores ese argumento no existe y el vector lleva nombres ("black", "orange"…). scale_fill_manual() intenta entonces emparejar esos nombres con los grupos ("A", "B"…), no encuentra ninguno y todos los puntos salen del mismo color, sin leyenda. Solo afecta a group sin group_custom.

El siguiente bloque solo actúa si vuestro R tiene la versión antigua de la función. Sustituye, dentro del entorno de hydrochem, palette.colors() por una versión que quita los nombres. No modifica nada en disco.

if (!"names" %in% names(formals(grDevices::palette.colors))) {
  imp <- parent.env(asNamespace("hydrochem"))   # entorno de importaciones
  unlockBinding("palette.colors", imp)
  assign("palette.colors",
         function(...) unname(grDevices::palette.colors(...)),
         envir = imp)
  lockBinding("palette.colors", imp)
  message("Parche de paleta aplicado (R ", getRversion(), ")")
}

3 Conceptos químicos de partida

Antes de usar el paquete conviene repasar cuatro ideas que aparecen en todas las funciones.

3.1 Iones mayoritarios

En la mayoría de las aguas naturales, más del 90 % de la materia disuelta corresponde a unos pocos iones, los iones mayoritarios:

• Cationes (carga positiva): calcio (Ca²⁺), magnesio (Mg²⁺), sodio (Na⁺) y potasio (K⁺).
• Aniones (carga negativa): bicarbonato (HCO₃⁻), carbonato (CO₃²⁻, solo relevante a pH > 8,3), cloruro (Cl⁻), sulfato (SO₄²⁻) y, en aguas afectadas por actividad agraria o urbana, nitrato (NO₃⁻).

El resto (Fe, Mn, Sr, Br, F, B, metales…) son iones minoritarios o traza: no pesan en el balance de cargas, pero sí en la calidad y en la interpretación del origen del agua. La sílice (SiO₂) se disuelve como especie neutra (H₄SiO₄) y no interviene en el balance de cargas, aunque sí en los sólidos disueltos.

3.2 Unidades: mg/L, mmol/L y meq/L

Los laboratorios suelen entregar las concentraciones en masa por volumen (mg/L). Pero para comparar iones entre sí hay que pasar a unidades químicas:

$$\text{mmol/L} = \frac{\text{mg/L}}{M} \qquad\qquad \text{meq/L} = \text{mmol/L} \times |z|$$

donde $M$ es la masa molar del ion (g/mol = mg/mmol) y $|z|$ el valor absoluto de su carga.

Por ejemplo, 80 mg/L de Ca²⁺ ($M$ = 40,078; $z$ = 2) son 80 / 40,078 = 2,00 mmol/L y 4,00 meq/L. Esos mismos 80 mg/L de Na⁺ ($M$ = 22,99; $z$ = 1) son 3,48 meq/L. En masa son iguales; en “capacidad de carga eléctrica” no.

Los miliequivalentes son la unidad natural de la hidroquímica porque reflejan la electroneutralidad: en una disolución la suma de cargas positivas iguala a la de cargas negativas. Por eso casi todos los índices y diagramas (Piper, Stiff, Schoeller, SAR…) se calculan en meq/L. hydrochem hace esa conversión internamente: basta con decirle en qué unidad vienen los datos mediante el argumento base_unit.

3.3 Formas elementales: nitrógeno nítrico, fósforo, azufre

Muchos laboratorios expresan algunos compuestos como elemento: nitrato como N (NO₃–N), fosfato como P (PO₄–P), amonio como N (NH₄–N). No es lo mismo 10 mg/L de NO₃ que 10 mg/L de NO₃–N: el factor es la relación de masas molares, 62,00 / 14,01 = 4,43. El paquete reconoce estas formas por su nombre de columna (NO3_N, PO4_P, NH4_N, NO2_N, H2S_HS_S) y puede convertirlas a la especie completa (convert_to_species = TRUE).

3.4 Datos censurados: valores bajo el límite de detección

Cuando la concentración es inferior a lo que el método analítico puede cuantificar, el laboratorio no da un número sino una desigualdad: "< 0,01". Es un dato censurado: sabemos que está entre 0 y el límite de detección (LD), pero no su valor. Para operar con él hay que sustituirlo por un valor numérico. Las convenciones habituales son LD/2, LD/√2, LD, 0 o dejarlo como ausente (NA). Ninguna es neutra: con muchos valores censurados la elección sesga medias y percentiles, y conviene documentarla.

4 Los datos de entrada

4.1 Reglas de formato

Para que hydrochem reconozca los datos, la tabla debe cumplir:

1. Formato ancho: una fila por muestra, una columna por parámetro.
2. Nombres de columna = símbolo químico exactamente como lo espera el paquete, respetando mayúsculas: Ca, Mg, Na, K, HCO3, CO3, Cl, SO4, NO3, SiO2, pH… Una columna llamada calcio, CA o Ca (mg/L) será ignorada sin aviso o provocará el error Missing obligatory ions.
3. Todas las concentraciones en la misma unidad (la que se indique en base_unit).
4. Los valores bajo el límite de detección pueden venir como texto ("<0.01", "< .5"); replace_bdl() los convierte.
5. Se admiten columnas adicionales (identificador, fecha, coordenadas, acuífero, grupo…). Algunas funciones las usan como agrupación (group), etiqueta (label) o identificador de panel (idcol).

El paquete trae una tabla interna con los 72 solutos que reconoce, su masa molar y su carga. No está exportada, pero se puede consultar con :::. Estos son los que intervienen en el balance iónico:

hydrochem:::ions_meta |>
  as_tibble() |>
  filter(BI %in% c("major", "minor")) |>
  select(columna = ion_r, nombre = ion, tipo = type,
         masa_molar = masse_molaire, carga = charge, balance = BI) |>
  arrange(desc(balance == "major"), tipo, columna) |>
  knitr::kable(digits = 3)

columna	nombre	tipo	masa_molar	carga	balance
CO3	Carbonate	anion	60.009	2	major
Cl	Chloride	anion	35.450	1	major
HCO3	Bicarbonate	anion	61.017	1	major
NO3	Nitrate	anion	62.005	1	major
NO3_N	Nitrate as nitrogen	anion	14.010	1	major
SO4	Sulfate	anion	96.060	2	major
Ca	Calcium	cation	40.078	2	major
K	Potassium	cation	39.098	1	major
Mg	Magnesium	cation	24.305	2	major
Na	Sodium	cation	22.990	1	major
Br	Bromide	anion	79.904	1	minor
F	Fluoride	anion	18.998	1	minor
I	Iodine	anion	126.904	1	minor
NO2	Nitrite	anion	46.005	1	minor
NO2_N	Nitrite as nitrogen	anion	14.010	1	minor
PO4	Phosphate	anion	94.971	3	minor
PO4_P	Phosphate as phosphorus	anion	30.970	3	minor
Al	Aluminium	cation	26.982	3	minor
Ba	Barium	cation	137.327	2	minor
Cu	Copper	cation	63.546	2	minor
Fe	Iron	cation	55.845	2	minor
Li	Lithium	cation	6.941	1	minor
Mn	Manganese	cation	54.938	2	minor
NH4	Ammonium	cation	18.039	1	minor
NH4_N	Ammonium as nitrogen	cation	14.010	1	minor
Ni	Nickel	cation	58.693	2	minor
Sr	Strontium	cation	87.620	2	minor
Zn	Zinc	cation	65.380	2	minor

4.2 Conjunto de ejemplo 1: hc_data

Datos de campo anonimizados: 59 muestras con los iones mayoritarios en mg/L.

glimpse(hc_data)

Rows: 59
Columns: 13
$ id   <chr> "id_1", "id_2", "id_3", "id_4", "id_5", "id_6", "id_7", "id_8", "…
$ date <dttm> 1981-01-17, 1981-01-24, 1981-01-31, 1981-02-07, 1981-02-14, 1981…
$ EC   <dbl> 380, 378, 381, 377, 374, 379, 374, 375, 375, 377, 368, 296, 302, …
$ HCO3 <dbl> 255.6, 257.1, 261.1, 267.8, 261.6, 260.5, 261.5, 263.5, 261.0, 25…
$ Cl   <dbl> 4.60, 4.80, 4.50, 3.60, 4.70, 4.70, 4.80, 5.00, 4.90, 3.90, 5.00,…
$ NO3  <dbl> 0.80, 0.87, 0.80, 0.78, 0.84, 0.83, 0.99, 1.00, 0.89, 0.91, 0.88,…
$ SO4  <dbl> 12.2, 12.6, 14.2, 14.1, 14.5, 14.7, 14.9, 15.9, 14.9, 15.0, 14.4,…
$ Ca   <dbl> 74.8, 74.7, 77.6, 77.4, 76.1, 79.3, 78.6, 78.8, 79.0, 78.8, 78.4,…
$ Mg   <dbl> 6.83, 6.82, 6.93, 6.90, 7.15, 7.24, 6.78, 6.97, 6.65, 6.53, 6.16,…
$ Na   <dbl> 1.59, 1.39, 1.40, 1.49, 1.47, 1.49, 1.48, 1.42, 1.56, 1.64, 1.51,…
$ K    <dbl> 0.48, 0.47, 0.48, 0.47, 0.49, 0.47, 0.48, 0.47, 0.46, 0.50, 0.45,…
$ SiO2 <dbl> 7.48, 7.28, 7.80, 7.48, 7.44, 7.44, 7.40, 7.44, 7.44, 7.48, 7.20,…
$ type <chr> "A", "A", "A", "A", "A", "A", "A", "A", "A", "A", "A", "A", "A", …

• id: identificador de la muestra (id_1, id_2…).
• date: fecha de muestreo (semanal, desde 1981).
• EC: conductividad eléctrica (el manual no indica la unidad; por su orden de magnitud, previsiblemente µS/cm).
• HCO3, Cl, NO3, SO4, Ca, Mg, Na, K, SiO2: concentraciones en mg/L. No hay CO3 ni pH, lo que provoca algunos avisos más adelante.
• type: variable de agrupación con tres niveles (A, B, C).

hc_data |> count(type)

4.3 Conjunto de ejemplo 2: data_qc

Datos reales de 146 muestras de aguas subterráneas del sistema acuífero en rocas sedimentarias paleozoicas de las Basses-Laurentides (Quebec, Canadá), procedentes de la tesis doctoral de V. Cloutier (2004, INRS). Es un caso más realista: 55 columnas, metadatos, iones minoritarios, valores censurados e isótopos.

glimpse(data_qc)

Rows: 146
Columns: 55
$ Sample_ID            <chr> "S1", "S2", "S3", "S4", "S5", "S6", "S7", "S8", "…
$ Sampling_date        <date> 2000-10-30, 2000-10-30, 2000-10-30, 2000-10-30, …
$ Easting_UTMNAD83Z18  <dbl> 569337.0, 569337.0, 575779.0, 575775.9, 579574.1,…
$ Northing_UTMNAD83Z18 <dbl> 5058046, 5058046, 5056577, 5056578, 5061306, 5061…
$ Elevation            <dbl> 74.180, 74.970, 79.990, 81.050, 69.890, 70.290, 7…
$ Geo_units            <chr> "Surface sediments", "Theresa F.", "Surface sedim…
$ Hydro_units          <chr> "Surface sediments", "Rock", "Surface sediments",…
$ Hydro_cond           <chr> "Surface sediments and springs", "Unconfined", "S…
$ Artesian_well        <chr> "No", "No", "No", "No", "No", "No", "No", "No", "…
$ Top_open_interval    <dbl> 4.89, 10.30, 2.76, 9.63, 8.85, 1.00, 13.10, 3.38,…
$ Bottom_open_interval <dbl> 6.80, 11.80, 4.60, 13.00, 10.67, 4.00, 14.63, 7.0…
$ Middle_open_interval <chr> "5.85", "11.05", "3.51", "11.13", "9.6", "2.5", "…
$ Cluster              <chr> "Cluster 7", "Cluster 7", "Cluster 7", "Cluster 7…
$ T                    <dbl> 8.24, 7.79, 10.87, 8.54, 8.11, 10.46, 9.10, 9.66,…
$ pH                   <dbl> 7.30, 7.38, 7.50, 7.25, 7.30, 7.51, 7.43, 7.14, 7…
$ EC                   <dbl> 718, 742, 409, 634, 741, 717, 692, 650, 435, 654,…
$ DO                   <dbl> 2.36, 0.98, 1.22, 1.63, 0.47, 0.26, 0.60, 1.96, N…
$ Eh                   <dbl> 187.032, 227.347, 294.191, 247.822, 247.123, 261.…
$ Ca                   <dbl> 81.0, 76.0, 53.0, 87.0, 83.0, 84.0, 72.0, 95.0, 4…
$ Mg                   <dbl> 36.0, 39.0, 23.0, 27.0, 47.0, 44.0, 28.0, 26.0, 2…
$ Na                   <dbl> 26.0, 21.0, 3.7, 49.0, 10.1, 9.4, 46.0, 7.2, 3.0,…
$ K                    <dbl> 3.80, 6.50, 0.76, 4.70, 4.80, 4.60, 5.00, 13.00, …
$ Cl                   <chr> "49", "72", "0.6", "66", "15", "0.8", "37", "15",…
$ SO4                  <chr> "22", "34", "49", "50", "80", "75", "130", "50", …
$ HCO3                 <dbl> 367.5482, 320.2608, 214.3431, 338.1995, 414.7502,…
$ CO3                  <dbl> 0.3028724, 0.3134761, 0.2860279, 0.2520686, 0.345…
$ Fe                   <chr> "1.4", "0.52", "0.34", "0.84", "0.48", "0.97", "2…
$ Mn                   <chr> "0.15", "0.069", "0.3", "0.08", "0.028", "0.13", …
$ Br                   <chr> "0.02", "0.04", "< .01", "0.04", "< .01", "< .01"…
$ Sr                   <chr> "0.24", "0.34", "0.107", "0.45", "0.32", "0.42", …
$ F                    <dbl> 0.20, 0.24, 0.04, 0.06, 0.20, 0.34, 0.22, 0.15, 0…
$ Ba                   <dbl> 0.22, 0.22, 0.06, 0.07, 0.17, 0.01, 0.07, 0.06, 0…
$ H2S_HS_S             <chr> "0.02", "< .02", "< .02", "< .02", "0.02", "0.02"…
$ SiO2                 <dbl> 15.7, 14.3, 10.7, 11.2, 13.7, 11.5, 9.4, 9.3, 11.…
$ B                    <chr> "< .002", "0.03", "< .002", "< .002", "0.02", "< …
$ Al                   <chr> "< .007", "< .007", "< .007", "0.7", "< .007", "<…
$ Li                   <chr> "0.005", "0.004", "0.001", "< .001", "0.006", "0.…
$ NO3_N                <chr> "< .02", "< .02", "< .02", "0.04", "< .02", "< .0…
$ NH4_N                <chr> "0.11", "0.03", "< .02", "0.02", "0.11", "0.06", …
$ PO4_P                <chr> "0.04", "0.01", "< .01", "< .01", "0.025", "0.025…
$ DOC                  <chr> "1.5", "1.5", "4", "1.9", "0.9", "1.4", "1.6", "1…
$ Ag                   <chr> "< .0003", "< .0003", "< .0003", "< .0003", NA, N…
$ As                   <chr> "0.001", "< .001", "0.001", "0.002", NA, NA, NA, …
$ Cd                   <chr> "0.0018", "0.0017", "0.0013", "0.0016", NA, NA, N…
$ Cr                   <chr> "0.003", "0.001", "0.002", "0.004", "0.006", "0.0…
$ Cu                   <chr> "< .001", "< .001", "< .001", "< .001", "< .001",…
$ Hg                   <chr> NA, NA, NA, NA, NA, NA, NA, NA, "< .0001", "< .00…
$ I                    <chr> "< .1", "< .1", "< .1", "< .1", NA, NA, NA, NA, "…
$ Ni                   <chr> "0.002", "0.002", "0.003", "0.004", "0.004", "0.0…
$ Pb                   <chr> "< .004", "< .004", "< .004", "< .004", NA, NA, N…
$ Se                   <chr> "< .001", "< .001", "< .001", "< .001", "< .001",…
$ Zn                   <chr> "0.078", "< .005", "< .005", "1.2", "< .005", "< …
$ U                    <chr> "< .005", "< .005", "< .005", "0.012", NA, NA, NA…
$ d2H                  <dbl> -69.3, -69.7, -55.8, -77.1, -71.4, -67.9, -77.4, …
$ d18O                 <dbl> -10.60, -10.65, -8.47, -11.57, -10.87, -10.43, -1…

Se distinguen varios grupos de columnas:

• Metadatos del punto: Sample_ID, Sampling_date, coordenadas UTM (NAD83, huso 18), Elevation, unidad geológica e hidrogeológica (Geo_units, Hydro_units), condición hidráulica (Hydro_cond: libre, confinado, semiconfinado, sedimentos superficiales/manantiales), pozo surgente (Artesian_well) e intervalo ranurado.
• Agrupación: Cluster, resultado de un análisis clúster (7 grupos).
• Parámetros físico-químicos de campo: temperatura T, pH, conductividad EC, oxígeno disuelto DO, potencial redox Eh.
• Mayoritarios (Ca, Mg, Na, K, Cl, SO4, HCO3, CO3) y minoritarios/traza (Fe, Mn, Br, Sr, F, Ba, B, Al, Li, NO3_N, NH4_N, PO4_P, metales…). Muchos están en formato texto porque contienen valores censurados ("< .01").
• Isótopos estables del agua: d2H y d18O (‰ respecto a VSMOW).

# Columnas de texto que en realidad son concentraciones con valores censurados
data_qc |>
  select(where(is.character)) |>
  select(-c(Sample_ID:Artesian_well, Middle_open_interval, Cluster)) |>
  names()

 [1] "Cl"       "SO4"      "Fe"       "Mn"       "Br"       "Sr"      
 [7] "H2S_HS_S" "B"        "Al"       "Li"       "NO3_N"    "NH4_N"   
[13] "PO4_P"    "DOC"      "Ag"       "As"       "Cd"       "Cr"      
[19] "Cu"       "Hg"       "I"        "Ni"       "Pb"       "Se"      
[25] "Zn"       "U"       

data_qc |> count(Hydro_cond)

4.4 Cómo adaptar vuestros propios datos

Los datos de laboratorio o de las bases de datos de seguimiento suelen venir en formato largo (punto, fecha, parámetro, valor) y con nombres de parámetro en castellano. El paso a la estructura que exige hydrochem es un pivot_wider() con una tabla de equivalencias:

lab_largo <- tribble(
  ~punto,  ~fecha,       ~parametro,     ~valor,
  "P-01",  "2025-03-12", "Calcio",       "85.2",
  "P-01",  "2025-03-12", "Magnesio",     "21.4",
  "P-01",  "2025-03-12", "Sodio",        "18.9",
  "P-01",  "2025-03-12", "Potasio",      "2.1",
  "P-01",  "2025-03-12", "Bicarbonatos", "290",
  "P-01",  "2025-03-12", "Cloruros",     "31.5",
  "P-01",  "2025-03-12", "Sulfatos",     "62.0",
  "P-01",  "2025-03-12", "Nitratos",     "<1",
  "P-02",  "2025-03-12", "Calcio",       "120.0",
  "P-02",  "2025-03-12", "Magnesio",     "35.0",
  "P-02",  "2025-03-12", "Sodio",        "95.0",
  "P-02",  "2025-03-12", "Potasio",      "4.5",
  "P-02",  "2025-03-12", "Bicarbonatos", "310",
  "P-02",  "2025-03-12", "Cloruros",     "160",
  "P-02",  "2025-03-12", "Sulfatos",     "180",
  "P-02",  "2025-03-12", "Nitratos",     "28"
)

equivalencias <- c(Calcio = "Ca", Magnesio = "Mg", Sodio = "Na", Potasio = "K",
                   Bicarbonatos = "HCO3", Cloruros = "Cl", Sulfatos = "SO4",
                   Nitratos = "NO3")

mis_datos <- lab_largo |>
  mutate(parametro = equivalencias[parametro]) |>
  pivot_wider(names_from = parametro, values_from = valor) |>
  replace_bdl(method = "DL/2") |>          # "<1" -> 0.5 y columna numérica
  mutate(across(Ca:SO4, as.numeric)) |>    # el resto de texto a número
  as_tibble()

BDL replacement applied to: NO3

mis_datos

Con eso mis_datos ya sirve para cualquier función del paquete. Más adelante se usa en ionic_balance() y water_type().

5 Preparación y control de calidad

5.1 replace_bdl(): valores bajo el límite de detección

Detecta en las columnas de texto los valores del tipo "<0.01", "< 0.005" o " < .5 ", los sustituye según el método elegido ("DL/2" por defecto, "DL/sqrt(2)", "DL", "0" o "NA") y convierte la columna a numérica. Solo toca las columnas que contienen al menos un valor censurado (o las indicadas en cols). Devuelve un data.table.

df <- data.frame(Ca = c(" <0.01", "5.2", "< 0.005"),
                 Cl = c("3.1", "<0.02", "2.5"))
replace_bdl(df, method = "DL/2")

BDL replacement applied to: Ca, Cl

Aplicado a data_qc con sustitución por cero:

data_cero <- replace_bdl(data_qc, method = "0")

BDL replacement applied to: Cl, SO4, Fe, Mn, Br, Sr, H2S_HS_S, B, Al, Li, NO3_N, NH4_N, PO4_P, DOC, Ag, As, Cd, Cr, Cu, Hg, I, Ni, Pb, Se, Zn, U

data_cero |> select(Sample_ID, Cl, SO4, Br, NO3_N, Zn) |> head()

A partir de aquí se usa data_qc con la sustitución por defecto (LD/2), como en los ejemplos del manual:

data <- replace_bdl(data_qc)

BDL replacement applied to: Cl, SO4, Fe, Mn, Br, Sr, H2S_HS_S, B, Al, Li, NO3_N, NH4_N, PO4_P, DOC, Ag, As, Cd, Cr, Cu, Hg, I, Ni, Pb, Se, Zn, U

Nota

Middle_open_interval sigue siendo texto: no contiene ningún "<", así que replace_bdl() no la procesa. Si se necesita numérica, hay que convertirla aparte.

5.2 convert_unit(): cambio de unidades

Convierte las columnas de iones indicadas (ions) entre mg/L, mmol/L y meq/L. Añade columnas nuevas con sufijo _mg, _mmol o _meq y conserva las originales.

hc_meq <- convert_unit(hc_data, base_unit = "mg/L", to_unit = "meq/L")
hc_meq |> select(id, Ca, Ca_meq, Mg, Mg_meq, HCO3, HCO3_meq) |> head()

Comprobación a mano para la primera muestra (74,8 mg/L de Ca):

74.8 / 40.078 * 2

[1] 3.732721

Con convert_to_species = TRUE las formas elementales pasan a especie. Por ejemplo, el nitrógeno nítrico de data_qc a mg/L de NO₃:

data |>
  convert_unit(ions = "NO3_N", base_unit = "mg/L", to_unit = "mg/L",
               convert_to_species = TRUE) |>
  select(Sample_ID, NO3_N, NO3_mg) |>
  filter(NO3_N > 0.1) |>
  head()

Ion 'NO3_N' converted to species form 'NO3'

5.3 ionic_balance(): error del balance iónico

Concepto. Por electroneutralidad, en un análisis perfecto la suma de cationes en meq/L ($T_c$) debería igualar a la de aniones ($T_a$). El error de balance mide la desviación:

$$EB\,(\%) = \frac{T_c - T_a}{T_c + T_a} \times 100$$

Un error grande indica un ion mal analizado, uno importante sin analizar (p. ej. CO₃²⁻ en aguas básicas, NO₃⁻ en zonas agrícolas, materia orgánica), errores de transcripción o unidades mal declaradas. Como criterio habitual se aceptan análisis con |EB| ≤ 5 %, relajándolo hasta ±10 % en aguas poco mineralizadas. Es el primer filtro antes de interpretar nada.

Con method = "major" (por defecto) solo entran los mayoritarios; con "major_minor" se añaden los minoritarios que haya en la tabla (Fe, Mn, Sr, NH₄, Br, F…).

hc_data$ib <- ionic_balance(hc_data, base_unit = "mg/L")

Warning: CO3 is missing

Ionic balance computed with: HCO3, Cl, NO3, SO4, Ca, Mg, Na, K

summary(hc_data$ib)

   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
-4.6239 -2.3248 -0.7509 -0.6631  0.3970  4.6926 

El aviso CO3 is missing es esperable: hc_data no tiene carbonato. En data_qc, comparando los dos métodos:

data <- data |>
  mutate(ib_major = ionic_balance(pick(everything())),
         ib_minor = ionic_balance(pick(everything()), method = "major_minor"))

data |>
  summarise(n = n(),
            fuera_5 = sum(abs(ib_major) > 5, na.rm = TRUE),
            fuera_10 = sum(abs(ib_major) > 10, na.rm = TRUE))

Una representación rápida (no incluida en el paquete) ayuda a ver la distribución de errores:

data |>
  ggplot(aes(ib_major)) +
  annotate("rect", xmin = -5, xmax = 5, ymin = -Inf, ymax = Inf, alpha = 0.15) +
  geom_histogram(binwidth = 1, fill = "steelblue", colour = "white") +
  labs(x = "Error de balance iónico (%)", y = "Nº de muestras") +
  theme_minimal()

Figura 1: Distribución del error de balance iónico en data_qc (iones mayoritarios). Banda gris: ±5 %.

Y con los datos propios del apartado anterior:

mis_datos |> mutate(eb = ionic_balance(pick(everything())))

5.4 stat_summary(): estadística descriptiva

Calcula de una vez un amplio resumen por parámetro: mínimo, percentiles (1, 5, 10, 25, 75, 90, 95, 99), mediana, máximo, media, media geométrica, moda, rango, rango intercuartílico, desviación típica, MAD, coeficiente de variación, número y porcentaje de datos y de NA. Se puede agrupar (group) y seleccionar estadísticos (stats).

Advertencia

hydrochem::stat_summary() tiene el mismo nombre que ggplot2::stat_summary(). Si cargáis el paquete después de ggplot2/tidyverse, el de hydrochem enmascara al de ggplot2. Para evitar sorpresas, llamad siempre a una de las dos con prefijo (hydrochem::stat_summary() o ggplot2::stat_summary()).

hydrochem::stat_summary(hc_data, params = c("Ca", "Mg", "Na"))

hydrochem::stat_summary(hc_data, group = "type", params = c("Ca", "Mg", "Na"),
                        stats = c("n", "min", "med", "mean", "max", "CV"))

6 Caracterización del agua

6.1 water_type(): facies hidroquímica

Concepto. La facies o tipo de agua resume la composición nombrando el catión y el anión dominantes en meq/L, en porcentaje sobre el total de su signo: un agua “bicarbonatada cálcica” (Ca-HCO₃) es típica de recarga en terrenos carbonatados; una “clorurada sódica” (Na-Cl) indica salinización, mezcla con agua marina, disolución de halita o un flujo largo y lento.

El paquete ofrece dos formatos:

• "short": solo el catión y el anión dominantes ("Ca-HCO3").
• "long" (por defecto): todos los iones que superan un umbral (threshold, 20 % por defecto), ordenados de mayor a menor ("Ca-Mg-HCO3-SO4"). Si ningún catión o ningún anión supera el umbral devuelve NA.

water_type(hc_data, base_unit = "mg/L", format = "short") |> table()

Ca-HCO3 
     59 

water_type(hc_data, base_unit = "mg/L", format = "long", threshold = 25) |> table()

Ca-HCO3 
     59 

En data_qc, el abanico de facies es mucho más amplio:

data <- data |> mutate(facies = water_type(pick(everything()), format = "short"))
data |> count(facies, sort = TRUE)

mis_datos |> mutate(facies = water_type(pick(everything())))

Nota

water_type() usa todos los cationes y aniones reconocidos presentes en la tabla, incluidos los minoritarios. Su peso en meq/L suele ser despreciable, pero si una columna de traza tiene un valor anómalo (p. ej. un error de unidades) puede alterar los porcentajes.

6.2 tds(): sólidos disueltos totales

Concepto. Los sólidos disueltos totales (TDS, total dissolved solids) son la masa total de sustancias disueltas por litro. tds() los estima sumando las concentraciones en mg/L de los iones:

• "major": Ca, Mg, Na, K, HCO₃, CO₃, Cl, SO₄, NO₃.
• "major_si": además, SiO₂ (si solo hay Si, lo convierte a SiO₂).
• "all_ions": todos los solutos reconocidos presentes.

tds(hc_data, base_unit = "mg/L") |> head()

Warning: CO3 is missing

TDS computed with: HCO3, Cl, NO3, SO4, Ca, Mg, Na, K

[1] 356.90 358.75 367.01 372.54 366.85 369.23

tds(hc_data, method = "major_si") |> head()

Warning: CO3 is missing

TDS computed with: HCO3, Cl, NO3, SO4, Ca, Mg, Na, K, SiO2

[1] 364.38 366.03 374.81 380.02 374.29 376.67

TipTDS calculado frente a residuo seco

La suma iónica no coincide con el residuo seco de laboratorio (evaporación a 180 °C): al evaporar, el bicarbonato pierde CO₂ y agua ($2\,\mathrm{HCO_3^-} \rightarrow \mathrm{CO_3^{2-}} + \mathrm{CO_2} + \mathrm{H_2O}$) y solo queda el 49 % de su masa. En aguas muy bicarbonatadas, el TDS por suma puede superar con claridad al residuo seco. Para compararlos, multiplicad el HCO₃ por 0,4917 antes de sumar.

6.3 total_hardness(): dureza total

Concepto. La dureza es la concentración de cationes divalentes, en la práctica Ca²⁺ + Mg²⁺, expresada como si fuera CaCO₃. Como el peso equivalente del CaCO₃ es 100,09 / 2 ≈ 50 mg/meq:

$$\text{Dureza (mg/L CaCO}_3) = 50 \times (\mathrm{Ca} + \mathrm{Mg})_{\text{meq/L}}$$

El argumento output permite expresarla también en grados franceses (°f, 1 °f = 10 mg/L CaCO₃), alemanes (°dH, 1 °dH = 17,848 mg/L) o ingleses (°e, 1 °e ≈ 14,3 mg/L). Una clasificación habitual (USGS) es: blanda < 60, moderadamente dura 60–120, dura 120–180 y muy dura > 180 mg/L CaCO₃.

total_hardness(hc_data, base_unit = "mg/L") |> head()

Total hardness computed with: Ca, Mg

Output is in mg/L_CaCO3

[1] 214.7373 214.4466 222.1351 221.5126 219.2975 227.6523

total_hardness(hc_data, output = "fH") |> head()

Total hardness computed with: Ca, Mg

Output is in fH

[1] 21.47373 21.44466 22.21351 22.15126 21.92975 22.76523

6.4 total_alkalinity(): alcalinidad total

Concepto. La alcalinidad es la capacidad del agua para neutralizar ácidos. En aguas naturales se debe casi por completo al sistema carbonatado:

$$\text{Alc} = [\mathrm{HCO_3^-}] + 2[\mathrm{CO_3^{2-}}] + [\mathrm{OH^-}] - [\mathrm{H^+}]$$

Se expresa en meq/L o en mg/L de CaCO₃ (×50). La función necesita HCO3 y pH, y suma CO3 si existe. Por eso el ejemplo del manual usa data_qc (que tiene pH) y no hc_data.

total_alkalinity(data_qc, base_unit = "mg/L") |> head()

Total alkalinity computed with: pH, HCO3, CO3

Output is in mg/L_CaCO3

[1] 301.6908 262.9589 176.1193 277.5564 340.4410 358.7769

total_alkalinity(data_qc, output = "meq/L") |> head()

Total alkalinity computed with: pH, HCO3, CO3

Output is in meq/L

[1] 6.033815 5.259179 3.522387 5.551128 6.808820 7.175538

Nota

En el rango de pH de las aguas naturales (6–9) los términos de OH⁻ y H⁺ son despreciables (< 0,01 meq/L). La alcalinidad resultante es, en la práctica, HCO₃ + CO₃ en meq/L. Revisando el código fuente, el término de OH⁻ se implementa como $10^{-14-\mathrm{pH}}$ en mol/L, sin pasar a meq/L, en lugar de $10^{\mathrm{pH}-14}$. A pH ≤ 9 el efecto es nulo a efectos prácticos, pero conviene saberlo si se trabaja con aguas muy alcalinas.

7 Índices hidrogeoquímicos y de aptitud para riego

Todos se calculan en meq/L. La función convierte desde base_unit y devuelve un vector. Los umbrales de interpretación que se indican son los valores orientativos habituales en la bibliografía, no forman parte del paquete. Antes de usarlos en un informe, contrastadlos con la fuente original o con la normativa aplicable.

7.1 Riesgo de sodificación del suelo

El exceso de sodio frente a calcio y magnesio en el agua de riego desplaza al Ca²⁺ y al Mg²⁺ del complejo de cambio del suelo. Las arcillas sódicas se dispersan, el suelo pierde estructura y cae su permeabilidad.

sar() – Relación de adsorción de sodio (SAR; Richards, 1954, USDA Handbook 60):

$$SAR = \frac{\mathrm{Na}}{\sqrt{(\mathrm{Ca} + \mathrm{Mg})/2}}$$

Clases habituales (sin considerar la salinidad): < 10 excelente, 10–18 buena, 18–26 dudosa, > 26 inadecuada. Se interpreta siempre junto con la conductividad (diagrama del USSL).

sar(hc_data, base_unit = "mg/L") |> head()

SAR computed with: Ca, Mg, Na

[1] 0.04719633 0.04128764 0.04085867 0.04354636 0.04317827 0.04295513

percent_sodium() – Porcentaje de sodio (%Na; clasificación de Wilcox, 1955):

$$\%\mathrm{Na} = \frac{\mathrm{Na}}{\mathrm{Ca} + \mathrm{Mg} + \mathrm{Na} + \mathrm{K}} \times 100$$

Clases habituales: < 20 excelente, 20–40 buena, 40–60 admisible, 60–80 dudosa, > 80 inadecuada.

percent_sodium(hc_data, base_unit = "mg/L") |> head()

Percent sodium computed with: Ca, Mg, Na, K

[1] 1.580398 1.386284 1.348500 1.437988 1.432885 1.399850

kelly_ratio() – Índice de Kelly:

$$KR = \frac{\mathrm{Na}}{\mathrm{Ca} + \mathrm{Mg}}$$

KR < 1 se considera apta para riego; KR > 1 indica exceso de sodio.

kelly_ratio(hc_data, base_unit = "mg/L") |> head()

Kelly ratio computed with: Ca, Mg, Na

[1] 0.01610366 0.01409712 0.01370710 0.01462927 0.01457869 0.01423473

rsc() – Carbonato sódico residual (RSC; Eaton, 1950):

$$RSC = (\mathrm{HCO_3} + \mathrm{CO_3}) - (\mathrm{Ca} + \mathrm{Mg})$$

Si el agua tiene más carbonatos que Ca + Mg, al concentrarse en el suelo precipitan carbonatos de calcio y magnesio y el sodio queda “libre” en proporción creciente. Clases habituales: < 1,25 meq/L segura, 1,25–2,5 marginal, > 2,5 inadecuada. Valores negativos indican que no hay carbonato residual.

rsc(hc_data, base_unit = "mg/L") |> head()

Warning: CO3 is missing

RSC computed with: HCO3, Ca, Mg

[1] -0.10573523 -0.07533869 -0.16355241 -0.04129742 -0.09860718 -0.28372949

7.2 Permeabilidad y magnesio

permeability_index() – Índice de permeabilidad (PI; Doneen, 1964):

$$PI = \frac{\mathrm{Na} + \sqrt{\mathrm{HCO_3}}}{\mathrm{Ca} + \mathrm{Mg} + \mathrm{Na}} \times 100$$

En la clasificación de Doneen: clase I (> 75 %) y clase II (25–75 %) aptas, clase III (< 25 %) no apta.

permeability_index(hc_data, base_unit = "mg/L") |> head()

Permeability index computed with: HCO3, Ca, Mg, Na

[1] 48.48565 48.58529 47.28457 48.04814 47.96816 46.14789

magnesium_hazard() – Riesgo de magnesio (MH; Szabolcs y Darab, 1964):

$$MH = \frac{\mathrm{Mg}}{\mathrm{Ca} + \mathrm{Mg}} \times 100$$

MH > 50 % se considera perjudicial: el exceso de Mg afecta a la estructura del suelo y a los rendimientos.

magnesium_hazard(hc_data, base_unit = "mg/L") |> head()

Magnesium hazard computed with: Ca, Mg

[1] 13.08632 13.08487 12.83573 12.81607 13.41457 13.08492

7.3 Índices cloro-alcalinos: intercambio iónico en el acuífero

Los índices cloro-alcalinos de Schoeller detectan si el agua ha intercambiado cationes con las arcillas o la materia orgánica del acuífero. Parten de que el cloruro es conservativo: en un agua “salina” pura, Na + K ≈ Cl.

$$CAI_1 = \frac{\mathrm{Cl} - (\mathrm{Na} + \mathrm{K})}{\mathrm{Cl}} \qquad CAI_2 = \frac{\mathrm{Cl} - (\mathrm{Na} + \mathrm{K})}{\mathrm{SO_4} + \mathrm{HCO_3} + \mathrm{CO_3} + \mathrm{NO_3}}$$

• CAI > 0: el agua ha cedido Na⁺ y K⁺ al terreno a cambio de Ca²⁺ y Mg²⁺ (intercambio inverso; típico de intrusión marina o de agua salina que entra en un acuífero de agua dulce).
• CAI < 0: el agua ha tomado Na⁺ y K⁺ del terreno y le ha cedido Ca²⁺ y Mg²⁺ (intercambio directo o ablandamiento natural; típico de aguas que avanzan por acuíferos con arcillas sódicas o de lavado de antiguas aguas marinas).

cai1(hc_data, base_unit = "mg/L") |> head()

CAI1 computed with: Cl, Na, K

[1] 0.3723974 0.4646860 0.4235568 0.2434149 0.4231914 0.4204880

cai2(hc_data, base_unit = "mg/L") |> head()

Warning: CO3 is missing

CAI2 computed with: HCO3, Cl, NO3, SO4, Na, K

[1] 0.01084453 0.01401335 0.01171960 0.00526488 0.01218983 0.01214898

Nota

cai2() exige NO3; en data_qc usa NO3_N y lo convierte automáticamente a NO₃.

7.4 Todos los índices a la vez (enfoque tidyverse)

Como cada índice es una función con la misma firma (data, base_unit), se puede aplicar la lista entera con purrr::map():

indices <- list(
  SAR  = sar,
  pNa  = percent_sodium,
  KR   = kelly_ratio,
  RSC  = rsc,
  PI   = permeability_index,
  MH   = magnesium_hazard,
  CAI1 = cai1,
  CAI2 = cai2,
  dureza_CaCO3 = total_hardness,
  TDS  = tds
)

hc_ind <- hc_data |>
  as_tibble() |>
  select(id, type) |>
  bind_cols(map(indices, \(f) f(hc_data)))

hc_ind

hc_ind |>
  pivot_longer(-c(id, type), names_to = "indice", values_to = "valor") |>
  ggplot(aes(type, valor, fill = type)) +
  geom_boxplot(show.legend = FALSE, outlier.size = 1) +
  facet_wrap(~indice, scales = "free_y", ncol = 5) +
  labs(x = "Grupo", y = NULL) +
  theme_minimal()

Figura 2: Distribución de los índices por grupo en hc_data.

8 Diagramas hidroquímicos

Los diagramas se presentan agrupados por su finalidad:

1. Clasificación de muchas muestras (diagramas triangulares): ternario, Piper y Durov.
2. Comparación de perfiles químicos en escala logarítmica: Schoeller.
3. Huella química de cada muestra: Stiff, Collins y Stabler.
4. Procesos que controlan la química: Gibbs y biplots.

Todos comparten una lógica de argumentos: base_unit, group (columna para colorear), group_custom (estilos manuales por grupo con vectores en fill, color, shape, size), label, plot_title, legend.position, legend.title y bg_color.

8.1 Diagramas triangulares

8.1.1 Fundamento: coordenadas ternarias

Si tres componentes se expresan en porcentaje de su suma ($a + b + c = 100$), cada muestra queda definida por dos de ellos y se puede situar en un triángulo equilátero: cada vértice representa el 100 % de un componente y el lado opuesto, el 0 %. Una muestra en el centro tiene un tercio de cada uno. El diagrama pierde la información de concentración absoluta y muestra solo proporciones relativas.

8.1.2 plot_ternary(): ternario libre

Representa tres solutos cualesquiera (first, second, third), sin restricciones de tipo (pueden mezclarse cationes, aniones o neutros). Se normalizan en mmol/L (por defecto) o mg/L (to_unit). Es útil para exploraciones a medida: Ca–Mg–Na, Cl–SO₄–HCO₃, Ca–Mg–Sr…

plot_ternary(hc_data, first = "Ca", second = "Mg", third = "Na")

Figura 3: Ternario Ca–Mg–Na básico.

plot_ternary(hc_data, first = "Ca", second = "Mg", third = "Na", group = "type")

Figura 4: Ternario con agrupación simple por type.

plot_ternary(
  hc_data, first = "Ca", second = "Mg", third = "Na",
  group = "type",
  group_custom = TRUE,
  fill  = c("steelblue", "tomato", "darkgreen"),
  color = c("black", "black", "black"),
  shape = c(21, 22, 24),
  size  = c(2.5, 3, 3.5)
)

Figura 5: Ternario con estilos personalizados por grupo.

Cuando un componente es muy minoritario, los puntos se amontonan contra un lado. La solución es transformarlo antes (manteniendo el nombre de columna) y reetiquetar el eje:

hc_mod <- hc_data
hc_mod$Mg <- hc_mod$Mg * 10
plot_ternary(hc_mod, first = "Ca", second = "Mg", third = "Na",
             first_label = "Ca", second_label = "10*Mg", third_label = "Na")

Figura 6: Ternario con la señal de Mg amplificada ×10.

8.1.3 plot_piper(): diagrama de Piper

Concepto (Piper, 1944). Es el diagrama de clasificación más extendido. Consta de:

• un triángulo de cationes (abajo a la izquierda): Ca²⁺, Mg²⁺ y Na⁺ + K⁺, en % meq/L;
• un triángulo de aniones (abajo a la derecha): HCO₃⁻ + CO₃²⁻, SO₄²⁻ y Cl⁻ (más NO₃⁻ si include_NO3 = TRUE), en % meq/L;
• un rombo central donde se proyecta cada muestra combinando ambos triángulos.

En el rombo, el vértice superior corresponde a aguas sulfatadas/cloruradas cálcico-magnésicas (dureza no carbonatada); el izquierdo, a bicarbonatadas cálcico-magnésicas (dureza carbonatada, aguas de recarga en carbonatos); el derecho, a cloruradas/sulfatadas sódicas (aguas salinas, marinas, evaporíticas); y el inferior, a bicarbonatadas sódicas (aguas que han sufrido intercambio iónico). Las mezclas entre dos aguas se alinean entre ambos extremos, y esa es una de sus mayores utilidades.

plot_piper(hc_data)

Figura 7: Piper básico sin agrupación.

plot_piper(hc_data, group = "type")

Figura 8: Piper agrupado por type.

plot_piper(
  hc_data,
  group = "type",
  group_custom = TRUE,
  fill  = c("steelblue", "tomato", "darkgreen"),
  color = c("black", "black", "black"),
  shape = c(21, 22, 24),
  size  = c(2.5, 3, 3.5)
)

Figura 9: Piper con estilos personalizados por grupo.

Con data_qc, que tiene facies mucho más variadas, el Piper es más informativo. Aquí se agrupa por condición hidráulica del acuífero:

plot_piper(data, group = "Hydro_cond", legend.title = "Condición")

Warning: Removed 4 rows containing missing values or values outside the scale range
(`geom_point()`).

Figura 10: Piper de data_qc agrupado por condición hidráulica.

El aviso Removed n rows containing missing values aparece porque dos muestras de data_qc no tienen sulfato (SO4 = NA): no se pueden situar en el triángulo de aniones ni en el rombo, y ggplot las omite. Lo mismo ocurre en otros diagramas con este conjunto de datos.

8.1.4 plot_durov(): diagrama de Durov

Concepto (Durov, 1948). Usa los mismos dos triángulos de cationes y aniones que el Piper, pero proyecta las muestras sobre un cuadrado (en vez de un rombo). Su ventaja es que admite paneles adicionales alineados con el cuadrado, con pH (abajo) y TDS (a la derecha). Así se ven a la vez la facies, la acidez y la mineralización. Es útil para identificar procesos de mezcla, intercambio iónico y disolución.

Argumentos propios: top_ions y left_ions (iones de cada triángulo), tds_method, plot_pH, plot_TDS y tds_log. Necesita la columna pH si plot_pH = TRUE, por eso los ejemplos usan data_qc.

plot_durov(data)

Ion 'NO3_N' converted to species form 'NO3'

TDS computed with: Ca, Mg, Na, K, Cl, SO4, HCO3, CO3, NO3

Warning: Removed 2 rows containing missing values or values outside the scale range
(`geom_point()`).
Removed 2 rows containing missing values or values outside the scale range
(`geom_point()`).
Removed 2 rows containing missing values or values outside the scale range
(`geom_point()`).

Figura 11: Durov básico con paneles de pH y TDS.

plot_durov(data, group = "Cluster")

Ion 'NO3_N' converted to species form 'NO3'

TDS computed with: Ca, Mg, Na, K, Cl, SO4, HCO3, CO3, NO3

Warning: Removed 2 rows containing missing values or values outside the scale range
(`geom_point()`).
Removed 2 rows containing missing values or values outside the scale range
(`geom_point()`).
Removed 2 rows containing missing values or values outside the scale range
(`geom_point()`).

Figura 12: Durov agrupado por clúster.

plot_durov(
  data,
  group = "Hydro_cond",
  group_custom = TRUE,
  fill  = c("steelblue", "tomato", "darkgreen", "yellow"),
  color = c("black", "black", "black", "black"),
  shape = c(21, 22, 24, 25),
  size  = c(2.5, 3, 3.5, 4)
)

Ion 'NO3_N' converted to species form 'NO3'

TDS computed with: Ca, Mg, Na, K, Cl, SO4, HCO3, CO3, NO3

Warning: Removed 2 rows containing missing values or values outside the scale range
(`geom_point()`).
Removed 2 rows containing missing values or values outside the scale range
(`geom_point()`).
Removed 2 rows containing missing values or values outside the scale range
(`geom_point()`).

Figura 13: Durov con estilos personalizados por condición hidráulica.

8.2 Comparación de perfiles: plot_schoeller()

Concepto (Schoeller, 1935). Cada muestra es una línea quebrada que une sus concentraciones en meq/L sobre un eje vertical logarítmico, con los iones en abscisas. La escala logarítmica permite comparar aguas de mineralización muy distinta:

• Líneas paralelas indican aguas con las mismas proporciones y distinta concentración (dilución, evaporación o mismo origen).
• Líneas que se cruzan indican cambios de facies.

El paquete dibuja además guías horizontales equivalentes en mg/L (mg_guide).

hc20 <- hc_data[1:20, ]   # primeras 20 muestras, como en el manual

plot_schoeller(hc20)

Figura 14: Schoeller básico (una línea por muestra).

plot_schoeller(hc20, legend.position = "none")

Figura 15: Schoeller sin leyenda.

plot_schoeller(hc20, mg_guide = "none")

Figura 16: Schoeller sin las guías de mg/L.

plot_schoeller(hc20, group = "id")

Figura 17: Schoeller agrupado por identificador de muestra.

plot_schoeller(hc20, group = "type")

Figura 18: Schoeller coloreado por type.

plot_schoeller(hc20, group = "type", group_custom = TRUE,
               line_color = c("red", "blue"))

Figura 19: Schoeller con colores personalizados por grupo (en las 20 primeras filas solo hay grupos A y B).

plot_schoeller(hc20, include_NO3 = TRUE)

Figura 20: Schoeller incluyendo el nitrato.

8.3 Huella química de cada muestra

Estos tres diagramas generan un panel por muestra. Sirven para comparar puntos concretos o para ponerlos sobre un mapa. Comparten los argumentos idcol (columna que identifica cada panel), nrow, ncol y facet_dir.

8.3.1 plot_stiff(): diagrama de Stiff

Concepto (Stiff, 1951). Cada muestra es un polígono: los cationes se representan a la izquierda de un eje vertical y los aniones a la derecha, en meq/L, en ejes horizontales a distintas alturas (en la disposición clásica, Na + K | Cl, Ca | HCO₃ + CO₃ y Mg | SO₄). La forma identifica la facies a simple vista y el tamaño, la mineralización. Por eso se usan mucho como símbolos en mapas hidroquímicos.

Con include_NO3 = "Cl" o "SO4" el nitrato se suma a ese anión. Con scales = "free_x" cada panel tiene su propia escala horizontal, algo útil si las mineralizaciones son muy distintas. Ojo: con escalas libres ya no se pueden comparar tamaños.

plot_stiff(hc20)

Figura 21: Stiff de las 20 primeras muestras (identificador automático).

plot_stiff(hc20, idcol = "id")

Figura 22: Stiff con identificador de muestra.

sel_data <- hc20[hc20$id %in% c("id_1", "id_5", "id_20"), ]
plot_stiff(sel_data)

Figura 23: Stiff de tres muestras seleccionadas.

plot_stiff(sel_data, group = "type")

Figura 24: Stiff de tres muestras coloreadas por grupo.

El orden de los paneles sigue el orden de aparición en los datos. Para controlarlo, basta con ordenar antes:

sel_data <- sel_data[order(sel_data$type, decreasing = TRUE), ]
plot_stiff(sel_data, idcol = "id", group = "type")

Figura 25: Stiff con orden de paneles controlado (grupo B primero).

plot_stiff(
  hc20,
  idcol = "id",
  group = "type",
  group_custom = TRUE,
  fill  = c("lightblue", "orange"),
  color = c("red", "red")
)

Figura 26: Stiff con relleno y contorno personalizados por grupo.

8.3.2 plot_collins(): diagrama de barras de Collins

Concepto (Collins, 1923). Cada muestra son dos columnas apiladas, una de cationes y otra de aniones, con la contribución de cada ion en % de meq/L. Si el análisis está equilibrado, las dos columnas suman lo mismo (100 %). La composición relativa se lee de forma inmediata. CO₃ se suma a HCO₃.

plot_collins(hc20, idcol = "id")

Figura 27: Diagramas de Collins de las 20 primeras muestras.

8.3.3 plot_stabler(): diagrama de Stabler

Concepto. Es la versión horizontal de la misma idea: para cada muestra, dos barras horizontales (cationes y aniones) divididas en segmentos proporcionales al % de meq/L de cada ion. Se usa sobre todo en informes como alternativa compacta al de Collins. El manual no lo menciona en la descripción del paquete, pero la función existe y está documentada.

plot_stabler(hc20, idcol = "id")

Figura 28: Diagramas de Stabler de las 20 primeras muestras.

8.4 Procesos que controlan la química

8.4.1 plot_gibbs(): diagrama de Gibbs

Concepto (Gibbs, 1970). Propuesto para aguas superficiales del mundo, relaciona el TDS (eje Y, logarítmico) con dos cocientes:

• izquierda: Na⁺ / (Na⁺ + Ca²⁺);
• derecha: Cl⁻ / (Cl⁻ + HCO₃⁻).

Las muestras se reparten en tres dominios, que el paquete marca con guías (show_guide) y, opcionalmente, con texto (show_text):

• Dominio de la precipitación (abajo a la derecha): TDS bajo y cocientes altos; aguas con química heredada de la lluvia.
• Dominio de la roca (centro a la izquierda): TDS intermedio y cocientes bajos; la química la controla la meteorización de los minerales (carbonatos, silicatos).
• Dominio de evaporación-cristalización (arriba a la derecha): TDS alto y cocientes altos; la evaporación concentra las sales y precipita la calcita, de modo que el agua se enriquece en Na y Cl.

Su aplicación a aguas subterráneas es muy frecuente, pero conviene recordar que el modelo original no contempla procesos como el intercambio iónico, la contaminación o la mezcla con aguas salinas.

plot_gibbs(hc_data)

Warning: CO3 is missing

TDS computed with: HCO3, Cl, NO3, SO4, Ca, Mg, Na, K

Figura 29: Gibbs básico.

plot_gibbs(hc_data, group = "type")

Warning: CO3 is missing

TDS computed with: HCO3, Cl, NO3, SO4, Ca, Mg, Na, K

Figura 30: Gibbs agrupado por type.

plot_gibbs(
  hc_data,
  group = "type",
  group_custom = TRUE,
  fill  = c("skyblue", "tomato", "green"),
  shape = c(21, 22, 23)
)

Warning: CO3 is missing

TDS computed with: HCO3, Cl, NO3, SO4, Ca, Mg, Na, K

Color vector length does not match number of groups. Using default 'black'.

Size vector length does not match number of groups. Using default 2.

Figura 31: Gibbs con relleno y forma personalizados.

plot_gibbs(hc_data, tds_method = "all_ions", show_text = TRUE)

Warning: CO3 is missing

TDS computed with: HCO3, Cl, NO3, SO4, Ca, Mg, Na, K, SiO2

Figura 32: Gibbs con el TDS calculado con todos los iones y etiquetas de los dominios.

Con data_qc la dispersión es mucho mayor y el diagrama resulta más ilustrativo:

plot_gibbs(data, group = "Hydro_cond", show_text = TRUE)

Ion 'NO3_N' converted to species form 'NO3'

TDS computed with: Ca, Mg, Na, K, Cl, SO4, HCO3, CO3, NO3

Figura 33: Gibbs de data_qc por condición hidráulica.

Nota

La documentación dice que los datos se convierten a meq/L, pero el código calcula los cocientes de Gibbs con concentraciones en mg/L, lo que es coherente con el trabajo original de Gibbs (que usó cocientes en peso). Si se comparan con cocientes en meq/L de otros trabajos, los valores no coincidirán.

8.4.2 plot_biplot(): relaciones iónicas

Concepto. Los diagramas binarios entre iones (o combinaciones de iones) en meq/L son la herramienta básica para identificar procesos geoquímicos, comparando la nube de puntos con la recta estequiométrica de cada reacción:

Relación (meq/L)	Recta de referencia	Proceso que sugiere
Na frente a Cl	1:1	Disolución de halita o aporte marino. Na > Cl: intercambio iónico directo o meteorización de silicatos; Na < Cl: intercambio inverso
Ca + Mg frente a HCO₃	1:1	Disolución de calcita/dolomita por CO₂
Ca frente a SO₄	1:1	Disolución de yeso/anhidrita
(Ca + Mg) − (HCO₃ + SO₄) frente a (Na + K) − Cl	pendiente −1	Intercambio catiónico
Ca frente a Mg	—	Relación calcita/dolomita

plot_biplot() admite expresiones aritméticas en la sintaxis "y ~ x" ("Ca + Mg ~ HCO3", "Na / Cl ~ SO4"). Cada elemento de pairs es un panel; el resultado es una rejilla de cowplot. Se le pueden pasar rectas de referencia (slope, intercept), escalas logarítmicas (log_x, log_y) y límites, ya sea comunes o por panel.

plot_biplot(data, pairs = c("Ca ~ Mg", "SO4 ~ Cl", "Cl + SO4 ~ NO3_N"))

Ion 'NO3_N' converted to species form 'NO3'

Warning: Removed 2 rows containing missing values or values outside the scale range
(`geom_point()`).

Warning: Removed 4 rows containing missing values or values outside the scale range
(`geom_point()`).

Figura 34: Biplots básicos.

plot_biplot(data, pairs = c("Ca ~ HCO3", "Na ~ Cl"), group = "Cluster")

Figura 35: Biplots agrupados por clúster.

plot_biplot(
  data,
  pairs = c("Ca ~ HCO3", "Na + K ~ Cl"),
  group = "Hydro_cond",
  group_custom = TRUE,
  fill  = c("steelblue", "tomato", "darkgreen", "yellow"),
  color = c("black", "black", "black", "black"),
  shape = c(21, 22, 24, 25),
  size  = c(2.5, 3, 3.5, 4),
  ncol  = 2
)

Figura 36: Biplots con estilos personalizados por condición hidráulica.

plot_biplot(data, pairs = c("Ca ~ Ca", "Na ~ Na"), slope = 1, intercept = 0)

Figura 37: Misma recta 1:1 en todos los paneles (ejemplo didáctico del manual: un ion frente a sí mismo cae siempre sobre la 1:1).

plot_biplot(
  data,
  pairs     = c("Ca ~ Ca", "Ca + Mg ~ HCO3", "Cl ~ Na"),
  slope     = c(NA, NA, 1),
  intercept = c(NA, NA, 0)
)

Figura 38: Recta de referencia solo en el tercer panel.

plot_biplot(data, pairs = c("Ca ~ Mg", "Cl ~ Na"), log_x = TRUE)

Figura 39: Eje X logarítmico en todos los paneles.

plot_biplot(
  data,
  pairs = c("Ca ~ Mg", "Ca + Mg ~ HCO3", "Cl ~ Na"),
  log_x = c(FALSE, TRUE, FALSE),
  log_y = c(FALSE, FALSE, TRUE)
)

Figura 40: Escala logarítmica solo en algunos paneles.

Un juego de biplots más “interpretativo”, con las rectas estequiométricas de las reacciones de la tabla:

plot_biplot(
  data,
  pairs = c("Na ~ Cl", "Ca + Mg ~ HCO3", "Ca ~ SO4"),
  slope = 1, intercept = 0,
  log_x = TRUE, log_y = TRUE,
  group = "Hydro_cond",
  ncol = 3
)

Warning: Removed 2 rows containing missing values or values outside the scale range
(`geom_point()`).

Figura 41: Relaciones iónicas de diagnóstico en data_qc, con rectas 1:1.

8.5 save_plot(): exportar los diagramas

Envoltorio de ggsave() que asigna un tamaño por defecto adecuado a cada tipo de diagrama (type: "piper", "ternary", "collins", "stabler", "durov", "gibbs", "schoeller", "stiff", "biplot"). Por defecto exporta a 300 ppp.

p <- plot_piper(hc_data)
save_plot(file.path(tempdir(), "piper.png"), type = "piper", plot = p)

p <- plot_stiff(hc_data[1:10, ], ncol = 5)
save_plot(file.path(tempdir(), "stiff.png"), type = "stiff", plot = p,
          width = 10, height = 4)

list.files(tempdir(), pattern = "png$")

[1] "piper.png" "stiff.png"

Como los diagramas son objetos ggplot, también se pueden retocar con la sintaxis habitual antes de guardarlos (títulos, temas…):

plot_gibbs(hc_data, group = "type") +
  labs(caption = "Fuente: hc_data (paquete hydrochem)")

Warning: CO3 is missing

TDS computed with: HCO3, Cl, NO3, SO4, Ca, Mg, Na, K

Figura 42: Un diagrama de Gibbs retocado con labs() de ggplot2.

9 Flujo de trabajo completo

Un esquema realista, de los datos brutos del laboratorio a la figura, en un único pipe:

qc_final <- data_qc |>
  replace_bdl(method = "DL/2") |>                 # 1. censurados
  as_tibble() |>
  mutate(
    eb     = ionic_balance(pick(everything())),   # 2. control de calidad
    facies = water_type(pick(everything()), format = "short"),
    TDS    = tds(pick(everything())),
    SAR    = sar(pick(everything())),
    CAI1   = cai1(pick(everything()))
  ) |>
  filter(abs(eb) <= 5)                            # 3. solo análisis válidos

qc_final |>
  summarise(n = n(),
            TDS_mediana = median(TDS),
            SAR_max = max(SAR),
            .by = Hydro_cond)

plot_piper(qc_final, group = "Hydro_cond", legend.title = "Condición hidráulica")

Warning: Removed 4 rows containing missing values or values outside the scale range
(`geom_point()`).

Figura 43: Piper de las muestras de data_qc con |EB| ≤ 5 %, agrupadas por condición hidráulica.

10 Observaciones tras revisar el código fuente

Es la versión 0.1.0: el paquete es muy útil, pero joven. Al contrastar el manual con el código he encontrado detalles que conviene conocer:

• Colores por grupo en R ≤ 4.3. Véase el aviso y el parche del apartado de instalación.
• Nombres de columna estrictos. Las funciones buscan exactamente Ca, HCO3, NO3_N… Una columna mal nombrada no siempre da error: en ionic_balance() o water_type() un minoritario mal nombrado simplemente no entra en el cálculo. Hay que revisar los mensajes "... computed with: ...".
• Una sola unidad para todos los iones. base_unit se aplica a todas las columnas. Si el laboratorio da unos parámetros en mg/L y otros en µg/L (metales traza), hay que homogeneizar antes.
• Gibbs en mg/L. La ayuda de plot_gibbs() dice que convierte a meq/L, pero los cocientes se calculan en mg/L (véase la nota en su apartado). Además, la ayuda describe tds_method = "major" como “solo Na, Ca, Cl, HCO₃”, pero en realidad llama a tds() y usa todos los mayoritarios.
• rsc() y CO₃. La ayuda dice que CO3 es obligatorio, pero el código solo emite un aviso y calcula con HCO₃ si falta.
• total_alkalinity(). El término de OH⁻ está mal planteado (véase su apartado), aunque es irrelevante en aguas naturales.
• plot_ternary(). La ayuda menciona meq/L como unidad destino, pero el argumento to_unit solo admite "mmol/L" y "mg/L".
• stat_summary() enmascara a ggplot2::stat_summary(). Conviene usar prefijo.
• Salidas data.table. convert_unit(), replace_bdl() y stat_summary() devuelven data.table. Con dplyr funcionan, pero se imprimen de otra forma; as_tibble() las homogeneiza.

11 Referencias

Paquete y datos

• Cinkus, G. y Bondu, R. (2026). hydrochem: Analysis and Visualization Tools for Hydrochemical Analysis. R package version 0.1.0. https://doi.org/10.32614/CRAN.package.hydrochem
• Cloutier, V. (2004). Origin and geochemical evolution of groundwater in the Paleozoic Basses-Laurentides sedimentary rock aquifer system, St. Lawrence Lowlands, Québec, Canada. Tesis doctoral, INRS-Eau, Terre & Environnement. https://espace.inrs.ca/id/eprint/1445/

Diagramas (citados en el propio paquete)

• Collins, W.D. (1923). Graphic representation of water analyses. Industrial & Engineering Chemistry, 15(4), 394. https://doi.org/10.1021/ie50160a030
• Durov, S.A. (1948). Natural waters and graphic representation of their composition. Doklady Akademii Nauk SSSR, 59, 87–90.
• Gibbs, R.J. (1970). Mechanisms controlling world water chemistry. Science, 170(3962), 1088–1090. https://doi.org/10.1126/science.170.3962.1088
• Piper, A.M. (1944). A graphic procedure in the geochemical interpretation of water-analyses. Eos, Transactions American Geophysical Union, 25(6), 914–928. https://doi.org/10.1029/TR025i006p00914
• Schoeller, H. (1935). Utilité de la notion des échanges de bases pour la comparaison des eaux souterraines. Bulletin de la Société Géologique de France.
• Stiff, H.A. (1951). The interpretation of chemical water analysis by means of patterns. Journal of Petroleum Technology, 3(10), 15–17. https://doi.org/10.2118/951376-G

Índices de riego (origen de los umbrales orientativos)

• Doneen, L.D. (1964). Notes on water quality in agriculture. Department of Water Science and Engineering, University of California, Davis.
• Eaton, F.M. (1950). Significance of carbonates in irrigation waters. Soil Science, 69(2), 123–133.
• Richards, L.A. (ed.) (1954). Diagnosis and improvement of saline and alkali soils. USDA Agriculture Handbook No. 60.
• Szabolcs, I. y Darab, C. (1964). The influence of irrigation water of high sodium carbonate content on soils. Proceedings of the 8th International Congress of ISSS, Transactions II, 803–812.
• Wilcox, L.V. (1955). Classification and use of irrigation waters. USDA Circular No. 969.

Importante

Los volúmenes y páginas de las referencias sin DOI (Durov, Eaton, Szabolcs y Darab, Doneen) y la referencia completa de Schoeller (1935) están citados de memoria. Las DOI proceden de la descripción oficial del paquete en CRAN. Verificadlos antes de citarlos en un documento oficial.

---

R4.6.0 | hydrochem 0.1.0 | ggplot2 4.0.3 | data.table 1.18.4