diff --git a/modules/emulator/NEWS.md b/modules/emulator/NEWS.md
index 0fd4703c018..19ca3675375 100644
--- a/modules/emulator/NEWS.md
+++ b/modules/emulator/NEWS.md
@@ -2,6 +2,10 @@
 
 * Added keywords and bug reporting URL to DESCRIPTION. No other files changed.
 
+## Fixed
+
+* `calcSpatialCov.matrix()`: removed a redundant nested loop that recomputed each matrix row `n` times unnecessarily, reducing complexity from O(n^3) to O(n^2) (#3964).
+
 
 
 # PEcAn.emulator 1.8.1
diff --git a/modules/emulator/R/calcSpatialCov.matrix.R b/modules/emulator/R/calcSpatialCov.matrix.R
index 7393ff98b9e..224073b76c2 100644
--- a/modules/emulator/R/calcSpatialCov.matrix.R
+++ b/modules/emulator/R/calcSpatialCov.matrix.R
@@ -10,13 +10,5 @@
 ##' @author Michael Dietze
 ##' @export
 calcSpatialCov.matrix <- function(d, psi, tau, ...) {
-  nl <- nrow(d)
-  H  <- matrix(0, nl, nl)
-  for (i in seq_len(nl)) {
-    # for(j in 1:nl){ H[i,j] <- tau*exp(-psi*d[i,j]) }
-    for (j in seq_len(nl)) {
-      H[i, ] <- tau * exp(-psi * d[i, ])
-    }
-  }
-  return(H)
+  tau * exp(-psi * d)
 }
diff --git a/modules/emulator/tests/testthat/test-calcSpatialCov.R b/modules/emulator/tests/testthat/test-calcSpatialCov.R
new file mode 100644
index 00000000000..5bedf75614f
--- /dev/null
+++ b/modules/emulator/tests/testthat/test-calcSpatialCov.R
@@ -0,0 +1,44 @@
+test_that("calcSpatialCov.matrix returns correct exponential covariance", {
+  d <- matrix(c(0, 1, 2,
+                1, 0, 1,
+                2, 1, 0), nrow = 3)
+  psi <- 0.5
+  tau <- 2.0
+
+  result <- calcSpatialCov.matrix(d, psi = psi, tau = tau)
+
+  expected <- tau * exp(-psi * d)
+  expect_equal(result, expected)
+})
+
+test_that("calcSpatialCov.matrix returns a matrix of the same dimensions as d", {
+  d <- matrix(runif(25), nrow = 5)
+  result <- calcSpatialCov.matrix(d, psi = 1, tau = 1)
+  expect_equal(dim(result), dim(d))
+})
+
+test_that("calcSpatialCov.matrix diagonal is tau when d has zeros on diagonal", {
+  n <- 4
+  d <- matrix(1, n, n)
+  diag(d) <- 0
+  tau <- 3.0
+  result <- calcSpatialCov.matrix(d, psi = 1, tau = tau)
+  expect_equal(diag(result), rep(tau, n))
+})
+
+test_that("calcSpatialCov.matrix is symmetric when d is symmetric", {
+  d <- as.matrix(dist(1:5))
+  result <- calcSpatialCov.matrix(d, psi = 0.3, tau = 1.5)
+  expect_equal(result, t(result))
+})
+
+test_that("calcSpatialCov.list and calcSpatialCov.matrix agree on single-component case", {
+  d <- as.matrix(dist(1:4))
+  psi <- 0.4
+  tau <- 2.0
+
+  result_matrix <- calcSpatialCov.matrix(d, psi = psi, tau = tau)
+  result_list   <- calcSpatialCov.list(list(d), psi = psi, tau = tau)
+
+  expect_equal(result_matrix, result_list)
+})