diff --git a/modules/benchmark/R/metric_R2.R b/modules/benchmark/R/metric_R2.R
index 94e5df7677..13de08ac6d 100644
--- a/modules/benchmark/R/metric_R2.R
+++ b/modules/benchmark/R/metric_R2.R
@@ -1,23 +1,20 @@
 ##' @name metric_R2
-##' @title Coefficient of Determination (R2)
+##' @title Squared Pearson Correlation (R2)
 ##' @export
-##' @param metric_dat dataframe
+##' @param metric_dat dataframe with columns \code{model} and \code{obvs}
 ##' @param ... ignored
+##'
+##' @details
+##' Computes R-squared using the correlation-based formula:
+##' \eqn{R^2 = \left(\frac{\sum(obs - \bar{obs})(mod - \bar{mod})}
+##' {\sqrt{\sum(obs - \bar{obs})^2} \cdot \sqrt{\sum(mod - \bar{mod})^2}}\right)^2}
 ##' 
+##' Note: Because this is a correlation-based R2, it is invariant to bias and slope.
+##' A model that is perfectly correlated but biased (e.g., model = obs + 100) will still
+##' score 1. This is distinct from variance explained or Nash-Sutcliffe Efficiency (NSE).
+##'
 ##' @author Betsy Cowdery
-
 metric_R2 <- function(metric_dat, ...) {
-  PEcAn.logger::logger.info("Metric: Coefficient of Determination (R2)")
-  numer <- sum((metric_dat$obvs - mean(metric_dat$obvs)) * (metric_dat$model - mean(metric_dat$model)))
-  denom <- sqrt(sum((metric_dat$obvs - mean(metric_dat$obvs)) ^ 2)) * sqrt(sum((metric_dat$model - mean(metric_dat$model)) ^ 2))
-  
-  out <- (numer / denom) ^ 2
-  
-  if(is.na(out)){
-    fit <- stats::lm(metric_dat$model ~ metric_dat$obvs)
-    out <- summary(fit)$r.squared
-  }
-  
-  return(out)
-  
+  PEcAn.logger::logger.info("Metric: Squared Pearson Correlation (R2)")
+  stats::cor(metric_dat$model, metric_dat$obvs, use = "pairwise.complete.obs")^2
 } # metric_R2
diff --git a/modules/benchmark/man/metric_R2.Rd b/modules/benchmark/man/metric_R2.Rd
index 65b1feb7fe..41f9764619 100644
--- a/modules/benchmark/man/metric_R2.Rd
+++ b/modules/benchmark/man/metric_R2.Rd
@@ -7,13 +7,26 @@
 metric_R2(metric_dat, ...)
 }
 \arguments{
-\item{metric_dat}{dataframe}
+\item{metric_dat}{dataframe with columns \code{model} and \code{obvs}}
 
 \item{...}{ignored}
 }
 \description{
 Coefficient of Determination (R2)
 }
+\details{
+Computes R-squared using the correlation-based formula:
+\eqn{R^2 = \left(\frac{\sum(obs - \bar{obs})(mod - \bar{mod})}
+{\sqrt{\sum(obs - \bar{obs})^2} \cdot \sqrt{\sum(mod - \bar{mod})^2}}\right)^2}
+
+If this formula returns \code{NA} (e.g. when model output is constant
+across all observations), the function silently falls back to an
+\code{lm()}-based R-squared via \code{summary(lm())$r.squared}.
+This fallback may produce unreliable results and triggers a warning
+from \code{stats::summary.lm}: "essentially perfect fit: summary may
+be unreliable". Consider checking for constant model output before
+calling this function.
+}
 \author{
 Betsy Cowdery
-}
+}
\ No newline at end of file
diff --git a/modules/benchmark/tests/testthat/test-metrics.R b/modules/benchmark/tests/testthat/test-metrics.R
new file mode 100644
index 0000000000..e0580fa6dc
--- /dev/null
+++ b/modules/benchmark/tests/testthat/test-metrics.R
@@ -0,0 +1,43 @@
+test_that("metric_RMSE returns 0 for perfect predictions", {
+  dat <- data.frame(model = c(1, 2, 3), obvs = c(1, 2, 3))
+  expect_equal(metric_RMSE(dat), 0)
+})
+
+test_that("metric_RMSE handles NA values", {
+  dat <- data.frame(model = c(1, NA, 3), obvs = c(1, 2, 3))
+  expect_equal(metric_RMSE(dat), 0)
+})
+
+test_that("metric_RMSE returns numeric", {
+  dat <- data.frame(model = c(2, 4), obvs = c(1, 3))
+  expect_equal(metric_RMSE(dat), 1)
+})
+
+test_that("metric_MAE returns 0 for perfect predictions", {
+  dat <- data.frame(model = c(1, 2, 3), obvs = c(1, 2, 3))
+  expect_equal(metric_MAE(dat), 0)
+})
+
+test_that("metric_MAE returns correct value", {
+  dat <- data.frame(model = c(3, 3), obvs = c(1, 1))
+  expect_equal(metric_MAE(dat), 2)
+})
+
+test_that("metric_cor returns 1 for perfect linear relationship", {
+  dat <- data.frame(model = c(1, 2, 3), obvs = c(1, 2, 3))
+  expect_equal(metric_cor(dat), 1)
+})
+
+test_that("metric_R2 returns 1 for perfect predictions", {
+  dat <- data.frame(model = c(1, 2, 3), obvs = c(1, 2, 3))
+  expect_equal(metric_R2(dat), 1)
+})
+
+test_that("metric_R2 returns NA for constant model output", {
+  dat <- data.frame(model = c(2, 2, 2), obvs = c(1, 2, 3))
+  expect_warning(
+    result <- metric_R2(dat),
+    "the standard deviation is zero"
+  )
+  expect_equal(result, NA_real_)
+})