diff --git a/workspace/R/2_prep.R b/workspace/R/2_prep.R
index ab879b287f55dde080d4160ae1d57059015c367b..b43ada7a15c6555e0e6bb09e354e10c79de77f50 100644
--- a/workspace/R/2_prep.R
+++ b/workspace/R/2_prep.R
@@ -394,8 +394,9 @@ finalize_gages <- function(SWIM, NHDPlusV2_gages, gfv11_gages, gages_ii_ref,
 #' @param nav_gpkg character path to geopackage where vpu base layers will be written
 #' @param vpu vpu that will be written into output list
 #' @param full_cat_table table of all catchments for each vpu
+#' @param rpu_vpu data.frame ontaining all rpu-vpu pairs
 #' @return a list with the vpu and three tables for flowlines, waterbodies, and catchments
-prepare_vpu_base_layers <- function(ref_gpkg, nav_gpkg, vpu, full_cat_table) {
+prepare_vpu_base_layers <- function(ref_gpkg, nav_gpkg, vpu, full_cat_table, rpu_vpu) {
 
   fline <- read_sf(ref_gpkg, "flowlines") 
   catchment <- read_sf(ref_gpkg, "catchments") 
@@ -439,5 +440,51 @@ prepare_vpu_base_layers <- function(ref_gpkg, nav_gpkg, vpu, full_cat_table) {
   write_sf(cat_network, nav_gpkg, nhd_network)
   write_sf(waterbodies, nav_gpkg, nhd_waterbody)
   
-  return(list(vpu = vpu, flowline = nhd, catchment = cat_network, waterbodies = waterbodies))
+  rpus <- rpu_vpu[rpu_vpu$vpuid == vpu,]$rpuid
+  
+  cats_rpu <- full_cat_table %>%
+    filter(RPUID %in% rpus)
+  
+  cats <- catchment %>% 
+    rename(FEATUREID = featureid, 
+           AREASQKM = areasqkm) %>%
+    mutate(VPUID = vpu) %>%
+    left_join(select(cats_rpu, FEATUREID, RPUID), by = c("FEATUREID")) %>%
+    filter(FEATUREID %in% 
+             unique(c(nhd$COMID, 
+                      full_cat_table$FEATUREID[full_cat_table$RPUID %in% rpus]))) %>%
+    mutate(hy_cats = ifelse(FEATUREID %in% nhd$COMID, 1, 0),
+           full_cats = ifelse(FEATUREID %in% cats_rpu$FEATUREID, 1, 0)) %>%
+    filter(full_cats == 1 | hy_cats == 1) %>%
+    st_sf()
+  
+  write_sf(cats, nav_gpkg, nhd_catchment)
+
+  rpu_list <- setNames(rep(list(list()), length(rpus)), rpus)
+  
+  for(rpu_code in rpus) {
+    
+    out_refac_gpkg <- file.path("cache", paste0(rpu_code, "_refactor.gpkg"))
+    
+    nhd_sub <- subset_rpu(nhd, rpu_code, run_make_standalone = TRUE) %>%
+      st_sf()
+    
+    write_sf(nhd_sub, out_refac_gpkg, nhd_flowline)
+    
+    cats_rpu <- full_cat_table %>%
+      filter(RPUID == rpu_code)
+    
+    catchments <- cats %>%
+      mutate(hy_cats = ifelse(FEATUREID %in% nhd_sub$COMID, 1, 0),
+             full_cats = ifelse(FEATUREID %in% cats_rpu$FEATUREID, 1, 0)) %>%
+      filter(full_cats == 1 | hy_cats == 1) %>%
+      st_sf()
+    
+    write_sf(catchments, out_refac_gpkg, nhd_catchment)
+    
+    rpu_list[[rpu_code]] <- list(flowline = nhd, catchment = catchments)
+    
+  }
+  
+  return(list(vpu = vpu, flowline = nhd, catchment = cats, catchment_network = cat_network, waterbodies = waterbodies, rpus = rpu_list))
 }
diff --git a/workspace/_targets.R b/workspace/_targets.R
index 8f9ff98e5f78b68705d1ef8505880021d291b7eb..cc9d049f9344c5834e9e1078c32b5f3d6aa4b165 100644
--- a/workspace/_targets.R
+++ b/workspace/_targets.R
@@ -80,6 +80,6 @@ list(
              pattern = map(vpu_codes)),
   
   # ~1GB of memory
-  tar_target(vpu_base, prepare_vpu_base_layers(ref_gpkg, nav_gpkg, vpu_codes, full_cat_table),
+  tar_target(vpu_base, prepare_vpu_base_layers(ref_gpkg, nav_gpkg, vpu_codes, full_cat_table, rpu_vpu),
              pattern = map(ref_gpkg, nav_gpkg, vpu_codes), deployment = "worker")
 )
\ No newline at end of file