diff --git a/CHANGELOG.md b/CHANGELOG.md
index 2280cfc967e..bbeba7cc989 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -45,6 +45,7 @@ For more information about this file see also [Keep a Changelog](http://keepacha
 - Added new feature of preparing initial conditions for MODIS LAI, AGB, ISCN SOC, and soil moisture across NA anchor sites.
 - Added GEDI AGB preparation workflow.
 - Added new feature of downloading datasets from the NASA DAAC ORNL database.
+- Added 'aggregate' as a new feature for downscaled data
 
 ### Fixed
 
diff --git a/CITATION.cff b/CITATION.cff
index e00ef7b29f6..7af92146298 100644
--- a/CITATION.cff
+++ b/CITATION.cff
@@ -120,7 +120,9 @@ authors:
   - given-names: Eric R. Scott
     affiliation: University of Arizona
     orcid: 'https://orcid.org/0000-0002-7430-7879'
-  
+  - given-names: Harunobu Ishii
+    affiliation: Boston University Software & Application Innovation Lab(SAIL)
+
 preferred-citation:
   type: article
   title: Facilitating feedbacks between field measurements and ecosystem models
diff --git a/docker/depends/pecan_package_dependencies.csv b/docker/depends/pecan_package_dependencies.csv
index c101b181f51..a258e3cde97 100644
--- a/docker/depends/pecan_package_dependencies.csv
+++ b/docker/depends/pecan_package_dependencies.csv
@@ -61,6 +61,7 @@
 "dplyr",">= 1.1.2","base/db","Imports",FALSE
 "ellipse","*","modules/assim.batch","Imports",FALSE
 "emdbook","*","modules/assim.sequential","Suggests",FALSE
+"exactextractr","*","modules/assim.sequential","Suggests",FALSE
 "foreach","*","base/remote","Imports",FALSE
 "foreach","*","modules/data.atmosphere","Suggests",FALSE
 "foreach","*","modules/data.remote","Imports",FALSE
diff --git a/modules/assim.sequential/DESCRIPTION b/modules/assim.sequential/DESCRIPTION
index 9b228e1cdee..8c9721d2483 100644
--- a/modules/assim.sequential/DESCRIPTION
+++ b/modules/assim.sequential/DESCRIPTION
@@ -32,6 +32,7 @@ Imports:
     stringr
 Suggests:
     corrplot,
+    exactextractr,
     ggrepel,
     emdbook,
     glue,
diff --git a/modules/assim.sequential/NAMESPACE b/modules/assim.sequential/NAMESPACE
index 8157eadc616..bb8aa415da4 100644
--- a/modules/assim.sequential/NAMESPACE
+++ b/modules/assim.sequential/NAMESPACE
@@ -21,6 +21,7 @@ export(SDA_control)
 export(SDA_remote_launcher)
 export(SDA_timeseries_plot)
 export(adj.ens)
+export(aggregate)
 export(alltocs)
 export(alr)
 export(assessParams)
diff --git a/modules/assim.sequential/R/aggregate.R b/modules/assim.sequential/R/aggregate.R
new file mode 100644
index 00000000000..220a85cdb43
--- /dev/null
+++ b/modules/assim.sequential/R/aggregate.R
@@ -0,0 +1,59 @@
+#' @title Aggregation Function
+#' @name aggregate
+#' @author Harunobu Ishii
+#'
+#' @param downscale_output  Raster file output from downscale_function.R. Read file in this way if stored locally: \code{downscale_output <- readRDS("xxx.rds")}
+#' @param polygon_data A spatial polygon object (e.g., an `sf` object) that defines the spatial units for aggregation. 
+#'                     This data should be in a coordinate reference system compatible with the raster data (e.g., "EPSG:4326").
+#' @param func A character string specifying the aggregation function to use (e.g., 'mean', 'sum').
+#' @details This function will aggregate previously downscaled carbon flux amount to a spatial unit of choice 
+#'
+#' @return It returns the `polygon_data` with added columns for mean and sum values of the aggregated raster data for each ensemble member.
+#' @import sf
+#' @import exactextractr
+#' @import raster
+#' @export
+#' @examples
+#'        \dontrun{
+#'        # Download a shapefile of U.S. (polygon data)
+#'        url <- "https://www2.census.gov/geo/tiger/GENZ2020/shp/cb_2020_us_state_20m.zip"
+#'        download.file(url, destfile = "polygon/us_states.zip")
+#'        
+#'        # Unzip the downloaded file and save locally
+#'        unzip("polygon/us_states.zip", exdir = "polygon/us_states")
+#'        us_states <- st_read("polygon/us_states/cb_2020_us_state_20m.shp")
+#'        saveRDS(us_states, "polygon/us_states.rds")
+#'        
+#'        # Load the saved polygon data with Massachusetts as an example
+#'        us_states <- readRDS("polygon/us_states.rds")
+#'        state <- "MA"
+#'        polygon_data <- st_transform(us_states[us_states$STUSPS == state, ], crs = "EPSG:4326")
+#'      
+#'        # Load the downscaled raster output
+#'        downscale_output <- readRDS("path/to/downscale_output.rds")
+#'        
+#'        # Slot in as argument to the aggregate function
+#'        result <- aggregate(downscale_output, polygon_data)
+#'        print(result)
+#'        }
+
+aggregate <- function(downscale_output, polygon_data, func = 'mean'){
+  grand_TTL <- 0
+  if (sf::st_crs(downscale_output$maps$ensemble1) != sf::st_crs(polygon_data)) {
+    stop("CRS of downscale_output and polygon_data must match.")
+  }
+
+  # Perform spatial operations on each raster
+  for (name in names(downscale_output$maps)) {
+    raster_data <- downscale_output$maps[[name]]
+    agg_values <- exactextractr::exact_extract(raster_data, polygon_data, fun = func)
+    
+    polygon_data[[paste0(name, "_", func)]] <- agg_values
+    grand_TTL = grand_TTL + agg_values
+  }
+  if(func == 'mean'){
+    grand_TTL = grand_TTL/length(downscale_output$maps)
+  }
+  polygon_data[[paste0("TTL_", func)]] <- grand_TTL
+  return (polygon_data)
+}
\ No newline at end of file
diff --git a/modules/assim.sequential/man/aggregate.Rd b/modules/assim.sequential/man/aggregate.Rd
new file mode 100644
index 00000000000..965dbdb145e
--- /dev/null
+++ b/modules/assim.sequential/man/aggregate.Rd
@@ -0,0 +1,52 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/aggregate.R
+\name{aggregate}
+\alias{aggregate}
+\title{Aggregation Function}
+\usage{
+aggregate(downscale_output, polygon_data, func = "mean")
+}
+\arguments{
+\item{downscale_output}{Raster file output from downscale_function.R. Read file in this way if stored locally: \code{downscale_output <- readRDS("xxx.rds")}}
+
+\item{polygon_data}{A spatial polygon object (e.g., an `sf` object) that defines the spatial units for aggregation. 
+This data should be in a coordinate reference system compatible with the raster data (e.g., "EPSG:4326").}
+
+\item{func}{A character string specifying the aggregation function to use (e.g., 'mean', 'sum').}
+}
+\value{
+It returns the `polygon_data` with added columns for mean and sum values of the aggregated raster data for each ensemble member.
+}
+\description{
+Aggregation Function
+}
+\details{
+This function will aggregate previously downscaled carbon flux amount to a spatial unit of choice
+}
+\examples{
+       \dontrun{
+       # Download a shapefile of U.S. (polygon data)
+       url <- "https://www2.census.gov/geo/tiger/GENZ2020/shp/cb_2020_us_state_20m.zip"
+       download.file(url, destfile = "polygon/us_states.zip")
+       
+       # Unzip the downloaded file and save locally
+       unzip("polygon/us_states.zip", exdir = "polygon/us_states")
+       us_states <- st_read("polygon/us_states/cb_2020_us_state_20m.shp")
+       saveRDS(us_states, "polygon/polygon_us_states.rds")
+       
+       # Load the saved polygon data with Massachusetts as an example
+       us_states <- readRDS("polygon/us_states.rds")
+       state <- "MA"
+       polygon_data <- st_transform(us_states[us_states$STUSPS == state, ], crs = "EPSG:4326")
+     
+       # Load the downscaled raster output
+       downscale_output <- readRDS("path/to/downscale_output.rds")
+       
+       # Slot in as argument to the aggregate function
+       result <- aggregate(downscale_output, polygon_data)
+       print(result)
+       }
+}
+\author{
+Harunobu Ishii
+}
diff --git a/modules/assim.sequential/tests/testthat/test_aggregation.R b/modules/assim.sequential/tests/testthat/test_aggregation.R
new file mode 100644
index 00000000000..8595aa8c05d
--- /dev/null
+++ b/modules/assim.sequential/tests/testthat/test_aggregation.R
@@ -0,0 +1,60 @@
+library(testthat)
+library(sf)
+library(raster)
+library(exactextractr)
+library(terra)
+source("../../R/aggregate.R")
+test_that("returns aggregated values for RI", {
+    # Load the saved polygon data with Massachusetts as an example
+    us_states <- readRDS("test_aggregation/us_states.rds")
+    state <- "RI"
+    polygon_data <- st_transform(us_states[us_states$STUSPS == state, ], crs = "EPSG:4326")
+
+    # Load the downscaled raster output
+    downscale_output <- list(
+      maps = list(
+        ensemble1 = "test_aggregation/ensemble1.tif",
+        ensemble2 = "test_aggregation/ensemble2.tif",
+        ensemble3 = "test_aggregation/ensemble3.tif"
+      )
+    )
+
+    read_raster <- function(file_path) {
+        rast(file_path)
+    }
+
+    downscale_output$maps <- lapply(downscale_output$maps, read_raster)
+    # Aggregate for RI
+    RI <- aggregate(downscale_output, polygon_data, func = 'mean')
+    comp <- RI$TTL_mean * 10^9
+    comparison_result <- (1.31 < comp & comp < 1.32)
+    expect_true(comparison_result)
+})
+
+test_that("returns error of unmatched CRS", {
+  # Load the saved polygon data with Massachusetts as an example
+  us_states <- readRDS("test_aggregation/us_states.rds")
+  state <- "RI"
+  polygon_data <- st_transform(us_states[us_states$STUSPS == state, ], crs = "EPSG:2222")
+  
+  # Load the downscaled raster output
+  downscale_output <- list(
+    maps = list(
+      ensemble1 = "test_aggregation/ensemble1.tif",
+      ensemble2 = "test_aggregation/ensemble2.tif",
+      ensemble3 = "test_aggregation/ensemble3.tif"
+    )
+  )
+  
+  read_raster <- function(file_path) {
+    rast(file_path)
+  }
+  
+  downscale_output$maps <- lapply(downscale_output$maps, read_raster)
+  expect_error(
+    aggregate(downscale_output, polygon_data, func = 'mean'),
+    "CRS of downscale_output and polygon_data must match."
+  )
+})
+
+
diff --git a/modules/assim.sequential/tests/testthat/test_aggregation/ensemble1.tif b/modules/assim.sequential/tests/testthat/test_aggregation/ensemble1.tif
new file mode 100644
index 00000000000..98e76a1be6d
Binary files /dev/null and b/modules/assim.sequential/tests/testthat/test_aggregation/ensemble1.tif differ
diff --git a/modules/assim.sequential/tests/testthat/test_aggregation/ensemble2.tif b/modules/assim.sequential/tests/testthat/test_aggregation/ensemble2.tif
new file mode 100644
index 00000000000..ae01689cdc3
Binary files /dev/null and b/modules/assim.sequential/tests/testthat/test_aggregation/ensemble2.tif differ
diff --git a/modules/assim.sequential/tests/testthat/test_aggregation/ensemble3.tif b/modules/assim.sequential/tests/testthat/test_aggregation/ensemble3.tif
new file mode 100644
index 00000000000..01478ca60d2
Binary files /dev/null and b/modules/assim.sequential/tests/testthat/test_aggregation/ensemble3.tif differ
diff --git a/modules/assim.sequential/tests/testthat/test_aggregation/us_states.rds b/modules/assim.sequential/tests/testthat/test_aggregation/us_states.rds
new file mode 100644
index 00000000000..dcecf5835a4
Binary files /dev/null and b/modules/assim.sequential/tests/testthat/test_aggregation/us_states.rds differ