diff --git a/rQSWATPlus/R/gwflow.R b/rQSWATPlus/R/gwflow.R
index f55b7e7..f1252e1 100644
--- a/rQSWATPlus/R/gwflow.R
+++ b/rQSWATPlus/R/gwflow.R
@@ -97,9 +97,9 @@ qswat_read_gwflow_config <- function(ini_file = NULL) {
 #' Set Up SWAT+ Groundwater Flow (gwflow) Tables
 #'
 #' Creates the gwflow tables in the project SQLite database and populates
-#' `gwflow_base` and `gwflow_solutes` from the supplied configuration.
-#' This is the R equivalent of the gwflow setup performed by the
-#' `GWFlow.createTables()` and related methods in the QGIS plugin's
+#' them from the supplied configuration and available spatial data.
+#' This is the R equivalent of the gwflow setup performed by
+#' `GWFlow.createTables()` and the GIS routines in the QGIS plugin's
 #' `gwflow.py`.
 #'
 #' The function mirrors the SWAT+ gwflow database schema defined in the
@@ -108,9 +108,11 @@ qswat_read_gwflow_config <- function(ini_file = NULL) {
 #' `gwflow_init_conc`, `gwflow_hrucell`, `gwflow_fpcell`,
 #' `gwflow_rivcell`, `gwflow_lsucell`, `gwflow_rescell`.
 #'
-#' After calling this function the database is ready to receive spatial
-#' cell data (zones, grid cells, river cells, etc.) that require GIS
-#' processing not yet automated in rQSWATPlus.
+#' When the project object contains spatial data (`lsu_sf`, `hru_sf`,
+#' `channels_sf` or `streams_sf`) the function also runs the GIS
+#' processing steps to populate the grid and cell connection tables.
+#' This happens automatically when called from [qswat_run()] with
+#' `use_gwflow = TRUE`.
 #'
 #' @param project A `qswat_project` object returned by
 #'   [qswat_write_database()].
@@ -119,25 +121,55 @@ qswat_read_gwflow_config <- function(ini_file = NULL) {
 #'   defaults are used.
 #' @param overwrite Logical.  If `TRUE`, drop and recreate existing gwflow
 #'   tables.  Default `FALSE`.
+#' @param aquifer_thickness Numeric (constant thickness in metres,
+#'   default `20`) **or** a file path to a raster (`.tif`) or polygon
+#'   shapefile (`.shp`) containing aquifer thickness values.  When a
+#'   shapefile is provided it must have a numeric attribute named
+#'   `"thickness"` or `"Avg_Thick"`.  When a raster is provided the
+#'   mean cell value (in metres) is extracted for each grid cell.
+#' @param conductivity_file Optional path to a polygon shapefile that
+#'   defines aquifer hydraulic-conductivity zones.  The shapefile must
+#'   contain a numeric column `"aquifer_k"` (K in m/day) or
+#'   `"logK_Ferr_"` (log10 intrinsic permeability × 100 as stored in
+#'   GLHYMPS).  If `NULL` a single default zone is created using
+#'   `default_aquifer_k`.
+#' @param default_aquifer_k Numeric.  Default hydraulic conductivity
+#'   (m/day) for the single zone created when `conductivity_file` is
+#'   `NULL`.  Default `1.0`.
 #'
 #' @return The `project` object with `project$use_gwflow` set to `TRUE`
 #'   (invisibly).
 #'
 #' @details
-#' The following tables are created:
+#' The following tables are created and, when spatial data are
+#' available, populated:
 #' \describe{
-#'   \item{gwflow_base}{Single-row table with global gwflow parameters.}
-#'   \item{gwflow_zone}{Aquifer hydraulic parameter zones.}
-#'   \item{gwflow_grid}{Active grid cells (populated by GIS workflow).}
-#'   \item{gwflow_out_days}{Optional output day pairs (year, jday).}
-#'   \item{gwflow_obs_locs}{Observation cell IDs.}
-#'   \item{gwflow_solutes}{Chemical solute definitions and initial concentrations.}
-#'   \item{gwflow_init_conc}{Per-cell initial solute concentrations.}
-#'   \item{gwflow_hrucell}{HRU-to-cell area mapping.}
-#'   \item{gwflow_fpcell}{Floodplain-cell connections.}
-#'   \item{gwflow_rivcell}{River-cell connections.}
+#'   \item{gwflow_base}{Single-row table with global gwflow parameters
+#'     (always populated; `row_count`/`col_count` updated after grid
+#'     creation).}
+#'   \item{gwflow_zone}{Aquifer hydraulic parameter zones (one default
+#'     zone, or one row per zone from `conductivity_file`).}
+#'   \item{gwflow_grid}{Active grid cells with elevation, aquifer
+#'     thickness, initial head, and zone assignment.}
+#'   \item{gwflow_out_days}{Optional output day pairs (not populated
+#'     automatically; left empty for user entry).}
+#'   \item{gwflow_obs_locs}{Observation cell IDs (not populated
+#'     automatically; left empty for user entry).}
+#'   \item{gwflow_solutes}{Chemical solute definitions and initial
+#'     concentrations (always populated with 10 default solutes).}
+#'   \item{gwflow_init_conc}{Per-cell initial solute concentrations
+#'     (not populated automatically).}
+#'   \item{gwflow_hrucell}{HRU-to-cell area mapping (populated when
+#'     HRU recharge mode is selected and `hru_sf` is available).}
+#'   \item{gwflow_fpcell}{Floodplain-cell connections (empty in
+#'     rQSWATPlus because LSUs are 1:1 with subbasins, with no
+#'     floodplain landscape units).}
+#'   \item{gwflow_rivcell}{River-cell connections (channel-cell
+#'     intersection lengths).}
 #'   \item{gwflow_lsucell}{LSU-cell area mapping.}
-#'   \item{gwflow_rescell}{Reservoir-cell connections.}
+#'   \item{gwflow_rescell}{Reservoir-cell connections (populated via
+#'     the `gis_water`–`gwflow_lsucell` join when water bodies are
+#'     present).}
 #' }
 #'
 #' @examples
@@ -149,7 +181,10 @@ qswat_read_gwflow_config <- function(ini_file = NULL) {
 #' @export
 qswat_setup_gwflow <- function(project,
                                 gwflow_config = NULL,
-                                overwrite = FALSE) {
+                                overwrite = FALSE,
+                                aquifer_thickness = 20.0,
+                                conductivity_file = NULL,
+                                default_aquifer_k = 1.0) {
 
   if (!inherits(project, "qswat_project")) {
     stop("'project' must be a qswat_project object.", call. = FALSE)
@@ -169,6 +204,25 @@ qswat_setup_gwflow <- function(project,
   .write_gwflow_base(con, gwflow_config)
   .write_gwflow_solutes(con, gwflow_config)
 
+  # Populate spatial tables (grid, zones, cell connections) when spatial
+  # data from qswat_create_hrus() / qswat_create_streams() is present.
+  if (!is.null(project$lsu_sf) && inherits(project$lsu_sf, "sf") &&
+      nrow(project$lsu_sf) > 0) {
+    .populate_gwflow_tables(
+      con               = con,
+      project           = project,
+      cfg               = gwflow_config,
+      aquifer_thickness = aquifer_thickness,
+      conductivity_file = conductivity_file,
+      default_aquifer_k = default_aquifer_k
+    )
+  } else {
+    message("Note: No spatial data in project object; gwflow_zone, ",
+            "gwflow_grid, gwflow_rivcell, gwflow_lsucell, and ",
+            "gwflow_hrucell remain empty.")
+    message("  Run qswat_setup_gwflow() after qswat_run() to populate them.")
+  }
+
   # Update project_config to flag gwflow as enabled
   DBI::dbExecute(con, "UPDATE project_config SET use_gwflow = 1")
 
@@ -400,3 +454,443 @@ qswat_setup_gwflow <- function(project,
   DBI::dbWriteTable(con, "gwflow_solutes", solutes,
                      append = TRUE, row.names = FALSE)
 }
+
+
+#' Populate gwflow Spatial Tables from Project GIS Data
+#'
+#' Runs the GIS processing steps that mirror the Python gwflow.py routines:
+#' fishnet grid creation, active-cell detection, elevation extraction,
+#' conductivity zone assignment, and spatial intersections with channels,
+#' LSUs, and (optionally) HRUs and water bodies.
+#'
+#' @param con  DBI connection to the project SQLite database.
+#' @param project `qswat_project` object with `lsu_sf`, `hru_sf`,
+#'   `channels_sf`/`streams_sf`, and `dem_file` populated.
+#' @param cfg  gwflow configuration list from [qswat_read_gwflow_config()].
+#' @param aquifer_thickness Numeric constant (metres) or file path.
+#' @param conductivity_file Path to conductivity shapefile, or `NULL`.
+#' @param default_aquifer_k Default K in m/day when no file is given.
+#' @noRd
+.populate_gwflow_tables <- function(con, project, cfg,
+                                     aquifer_thickness = 20.0,
+                                     conductivity_file = NULL,
+                                     default_aquifer_k = 1.0) {
+
+  watershed_sf <- project$lsu_sf
+  crs_obj      <- sf::st_crs(watershed_sf)
+  cell_size    <- cfg$cell_size
+
+  # ------------------------------------------------------------------
+  # 1. Build fishnet grid over the watershed bounding box
+  # ------------------------------------------------------------------
+  message("  [gwflow] Creating fishnet grid (cell_size = ", cell_size, " m)...")
+  bbox_obj <- sf::st_bbox(watershed_sf)
+
+  grid_geom <- sf::st_make_grid(
+    watershed_sf,
+    cellsize = cell_size,
+    what     = "polygons"
+  )
+  grid_sf <- sf::st_sf(
+    cell_id  = seq_along(grid_geom),
+    geometry = grid_geom,
+    crs      = crs_obj
+  )
+
+  # Compute grid dimensions for gwflow_base
+  n_cols <- ceiling((bbox_obj["xmax"] - bbox_obj["xmin"]) / cell_size)
+  n_rows <- ceiling((bbox_obj["ymax"] - bbox_obj["ymin"]) / cell_size)
+
+  # ------------------------------------------------------------------
+  # 2. Identify active cells (those that intersect the watershed)
+  # ------------------------------------------------------------------
+  message("  [gwflow] Identifying active cells...")
+  watershed_union <- sf::st_union(watershed_sf)
+  hits <- which(lengths(sf::st_intersects(grid_sf, watershed_union)) > 0)
+  if (length(hits) == 0) {
+    message("  [gwflow] Warning: no grid cells intersect the watershed.")
+    return(invisible(NULL))
+  }
+  grid_active <- grid_sf[hits, , drop = FALSE]
+
+  # ------------------------------------------------------------------
+  # 3. Assign elevation from DEM
+  # ------------------------------------------------------------------
+  message("  [gwflow] Extracting DEM elevation per cell...")
+  cell_elev <- rep(0.0, nrow(grid_active))
+  if (!is.null(project$dem_file) && file.exists(project$dem_file)) {
+    dem_rast <- terra::rast(project$dem_file)
+    # Project grid to DEM CRS for extraction
+    grid_v_dem <- terra::vect(
+      sf::st_transform(grid_active, crs = terra::crs(dem_rast))
+    )
+    ex <- terra::extract(dem_rast, grid_v_dem, fun = "mean", na.rm = TRUE)
+    if (ncol(ex) >= 2) {
+      vals <- ex[, 2]
+      vals[is.na(vals)] <- 0.0
+      cell_elev <- as.numeric(vals)
+    }
+  }
+
+  # ------------------------------------------------------------------
+  # 4. Assign aquifer thickness
+  # ------------------------------------------------------------------
+  message("  [gwflow] Assigning aquifer thickness...")
+  # Resolve the default/fallback thickness once (avoids repeated numeric coercion)
+  default_thick <- if (is.numeric(aquifer_thickness)) {
+    aquifer_thickness[1]
+  } else {
+    20.0  # numeric default when a file path is given
+  }
+  cell_thick <- rep(default_thick, nrow(grid_active))
+
+  if (is.character(aquifer_thickness) && file.exists(aquifer_thickness)) {
+    ext <- tolower(tools::file_ext(aquifer_thickness))
+    if (ext %in% c("tif", "tiff", "img")) {
+      # Raster source
+      thick_rast <- terra::rast(aquifer_thickness)
+      grid_v_aq  <- terra::vect(
+        sf::st_transform(grid_active, crs = terra::crs(thick_rast))
+      )
+      ex_t <- terra::extract(thick_rast, grid_v_aq, fun = "mean", na.rm = TRUE)
+      if (ncol(ex_t) >= 2) {
+        vals <- ex_t[, 2]
+        vals[is.na(vals) | vals <= 0] <- default_thick
+        cell_thick <- as.numeric(vals)
+      }
+    } else if (ext %in% c("shp", "gpkg", "geojson")) {
+      # Vector source: spatial join and average
+      thick_sf <- sf::st_read(aquifer_thickness, quiet = TRUE)
+      thick_sf  <- sf::st_transform(thick_sf, crs = crs_obj)
+      # Detect the thickness column
+      t_col <- intersect(c("thickness", "Avg_Thick"),
+                         names(thick_sf))[1]
+      if (!is.na(t_col)) {
+        grid_j <- sf::st_join(
+          grid_active[, "cell_id"],
+          thick_sf[, c(t_col, "geometry")],
+          join = sf::st_intersects
+        )
+        t_agg <- stats::aggregate(
+          grid_j[[t_col]],
+          by  = list(cell_id = grid_j$cell_id),
+          FUN = function(x) mean(x[x > 0], na.rm = TRUE)
+        )
+        idx <- match(grid_active$cell_id, t_agg$cell_id)
+        vals <- t_agg$x[idx]
+        vals[is.na(vals) | vals <= 0] <- default_thick
+        cell_thick <- as.numeric(vals)
+      }
+    }
+  }
+
+  # ------------------------------------------------------------------
+  # 5. Assign conductivity zones
+  # ------------------------------------------------------------------
+  message("  [gwflow] Assigning conductivity zones...")
+  cell_zone <- rep(1L, nrow(grid_active))
+
+  if (!is.null(conductivity_file) && file.exists(conductivity_file)) {
+    k_sf <- sf::st_read(conductivity_file, quiet = TRUE)
+    k_sf <- sf::st_transform(k_sf, crs = crs_obj)
+
+    # Detect the K column
+    if ("aquifer_k" %in% names(k_sf)) {
+      k_col <- "aquifer_k"
+    } else if ("logK_Ferr_" %in% names(k_sf)) {
+      # GLHYMPS convention: logK_Ferr_ is log10(K_intrinsic) × 100 (m²).
+      # Convert to hydraulic conductivity K [m/day]:
+      #   K_intrinsic [m²] → K_saturated [m/s] via:  K_sat = K_intr * ρg/μ
+      #     ρ = 1000 kg/m³ (water density)
+      #     g = 9.81 m/s²  (gravitational acceleration)
+      #     μ = 0.001 Pa·s (dynamic viscosity at ~20 °C)
+      #   Then convert [m/s] → [m/day]:  × 86400 s/day
+      k_sf$aquifer_k <- (10^(k_sf$logK_Ferr_ / 100)) * 1000 * 9.81 / 0.001 * 86400
+      k_col <- "aquifer_k"
+    } else {
+      message("  [gwflow] Warning: conductivity shapefile has no 'aquifer_k' ",
+              "or 'logK_Ferr_' column; using single default zone.")
+      k_col <- NULL
+    }
+
+    if (!is.null(k_col)) {
+      # Assign sequential zone IDs to conductivity polygons
+      k_sf$zone_id <- seq_len(nrow(k_sf))
+
+      # Spatial join: assign zone to each active cell (take first match)
+      grid_j <- sf::st_join(
+        grid_active[, "cell_id"],
+        k_sf[, c("zone_id", k_col, "geometry")],
+        join = sf::st_intersects
+      )
+      # If a cell hits multiple zones take the one with the smallest zone_id
+      z_agg <- stats::aggregate(
+        grid_j$zone_id,
+        by  = list(cell_id = grid_j$cell_id),
+        FUN = function(x) min(x, na.rm = TRUE)
+      )
+      idx <- match(grid_active$cell_id, z_agg$cell_id)
+      zone_vals <- z_agg$x[idx]
+      zone_vals[is.na(zone_vals)] <- 1L
+      cell_zone <- as.integer(zone_vals)
+
+      # Build zone table
+      k_agg <- stats::aggregate(
+        k_sf[[k_col]],
+        by  = list(zone_id = k_sf$zone_id),
+        FUN = function(x) mean(x, na.rm = TRUE)
+      )
+      zone_df <- data.frame(
+        zone_id             = k_agg$zone_id,
+        aquifer_k           = k_agg$x,
+        specific_yield      = cfg$init_sy,
+        streambed_k         = cfg$streambed_k,
+        streambed_thickness = cfg$streambed_thick,
+        stringsAsFactors    = FALSE
+      )
+    } else {
+      # No usable K column → fall through to single default zone below
+      conductivity_file <- NULL
+    }
+  }
+
+  if (is.null(conductivity_file)) {
+    # Single default zone
+    zone_df <- data.frame(
+      zone_id            = 1L,
+      aquifer_k          = default_aquifer_k,
+      specific_yield     = cfg$init_sy,
+      streambed_k        = cfg$streambed_k,
+      streambed_thickness = cfg$streambed_thick,
+      stringsAsFactors   = FALSE
+    )
+  }
+
+  # ------------------------------------------------------------------
+  # 6. Detect boundary cells
+  # ------------------------------------------------------------------
+  message("  [gwflow] Detecting boundary cells...")
+  boundary_line <- sf::st_cast(
+    sf::st_boundary(watershed_union), "MULTILINESTRING"
+  )
+  is_boundary <- as.integer(
+    lengths(sf::st_intersects(grid_active, boundary_line)) > 0
+  )
+
+  # status = 1 (active interior) or 2 (active + boundary)
+  cell_status <- 1L + is_boundary
+
+  # ------------------------------------------------------------------
+  # 7. Write gwflow_zone
+  # ------------------------------------------------------------------
+  message("  [gwflow] Writing gwflow_zone (", nrow(zone_df), " zone(s))...")
+  DBI::dbWriteTable(con, "gwflow_zone", zone_df,
+                    append = TRUE, row.names = FALSE)
+
+  # ------------------------------------------------------------------
+  # 8. Write gwflow_grid (active cells only, mirroring Python status logic)
+  # ------------------------------------------------------------------
+  message("  [gwflow] Writing gwflow_grid (", nrow(grid_active),
+          " active cells)...")
+  grid_df <- data.frame(
+    cell_id           = grid_active$cell_id,
+    status            = cell_status,
+    zone              = cell_zone,
+    elevation         = cell_elev,
+    aquifer_thickness = cell_thick,
+    extinction_depth  = cfg$exdp,
+    initial_head      = cell_elev - cfg$wt_depth,
+    tile              = 0L,
+    stringsAsFactors  = FALSE
+  )
+  DBI::dbWriteTable(con, "gwflow_grid", grid_df,
+                    append = TRUE, row.names = FALSE)
+
+  # Update gwflow_base row/col counts
+  DBI::dbExecute(
+    con,
+    sprintf("UPDATE gwflow_base SET row_count = %d, col_count = %d",
+            as.integer(n_rows), as.integer(n_cols))
+  )
+
+  # ------------------------------------------------------------------
+  # 9. gwflow_rivcell: channel–cell intersections
+  # ------------------------------------------------------------------
+  channels_sf <- project$channels_sf
+  if (is.null(channels_sf) && !is.null(project$streams_sf)) {
+    channels_sf <- project$streams_sf
+  }
+
+  if (!is.null(channels_sf) && inherits(channels_sf, "sf") &&
+      nrow(channels_sf) > 0) {
+    message("  [gwflow] Computing channel–cell intersections...")
+
+    tryCatch({
+      # Ensure matching CRS
+      channels_sf <- sf::st_transform(channels_sf, crs = crs_obj)
+
+      # Read gis_channels id / subbasin mapping from DB
+      gis_ch <- DBI::dbGetQuery(con, "SELECT id, subbasin FROM gis_channels")
+
+      # Intersect: keep only grid cells that a channel runs through
+      # Build a minimal sf with just WSNO (= subbasin ID) and geometry
+      wsno_vals <- if ("WSNO" %in% names(channels_sf)) channels_sf$WSNO
+                   else if ("LINKNO" %in% names(channels_sf)) channels_sf$LINKNO
+                   else seq_len(nrow(channels_sf))
+
+      channels_work <- sf::st_sf(
+        WSNO     = wsno_vals,
+        geometry = sf::st_geometry(channels_sf),
+        crs      = crs_obj
+      )
+
+      riv_inter <- sf::st_intersection(
+        channels_work,
+        grid_active[, c("cell_id", "geometry")]
+      )
+
+      if (nrow(riv_inter) > 0) {
+        # Only keep line/multiline geometries (the channel segments)
+        keep <- sf::st_geometry_type(riv_inter) %in%
+          c("LINESTRING", "MULTILINESTRING")
+        riv_inter <- riv_inter[keep, , drop = FALSE]
+      }
+
+      if (nrow(riv_inter) > 0) {
+        riv_inter$length_m <- as.numeric(sf::st_length(riv_inter))
+
+        # Map WSNO → gis_channels.id
+        ch_id <- gis_ch$id[match(riv_inter$WSNO, gis_ch$subbasin)]
+
+        riv_df <- data.frame(
+          cell_id  = riv_inter$cell_id,
+          channel  = ch_id,
+          length_m = riv_inter$length_m,
+          stringsAsFactors = FALSE
+        )
+        riv_df <- riv_df[!is.na(riv_df$channel) & riv_df$length_m > 0, ]
+
+        if (nrow(riv_df) > 0) {
+          message("  [gwflow] Writing gwflow_rivcell (",
+                  nrow(riv_df), " row(s))...")
+          DBI::dbWriteTable(con, "gwflow_rivcell", riv_df,
+                            append = TRUE, row.names = FALSE)
+        }
+      }
+    }, error = function(e) {
+      message("  [gwflow] Warning: could not compute channel–cell ",
+              "intersections: ", conditionMessage(e))
+    })
+  }
+
+  # ------------------------------------------------------------------
+  # 10. gwflow_lsucell: LSU–cell area intersections
+  # ------------------------------------------------------------------
+  message("  [gwflow] Computing LSU–cell intersections...")
+  tryCatch({
+    gis_lsus <- DBI::dbGetQuery(con, "SELECT id, subbasin FROM gis_lsus")
+
+    lsu_inter <- sf::st_intersection(
+      watershed_sf[, c("subbasin", "geometry")],
+      grid_active[, c("cell_id", "geometry")]
+    )
+
+    if (nrow(lsu_inter) > 0) {
+      lsu_inter$area_m2 <- as.numeric(sf::st_area(lsu_inter))
+      lsu_id <- gis_lsus$id[match(lsu_inter$subbasin, gis_lsus$subbasin)]
+
+      lsu_df <- data.frame(
+        cell_id = lsu_inter$cell_id,
+        lsu     = lsu_id,
+        area_m2 = lsu_inter$area_m2,
+        stringsAsFactors = FALSE
+      )
+      lsu_df <- lsu_df[!is.na(lsu_df$lsu) & lsu_df$area_m2 > 0, ]
+
+      if (nrow(lsu_df) > 0) {
+        message("  [gwflow] Writing gwflow_lsucell (",
+                nrow(lsu_df), " row(s))...")
+        DBI::dbWriteTable(con, "gwflow_lsucell", lsu_df,
+                          append = TRUE, row.names = FALSE)
+      }
+    }
+  }, error = function(e) {
+    message("  [gwflow] Warning: could not compute LSU–cell intersections: ",
+            conditionMessage(e))
+  })
+
+  # ------------------------------------------------------------------
+  # 11. gwflow_hrucell: HRU–cell area intersections (HRU recharge mode)
+  # ------------------------------------------------------------------
+  hru_recharge <- cfg$hruorlsu_recharge == 1L || cfg$hruorlsu_recharge == 3L
+
+  if (hru_recharge &&
+      !is.null(project$hru_sf) && inherits(project$hru_sf, "sf") &&
+      nrow(project$hru_sf) > 0) {
+    message("  [gwflow] Computing HRU–cell intersections...")
+    tryCatch({
+      hru_inter <- sf::st_intersection(
+        project$hru_sf[, c("hru_id", "geometry")],
+        grid_active[, c("cell_id", "geometry")]
+      )
+
+      if (nrow(hru_inter) > 0) {
+        hru_inter$area_m2 <- as.numeric(sf::st_area(hru_inter))
+
+        hru_df <- data.frame(
+          cell_id = hru_inter$cell_id,
+          hru     = hru_inter$hru_id,
+          area_m2 = hru_inter$area_m2,
+          stringsAsFactors = FALSE
+        )
+        hru_df <- hru_df[!is.na(hru_df$hru) & hru_df$area_m2 > 0, ]
+
+        if (nrow(hru_df) > 0) {
+          message("  [gwflow] Writing gwflow_hrucell (",
+                  nrow(hru_df), " row(s))...")
+          DBI::dbWriteTable(con, "gwflow_hrucell", hru_df,
+                            append = TRUE, row.names = FALSE)
+        }
+      }
+    }, error = function(e) {
+      message("  [gwflow] Warning: could not compute HRU–cell intersections: ",
+              conditionMessage(e))
+    })
+  }
+
+  # ------------------------------------------------------------------
+  # 12. gwflow_rescell: reservoir–cell connections via gis_water join
+  #     (mirrors the Python SQL: SELECT gwflow_lsucell.cell_id,
+  #      gis_water.id, gis_water.elev FROM gis_water
+  #      INNER JOIN gwflow_lsucell USING (lsu))
+  # ------------------------------------------------------------------
+  tryCatch({
+    gis_water_n <- DBI::dbGetQuery(
+      con, "SELECT COUNT(*) AS n FROM gis_water"
+    )$n
+    lsucell_n <- DBI::dbGetQuery(
+      con, "SELECT COUNT(*) AS n FROM gwflow_lsucell"
+    )$n
+
+    if (gis_water_n > 0 && lsucell_n > 0) {
+      res_rows <- DBI::dbGetQuery(con, "
+        SELECT gwflow_lsucell.cell_id AS cell_id,
+               gis_water.id          AS res,
+               gis_water.elev        AS res_stage
+        FROM gis_water
+        INNER JOIN gwflow_lsucell ON gis_water.lsu = gwflow_lsucell.lsu
+      ")
+      if (nrow(res_rows) > 0) {
+        message("  [gwflow] Writing gwflow_rescell (",
+                nrow(res_rows), " row(s))...")
+        DBI::dbWriteTable(con, "gwflow_rescell", res_rows,
+                          append = TRUE, row.names = FALSE)
+      }
+    }
+  }, error = function(e) {
+    message("  [gwflow] Note: could not populate gwflow_rescell: ",
+            conditionMessage(e))
+  })
+
+  invisible(NULL)
+}
diff --git a/rQSWATPlus/R/validate_database.R b/rQSWATPlus/R/validate_database.R
index de51294..e0b407e 100644
--- a/rQSWATPlus/R/validate_database.R
+++ b/rQSWATPlus/R/validate_database.R
@@ -40,6 +40,13 @@
 #'   \item **HAWQS tables**: HAWQS-specific tables (plant_HAWQS,
 #'     urban_HAWQS, CDL landuse field tables, statsgo_ssurgo_lkey)
 #'     from QSWATPlusProjHAWQS.sqlite are present
+#'   \item **gwflow tables**: when `use_gwflow = 1` in project_config,
+#'     all twelve gwflow tables exist and the core tables
+#'     (gwflow_base, gwflow_zone, gwflow_grid, gwflow_solutes,
+#'     gwflow_rivcell) are populated; recharge-specific cell mapping
+#'     tables (gwflow_hrucell for HRU recharge, gwflow_lsucell for
+#'     LSU recharge) are non-empty according to the recharge mode
+#'     stored in gwflow_base
 #' }
 #'
 #' @examples
@@ -394,6 +401,143 @@ qswat_check_database <- function(db_file, verbose = TRUE) {
                 paste(missing_hawqs, collapse = ", ")))
   }
 
+  # ---- 10. gwflow tables (if gwflow is enabled) ----
+  use_gwflow <- tryCatch({
+    pc <- DBI::dbGetQuery(
+      con, "SELECT use_gwflow FROM project_config LIMIT 1"
+    )
+    if (nrow(pc) > 0 && !is.na(pc$use_gwflow)) as.integer(pc$use_gwflow)
+    else 0L
+  }, error = function(e) 0L)
+
+  if (use_gwflow == 1L) {
+    gwflow_all_tables <- c(
+      "gwflow_base", "gwflow_zone", "gwflow_grid",
+      "gwflow_out_days", "gwflow_obs_locs", "gwflow_solutes",
+      "gwflow_init_conc", "gwflow_hrucell", "gwflow_fpcell",
+      "gwflow_rivcell", "gwflow_lsucell", "gwflow_rescell"
+    )
+    missing_gw <- gwflow_all_tables[!gwflow_all_tables %in% existing_tables]
+    ok <- length(missing_gw) == 0
+    record(
+      "gwflow:tables_exist",
+      ok,
+      if (ok) "All gwflow tables exist"
+      else paste0("gwflow is enabled but ", length(missing_gw),
+                  " table(s) are missing: ",
+                  paste(missing_gw, collapse = ", "))
+    )
+
+    if (ok) {
+      # gwflow_base: exactly 1 row
+      n_base <- DBI::dbGetQuery(
+        con, "SELECT COUNT(*) AS n FROM gwflow_base"
+      )$n
+      ok_base <- n_base == 1L
+      record(
+        "gwflow:base_populated",
+        ok_base,
+        if (ok_base) "gwflow_base has 1 configuration row"
+        else paste0("gwflow_base must have exactly 1 row (found ", n_base, ")")
+      )
+
+      # gwflow_zone: at least 1 row
+      n_zone <- DBI::dbGetQuery(
+        con, "SELECT COUNT(*) AS n FROM gwflow_zone"
+      )$n
+      ok_zone <- n_zone > 0
+      record(
+        "gwflow:zone_populated",
+        ok_zone,
+        if (ok_zone) paste0("gwflow_zone has ", n_zone, " zone(s)")
+        else "gwflow_zone is empty (no aquifer zones defined)"
+      )
+
+      # gwflow_grid: at least 1 row
+      n_grid <- DBI::dbGetQuery(
+        con, "SELECT COUNT(*) AS n FROM gwflow_grid"
+      )$n
+      ok_grid <- n_grid > 0
+      record(
+        "gwflow:grid_populated",
+        ok_grid,
+        if (ok_grid) paste0("gwflow_grid has ", n_grid, " cell(s)")
+        else "gwflow_grid is empty (no grid cells defined)"
+      )
+
+      # gwflow_solutes: at least 1 row (Python writes 10 default solutes)
+      n_sol <- DBI::dbGetQuery(
+        con, "SELECT COUNT(*) AS n FROM gwflow_solutes"
+      )$n
+      ok_sol <- n_sol > 0
+      record(
+        "gwflow:solutes_populated",
+        ok_sol,
+        if (ok_sol) paste0("gwflow_solutes has ", n_sol, " solute(s)")
+        else "gwflow_solutes is empty (no solutes defined)"
+      )
+
+      # gwflow_rivcell: at least 1 row
+      n_riv <- DBI::dbGetQuery(
+        con, "SELECT COUNT(*) AS n FROM gwflow_rivcell"
+      )$n
+      ok_riv <- n_riv > 0
+      record(
+        "gwflow:rivcell_populated",
+        ok_riv,
+        if (ok_riv) paste0("gwflow_rivcell has ", n_riv, " river cell(s)")
+        else "gwflow_rivcell is empty (no river-cell connections defined)"
+      )
+
+      # Recharge-specific checks (mirrors gwflow.py HRUorLSU_recharge logic):
+      #   recharge = 1 → HRU-cell only   → gwflow_hrucell must be populated
+      #   recharge = 2 → LSU-cell only   → gwflow_lsucell must be populated
+      #   recharge = 3 → both HRU + LSU  → both must be populated
+      # Only evaluate when gwflow_base has exactly 1 row (ok_base is TRUE);
+      # otherwise the table is already flagged as empty above.
+      if (ok_base) {
+        recharge <- as.integer(
+          DBI::dbGetQuery(con,
+            "SELECT recharge FROM gwflow_base LIMIT 1"
+          )$recharge
+        )
+
+        hru_recharge <- recharge == 1L || recharge == 3L
+        lsu_recharge <- recharge == 2L || recharge == 3L
+
+        if (hru_recharge) {
+          n_hru <- DBI::dbGetQuery(
+            con, "SELECT COUNT(*) AS n FROM gwflow_hrucell"
+          )$n
+          ok_hru <- n_hru > 0
+          record(
+            "gwflow:hrucell_populated",
+            ok_hru,
+            if (ok_hru) paste0("gwflow_hrucell has ", n_hru,
+                                " HRU-cell mapping(s)")
+            else paste0("gwflow_hrucell is empty but recharge mode (", recharge,
+                        ") requires HRU-cell connections")
+          )
+        }
+
+        if (lsu_recharge) {
+          n_lsu <- DBI::dbGetQuery(
+            con, "SELECT COUNT(*) AS n FROM gwflow_lsucell"
+          )$n
+          ok_lsu <- n_lsu > 0
+          record(
+            "gwflow:lsucell_populated",
+            ok_lsu,
+            if (ok_lsu) paste0("gwflow_lsucell has ", n_lsu,
+                                " LSU-cell mapping(s)")
+            else paste0("gwflow_lsucell is empty but recharge mode (", recharge,
+                        ") requires LSU-cell connections")
+          )
+        }
+      }
+    }
+  }
+
   result <- .build_check_result(checks, warnings_list, errors_list)
   if (verbose) .print_check_result(result)
   return(invisible(result))
diff --git a/rQSWATPlus/R/workflow.R b/rQSWATPlus/R/workflow.R
index 1687d54..46d2d3a 100644
--- a/rQSWATPlus/R/workflow.R
+++ b/rQSWATPlus/R/workflow.R
@@ -40,6 +40,15 @@
 #' @param gwflow_config A named list of gwflow settings as returned by
 #'   [qswat_read_gwflow_config()]. Only used when `use_gwflow = TRUE`.
 #'   If `NULL` the bundled `gwflow.ini` defaults are used.
+#' @param aquifer_thickness Numeric constant (aquifer thickness in metres,
+#'   default `20`) or a file path to a raster or polygon shapefile
+#'   containing aquifer thickness values. Passed to [qswat_setup_gwflow()]
+#'   when `use_gwflow = TRUE`.
+#' @param conductivity_file Optional path to a polygon shapefile defining
+#'   aquifer hydraulic-conductivity zones. Passed to [qswat_setup_gwflow()]
+#'   when `use_gwflow = TRUE`. If `NULL` a single default zone is used.
+#' @param default_aquifer_k Default hydraulic conductivity (m/day) for the
+#'   single zone created when `conductivity_file = NULL`. Default `1.0`.
 #' @param use_aquifers Logical. Create SWAT+ aquifer objects. Default
 #'   TRUE.
 #' @param n_processes Integer. Number of MPI processes. Default 1.
@@ -97,6 +106,9 @@ qswat_run <- function(project_dir,
                       target_hrus = NULL,
                       use_gwflow = FALSE,
                       gwflow_config = NULL,
+                      aquifer_thickness = 20.0,
+                      conductivity_file = NULL,
+                      default_aquifer_k = 1.0,
                       use_aquifers = TRUE,
                       n_processes = 1L,
                       quiet = FALSE,
@@ -160,8 +172,14 @@ qswat_run <- function(project_dir,
   # Optional Step 6: Set up gwflow tables
   if (isTRUE(use_gwflow)) {
     if (!quiet) message("\n=== Step 6/6: Setting up gwflow tables ===")
-    project <- qswat_setup_gwflow(project, gwflow_config = gwflow_config,
-                                  overwrite = TRUE)
+    project <- qswat_setup_gwflow(
+      project,
+      gwflow_config     = gwflow_config,
+      overwrite         = TRUE,
+      aquifer_thickness = aquifer_thickness,
+      conductivity_file = conductivity_file,
+      default_aquifer_k = default_aquifer_k
+    )
   }
 
   if (!quiet) message("\n=== QSWATPlus workflow complete! ===")
diff --git a/rQSWATPlus/man/qswat_check_database.Rd b/rQSWATPlus/man/qswat_check_database.Rd
index 95637d4..de39c04 100644
--- a/rQSWATPlus/man/qswat_check_database.Rd
+++ b/rQSWATPlus/man/qswat_check_database.Rd
@@ -50,6 +50,13 @@ present and non-empty
 \item \strong{HAWQS tables}: HAWQS-specific tables (plant_HAWQS,
 urban_HAWQS, CDL landuse field tables, statsgo_ssurgo_lkey)
 from QSWATPlusProjHAWQS.sqlite are present
+\item \strong{gwflow tables}: when \code{use_gwflow = 1} in project_config,
+all twelve gwflow tables exist and the core tables
+(gwflow_base, gwflow_zone, gwflow_grid, gwflow_solutes,
+gwflow_rivcell) are populated; recharge-specific cell mapping
+tables (gwflow_hrucell for HRU recharge, gwflow_lsucell for
+LSU recharge) are non-empty according to the recharge mode
+stored in gwflow_base
 }
 }
 \examples{
diff --git a/rQSWATPlus/man/qswat_run.Rd b/rQSWATPlus/man/qswat_run.Rd
index bdce0ae..f836849 100644
--- a/rQSWATPlus/man/qswat_run.Rd
+++ b/rQSWATPlus/man/qswat_run.Rd
@@ -24,6 +24,9 @@ qswat_run(
   target_hrus = NULL,
   use_gwflow = FALSE,
   gwflow_config = NULL,
+  aquifer_thickness = 20,
+  conductivity_file = NULL,
+  default_aquifer_k = 1,
   use_aquifers = TRUE,
   n_processes = 1L,
   quiet = FALSE,
@@ -87,6 +90,18 @@ project database. Default FALSE.}
 \code{\link[=qswat_read_gwflow_config]{qswat_read_gwflow_config()}}. Only used when \code{use_gwflow = TRUE}.
 If \code{NULL} the bundled \code{gwflow.ini} defaults are used.}
 
+\item{aquifer_thickness}{Numeric constant (aquifer thickness in metres,
+default \code{20}) or a file path to a raster or polygon shapefile
+containing aquifer thickness values. Passed to \code{\link[=qswat_setup_gwflow]{qswat_setup_gwflow()}}
+when \code{use_gwflow = TRUE}.}
+
+\item{conductivity_file}{Optional path to a polygon shapefile defining
+aquifer hydraulic-conductivity zones. Passed to \code{\link[=qswat_setup_gwflow]{qswat_setup_gwflow()}}
+when \code{use_gwflow = TRUE}. If \code{NULL} a single default zone is used.}
+
+\item{default_aquifer_k}{Default hydraulic conductivity (m/day) for the
+single zone created when \code{conductivity_file = NULL}. Default \code{1.0}.}
+
 \item{use_aquifers}{Logical. Create SWAT+ aquifer objects. Default
 TRUE.}
 
diff --git a/rQSWATPlus/man/qswat_setup_gwflow.Rd b/rQSWATPlus/man/qswat_setup_gwflow.Rd
index 0dc79e6..102a1d2 100644
--- a/rQSWATPlus/man/qswat_setup_gwflow.Rd
+++ b/rQSWATPlus/man/qswat_setup_gwflow.Rd
@@ -4,7 +4,14 @@
 \alias{qswat_setup_gwflow}
 \title{Set Up SWAT+ Groundwater Flow (gwflow) Tables}
 \usage{
-qswat_setup_gwflow(project, gwflow_config = NULL, overwrite = FALSE)
+qswat_setup_gwflow(
+  project,
+  gwflow_config = NULL,
+  overwrite = FALSE,
+  aquifer_thickness = 20,
+  conductivity_file = NULL,
+  default_aquifer_k = 1
+)
 }
 \arguments{
 \item{project}{A \code{qswat_project} object returned by
@@ -16,6 +23,24 @@ defaults are used.}
 
 \item{overwrite}{Logical.  If \code{TRUE}, drop and recreate existing gwflow
 tables.  Default \code{FALSE}.}
+
+\item{aquifer_thickness}{Numeric (constant thickness in metres,
+default \code{20}) \strong{or} a file path to a raster (\code{.tif}) or polygon
+shapefile (\code{.shp}) containing aquifer thickness values.  When a
+shapefile is provided it must have a numeric attribute named
+\code{"thickness"} or \code{"Avg_Thick"}.  When a raster is provided the
+mean cell value (in metres) is extracted for each grid cell.}
+
+\item{conductivity_file}{Optional path to a polygon shapefile that
+defines aquifer hydraulic-conductivity zones.  The shapefile must
+contain a numeric column \code{"aquifer_k"} (K in m/day) or
+\code{"logK_Ferr_"} (log10 intrinsic permeability × 100 as stored in
+GLHYMPS).  If \code{NULL} a single default zone is created using
+\code{default_aquifer_k}.}
+
+\item{default_aquifer_k}{Numeric.  Default hydraulic conductivity
+(m/day) for the single zone created when \code{conductivity_file} is
+\code{NULL}.  Default \code{1.0}.}
 }
 \value{
 The \code{project} object with \code{project$use_gwflow} set to \code{TRUE}
@@ -23,9 +48,9 @@ The \code{project} object with \code{project$use_gwflow} set to \code{TRUE}
 }
 \description{
 Creates the gwflow tables in the project SQLite database and populates
-\code{gwflow_base} and \code{gwflow_solutes} from the supplied configuration.
-This is the R equivalent of the gwflow setup performed by the
-\code{GWFlow.createTables()} and related methods in the QGIS plugin's
+them from the supplied configuration and available spatial data.
+This is the R equivalent of the gwflow setup performed by
+\code{GWFlow.createTables()} and the GIS routines in the QGIS plugin's
 \code{gwflow.py}.
 }
 \details{
@@ -35,24 +60,41 @@ Python plugin: \code{gwflow_base}, \code{gwflow_zone}, \code{gwflow_grid},
 \code{gwflow_init_conc}, \code{gwflow_hrucell}, \code{gwflow_fpcell},
 \code{gwflow_rivcell}, \code{gwflow_lsucell}, \code{gwflow_rescell}.
 
-After calling this function the database is ready to receive spatial
-cell data (zones, grid cells, river cells, etc.) that require GIS
-processing not yet automated in rQSWATPlus.
+When the project object contains spatial data (\code{lsu_sf}, \code{hru_sf},
+\code{channels_sf} or \code{streams_sf}) the function also runs the GIS
+processing steps to populate the grid and cell connection tables.
+This happens automatically when called from \code{\link[=qswat_run]{qswat_run()}} with
+\code{use_gwflow = TRUE}.
 
-The following tables are created:
+The following tables are created and, when spatial data are
+available, populated:
 \describe{
-\item{gwflow_base}{Single-row table with global gwflow parameters.}
-\item{gwflow_zone}{Aquifer hydraulic parameter zones.}
-\item{gwflow_grid}{Active grid cells (populated by GIS workflow).}
-\item{gwflow_out_days}{Optional output day pairs (year, jday).}
-\item{gwflow_obs_locs}{Observation cell IDs.}
-\item{gwflow_solutes}{Chemical solute definitions and initial concentrations.}
-\item{gwflow_init_conc}{Per-cell initial solute concentrations.}
-\item{gwflow_hrucell}{HRU-to-cell area mapping.}
-\item{gwflow_fpcell}{Floodplain-cell connections.}
-\item{gwflow_rivcell}{River-cell connections.}
+\item{gwflow_base}{Single-row table with global gwflow parameters
+(always populated; \code{row_count}/\code{col_count} updated after grid
+creation).}
+\item{gwflow_zone}{Aquifer hydraulic parameter zones (one default
+zone, or one row per zone from \code{conductivity_file}).}
+\item{gwflow_grid}{Active grid cells with elevation, aquifer
+thickness, initial head, and zone assignment.}
+\item{gwflow_out_days}{Optional output day pairs (not populated
+automatically; left empty for user entry).}
+\item{gwflow_obs_locs}{Observation cell IDs (not populated
+automatically; left empty for user entry).}
+\item{gwflow_solutes}{Chemical solute definitions and initial
+concentrations (always populated with 10 default solutes).}
+\item{gwflow_init_conc}{Per-cell initial solute concentrations
+(not populated automatically).}
+\item{gwflow_hrucell}{HRU-to-cell area mapping (populated when
+HRU recharge mode is selected and \code{hru_sf} is available).}
+\item{gwflow_fpcell}{Floodplain-cell connections (empty in
+rQSWATPlus because LSUs are 1:1 with subbasins, with no
+floodplain landscape units).}
+\item{gwflow_rivcell}{River-cell connections (channel-cell
+intersection lengths).}
 \item{gwflow_lsucell}{LSU-cell area mapping.}
-\item{gwflow_rescell}{Reservoir-cell connections.}
+\item{gwflow_rescell}{Reservoir-cell connections (populated via
+the \code{gis_water}–\code{gwflow_lsucell} join when water bodies are
+present).}
 }
 }
 \examples{
diff --git a/rQSWATPlus/tests/testthat/test-validate-database.R b/rQSWATPlus/tests/testthat/test-validate-database.R
index 519339e..a2d24de 100644
--- a/rQSWATPlus/tests/testthat/test-validate-database.R
+++ b/rQSWATPlus/tests/testthat/test-validate-database.R
@@ -309,3 +309,238 @@ test_that("qswat_check_database warns when reference table exists but is empty",
               info = paste("Expected warning about empty plants_plt.",
                            "Got warnings:", paste(result$warnings, collapse = "; ")))
 })
+
+# ---- gwflow validation tests ----
+
+# Helper to build a minimal valid database with optional gwflow setup
+.make_minimal_db <- function(db_file, use_gwflow = 0L) {
+  con <- DBI::dbConnect(RSQLite::SQLite(), db_file)
+  rQSWATPlus:::.create_db_tables(con)
+  DBI::dbExecute(con, "INSERT INTO gis_subbasins VALUES (1,100,5,100,5,0,0,500,400,600,0)")
+  DBI::dbExecute(con, "INSERT INTO gis_hrus VALUES (1,1,'AGRL','TX047',1,10,3,0,0,500,50,50,50,100)")
+  DBI::dbExecute(con, "INSERT INTO gis_channels VALUES (1,1,1,1,1000,0.01,1.0,0.5,0,0,0,0)")
+  DBI::dbExecute(con, "INSERT INTO gis_lsus VALUES (1,0,1,1,100,5,100,5,1.0,0.5,0,0,500)")
+  DBI::dbExecute(con, "INSERT INTO gis_routing VALUES (1,'sub','tot',0,'outlet',100)")
+  DBI::dbExecute(con,
+    paste0("INSERT INTO project_config (id, project_name, delineation_done, ",
+           "hrus_done, use_gwflow) VALUES (1,'test',1,1,", use_gwflow, ")"))
+  con
+}
+
+test_that("qswat_check_database skips gwflow checks when use_gwflow=0", {
+  skip_if_not_installed("RSQLite")
+  skip_if_not_installed("DBI")
+
+  db_file <- tempfile(fileext = ".sqlite")
+  on.exit(unlink(db_file), add = TRUE)
+
+  con <- .make_minimal_db(db_file, use_gwflow = 0L)
+  DBI::dbDisconnect(con)
+
+  result <- qswat_check_database(db_file, verbose = FALSE)
+
+  # No gwflow checks should be present
+  gwflow_checks <- result$checks[grepl("^gwflow:", result$checks$check), ]
+  expect_equal(nrow(gwflow_checks), 0,
+               label = "no gwflow check rows when use_gwflow=0")
+  expect_true(result$passed,
+              info = paste("Errors:", paste(result$errors, collapse = "; ")))
+})
+
+test_that("qswat_check_database fails when gwflow enabled but tables missing", {
+  skip_if_not_installed("RSQLite")
+  skip_if_not_installed("DBI")
+
+  db_file <- tempfile(fileext = ".sqlite")
+  on.exit(unlink(db_file), add = TRUE)
+
+  # Enable gwflow but do NOT create the gwflow tables
+  con <- .make_minimal_db(db_file, use_gwflow = 1L)
+  DBI::dbDisconnect(con)
+
+  result <- qswat_check_database(db_file, verbose = FALSE)
+
+  expect_false(result$passed)
+  expect_true(any(grepl("gwflow is enabled but", result$errors)),
+              info = paste("Expected missing-tables error. Got:",
+                           paste(result$errors, collapse = "; ")))
+})
+
+test_that("qswat_check_database fails when gwflow core tables are empty", {
+  skip_if_not_installed("RSQLite")
+  skip_if_not_installed("DBI")
+
+  db_file <- tempfile(fileext = ".sqlite")
+  on.exit(unlink(db_file), add = TRUE)
+
+  # Create all gwflow tables but leave them empty
+  con <- .make_minimal_db(db_file, use_gwflow = 1L)
+  rQSWATPlus:::.create_gwflow_tables(con)
+  DBI::dbDisconnect(con)
+
+  result <- qswat_check_database(db_file, verbose = FALSE)
+
+  expect_false(result$passed)
+  gwflow_fails <- result$errors[grepl("^gwflow", result$checks$check[
+    result$checks$status == "FAIL"])]
+  # At minimum gwflow_base, gwflow_zone, gwflow_grid, gwflow_solutes,
+  # gwflow_rivcell should all fail the non-empty check
+  fail_checks <- result$checks$check[result$checks$status == "FAIL"]
+  expect_true(any(grepl("gwflow:", fail_checks)),
+              info = paste("Expected gwflow check failures. Failing checks:",
+                           paste(fail_checks, collapse = ", ")))
+})
+
+test_that("qswat_check_database passes gwflow checks with fully populated tables (LSU recharge)", {
+  skip_if_not_installed("RSQLite")
+  skip_if_not_installed("DBI")
+
+  db_file <- tempfile(fileext = ".sqlite")
+  on.exit(unlink(db_file), add = TRUE)
+
+  con <- .make_minimal_db(db_file, use_gwflow = 1L)
+  rQSWATPlus:::.create_gwflow_tables(con)
+
+  # Populate gwflow_zone (1 zone)
+  DBI::dbExecute(con,
+    "INSERT INTO gwflow_zone VALUES (1, 10.0, 0.2, 0.005, 0.5)")
+
+  # Populate gwflow_grid (1 active cell)
+  DBI::dbExecute(con,
+    "INSERT INTO gwflow_grid VALUES (1, 1, 1, 500.0, 20.0, 1.0, 495.0, 0)")
+
+  # Populate gwflow_base (recharge=2 → LSU recharge).
+  # Column order: cell_size, row_count, col_count, boundary_conditions,
+  #   recharge, soil_transfer, saturation_excess, external_pumping,
+  #   tile_drainage, reservoir_exchange, wetland_exchange, floodplain_exchange,
+  #   canal_seepage, solute_transport, transport_steps, disp_coef,
+  #   recharge_delay, et_extinction_depth, water_table_depth, river_depth,
+  #   tile_depth, tile_area, tile_k, tile_groups, resbed_thickness, resbed_k,
+  #   wet_thickness, daily_output, annual_output, aa_output,
+  #   daily_output_row, daily_output_col, timestep_balance
+  DBI::dbExecute(con, paste0(
+    "INSERT INTO gwflow_base VALUES (",
+    "200,1,1,1,2,1,1,0,0,1,1,1,0,1,1,5.0,0,1.0,5.0,5.0,1.22,50.0,5.0,",
+    "0,2.0,9.99e-6,0.25,1,1,1,0,0,1.0)"))
+
+  # Populate gwflow_solutes (10 default solutes)
+  solutes <- list(
+    c("no3-n", 1, -0.0001, 3,    "single", 3.0),
+    c("p",     1, 0,       0.05, "single", 0.05),
+    c("so4",   1, 0,       0,    "single", 100),
+    c("ca",    1, 0,       0,    "single", 50),
+    c("mg",    1, 0,       0,    "single", 30),
+    c("na",    1, 0,       0,    "single", 40),
+    c("k",     1, 0,       0,    "single", 1),
+    c("cl",    1, 0,       0,    "single", 25),
+    c("co3",   1, 0,       0,    "single", 1),
+    c("hco3",  1, 0,       0,    "single", 80)
+  )
+  for (s in solutes) {
+    DBI::dbExecute(con, sprintf(
+      "INSERT INTO gwflow_solutes VALUES ('%s',%s,%s,%s,'%s',%s)",
+      s[[1]], s[[2]], s[[3]], s[[4]], s[[5]], s[[6]]))
+  }
+
+  # Populate gwflow_rivcell (1 river-cell connection)
+  DBI::dbExecute(con, "INSERT INTO gwflow_rivcell VALUES (1, 1, 150.0)")
+
+  # Populate gwflow_lsucell (LSU recharge mode)
+  DBI::dbExecute(con, "INSERT INTO gwflow_lsucell VALUES (1, 1, 40000.0)")
+
+  DBI::dbDisconnect(con)
+
+  result <- qswat_check_database(db_file, verbose = FALSE)
+
+  gwflow_fail_checks <- result$checks$check[
+    grepl("^gwflow:", result$checks$check) &
+      result$checks$status == "FAIL"]
+  expect_equal(length(gwflow_fail_checks), 0L,
+               label = paste("gwflow checks should all pass. Failed:",
+                             paste(gwflow_fail_checks, collapse = ", ")))
+})
+
+test_that("qswat_check_database fails when gwflow_lsucell empty in LSU recharge mode", {
+  skip_if_not_installed("RSQLite")
+  skip_if_not_installed("DBI")
+
+  db_file <- tempfile(fileext = ".sqlite")
+  on.exit(unlink(db_file), add = TRUE)
+
+  con <- .make_minimal_db(db_file, use_gwflow = 1L)
+  rQSWATPlus:::.create_gwflow_tables(con)
+
+  DBI::dbExecute(con,
+    "INSERT INTO gwflow_zone VALUES (1, 10.0, 0.2, 0.005, 0.5)")
+  DBI::dbExecute(con,
+    "INSERT INTO gwflow_grid VALUES (1, 1, 1, 500.0, 20.0, 1.0, 495.0, 0)")
+  # recharge=2 → LSU recharge (column 5 of gwflow_base)
+  DBI::dbExecute(con, paste0(
+    "INSERT INTO gwflow_base VALUES (",
+    "200,1,1,1,2,1,1,0,0,1,1,1,0,1,1,5.0,0,1.0,5.0,5.0,1.22,50.0,5.0,",
+    "0,2.0,9.99e-6,0.25,1,1,1,0,0,1.0)"))
+  for (s in list(
+    c("no3-n","single"), c("p","single"), c("so4","single"),
+    c("ca","single"),   c("mg","single"), c("na","single"),
+    c("k","single"),    c("cl","single"), c("co3","single"),
+    c("hco3","single")
+  )) {
+    DBI::dbExecute(con, sprintf(
+      "INSERT INTO gwflow_solutes VALUES ('%s',1,0,0,'%s',1)",
+      s[[1]], s[[2]]))
+  }
+  DBI::dbExecute(con, "INSERT INTO gwflow_rivcell VALUES (1, 1, 150.0)")
+  # gwflow_lsucell intentionally left empty
+
+  DBI::dbDisconnect(con)
+
+  result <- qswat_check_database(db_file, verbose = FALSE)
+
+  expect_false(result$passed)
+  expect_true(any(grepl("gwflow_lsucell is empty", result$errors)),
+              info = paste("Expected lsucell error. Got:",
+                           paste(result$errors, collapse = "; ")))
+})
+
+test_that("qswat_check_database fails when gwflow_hrucell empty in HRU recharge mode", {
+  skip_if_not_installed("RSQLite")
+  skip_if_not_installed("DBI")
+
+  db_file <- tempfile(fileext = ".sqlite")
+  on.exit(unlink(db_file), add = TRUE)
+
+  con <- .make_minimal_db(db_file, use_gwflow = 1L)
+  rQSWATPlus:::.create_gwflow_tables(con)
+
+  DBI::dbExecute(con,
+    "INSERT INTO gwflow_zone VALUES (1, 10.0, 0.2, 0.005, 0.5)")
+  DBI::dbExecute(con,
+    "INSERT INTO gwflow_grid VALUES (1, 1, 1, 500.0, 20.0, 1.0, 495.0, 0)")
+  # recharge=1 → HRU recharge (column 5 of gwflow_base)
+  DBI::dbExecute(con, paste0(
+    "INSERT INTO gwflow_base VALUES (",
+    "200,1,1,1,1,1,1,0,0,1,1,1,0,1,1,5.0,0,1.0,5.0,5.0,1.22,50.0,5.0,",
+    "0,2.0,9.99e-6,0.25,1,1,1,0,0,1.0)"))
+  for (s in list(
+    c("no3-n","single"), c("p","single"), c("so4","single"),
+    c("ca","single"),   c("mg","single"), c("na","single"),
+    c("k","single"),    c("cl","single"), c("co3","single"),
+    c("hco3","single")
+  )) {
+    DBI::dbExecute(con, sprintf(
+      "INSERT INTO gwflow_solutes VALUES ('%s',1,0,0,'%s',1)",
+      s[[1]], s[[2]]))
+  }
+  DBI::dbExecute(con, "INSERT INTO gwflow_rivcell VALUES (1, 1, 150.0)")
+  # gwflow_hrucell intentionally left empty (gwflow_lsucell also empty - not
+  # required for recharge=1)
+
+  DBI::dbDisconnect(con)
+
+  result <- qswat_check_database(db_file, verbose = FALSE)
+
+  expect_false(result$passed)
+  expect_true(any(grepl("gwflow_hrucell is empty", result$errors)),
+              info = paste("Expected hrucell error. Got:",
+                           paste(result$errors, collapse = "; ")))
+})
diff --git a/rQSWATPlus/tests/testthat/test-workflow.R b/rQSWATPlus/tests/testthat/test-workflow.R
index 67191bd..0fb23b9 100644
--- a/rQSWATPlus/tests/testthat/test-workflow.R
+++ b/rQSWATPlus/tests/testthat/test-workflow.R
@@ -159,8 +159,28 @@ test_that("qswat_run with use_gwflow=TRUE initialises gwflow tables", {
   for (tbl in gwflow_tables) {
     expect_true(tbl %in% tbls, label = paste0("table '", tbl, "' exists"))
   }
+
+  # Verify that spatial tables are populated (not empty)
+  base_row <- DBI::dbGetQuery(con, "SELECT row_count, col_count FROM gwflow_base")
+  expect_gt(base_row$row_count, 0L,
+            label = "gwflow_base.row_count updated after grid creation")
+  expect_gt(base_row$col_count, 0L,
+            label = "gwflow_base.col_count updated after grid creation")
+
+  n_zone <- DBI::dbGetQuery(con, "SELECT COUNT(*) AS n FROM gwflow_zone")$n
+  expect_gt(n_zone, 0L, label = "gwflow_zone is populated")
+
+  n_grid <- DBI::dbGetQuery(con, "SELECT COUNT(*) AS n FROM gwflow_grid")$n
+  expect_gt(n_grid, 0L, label = "gwflow_grid is populated")
+
+  n_riv <- DBI::dbGetQuery(con, "SELECT COUNT(*) AS n FROM gwflow_rivcell")$n
+  expect_gt(n_riv, 0L, label = "gwflow_rivcell is populated")
+
+  n_lsu <- DBI::dbGetQuery(con, "SELECT COUNT(*) AS n FROM gwflow_lsucell")$n
+  expect_gt(n_lsu, 0L, label = "gwflow_lsucell is populated")
 })
 
+
 # Full end-to-end integration test: build project and verify SWAT+ Editor readiness
 test_that("example dataset produces a SWAT+ Editor-ready database", {
   skip_if_no_taudem()