diff --git a/R/get-coverage.R b/R/get-coverage.R
index 66974e41c..d7cf58285 100644
--- a/R/get-coverage.R
+++ b/R/get-coverage.R
@@ -68,7 +68,7 @@ get_coverage <- function(forecast, by = "model") {
   # convert to wide interval format and compute interval coverage --------------
   interval_forecast <- quantile_to_interval(forecast, format = "wide")
   interval_forecast[,
-    interval_coverage := (observed <= upper) & (observed >= lower)
+    interval_coverage := check_interval_coverage(observed, lower, upper)
   ][, c("lower", "upper", "observed") := NULL]
   interval_forecast[, interval_coverage_deviation :=
                       interval_coverage - interval_range / 100]
diff --git a/R/helper-quantile-interval-range.R b/R/helper-quantile-interval-range.R
index cba7b5ff7..7ada9cab9 100644
--- a/R/helper-quantile-interval-range.R
+++ b/R/helper-quantile-interval-range.R
@@ -191,3 +191,19 @@ get_range_from_quantile <- function(quantile_level) {
   )
   return(interval_range)
 }
+
+
+#' Check whether observed values fall inside a prediction interval
+#' @description
+#' Internal helper that computes whether each observed value falls within the
+#' bounds defined by `lower` and `upper`. Used by both [get_coverage()] and
+#' [interval_coverage()] to avoid duplicating the bounds-check logic.
+#' @param observed Numeric vector of observed values.
+#' @param lower Numeric vector of lower interval bounds.
+#' @param upper Numeric vector of upper interval bounds.
+#' @returns A logical vector indicating whether each observed value falls
+#'   within the corresponding interval (inclusive on both bounds).
+#' @keywords internal
+check_interval_coverage <- function(observed, lower, upper) {
+  (observed >= lower) & (observed <= upper)
+}
diff --git a/R/metrics-quantile.R b/R/metrics-quantile.R
index fd5775076..bea365749 100644
--- a/R/metrics-quantile.R
+++ b/R/metrics-quantile.R
@@ -360,7 +360,7 @@ interval_coverage <- function(observed, predicted,
   r <- interval_range
   reformatted <- quantile_to_interval(observed, predicted, quantile_level)
   reformatted <- reformatted[interval_range %in% r]
-  reformatted[, interval_coverage := (observed >= lower) & (observed <= upper)]
+  reformatted[, interval_coverage := check_interval_coverage(observed, lower, upper)]
   return(reformatted$interval_coverage)
 }
 
diff --git a/man/check_interval_coverage.Rd b/man/check_interval_coverage.Rd
new file mode 100644
index 000000000..c8363d0c7
--- /dev/null
+++ b/man/check_interval_coverage.Rd
@@ -0,0 +1,25 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/helper-quantile-interval-range.R
+\name{check_interval_coverage}
+\alias{check_interval_coverage}
+\title{Check whether observed values fall inside a prediction interval}
+\usage{
+check_interval_coverage(observed, lower, upper)
+}
+\arguments{
+\item{observed}{Numeric vector of observed values.}
+
+\item{lower}{Numeric vector of lower interval bounds.}
+
+\item{upper}{Numeric vector of upper interval bounds.}
+}
+\value{
+A logical vector indicating whether each observed value falls
+within the corresponding interval (inclusive on both bounds).
+}
+\description{
+Internal helper that computes whether each observed value falls within the
+bounds defined by \code{lower} and \code{upper}. Used by both \code{\link[=get_coverage]{get_coverage()}} and
+\code{\link[=interval_coverage]{interval_coverage()}} to avoid duplicating the bounds-check logic.
+}
+\keyword{internal}
diff --git a/tests/testthat/test-get-coverage.R b/tests/testthat/test-get-coverage.R
index 13d992e06..f23d87743 100644
--- a/tests/testthat/test-get-coverage.R
+++ b/tests/testthat/test-get-coverage.R
@@ -29,6 +29,118 @@ test_that("get_coverage() outputs an object of class c('data.table', 'data.frame
   expect_s3_class(cov, c("data.table", "data.frame"), exact = TRUE)
 })
 
+test_that("get_coverage() interval coverage matches interval_coverage() for same data", {
+  # Regression guard: both functions independently compute the same bounds check
+  fc <- data.table::copy(example_quantile[model == "EuroCOVIDhub-ensemble"])
+  fc <- fc[!is.na(predicted)]
+  fc_obj <- as_forecast_quantile(fc)
+
+  cov <- get_coverage(fc_obj, by = get_forecast_unit(fc_obj))
+
+  # Compare for 50% interval — get_coverage returns multiple rows per forecast
+  # (one per quantile_level), but interval_coverage is the same for all rows
+  # with the same interval_range. Take unique per forecast unit + interval_range.
+  cov_50 <- unique(cov[interval_range == 50, c(get_forecast_unit(fc_obj),
+                                                "interval_range",
+                                                "interval_coverage"),
+                        with = FALSE])
+
+  # Get matching numeric data for interval_coverage()
+  obs <- fc[quantile_level == 0.5]$observed
+  pred_mat <- as.matrix(
+    data.table::dcast(
+      fc, ... ~ quantile_level, value.var = "predicted"
+    )[, .SD, .SDcols = as.character(sort(unique(fc$quantile_level)))]
+  )
+  ql <- sort(unique(fc$quantile_level))
+
+  ic_50 <- interval_coverage(obs, pred_mat, ql, interval_range = 50)
+  expect_equal(cov_50$interval_coverage, as.numeric(ic_50))
+})
+
+test_that("get_coverage() produces correct interval_coverage for known inputs", {
+  # Hand-crafted data with known expected coverage
+  dt1 <- data.table::data.table(
+    observed = rep(5, 3),
+    model = "m1", target_type = "t1",
+    target_end_date = as.Date("2020-01-01"), location = "loc1",
+    quantile_level = c(0.25, 0.5, 0.75),
+    predicted = c(3, 5, 7)
+  )
+  dt2 <- data.table::data.table(
+    observed = rep(10, 3),
+    model = "m1", target_type = "t1",
+    target_end_date = as.Date("2020-01-02"), location = "loc1",
+    quantile_level = c(0.25, 0.5, 0.75),
+    predicted = c(3, 5, 7)
+  )
+  dt <- rbind(dt1, dt2)
+  fc <- as_forecast_quantile(dt)
+
+  cov <- get_coverage(fc, by = get_forecast_unit(fc))
+
+  # For observed=5, 50% interval [3,7]: TRUE (5 >= 3 and 5 <= 7)
+  cov_50_obs5 <- cov[target_end_date == as.Date("2020-01-01") &
+                        interval_range == 50]
+  # interval_coverage is the same for all quantile_levels in this range
+  expect_true(all(cov_50_obs5$interval_coverage == TRUE))
+
+  # For observed=10, 50% interval [3,7]: FALSE (10 > 7)
+  cov_50_obs10 <- cov[target_end_date == as.Date("2020-01-02") &
+                         interval_range == 50]
+  expect_true(all(cov_50_obs10$interval_coverage == FALSE))
+
+  # Quantile coverage for quantile_level=0.5: TRUE for observed=5, FALSE for observed=10
+  qcov_obs5 <- cov[target_end_date == as.Date("2020-01-01") &
+                      quantile_level == 0.5]
+  expect_equal(nrow(qcov_obs5), 1)
+  expect_true(as.logical(qcov_obs5$quantile_coverage))
+
+  qcov_obs10 <- cov[target_end_date == as.Date("2020-01-02") &
+                       quantile_level == 0.5]
+  expect_equal(nrow(qcov_obs10), 1)
+  expect_false(as.logical(qcov_obs10$quantile_coverage))
+})
+
+test_that("get_coverage() and interval_coverage() agree when observation outside all intervals", {
+  dt <- data.table::data.table(
+    observed = rep(100, 5),
+    model = "m1", target_type = "t1",
+    target_end_date = as.Date("2020-01-01"), location = "loc1",
+    quantile_level = c(0.1, 0.25, 0.5, 0.75, 0.9),
+    predicted = c(1, 3, 5, 7, 9)
+  )
+  fc <- as_forecast_quantile(dt)
+  cov <- get_coverage(fc, by = get_forecast_unit(fc))
+
+  # All interval_coverage should be FALSE
+  expect_true(all(cov$interval_coverage == FALSE))
+
+  # interval_coverage() should agree
+  pred_mat <- matrix(c(1, 3, 5, 7, 9), nrow = 1)
+  ql <- c(0.1, 0.25, 0.5, 0.75, 0.9)
+  expect_false(interval_coverage(100, pred_mat, ql, interval_range = 50))
+  expect_false(interval_coverage(100, pred_mat, ql, interval_range = 80))
+})
+
+test_that("refactored interval coverage produces identical output to original", {
+  # Comprehensive regression guard using full example dataset
+  cov <- get_coverage(example_quantile, by = get_forecast_unit(example_quantile))
+  scores <- score(example_quantile)
+
+  # Compare interval_coverage from get_coverage() for range=50 with score()'s interval_coverage_50
+  cov_50 <- cov[interval_range == 50]
+  # Merge on forecast unit to compare
+  fu <- get_forecast_unit(example_quantile)
+  merged <- merge(cov_50, scores, by = fu)
+  expect_equal(merged$interval_coverage, as.numeric(merged$interval_coverage_50))
+
+  # Same for range=90
+  cov_90 <- cov[interval_range == 90]
+  merged_90 <- merge(cov_90, scores, by = fu)
+  expect_equal(merged_90$interval_coverage, as.numeric(merged_90$interval_coverage_90))
+})
+
 test_that("get_coverage() can deal with non-symmetric prediction intervals", {
   # the expected result is that `get_coverage()` just works. However,
   # all interval coverages with missing values should just be `NA`
diff --git a/tests/testthat/test-metrics-quantile.R b/tests/testthat/test-metrics-quantile.R
index b7af908e1..e5a2846b7 100644
--- a/tests/testthat/test-metrics-quantile.R
+++ b/tests/testthat/test-metrics-quantile.R
@@ -687,6 +687,25 @@ test_that("interval_coverage rejects wrong inputs", {
   )
 })
 
+test_that("interval_coverage() produces correct results for boundary cases", {
+  # Observation exactly on lower bound, upper bound, and inside
+  obs <- c(3, 7, 5)
+  pred <- matrix(c(3, 5, 7), nrow = 3, ncol = 3, byrow = TRUE)
+  ql <- c(0.25, 0.5, 0.75)
+  result <- interval_coverage(obs, pred, ql, interval_range = 50)
+  expect_equal(result, c(TRUE, TRUE, TRUE))
+})
+
+test_that("interval_coverage() handles multiple interval ranges correctly", {
+  obs <- c(5)
+  pred <- matrix(c(1, 3, 5, 7, 9), nrow = 1)
+  ql <- c(0.1, 0.25, 0.5, 0.75, 0.9)
+  # 50% interval: [3, 7]
+  expect_equal(interval_coverage(obs, pred, ql, interval_range = 50), TRUE)
+  # 80% interval: [1, 9]
+  expect_equal(interval_coverage(obs, pred, ql, interval_range = 80), TRUE)
+})
+
 test_that("interval_coverage_quantile throws a warning when a required quantile is not available", {
   dropped_quantile_pred <- predicted[, -4]
   dropped_quantiles <- quantile_level[-4]