GHCNd Station Drought Monitoring Pipeline: Methods

Zachary H. Hoylman — Montana Climate Office, University of Montana

Updated: Version 1.1

Abstract The Montana Climate Office (MCO) GHCNd Station Drought Monitoring Pipeline is an operational system that computes standardized drought indices on a daily basis from individual weather station observations across the contiguous United States (CONUS). The pipeline ingests daily station data from the Global Historical Climatology Network – Daily (GHCNd) [1], a quality-controlled archive of in situ observations from over 100,000 stations worldwide, updated daily by NOAA’s National Centers for Environmental Information (NCEI).

From the core GHCNd variables — precipitation (PRCP), maximum temperature (TMAX), and minimum temperature (TMIN) — the pipeline derives three standardized drought index families: the Standardized Precipitation Index (SPI), the Standardized Precipitation-Evapotranspiration Index (SPEI), and the Evaporative Demand Drought Index (EDDI). Reference evapotranspiration (ET0) is estimated using the Hargreaves-Samani method [2] with extraterrestrial radiation computed from FAO-56 equations [3]. Ancillary products include precipitation accumulations (mm), percent of normal, departure from normal, and empirical percentile ranks. Each metric is computed across multiple aggregation timescales ranging from 15 to 730 days, as well as water-year and year-to-date accumulations.

Unlike the companion MCO CONUS Drought Monitoring Pipeline, which operates on the ~4 km gridMET dataset [14], this pipeline produces point-based station estimates. Station-level drought indices are complementary to the gridded products: they reflect actual observed conditions rather than interpolated surfaces, and they provide quasi-independent validation points for gridded analyses (many GHCNd stations also inform the surface-creation process for gridded products such as PRISM, so the two are not fully independent). To keep the two pipelines directly comparable, the station pipeline’s period of record begins in 1979 (START_YEAR=1979), deliberately matching the gridMET period of record so station-based and gridded indices are computed over equivalent historical baselines. Only stations that are currently reporting and have at least 30 years of historical data are included, ensuring robust climatological fitting for operational drought assessment.

2. Data

2.1 GHCNd Station Data

All meteorological inputs are obtained from the Global Historical Climatology Network – Daily (GHCNd) [1], an integrated database of daily climate summaries from land surface stations maintained by NOAA NCEI. GHCNd is subject to a suite of automated quality assurance procedures including duplicate checks, climatological outlier detection, spatial consistency tests, and temporal gap checks [4]. Observations that fail any quality check are flagged and excluded from this pipeline.

The following GHCNd elements are used:

Variable GHCNd element Native units Converted units Used for
PrecipitationPRCPtenths of mmmmSPI, SPEI, precip accumulations
Max temperatureTMAXtenths of °C°CHargreaves-Samani ET0
Min temperatureTMINtenths of °C°CHargreaves-Samani ET0

Data Access

GHCNd data are obtained from the NCEI public archive as by-year CSV files (YYYY.csv.gz) from ncei.noaa.gov/pub/data/ghcn/daily/by_year/. Each file contains all station observations for a single calendar year (~170 MB compressed for recent years). The pipeline downloads the full historical archive on first run (cold start) and subsequently re-downloads only the current-year and prior-year files on each operational run (warm start) — the prior year is refreshed only when the local copy is more than 7 days old — since older historical years are effectively immutable.

Quality Control

Each GHCNd observation carries a quality flag (QFLAG) from NCEI’s automated QC suite. Any observation with a non-blank quality flag is excluded prior to analysis. This removes values that failed duplicate checks, climatological range tests, spatial consistency checks, or other QC procedures described in Durre et al. (2010) [4].

2.2 Station Selection Criteria

The pipeline targets operational drought assessment and therefore restricts analysis to stations that are both currently reporting and have sufficient historical depth for robust climatological fitting. Station eligibility is determined from the ghcnd-inventory.txt metadata file, which records the first and last year of data for each station-element combination.

Filtering Criteria

  • CONUS bounding box — Only stations with US country prefix and coordinates within 24.5–49.5°N, 125–66.5°W are included, excluding Alaska, Hawaii, and US territories.
  • Recency — The station must have reported PRCP within the last calendar year (last year in inventory ≥ current year − 1).
  • Record length — The station must have at least 30 years of PRCP data (last year − first year + 1 ≥ 30).
  • SPEI/EDDI eligibility — Stations additionally require TMAX and TMIN with the same recency and record length thresholds. Approximately 75–80% of PRCP-eligible stations also report temperature.
  • Reporting latency — At computation time, stations whose most recent observation is older than the configured latency threshold (MAX_REPORTING_LATENCY, operational default: 3 days) are excluded from drought index computation for that run.

Station freshness is therefore enforced at three levels: (1) an inventory-level recency screen (reported within the last calendar year), (2) a run-time reporting-latency screen (most recent observation within 3 days of the pipeline run), and (3) a display-level filter in the D³ dashboard, which additionally hides stations whose most recent observation lags the currently displayed raster’s valid date by more than 4 days (GHCN_MAX_LAG_DAYS), so that lagged station accumulations are never compared against a fresher gridded field.

Current station counts (CONUS, April 2026): Approximately 5,277 stations are eligible for SPI (precipitation only), and approximately 4,096 are additionally eligible for SPEI and EDDI (precipitation + temperature).

3. Climatological Reference Periods

Scientific Motivation

Standardized drought indices measure anomalies relative to a baseline climatology. The choice of baseline profoundly affects index values: for example, the same temperature, reference ET, or precipitation value could indicate moderate drought against a 1950–1980 baseline but near-normal conditions against a 1991–2020 WMO normal. In other words, as the climate warms, the statistical properties of temperature, PET, and water balance may trend systematically in any of the distribution’s moments (mean, standard deviation, skewness, kurtosis, etc.), violating the stationarity assumption embedded in fixed-period normals. Hoylman et al. (2022) [5] demonstrate that non-stationarity in baseline climate substantially alters drought frequency and severity assessments, motivating the use of flexible, contemporary reference periods for operational monitoring.

Supported Modes

Spec syntax Mode Description
rolling:30 Rolling 30-year window Rolling 30-year window ending at the current year; tracks contemporary climatology (default). rolling:N generalizes to any window length.
fixed:1991:2020 Fixed interval WMO 1991–2020 standard normal
full Full period of record All years from START_YEAR to present; maximizes sample size

For each timescale and reference period combination, a vector of one aggregated value per year is assembled by extracting the rolling window sum (or mean) ending on the same calendar month-day as the most recent observation. This vector serves as input to the distribution fitting procedure for SPI, SPEI, and EDDI.

Minimum data requirement: At least 10 valid years must be present in the climatological vector for a distribution fit to proceed — matching the 10-anchor-year requirement of the companion gridded CONUS pipeline. For the Gamma fit (SPI), at least 3 positive (non-zero) values with non-zero variance are additionally required. Stations or timescales that do not meet these thresholds receive NA.

This floor is not a relaxation of the 30-year station screen in Section 2.2. The 30-year criterion is a station-level test on inventory metadata (the total span of years a station has reported an element). The fitting sample for a specific timescale × reference period combination can be much smaller than the station’s record span: with MIN_OBS_FRACTION=1 (the operational setting), any year with even a single missing day in the aggregation window is excluded from the sample; a fixed:1991:2020 period caps the sample at 30 years regardless of record length; and long windows (e.g. 730 days) consume years at the start of the record. The 10-valid-year floor guards against poorly constrained fits for a specific window at stations with gappy records — it is not the expected sample size, which for most stations is close to the full reference period.

4. Reference Evapotranspiration (ET0) Estimation

SPEI and EDDI require an estimate of atmospheric evaporative demand. GHCNd provides only temperature data (TMAX, TMIN) — wind speed, humidity, and solar radiation are not available from the standard GHCNd archive. The Penman-Monteith method is therefore not feasible. The Thornthwaite method was rejected because it operates only at monthly timesteps and performs poorly in arid, semiarid, and high-elevation climates [6].

4.1 Hargreaves-Samani Method

Daily reference evapotranspiration is estimated using the Hargreaves-Samani (HS) equation [2], which requires only daily maximum and minimum temperature plus extraterrestrial radiation (computed from latitude and day of year):

ET0 = 0.0023 ⋅ Ra* ⋅ (Tmean + 17.8) ⋅ (Tmax − Tmin)0.5

where:

  • ET0 is daily reference evapotranspiration (mm/day)
  • Ra* is extraterrestrial radiation already converted to equivalent evaporation depth (mm/day). That is, Ra* = 0.408 ⋅ Ra where Ra is in MJ m−2 day−1 (FAO-56 Eq. 20; see Section 4.2). The 0.408 factor is applied once when computing Ra* and is not reapplied here.
  • Tmean = (Tmax + Tmin) / 2 (°C)
  • Tmax − Tmin is the diurnal temperature range (°C), clamped ≥ 0
  • ET0 is clamped ≥ 0 (negative values are not physically meaningful)

The HS method has been shown to provide reasonable estimates of ET0 across a wide range of climates when compared to Penman-Monteith, particularly for drought monitoring applications where relative anomalies matter more than absolute accuracy [7].

Limitation: the use of Hargreaves-Samani is a constraint imposed by the variables available in GHCNd, and should be treated as a limitation of the station products. HS is a temperature-based proxy: it does not capture day-to-day variability in wind speed, humidity, or solar radiation. Systematic differences between HS-based ET0 and Penman-Monteith-based ET0 — and therefore between this pipeline’s SPEI and EDDI and those of the companion gridded CONUS pipeline, which uses gridMET’s Penman-Monteith reference ET — should be expected, particularly under unusually windy, humid, or cloudy conditions.

4.2 Extraterrestrial Radiation (FAO-56)

Extraterrestrial radiation (Ra) is computed from station latitude and day of year using the standard FAO-56 equations [3]:

  1. Solar declination (Eq. 24):
    δ = 0.409 ⋅ sin(2π/365 ⋅ J − 1.39)
    where J is the day of year (1–366).
  2. Inverse relative Earth-Sun distance (Eq. 23):
    dr = 1 + 0.033 ⋅ cos(2π/365 ⋅ J)
  3. Sunset hour angle (Eq. 25):
    ωs = arccos[ −tan(φ) ⋅ tan(δ) ]
    where φ is station latitude in radians. For polar conditions where |tan(φ) ⋅ tan(δ)| ≥ 1, ωs is set to π (24-hour daylight) or 0 (no daylight).
  4. Extraterrestrial radiation (Eq. 21):
    Ra = (24 ⋅ 60 / π) ⋅ Gsc ⋅ dr ⋅ [ ωs sin(φ) sin(δ) + cos(φ) cos(δ) sin(ωs) ]
    where Gsc = 0.0820 MJ m−2 min−1 is the solar constant.
  5. Conversion to equivalent evaporation (Eq. 20):
    Ra (mm/day) = 0.408 ⋅ Ra (MJ m−2 day−1)
    where 0.408 = 1/λ and λ = 2.45 MJ kg−1 is the latent heat of vaporization.

5. Drought Index Methods

All indices are computed independently at each station. For a given aggregation window of W days, the pipeline computes two distinct quantities:

  1. The current observation, XW(now) — the rolling window of W days ending at the most recent valid observation date. This value is independent of the climatological reference period.
  2. The reference distribution, {XW(y)} — one aggregated value per reference year y, anchored on the same calendar month-day as the current observation. The set of years y is determined by CLIM_PERIODS.

A probability distribution (or empirical CDF) is then fitted to the reference distribution, and the current observation is evaluated against that fit to produce a standardized anomaly. For rolling:N and full modes the current year is naturally part of the reference distribution, so the current observation also appears as the last element of the reference set. For fixed:YYYY:YYYY modes where the current year is outside the fixed range, the current observation is external to the reference set and is evaluated against (but not included in) the fitted distribution. Treating the current value and the reference distribution as separate inputs ensures correct behavior in all modes.

Aggregation Timescales

5.1 Standardized Precipitation Index (SPI)

Motivation

The SPI [8] expresses precipitation anomalies in units of standard deviations relative to a fitted climatological distribution, making it comparable across locations and timescales with differing precipitation regimes. It requires only precipitation data, making it computable for all stations in the network.

Input

Daily precipitation (PRCP), in millimeters. Rolling window sums.

Statistical Method

Precipitation is inherently bounded at zero and right-skewed; a two-parameter Gamma distribution (shape α, scale β) is therefore fitted to the climatological sample for each station and aggregation window. Parameter estimation follows the L-moment / probability-weighted moment (PWM) approach of Hosking (1990) [9], as implemented in the R package lmomco:

  1. Zero handling (Stagge et al., 2015) — Following the methodology of Stagge et al. (2015) [10], zero (or below-threshold) precipitation values are handled via a mixed distribution approach using the Weibull plotting position. The probability of zero precipitation is estimated as p0 = m / (n + 1), where m is the number of zero values and n is the total sample size. Zero-valued observations are assigned the center of probability mass below the threshold:
    p̄0 = (m + 1) / [ 2(n + 1) ]
    Non-zero values are mapped to the upper portion of the probability space:
    H(x) = p0 + (1 − p0) ⋅ G(x; α, β)
    where G is the Gamma CDF fitted only to positive values. This approach correctly represents the discrete mass at zero without the artifacts introduced by replacing zeros with small positive constants.
  2. PWM estimation — Unbiased probability-weighted moments β0, β1 are computed from the positive (non-zero) values in the climatological sample.
  3. L-moment conversion — PWMs are converted to L-moments (λ1, λ2, τ3).
  4. Parameter estimation — lmomco::pargam() returns (α, β) from L-moments of the positive values only.
  5. CDF evaluation — For non-zero observations, the mixed CDF H(x) is evaluated to obtain a non-exceedance probability p. For zero-valued observations, p = p̄0 (center of probability mass).
  6. Normal quantile transform — SPI = Φ−1(p), where Φ−1 is the standard normal quantile function. No clamping is applied to the resulting SPI values; in the rare case that the fitted CDF saturates numerically (p = 0 or 1), the resulting non-finite value is set to NA.
SPI = Φ−1[ H( XW ) ]   (exact zeros take p = p̄0)

Interpretation

Category Description Percentile Range SPI / SPEI Values
Normal or wet conditionsx > 30x > −0.524
D0Abnormally Dry30 ≥ x > 20−0.524 ≥ x > −0.842
D1Moderate Drought20 ≥ x > 10−0.842 ≥ x > −1.282
D2Severe Drought10 ≥ x > 5−1.282 ≥ x > −1.645
D3Extreme Drought5 ≥ x > 2−1.645 ≥ x > −2.054
D4Exceptional Droughtx < 2x < −2.054

5.2 Ancillary Precipitation Metrics

Metric Formula Units Description
Percent of Normal (PON) ( XW / μclim ) × 100 % Current accumulation as percentage of climatological mean
Departure XW − μclim mm Absolute anomaly relative to climatological mean
Percentile F̂n( XW ) 0–1 Empirical ECDF rank within climatological sample
Accumulation XW mm Raw rolling precipitation sum

5.3 Standardized Precipitation-Evapotranspiration Index (SPEI)

Motivation

The SPEI [11] extends SPI by incorporating atmospheric evaporative demand, making it sensitive to warming-driven increases in ET0 even when precipitation is unchanged. In a warming climate, SPEI captures atmospheric-demand-driven drought stress that SPI alone would miss.

Input

Daily water balance = precipitation (PRCP) − reference evapotranspiration (ET0 from Hargreaves-Samani), in mm/day. Negative values indicate a moisture deficit; positive values indicate a moisture surplus.

Aggregation

Rolling window sums of the daily water balance over the same timescales as SPI; equivalently, the rolling sum of precipitation minus the rolling sum of ET0, both accumulated over the same window. Units are millimeters of accumulated moisture surplus (positive) or deficit (negative). The Generalized Logistic (GLO) distribution is then fitted to this aggregated quantity for each calendar day-of-year within the climatological reference period, and the current observation is transformed via the same fitted distribution.

Statistical Method

The water balance variable is unbounded below (large ET0 on dry days) and has heavier distributional tails than precipitation alone. The GLO distribution, which accommodates both heavy tails and asymmetry, is therefore used in place of the Gamma:

  1. L-moments estimated from the climatological sample via unbiased PWMs.
  2. lmomco::parglo() fits the three-parameter GLO (ξ, α, k).
  3. GLO CDF evaluated at current observation → probability p.
  4. SPEI = Φ−1(p).
SPEI = Φ−1[ FGLO(ξ,α,k)( P − ET0 )W ]

Interpretation

SPEI values follow the same classification thresholds as SPI. Because it incorporates evaporative demand, SPEI typically shows stronger drying trends over recent decades than SPI alone, consistent with the non-stationarity findings of [5]. Only stations with both PRCP and temperature data (TMAX + TMIN) can produce SPEI, as Hargreaves-Samani ET0 is required for the water balance computation.

5.4 Evaporative Demand Drought Index (EDDI)

Motivation

EDDI [12] isolates the atmospheric demand component of the hydrological cycle. Because it is derived solely from ET0 — without precipitation — it captures early warning signals of drought stress driven by high temperature and diurnal range before soil moisture depletion has occurred.

Input

Hargreaves-Samani ET0 (mm/day). Rolling window sums (cumulative evaporative demand).

Statistical Method

EDDI uses a nonparametric rank-based approach following Hobbins et al. (2016) [12]:

  1. Rank — Observations within the climatological sample are ranked from highest (rank 1) to lowest. A rank of 1 therefore corresponds to maximum evaporative demand (drought stress).
  2. Tukey plotting position — Following Hobbins et al. (2016) [12], the empirical probability is estimated via the Tukey plotting position:
    pi = (i − 0.33) / (n + 0.33)
    where i is the rank of the aggregated ET0 in the historical time series (i = 1 for maximum ET0) and n is the number of observations being ranked. When the reference period is rolling or full-record, the current observation is the final member of the climatological sample (n = sample size); under a fixed reference period that excludes the current year, the current observation is appended to the sample and ranked within the resulting n + 1 values. Because rank 1 = highest demand, p is small for drought conditions.
  3. Inverse normal approximation — EDDI is derived via the rational approximation of Abramowitz and Stegun (1965) [13], following the formulation in Hobbins et al. (2016) [12] and Vicente-Serrano et al. (2010) [15]. For p ≤ 0.5, W = √(−2⋅ln(p)); for p > 0.5, p is replaced with 1−p and the sign of EDDI is reversed:
    EDDI = W − ( C0 + C1W + C2W2 ) / ( 1 + d1W + d2W2 + d3W3 )
    with constants C0 = 2.515517, C1 = 0.802853, C2 = 0.010328, d1 = 1.432788, d2 = 0.189269, and d3 = 0.001308, as tabulated in Hobbins et al. (2016).
  4. Sign convention — High ET0 (low p) produces positive EDDI values, indicating drier-than-normal conditions. Negative values indicate wet anomalies. A zero EDDI indicates that accumulated ET0 equals the climatological median.
EDDI = −Φ−1[ pTukey( ET0,W ) ]   (rank 1 = highest ET0; the leading minus enforces positive EDDI for drought)

Interpretation

Positive EDDI indicates above-normal atmospheric demand (drought precursor or ongoing drought amplifier). EDDI provides a pure atmospheric signal unconfounded by antecedent precipitation and is particularly valuable for flash drought early warning.

6. Processing Architecture

Station-Based Parallelism

Unlike the companion gridded pipeline, which decomposes the CONUS domain into spatial tiles, this pipeline parallelizes across stations. Each station’s time series is processed independently via pbmcapply::pbmclapply() (fork-based parallel apply). This is naturally embarrassingly parallel — stations share no state.

Cold Start vs. Warm Start

On first run (cold start), all 48 yearly CSV files (START_YEAR=1979 through the current year — the start year deliberately matches the gridMET period of record used by the companion gridded pipeline) are downloaded from NCEI and parsed into per-station RDS files. Each station’s RDS contains a complete daily time series with columns for date, PRCP, TMAX, and TMIN. On subsequent operational runs (warm start), only the current year and prior year CSV files are re-downloaded and merged into existing per-station RDS files, reducing I/O by approximately 200×.

S3 Caching (AWS Fargate)

When deployed on AWS Fargate, all per-station RDS files, raw CSVs, and derived outputs are synced to an S3 bucket at the end of each run. The next Fargate run restores this cache before processing, avoiding a full cold start on ephemeral compute infrastructure. Derived outputs are published twice: to a date-stamped archive (derived/ghcnd_drought/{YYYY-MM-DD}/) and to an operational latest/ copy that always reflects the most recent run. Intermediate syncs after each pipeline step keep latest/ current, so a failed run can resume without recomputing completed steps. Local development runs bypass S3 sync entirely.

Pipeline Steps

Step 0: Filter active CONUS stations from GHCNd inventory
Step 1: Download/parse GHCNd by-year CSVs → per-station RDS files
Step 2: Compute Hargreaves-Samani daily ET₀ for SPEI-eligible stations
Step 3: Compute SPI (all precip stations)
Step 4: Compute SPEI (stations with PET)
Step 5: Compute EDDI (stations with PET)
Step 6: Compute precipitation accumulations + ancillary metrics
Step 7: Assemble outputs (GeoJSON, CSV, per-station JSON, manifest)

7. Output Products

File Format Description
GHCNd_drought_current.geojson.gz GeoJSON (gzipped) All stations with all drought indices as Feature properties; ~2 MB compressed
all_stations.csv CSV Wide-format table: one row per station, all metrics as columns
stations/{ID}.json JSON Per-station detail files with all indices, accumulations, and metadata
station_catalog.csv CSV Station metadata, record dates, and available index flags
manifest.csv CSV Per-dataset date manifest for tracking data freshness
stations_base.geojson.gz GeoJSON (gzipped) Lazy-load base layer: station geometry and metadata only (id, name, state, elevation, data date)
slices/{var}_{ts}_{period}.json.gz JSON (gzipped) Lazy-load value slices: one small file per variable × timescale × reference period, mapping station id → value
manifest.json.gz JSON (gzipped) Index of all available slices; clients load this first, then fetch only the slices they display

The base + slices + manifest layout allows web clients to render the station map without downloading the full GHCNd_drought_current.geojson.gz: the base geometry loads once, and individual index/timescale/period combinations are fetched on demand.

Index Summary

Index Stations Units / Range Distribution
SPIAll with PRCP (~5,277)Dimensionless z-scoreGamma (L-moments)
SPEIPRCP + TMAX + TMIN (~4,096)Dimensionless z-scoreGeneralized Logistic (L-moments)
EDDITMAX + TMIN (~4,096)Dimensionless z-scoreNonparametric (Tukey)
Precip accumulationAll with PRCPmmRaw sum
Percent of NormalAll with PRCP%Climatological mean
DepartureAll with PRCPmmClimatological mean
PercentileAll with PRCP0–1Empirical ECDF

8. Software & Reproducibility

Runtime Environment

The pipeline runs inside a Docker container based on rocker/r-ver:4.4.1 (R 4.4.1). The container includes the AWS CLI for S3 synchronization and minimal spatial libraries (GDAL/GEOS/PROJ for GeoJSON output via sf).

Key R Packages

PackageRole
data.tableFast CSV parsing and in-memory data manipulation
lmomcoL-moment parameter estimation; Gamma, GLO distributions
sfGeoJSON output generation
jsonlitePer-station JSON and compact GeoJSON output
pbmcapplyFork-based parallel apply with progress bar
purrrFunctional programming utilities

Source Code Structure

mco-GHCNd/
├── R/
│   ├── drought-functions.R        # Hargreaves-Samani, SPI, SPEI, EDDI fitting
│   ├── pipeline-common.R          # Shared infrastructure
│   ├── 0_ghcnd-station-filter.R   # Station selection
│   ├── 1_ghcnd-cache.R            # Data download and parsing
│   ├── 2_compute-pet.R            # Hargreaves-Samani ET₀
│   ├── 3_metrics-spi.R            # SPI computation
│   ├── 4_metrics-spei.R           # SPEI computation
│   ├── 5_metrics-eddi.R           # EDDI computation
│   ├── 6_metrics-precip-accum.R   # Precipitation accumulations
│   └── 7_output-assembly.R        # Output generation
├── pipeline/
│   ├── run_once.sh                # Orchestration script
│   └── run_test.sh                # Test runner
├── Dockerfile
├── docker-compose.yml
└── terraform/                     # AWS infrastructure (Fargate, S3, ECR)

Deployment

The pipeline is deployed on AWS Fargate via Amazon ECS, triggered nightly at 10:00 PM Mountain Time by an EventBridge Scheduler rule. Outputs are written to a public S3 bucket with CORS enabled for direct web consumption. Infrastructure is managed via Terraform.

9. References

  1. Menne MJ, Durre I, Vose RS, Gleason BE, Houston TG. An overview of the Global Historical Climatology Network-Daily database. Journal of Atmospheric and Oceanic Technology. 2012;29(7):897–910. doi:10.1175/JTECH-D-11-00103.1
  2. Hargreaves GH, Samani ZA. Reference crop evapotranspiration from temperature. Applied Engineering in Agriculture. 1985;1(2):96–99. doi:10.13031/2013.26773
  3. Allen RG, Pereira LS, Raes D, Smith M. Crop evapotranspiration: guidelines for computing crop water requirements. FAO Irrigation and Drainage Paper 56. Rome: Food and Agriculture Organization of the United Nations; 1998.
  4. Durre I, Menne MJ, Gleason BE, Houston TG, Vose RS. Comprehensive automated quality assurance of daily surface observations. Journal of Applied Meteorology and Climatology. 2010;49(8):1615–1633. doi:10.1175/2010JAMC2375.1
  5. Hoylman ZH, Bocinsky RK, Jencso KG. Drought assessment has been outpaced by climate change: empirical arguments for a paradigm shift. Nature Communications. 2022;13:2715. doi:10.1038/s41467-022-30316-5
  6. Droogers P, Allen RG. Estimating reference evapotranspiration under inaccurate data conditions. Irrigation and Drainage Systems. 2002;16(1):33–45. doi:10.1023/A:1015508322413
  7. Begueria S, Vicente-Serrano SM, Angulo-Martinez M. A multiscalar global drought dataset: the SPEIbase. Bulletin of the American Meteorological Society. 2010;91(10):1351–1356. doi:10.1175/2010BAMS2988.1
  8. McKee TB, Doesken NJ, Kleist J. The relationship of drought frequency and duration to time scales. Proceedings of the 8th Conference on Applied Climatology. 1993;17(22):179–183. American Meteorological Society.
  9. Hosking JRM. L-moments: analysis and estimation of distributions using linear combinations of order statistics. Journal of the Royal Statistical Society, Series B. 1990;52(1):105–124. doi:10.1111/j.2517-6161.1990.tb01775.x
  10. Stagge JH, Tallaksen LM, Gudmundsson L, Van Loon AF, Stahl K. Candidate distributions for climatological drought indices (SPI and SPEI). International Journal of Climatology. 2015;35(13):4027–4040. doi:10.1002/joc.4267
  11. Begueria S, Vicente-Serrano SM, Reig F, Latorre B. Standardized precipitation evapotranspiration index (SPEI) revisited: parameter fitting, evapotranspiration models, tools, datasets and related applications. International Journal of Climatology. 2014;34(10):3001–3023. doi:10.1002/joc.3887
  12. Hobbins MT, Wood A, McEvoy DJ, Huntington JL, Morton C, Anderson M, Hain C. The Evaporative Demand Drought Index. Part I: Linking drought evolution to variations in evaporative demand. Journal of Hydrometeorology. 2016;17(6):1745–1761. doi:10.1175/JHM-D-15-0121.1
  13. Abramowitz M, Stegun IA, eds. Handbook of Mathematical Functions with Formulas, Graphs, and Mathematical Tables. National Bureau of Standards Applied Mathematics Series 55. Washington, DC: U.S. Government Printing Office; 1965.
  14. Abatzoglou JT. Development of gridded surface meteorological data for ecological applications and modelling. International Journal of Climatology. 2013;33(1):121–131. doi:10.1002/joc.3413
  15. Vicente-Serrano SM, Begueria S, Lopez-Moreno JI. A multiscalar drought index sensitive to global warming: the Standardized Precipitation Evapotranspiration Index. Journal of Climate. 2010;23(7):1696–1718. doi:10.1175/2009JCLI2909.1