diff --git a/gyrokinetic/apps/gk_field.c b/gyrokinetic/apps/gk_field.c
index 8e81b86f67..108c1c82e8 100644
--- a/gyrokinetic/apps/gk_field.c
+++ b/gyrokinetic/apps/gk_field.c
@@ -306,9 +306,6 @@ gk_field_new(struct gkyl_gk *gk, struct gkyl_gyrokinetic_app *app)
     }
   }
 
-  // Initialize biased walls.
-  gk_field_biased_wall_new(app, f);
-
   return f;
 }
 
@@ -432,8 +429,5 @@ gk_field_release(const gkyl_gyrokinetic_app* app, struct gk_field *f)
   // Release energy diagnostics.
   gk_field_energy_release(app, f);
 
-  // Release biased walls.
-  gk_field_biased_wall_release(app, f);
-
   gkyl_free(f);
 }
diff --git a/gyrokinetic/apps/gk_field_biased_wall.c b/gyrokinetic/apps/gk_field_biased_wall.c
deleted file mode 100644
index 3658b29a71..0000000000
--- a/gyrokinetic/apps/gk_field_biased_wall.c
+++ /dev/null
@@ -1,96 +0,0 @@
-#include <gkyl_alloc.h>
-#include <gkyl_array_ops.h>
-#include <gkyl_gyrokinetic_priv.h>
-#include <gkyl_gk_field_priv.h>
-
-void
-gk_field_biased_wall_new(struct gkyl_gyrokinetic_app *app, struct gk_field *f)
-{
-  f->phi_wall_lo = mkarr(app->use_gpu, app->basis.num_basis, app->local_ext.volume);
-  f->has_phi_wall_lo = false;
-  f->phi_wall_lo_evolve = false;
-  if (f->info.phi_wall_lo)
-  {
-    f->has_phi_wall_lo = true;
-    if (f->info.phi_wall_lo_evolve) {
-      f->phi_wall_lo_evolve = f->info.phi_wall_lo_evolve;
-    }
-
-    f->phi_wall_lo_host = f->phi_wall_lo;
-    if (app->use_gpu) {
-      f->phi_wall_lo_host = mkarr(false, f->phi_wall_lo->ncomp, f->phi_wall_lo->size);
-    }
-
-    f->phi_wall_lo_proj = gkyl_eval_on_nodes_new(&app->grid, &app->basis,
-                                                 1, f->info.phi_wall_lo, f->info.phi_wall_lo_ctx);
-
-    // Compute phi_wall_lo at t = 0
-    gkyl_eval_on_nodes_advance(f->phi_wall_lo_proj, 0.0, &app->local_ext, f->phi_wall_lo_host);
-    if (app->use_gpu) { // note: phi_wall_lo_host is same as phi_wall_lo when not on GPUs
-      gkyl_array_copy(f->phi_wall_lo, f->phi_wall_lo_host);
-    }
-  }
-
-  // Set up biased upper wall (same size as electrostatic potential), by default is 0.0
-  f->phi_wall_up = mkarr(app->use_gpu, app->basis.num_basis, app->local_ext.volume);
-  f->has_phi_wall_up = false;
-  f->phi_wall_up_evolve = false;
-  if (f->info.phi_wall_up)
-  {
-    f->has_phi_wall_up = true;
-    if (f->info.phi_wall_up_evolve) {
-      f->phi_wall_up_evolve = f->info.phi_wall_up_evolve;
-    }
-
-    f->phi_wall_up_host = f->phi_wall_up;
-    if (app->use_gpu) {
-      f->phi_wall_up_host = mkarr(false, f->phi_wall_up->ncomp, f->phi_wall_up->size);
-    }
-
-    f->phi_wall_up_proj = gkyl_eval_on_nodes_new(&app->grid, &app->basis,
-                                                 1, f->info.phi_wall_up, f->info.phi_wall_up_ctx);
-
-    // Compute phi_wall_up at t = 0.
-    gkyl_eval_on_nodes_advance(f->phi_wall_up_proj, 0.0, &app->local_ext, f->phi_wall_up_host);
-    if (app->use_gpu) { // Note: phi_wall_up_host is same as phi_wall_up when not on GPUs.
-      gkyl_array_copy(f->phi_wall_up, f->phi_wall_up_host);
-    }
-  }
-}
-
-void
-gk_field_calc_phi_wall(gkyl_gyrokinetic_app *app, struct gk_field *field, double tm)
-{
-  if (field->has_phi_wall_lo && field->phi_wall_lo_evolve) {
-    gkyl_eval_on_nodes_advance(field->phi_wall_lo_proj, tm, &app->local_ext, field->phi_wall_lo_host);
-    if (app->use_gpu) {
-      gkyl_array_copy(field->phi_wall_lo, field->phi_wall_lo_host);
-    }
-  }
-  if (field->has_phi_wall_up && field->phi_wall_up_evolve) {
-    gkyl_eval_on_nodes_advance(field->phi_wall_up_proj, tm, &app->local_ext, field->phi_wall_up_host);
-    if (app->use_gpu) {
-      gkyl_array_copy(field->phi_wall_up, field->phi_wall_up_host);
-    }
-  }
-}
-
-void
-gk_field_biased_wall_release(const struct gkyl_gyrokinetic_app *app, struct gk_field *f)
-{
-  gkyl_array_release(f->phi_wall_lo);
-  if (f->has_phi_wall_lo) {
-    gkyl_eval_on_nodes_release(f->phi_wall_lo_proj);
-    if (app->use_gpu) {
-      gkyl_array_release(f->phi_wall_lo_host);
-    }
-  }
-
-  gkyl_array_release(f->phi_wall_up);
-  if (f->has_phi_wall_up) { 
-    gkyl_eval_on_nodes_release(f->phi_wall_up_proj);
-    if (app->use_gpu) {
-      gkyl_array_release(f->phi_wall_up_host);
-    }
-  }
-}
diff --git a/gyrokinetic/apps/gk_species.c b/gyrokinetic/apps/gk_species.c
index 2b93be63cf..1c1b68d641 100644
--- a/gyrokinetic/apps/gk_species.c
+++ b/gyrokinetic/apps/gk_species.c
@@ -197,7 +197,8 @@ gk_species_rhs_implicit_static(gkyl_gyrokinetic_app *app, struct gk_species *spe
 }
 
 static void
-gk_species_apply_bc_dynamic(gkyl_gyrokinetic_app *app, const struct gk_species *species, struct gkyl_array *f)
+gk_species_apply_bc_dynamic(gkyl_gyrokinetic_app *app, const struct gk_species *species,
+  double tm, struct gkyl_array *f)
 {
   struct timespec wst = gkyl_wall_clock();
   
@@ -210,8 +211,9 @@ gk_species_apply_bc_dynamic(gkyl_gyrokinetic_app *app, const struct gk_species *
 
       switch (species->lower_bc[d].type) {
         case GKYL_BC_GK_SPECIES_SHEATH:
+          gk_species_phi_wall_advance(app, &species->phi_wall_lo, tm);
           gkyl_bc_sheath_gyrokinetic_advance(species->bc_sheath_lo, app->field->phi_smooth, 
-            app->field->phi_wall_lo, f, &app->local);
+            species->phi_wall_lo.phi, f, &app->local);
           break;
         case GKYL_BC_GK_SPECIES_TWISTSHIFT:
           gkyl_bc_twistshift_advance(species->bc_ts_lo, f, f);
@@ -232,8 +234,9 @@ gk_species_apply_bc_dynamic(gkyl_gyrokinetic_app *app, const struct gk_species *
 
       switch (species->upper_bc[d].type) {
         case GKYL_BC_GK_SPECIES_SHEATH:
+          gk_species_phi_wall_advance(app, &species->phi_wall_up, tm);
           gkyl_bc_sheath_gyrokinetic_advance(species->bc_sheath_up, app->field->phi_smooth, 
-            app->field->phi_wall_up, f, &app->local);
+            species->phi_wall_up.phi, f, &app->local);
           break;
         case GKYL_BC_GK_SPECIES_TWISTSHIFT:
           gkyl_bc_twistshift_advance(species->bc_ts_up, f, f);
@@ -266,7 +269,8 @@ gk_species_apply_bc_dynamic(gkyl_gyrokinetic_app *app, const struct gk_species *
 }
 
 static void
-gk_species_apply_bc_static(gkyl_gyrokinetic_app *app, const struct gk_species *species, struct gkyl_array *f)
+gk_species_apply_bc_static(gkyl_gyrokinetic_app *app, const struct gk_species *species,
+  double tm, struct gkyl_array *f)
 {
   // do nothing
 }
@@ -677,6 +681,7 @@ gk_species_release_dynamic(const gkyl_gyrokinetic_app* app, const struct gk_spec
   for (int d=0; d<app->cdim; ++d) {
     if (s->lower_bc[d].type == GKYL_BC_GK_SPECIES_SHEATH) {
       gkyl_bc_sheath_gyrokinetic_release(s->bc_sheath_lo);
+      gk_species_phi_wall_release(app, &s->phi_wall_lo);
     }
     else if (s->lower_bc[d].type == GKYL_BC_GK_SPECIES_TWISTSHIFT) { 
       gkyl_bc_twistshift_release(s->bc_ts_lo);
@@ -690,6 +695,7 @@ gk_species_release_dynamic(const gkyl_gyrokinetic_app* app, const struct gk_spec
     
     if (s->upper_bc[d].type == GKYL_BC_GK_SPECIES_SHEATH) {
       gkyl_bc_sheath_gyrokinetic_release(s->bc_sheath_up);
+      gk_species_phi_wall_release(app, &s->phi_wall_up);
     }
     else if (s->upper_bc[d].type == GKYL_BC_GK_SPECIES_TWISTSHIFT) {
       gkyl_bc_twistshift_release(s->bc_ts_up);
@@ -877,6 +883,7 @@ gk_species_init_dynamic(struct gkyl_gk *gk_app_inp, struct gkyl_gyrokinetic_app
       struct gkyl_range *sol_ghost = gk_app_inp->geometry.has_LCFS? &gks->local_lower_ghost_par_sol : &gks->local_lower_ghost[d];
       gks->bc_sheath_lo = gkyl_bc_sheath_gyrokinetic_new(d, GKYL_LOWER_EDGE, gks->basis_on_dev, 
         sol_skin, sol_ghost, gks->vel_map, cdim, 2.0*(gks->info.charge/gks->info.mass), app->use_gpu);
+      gk_species_phi_wall_init(app, &gks->lower_bc[d], &gks->phi_wall_lo);
     }
     else if (gks->lower_bc[d].type == GKYL_BC_GK_SPECIES_TWISTSHIFT) {
       assert(cdim == 3);
@@ -937,6 +944,7 @@ gk_species_init_dynamic(struct gkyl_gk *gk_app_inp, struct gkyl_gyrokinetic_app
       struct gkyl_range *sol_ghost = gk_app_inp->geometry.has_LCFS? &gks->local_upper_ghost_par_sol : &gks->local_upper_ghost[d];
       gks->bc_sheath_up = gkyl_bc_sheath_gyrokinetic_new(d, GKYL_UPPER_EDGE, gks->basis_on_dev, 
         sol_skin, sol_ghost, gks->vel_map, cdim, 2.0*(gks->info.charge/gks->info.mass), app->use_gpu);
+      gk_species_phi_wall_init(app, &gks->upper_bc[d], &gks->phi_wall_up);
     }
     else if (gks->upper_bc[d].type == GKYL_BC_GK_SPECIES_TWISTSHIFT) {
       assert(cdim == 3);
@@ -1497,7 +1505,7 @@ gk_species_init(struct gkyl_gk *gk_app_inp, struct gkyl_gyrokinetic_app *app, st
       gks->upper_bc[d].type = GKYL_BC_GK_SPECIES_PERIODIC;
     }
   }
- 
+
   // Species properties metadata.
   struct gkyl_msgpack_map_elem io_meta_sprop[] = {
     { .key = "mass", .elem_type = GKYL_MP_DOUBLE, .dval = gks->info.mass },
@@ -1935,9 +1943,10 @@ gk_species_copy_range(struct gk_species *species, struct gkyl_array *out,
 }
 
 void
-gk_species_apply_bc(gkyl_gyrokinetic_app *app, const struct gk_species *species, struct gkyl_array *f)
+gk_species_apply_bc(gkyl_gyrokinetic_app *app, const struct gk_species *species,
+  double tm, struct gkyl_array *f)
 {
-  species->bc_func(app, species, f);
+  species->bc_func(app, species, tm, f);
 }
 
 void
diff --git a/gyrokinetic/apps/gk_species_fdot_multiplier.c b/gyrokinetic/apps/gk_species_fdot_multiplier.c
index 7992e1f6d1..cc9d74a13d 100644
--- a/gyrokinetic/apps/gk_species_fdot_multiplier.c
+++ b/gyrokinetic/apps/gk_species_fdot_multiplier.c
@@ -144,11 +144,11 @@ gk_species_fdot_multiplier_advance_loss_cone_mult(gkyl_gyrokinetic_app *app,
   int zdim = app->cdim - 1;
   if (gks->lower_bc[zdim].type == GKYL_BC_GK_SPECIES_SHEATH) {
     gkyl_comm_array_allgather(app->comm, &app->local, &app->global,
-      app->field->phi_wall_lo, fdmul->phi_wall_lo_global);
+      gks->phi_wall_lo.phi, fdmul->phi_wall_lo_global);
   }
   if (gks->upper_bc[zdim].type == GKYL_BC_GK_SPECIES_SHEATH) {
     gkyl_comm_array_allgather(app->comm, &app->local, &app->global,
-      app->field->phi_wall_up, fdmul->phi_wall_up_global);
+      gks->phi_wall_up.phi, fdmul->phi_wall_up_global);
   }
   gkyl_loss_cone_mask_gyrokinetic_advance(fdmul->lcm_proj_op, &gks->local, &app->global,
     fdmul->bmag_global, fdmul->phi_global, fdmul->phi_wall_lo_global,
@@ -338,6 +338,8 @@ gk_species_fdot_multiplier_init_comp(gkyl_gyrokinetic_app *app, struct gk_specie
     }
     else if (fdmul->type == GKYL_GK_FDOT_MULTIPLIER_LOSS_CONE) {
       // Operator that projects the loss cone mask.
+      fdmul->phi_wall_lo_global = 0;
+      fdmul->phi_wall_up_global = 0;
       int zdim = app->cdim - 1;
       struct gkyl_loss_cone_mask_gyrokinetic_inp inp_proj = {
         .conf_basis = &app->basis,
@@ -354,10 +356,12 @@ gk_species_fdot_multiplier_init_comp(gkyl_gyrokinetic_app *app, struct gk_specie
       fdmul->bmag_global = mkarr(app->use_gpu, app->gk_geom->geo_corn.bmag->ncomp,
         app->global_ext.volume);
       fdmul->phi_global = mkarr(app->use_gpu, app->basis.num_basis, app->global_ext.volume);
-      fdmul->phi_wall_lo_global = mkarr(app->use_gpu, app->basis.num_basis,
-        app->global_ext.volume);
-      fdmul->phi_wall_up_global = mkarr(app->use_gpu, app->basis.num_basis,
-        app->global_ext.volume);
+      if (inp_proj.lower_boundary == GKYL_LOSS_CONE_BC_SHEATH)
+        fdmul->phi_wall_lo_global = mkarr(app->use_gpu, app->basis.num_basis,
+          app->global_ext.volume);
+      if (inp_proj.upper_boundary == GKYL_LOSS_CONE_BC_SHEATH)
+        fdmul->phi_wall_up_global = mkarr(app->use_gpu, app->basis.num_basis,
+          app->global_ext.volume);
 
       gkyl_comm_array_allgather(app->comm, &app->local, &app->global, app->gk_geom->geo_corn.bmag,
         fdmul->bmag_global);
@@ -540,8 +544,10 @@ gk_species_fdot_multiplier_release_comp(const struct gkyl_gyrokinetic_app *app,
     gkyl_array_release(fdmul->buffer);
     gkyl_array_release(fdmul->bmag_global);
     gkyl_array_release(fdmul->phi_global);
-    gkyl_array_release(fdmul->phi_wall_lo_global);
-    gkyl_array_release(fdmul->phi_wall_up_global);
+    if (fdmul->phi_wall_lo_global)
+      gkyl_array_release(fdmul->phi_wall_lo_global);
+    if (fdmul->phi_wall_up_global)
+      gkyl_array_release(fdmul->phi_wall_up_global);
     gkyl_loss_cone_mask_gyrokinetic_release(fdmul->lcm_proj_op);
   }
   else if (fdmul->type == GKYL_GK_FDOT_MULTIPLIER_CONSTANT) {
diff --git a/gyrokinetic/apps/gk_species_phi_wall.c b/gyrokinetic/apps/gk_species_phi_wall.c
new file mode 100644
index 0000000000..c5f211947b
--- /dev/null
+++ b/gyrokinetic/apps/gk_species_phi_wall.c
@@ -0,0 +1,64 @@
+#include <gkyl_array_ops.h>
+#include <gkyl_gyrokinetic_priv.h>
+
+static void
+gk_species_phi_wall_advance_disabled(gkyl_gyrokinetic_app *app,
+  const struct gk_species_wall_potential *wall, double tm)
+{
+}
+
+static void
+gk_species_phi_wall_advance_enabled(gkyl_gyrokinetic_app *app,
+  const struct gk_species_wall_potential *wall, double tm)
+{
+  gkyl_eval_on_nodes_advance(wall->projector, tm, &app->local_ext, wall->phi_host);
+  if (app->use_gpu)
+    gkyl_array_copy(wall->phi, wall->phi_host);
+}
+
+void
+gk_species_phi_wall_init(gkyl_gyrokinetic_app *app,
+  const struct gkyl_gyrokinetic_bc *bc, struct gk_species_wall_potential *wall)
+{
+  *wall = (struct gk_species_wall_potential) {
+    .advance_func = gk_species_phi_wall_advance_disabled,
+  };
+  if (bc->type != GKYL_BC_GK_SPECIES_SHEATH)
+    return;
+
+  wall->phi = mkarr(app->use_gpu, app->basis.num_basis, app->local_ext.volume);
+  gkyl_array_clear(wall->phi, 0.0);
+
+  wall->phi_host = wall->phi;
+  if (bc->aux_profile) {
+    if (app->use_gpu)
+      wall->phi_host = mkarr(false, wall->phi->ncomp, wall->phi->size);
+
+    wall->projector = gkyl_eval_on_nodes_new(&app->grid, &app->basis,
+      1, bc->aux_profile, bc->aux_ctx);
+    wall->advance_func = gk_species_phi_wall_advance_enabled;
+    gk_species_phi_wall_advance(app, wall, 0.0);
+  }
+}
+
+void
+gk_species_phi_wall_advance(gkyl_gyrokinetic_app *app,
+  const struct gk_species_wall_potential *wall, double tm)
+{
+  wall->advance_func(app, wall, tm);
+}
+
+void
+gk_species_phi_wall_release(const gkyl_gyrokinetic_app *app,
+  const struct gk_species_wall_potential *wall)
+{
+  if (!wall->phi)
+    return;
+
+  gkyl_array_release(wall->phi);
+  if (wall->projector) {
+    gkyl_eval_on_nodes_release(wall->projector);
+    if (app->use_gpu)
+      gkyl_array_release(wall->phi_host);
+  }
+}
diff --git a/gyrokinetic/apps/gkyl_gk_field_priv.h b/gyrokinetic/apps/gkyl_gk_field_priv.h
index da502d8cd9..eee5b8281b 100644
--- a/gyrokinetic/apps/gkyl_gk_field_priv.h
+++ b/gyrokinetic/apps/gkyl_gk_field_priv.h
@@ -10,7 +10,7 @@
  * This header contains forward declarations for internal field-related
  * functions used by the gyrokinetic application. These functions handle
  * FEM projections, field solves, boundary conditions, energy diagnostics,
- * FLR corrections, polarization potentials, and biased wall configurations.
+ * FLR corrections and polarization potentials.
  */
 
 /** FEM Initialization Functions **/
@@ -147,35 +147,3 @@ gk_field_invert_flr(gkyl_gyrokinetic_app *app, struct gk_field *field,
 void
 gk_field_invert_flr_none(gkyl_gyrokinetic_app *app, struct gk_field *field,
   struct gkyl_array *phi);
-
-/** Biased Wall Functions **/
-
-/**
- * Initialize biased wall boundary condition objects.
- * Sets up machinery for applying time-dependent wall potentials.
- *
- * @param app Gyrokinetic application object.
- * @param f Field object to initialize biased wall for.
- */
-void
-gk_field_biased_wall_new(struct gkyl_gyrokinetic_app *app, struct gk_field *f);
-
-/**
- * Release biased wall resources.
- *
- * @param app Gyrokinetic application object.
- * @param f Field object whose biased wall resources are to be released.
- */
-void
-gk_field_biased_wall_release(const struct gkyl_gyrokinetic_app *app, struct gk_field *f);
-
-/**
- * Calculate and apply the wall potential at the current simulation time.
- * Updates phi at the wall boundaries based on the biased wall configuration.
- *
- * @param app Gyrokinetic application object.
- * @param field Field object containing the potential.
- * @param tm Current simulation time.
- */
-void
-gk_field_calc_phi_wall(gkyl_gyrokinetic_app *app, struct gk_field *field, double tm);
diff --git a/gyrokinetic/apps/gkyl_gyrokinetic.h b/gyrokinetic/apps/gkyl_gyrokinetic.h
index 65c8da144b..31938b0e63 100644
--- a/gyrokinetic/apps/gkyl_gyrokinetic.h
+++ b/gyrokinetic/apps/gkyl_gyrokinetic.h
@@ -213,7 +213,8 @@ struct gkyl_gyrokinetic_bc {
   enum gkyl_edge_loc edge; // Which edge this BC is for.
   enum gkyl_gyrokinetic_bc_type type; // BC type flag.
   double value[3]; // Meaning depends on type.
-  void (*aux_profile)(double t, const double *xn, double *fout, void *ctx); // Auxiliary function (e.g. wall potential).
+  // Auxiliary profile. For a species sheath BC, this is the material-wall potential.
+  void (*aux_profile)(double t, const double *xn, double *fout, void *ctx);
   void *aux_ctx; // Context for aux_profile.
   struct gkyl_gyrokinetic_projection projection; // Projection object input (e.g. for FIXED_FUNC).
   struct gkyl_gyrokinetic_emission_inp emission; 
@@ -569,16 +570,6 @@ struct gkyl_gyrokinetic_field {
   void (*init_field_profile)(double t, const double *xn, double *out, void *ctx);
   void *init_field_profile_ctx;
 
-  void *phi_wall_lo_ctx; // context for biased wall potential on lower wall
-  // pointer to biased wall potential on lower wall function
-  void (*phi_wall_lo)(double t, const double *xn, double *phi_wall_lo_out, void *ctx);
-  bool phi_wall_lo_evolve; // set to true if biased wall potential on lower wall function is time dependent  
-
-  void *phi_wall_up_ctx; // context for biased wall potential on upper wall
-  // pointer to biased wall potential on upper wall function
-  void (*phi_wall_up)(double t, const double *xn, double *phi_wall_up_out, void *ctx);
-  bool phi_wall_up_evolve; // set to true if biased wall potential on upper wall function is time dependent  
-
   struct gkyl_poisson_bias_line_list *bias_line_list; // Biased lines constraining the solution.
 };
 
diff --git a/gyrokinetic/apps/gkyl_gyrokinetic_multib.h b/gyrokinetic/apps/gkyl_gyrokinetic_multib.h
index e47288d6f1..3706894eab 100644
--- a/gyrokinetic/apps/gkyl_gyrokinetic_multib.h
+++ b/gyrokinetic/apps/gkyl_gyrokinetic_multib.h
@@ -110,15 +110,6 @@ struct gkyl_gyrokinetic_multib_field_pb {
 
   bool time_rate_diagnostics; // Writes the time rate of change of field energy.
 
-  void *phi_wall_lo_ctx; // context for biased wall potential on lower wall
-  // pointer to biased wall potential on lower wall function
-  void (*phi_wall_lo)(double t, const double *xn, double *phi_wall_lo_out, void *ctx);
-  bool phi_wall_lo_evolve; // set to true if biased wall potential on lower wall function is time dependent  
-
-  void *phi_wall_up_ctx; // context for biased wall potential on upper wall
-  // pointer to biased wall potential on upper wall function
-  void (*phi_wall_up)(double t, const double *xn, double *phi_wall_up_out, void *ctx);
-  bool phi_wall_up_evolve; // set to true if biased wall potential on upper wall function is time dependent  
 };
 
 // Field input 
diff --git a/gyrokinetic/apps/gkyl_gyrokinetic_priv.h b/gyrokinetic/apps/gkyl_gyrokinetic_priv.h
index 4621a1f499..fb15362ffe 100644
--- a/gyrokinetic/apps/gkyl_gyrokinetic_priv.h
+++ b/gyrokinetic/apps/gkyl_gyrokinetic_priv.h
@@ -999,6 +999,15 @@ struct gk_positivity {
     struct gk_positivity *pos);
 };
 
+// Projected material-wall potential associated with a species sheath BC.
+struct gk_species_wall_potential {
+  struct gkyl_array *phi;
+  struct gkyl_array *phi_host;
+  gkyl_eval_on_nodes *projector;
+  void (*advance_func)(gkyl_gyrokinetic_app *app,
+    const struct gk_species_wall_potential *wall, double tm);
+};
+
 // Species data.
 struct gk_species {
   struct gkyl_gyrokinetic_species info; // Input data.
@@ -1061,6 +1070,8 @@ struct gk_species {
 
   // Boundary conditions on lower/upper edges in each direction.
   struct gkyl_gyrokinetic_bc lower_bc[GKYL_MAX_CDIM], upper_bc[GKYL_MAX_CDIM];
+  // Wall potentials supplied by the parallel sheath BCs.
+  struct gk_species_wall_potential phi_wall_lo, phi_wall_up;
   // gyrokinetic sheath boundary conditions
   struct gkyl_bc_sheath_gyrokinetic *bc_sheath_lo;
   struct gkyl_bc_sheath_gyrokinetic *bc_sheath_up;
@@ -1126,7 +1137,7 @@ struct gk_species {
   double (*rhs_implicit_func)(gkyl_gyrokinetic_app *app, struct gk_species *species,
     const struct gkyl_array *fin, struct gkyl_array *rhs, struct gkyl_array **bflux_moms, double dt);
   void (*bc_func)(gkyl_gyrokinetic_app *app, const struct gk_species *species,
-    struct gkyl_array *f);
+    double tm, struct gkyl_array *f);
   void (*release_func)(const gkyl_gyrokinetic_app* app, const struct gk_species *s);
   void (*step_f_func)(struct gkyl_array* out, double dt, const struct gkyl_array* inp); 
   void (*combine_func)(struct gkyl_array *out, double c1,
@@ -1367,18 +1378,6 @@ struct gk_field {
   gkyl_dynvec integ_energy_dot; // d/dt of integrated energy components.
   bool is_first_energy_dot_write_call; // flag for d(energy)/dt dynvec written first time
 
-  bool has_phi_wall_lo; // flag to indicate there is biased wall potential on lower wall
-  bool phi_wall_lo_evolve; // flag to indicate biased wall potential on lower wall is time dependent
-  struct gkyl_array *phi_wall_lo; // biased wall potential on lower wall
-  struct gkyl_array *phi_wall_lo_host; // host copy for use in IO and projecting
-  gkyl_eval_on_nodes *phi_wall_lo_proj; // projector for biased wall potential on lower wall 
-
-  bool has_phi_wall_up; // flag to indicate there is biased wall potential on upper wall
-  bool phi_wall_up_evolve; // flag to indicate biased wall potential on upper wall is time dependent
-  struct gkyl_array *phi_wall_up; // biased wall potential on upper wall
-  struct gkyl_array *phi_wall_up_host; // host copy for use in IO and projecting
-  gkyl_eval_on_nodes *phi_wall_up_proj; // projector for biased wall potential on upper wall 
-
   // Pointer to function that computes the time rate of change of the energy.
   void (*calc_energy_dt_func)(gkyl_gyrokinetic_app *app, const struct gk_field *field, double dt, double *energy_reduced);
 
@@ -3243,9 +3242,11 @@ void gk_species_apply_pos_shift(gkyl_gyrokinetic_app* app, struct gk_species *gk
  *
  * @param app gyrokinetic app object.
  * @param species Pointer to species.
+ * @param tm Time at which to apply the BCs.
  * @param f Field to apply BCs.
  */
-void gk_species_apply_bc(gkyl_gyrokinetic_app *app, const struct gk_species *species, struct gkyl_array *f);
+void gk_species_apply_bc(gkyl_gyrokinetic_app *app, const struct gk_species *species,
+  double tm, struct gkyl_array *f);
 
 /**
  * Fill stat object in app with collision timers.
@@ -3339,6 +3340,35 @@ gk_species_calc_int_mom_dt(gkyl_gyrokinetic_app* app, struct gk_species *gks, do
  */
 void gk_species_release(const gkyl_gyrokinetic_app* app, const struct gk_species *s);
 
+/**
+ * Initialize a wall potential owned by a species sheath BC.
+ *
+ * @param app Gyrokinetic app object.
+ * @param bc Species sheath BC input.
+ * @param wall Wall-potential object.
+ */
+void gk_species_phi_wall_init(gkyl_gyrokinetic_app *app,
+  const struct gkyl_gyrokinetic_bc *bc, struct gk_species_wall_potential *wall);
+
+/**
+ * Evaluate a species sheath wall-potential profile at the requested time.
+ *
+ * @param app Gyrokinetic app object.
+ * @param wall Wall-potential object.
+ * @param tm Time at which to evaluate the profiles.
+ */
+void gk_species_phi_wall_advance(gkyl_gyrokinetic_app *app,
+  const struct gk_species_wall_potential *wall, double tm);
+
+/**
+ * Release a species sheath wall-potential resource.
+ *
+ * @param app Gyrokinetic app object.
+ * @param wall Wall-potential object.
+ */
+void gk_species_phi_wall_release(const gkyl_gyrokinetic_app *app,
+  const struct gk_species_wall_potential *wall);
+
 /** gk_neut_species_moment API */
 
 /**
@@ -4134,15 +4164,6 @@ void gk_neut_species_release(const gkyl_gyrokinetic_app* app, const struct gk_ne
  */
 struct gk_field* gk_field_new(struct gkyl_gk *gk, struct gkyl_gyrokinetic_app *app);
 
-/**
- * Compute biased wall potentials.
- *
- * @param app gyrokinetic app object.
- * @param field Pointer to field.
- * @param tm Time to compute biased wall potentials at.
- */
-void gk_field_calc_phi_wall(gkyl_gyrokinetic_app *app, struct gk_field *field, double tm);
-
 /**
  * Accumulate charge density for Poisson solve.
  *
diff --git a/gyrokinetic/apps/gyrokinetic.c b/gyrokinetic/apps/gyrokinetic.c
index 526aeebc67..4edc456b48 100644
--- a/gyrokinetic/apps/gyrokinetic.c
+++ b/gyrokinetic/apps/gyrokinetic.c
@@ -690,12 +690,6 @@ gyrokinetic_calc_field_enabled(gkyl_gyrokinetic_app* app, double tcurr,
   // Compute electrostatic potential from gyrokinetic Poisson's equation.
   gk_field_accumulate_rho_c(app, app->field, fin, bflux);
 
-  // Compute biased wall potential if present and time-dependent.
-  // Note: biased wall potential use eval_on_nodes. 
-  // so does copy to GPU every call if app->use_gpu = true.
-  if (app->field->phi_wall_lo_evolve || app->field->phi_wall_up_evolve)
-    gk_field_calc_phi_wall(app, app->field, tcurr);
-
   // Solve the field equation.
   gk_field_rhs(app, app->field);
   app->stat.field_tm += gkyl_time_diff_now_sec(wtm);
@@ -999,7 +993,7 @@ gyrokinetic_calc_field_and_apply_bc(gkyl_gyrokinetic_app* app, double tcurr,
   // Apply boundary conditions.
   struct timespec wst = gkyl_wall_clock();
   for (int i=0; i<app->num_species; ++i) {
-    gk_species_apply_bc(app, &app->species[i], distf[i]);
+    gk_species_apply_bc(app, &app->species[i], tcurr, distf[i]);
   }
   for (int i=0; i<app->num_neut_species; ++i) {
     gk_neut_species_apply_bc(app, &app->neut_species[i], distf_neut[i]);
@@ -1115,7 +1109,7 @@ gkyl_gyrokinetic_app_apply_ic(gkyl_gyrokinetic_app* app, double t0)
 
   // Apply boundary conditions.
   for (int i=0; i<app->num_species; ++i) {
-    gk_species_apply_bc(app, &app->species[i], distf[i]);
+    gk_species_apply_bc(app, &app->species[i], t0, distf[i]);
   }
   for (int i=0; i<app->num_neut_species; ++i) {
     if (!app->neut_species[i].info.is_static) {
@@ -3506,7 +3500,7 @@ gkyl_gyrokinetic_app_read_from_frame(gkyl_gyrokinetic_app *app, int frame)
 
     // Apply boundary conditions.
     for (int i=0; i<app->num_species; ++i) {
-      gk_species_apply_bc(app, &app->species[i], distf[i]);
+      gk_species_apply_bc(app, &app->species[i], rstat.stime, distf[i]);
     }
     for (int i=0; i<app->num_neut_species; ++i) {
       gk_neut_species_apply_bc(app, &app->neut_species[i], distf_neut[i]);
diff --git a/gyrokinetic/apps/gyrokinetic_multib.c b/gyrokinetic/apps/gyrokinetic_multib.c
index de0352dee2..792032fb92 100644
--- a/gyrokinetic/apps/gyrokinetic_multib.c
+++ b/gyrokinetic/apps/gyrokinetic_multib.c
@@ -193,10 +193,8 @@ singleb_app_new_solver(const struct gkyl_gyrokinetic_multib *mbinp, int bid,
     species_inp.polarization_density = sp_pb->polarization_density;  
 
     // By default, skip BCs altogether.
-    for (int i=0; i<2*GKYL_MAX_CDIM; i++) {
+    for (int i=0; i<2*GKYL_MAX_CDIM; i++)
       species_inp.bcs[i].type = GKYL_BC_GK_SKIP;
-      species_inp.bcs[i].type = GKYL_BC_GK_SKIP;
-    }
 
     int pardir = cdim-1;
     int num_below = gkyl_multib_conn_get_num_connected(mbapp->block_topo,
@@ -285,10 +283,8 @@ singleb_app_new_solver(const struct gkyl_gyrokinetic_multib *mbinp, int bid,
     neut_species_inp.source = nsp_pb->source;
 
     // By default, skip BCs altogether.
-    for (int i=0; i<2*GKYL_MAX_CDIM; i++) {
-      neut_species_inp.bcs[i].type = GKYL_BC_GK_SKIP;
+    for (int i=0; i<2*GKYL_MAX_CDIM; i++)
       neut_species_inp.bcs[i].type = GKYL_BC_GK_SKIP;
-    }
 
     // Set species physical BCs.
     int bc_count_nsp[num_blocks];
@@ -339,7 +335,7 @@ singleb_app_new_solver(const struct gkyl_gyrokinetic_multib *mbinp, int bid,
         field_inp.poisson_bcs[i].value[k] = -1.0e3;
     }
   }
-  
+
   for (int i=0; i<fld->num_physical_bcs; i++) { 
     if (bid == fld->bcs[i].bidx) {
       struct gkyl_gyrokinetic_bc *bc_curr = gk_fetch_bc_with_dir_edge(field_inp.poisson_bcs, 2*cdim, fld->bcs[i].dir, fld->bcs[i].edge);
@@ -361,7 +357,7 @@ singleb_app_new_solver(const struct gkyl_gyrokinetic_multib *mbinp, int bid,
       }
     }
   }
-  
+
   if (!fld->duplicate_across_blocks) {
     for (int i=0; i<num_blocks; ++i) {
       if (bid == fld->blocks[i].block_id) {
@@ -374,15 +370,6 @@ singleb_app_new_solver(const struct gkyl_gyrokinetic_multib *mbinp, int bid,
   field_inp.polarization_bmag = fld_pb->polarization_bmag ? fld_pb->polarization_bmag : mbapp->bmag_ref;
   field_inp.kperpSq = fld_pb->kperpSq;
   field_inp.time_rate_diagnostics = fld_pb->time_rate_diagnostics; 
-  
-  field_inp.phi_wall_lo_ctx = fld_pb->phi_wall_lo_ctx; 
-  field_inp.phi_wall_lo = fld_pb->phi_wall_lo; 
-  field_inp.phi_wall_lo_evolve = fld_pb->phi_wall_lo_evolve; 
-  
-  field_inp.phi_wall_up_ctx = fld_pb->phi_wall_up_ctx; 
-  field_inp.phi_wall_up = fld_pb->phi_wall_up; 
-  field_inp.phi_wall_up_evolve = fld_pb->phi_wall_up_evolve;   
-  
 
   // Copy field input into app input.
   memcpy(&app_inp.field, &field_inp, sizeof(struct gkyl_gyrokinetic_field));  
@@ -932,7 +919,7 @@ gyrokinetic_multib_apply_bc(struct gkyl_gyrokinetic_multib_app* app, double tcur
     int li_charged = b * app->num_species;
     int li_neut = b * app->num_neut_species;
     for (int i=0; i<app->num_species; ++i) {
-      gk_species_apply_bc(sbapp, &sbapp->species[i], distf[li_charged+i]);
+      gk_species_apply_bc(sbapp, &sbapp->species[i], tcurr, distf[li_charged+i]);
     }
     for (int i=0; i<app->num_neut_species; ++i) {
       gk_neut_species_apply_bc(sbapp, &sbapp->neut_species[i], distf_neut[li_neut+i]);
diff --git a/gyrokinetic/unit/ctest_gk_species_wall_potential.c b/gyrokinetic/unit/ctest_gk_species_wall_potential.c
new file mode 100644
index 0000000000..7b10f8abdb
--- /dev/null
+++ b/gyrokinetic/unit/ctest_gk_species_wall_potential.c
@@ -0,0 +1,151 @@
+#include <acutest.h>
+
+#include <math.h>
+
+#include <gkyl_gyrokinetic_priv.h>
+
+struct wall_profile_ctx {
+  double offset;
+  double time_slope;
+  double space_slope;
+};
+
+static void
+wall_profile(double t, const double *xn, double *fout, void *ctx)
+{
+  const struct wall_profile_ctx *profile = ctx;
+  fout[0] = profile->offset + profile->time_slope*t + profile->space_slope*xn[0];
+}
+
+static gkyl_gyrokinetic_app
+make_app(void)
+{
+  gkyl_gyrokinetic_app app = {
+    .cdim = 1,
+    .use_gpu = false,
+  };
+  double lower[] = { -1.0 }, upper[] = { 1.0 };
+  int cells[] = { 4 }, ghost[] = { 1 };
+
+  gkyl_rect_grid_init(&app.grid, app.cdim, lower, upper, cells);
+  gkyl_cart_modal_serendip(&app.basis, app.cdim, 1);
+  gkyl_create_grid_ranges(&app.grid, ghost, &app.local_ext, &app.local);
+
+  return app;
+}
+
+static void
+check_profile(const gkyl_gyrokinetic_app *app, const struct gkyl_array *phi,
+  const struct wall_profile_ctx *profile, double tm)
+{
+  struct gkyl_range_iter iter;
+  gkyl_range_iter_init(&iter, &app->local_ext);
+  while (gkyl_range_iter_next(&iter)) {
+    long lidx = gkyl_range_idx(&app->local_ext, iter.idx);
+    const double *phi_c = gkyl_array_cfetch(phi, lidx);
+    double xc[GKYL_MAX_DIM];
+    gkyl_rect_grid_cell_center(&app->grid, iter.idx, xc);
+
+    const double eta[] = { -1.0, 0.0, 1.0 };
+    for (int i=0; i<3; ++i) {
+      double x = xc[0] + 0.5*app->grid.dx[0]*eta[i];
+      double expected = profile->offset + profile->time_slope*tm + profile->space_slope*x;
+      TEST_CHECK(fabs(app->basis.eval_expand(&eta[i], phi_c)-expected) < 1e-12);
+    }
+  }
+}
+
+static void
+test_species_wall_profiles(void)
+{
+  gkyl_gyrokinetic_app app = make_app();
+  struct wall_profile_ctx lower_ctx = {
+    .offset = 1.5,
+    .time_slope = 2.0,
+    .space_slope = -0.25,
+  };
+  struct wall_profile_ctx upper_ctx = {
+    .offset = -0.5,
+    .time_slope = -1.0,
+    .space_slope = 0.75,
+  };
+  struct gk_species species = {
+    .lower_bc[0] = {
+      .type = GKYL_BC_GK_SPECIES_SHEATH,
+      .aux_profile = wall_profile,
+      .aux_ctx = &lower_ctx,
+    },
+    .upper_bc[0] = {
+      .type = GKYL_BC_GK_SPECIES_SHEATH,
+      .aux_profile = wall_profile,
+      .aux_ctx = &upper_ctx,
+    },
+  };
+
+  gk_species_phi_wall_init(&app, &species.lower_bc[0], &species.phi_wall_lo);
+  gk_species_phi_wall_init(&app, &species.upper_bc[0], &species.phi_wall_up);
+  TEST_ASSERT(species.phi_wall_lo.phi != 0);
+  TEST_ASSERT(species.phi_wall_up.phi != 0);
+  TEST_ASSERT(species.phi_wall_lo.projector != 0);
+  TEST_ASSERT(species.phi_wall_up.projector != 0);
+  TEST_ASSERT(species.phi_wall_lo.advance_func != 0);
+  TEST_ASSERT(species.phi_wall_up.advance_func != 0);
+  check_profile(&app, species.phi_wall_lo.phi, &lower_ctx, 0.0);
+  check_profile(&app, species.phi_wall_up.phi, &upper_ctx, 0.0);
+
+  double tm = 0.625;
+  gk_species_phi_wall_advance(&app, &species.phi_wall_lo, tm);
+  check_profile(&app, species.phi_wall_lo.phi, &lower_ctx, tm);
+  check_profile(&app, species.phi_wall_up.phi, &upper_ctx, 0.0);
+  gk_species_phi_wall_advance(&app, &species.phi_wall_up, tm);
+  check_profile(&app, species.phi_wall_up.phi, &upper_ctx, tm);
+
+  gk_species_phi_wall_release(&app, &species.phi_wall_lo);
+  gk_species_phi_wall_release(&app, &species.phi_wall_up);
+}
+
+static void
+test_grounded_and_non_sheath_walls(void)
+{
+  gkyl_gyrokinetic_app app = make_app();
+  struct gk_species species = {
+    .lower_bc[0] = {
+      .type = GKYL_BC_GK_SPECIES_SHEATH,
+    },
+    .upper_bc[0] = {
+      .type = GKYL_BC_GK_SPECIES_REFLECT,
+    },
+  };
+
+  gk_species_phi_wall_init(&app, &species.lower_bc[0], &species.phi_wall_lo);
+  gk_species_phi_wall_init(&app, &species.upper_bc[0], &species.phi_wall_up);
+  TEST_ASSERT(species.phi_wall_lo.phi != 0);
+  TEST_CHECK(species.phi_wall_lo.phi_host == species.phi_wall_lo.phi);
+  TEST_CHECK(species.phi_wall_lo.projector == 0);
+  TEST_ASSERT(species.phi_wall_lo.advance_func != 0);
+  TEST_CHECK(species.phi_wall_up.phi == 0);
+  TEST_CHECK(species.phi_wall_up.phi_host == 0);
+  TEST_CHECK(species.phi_wall_up.projector == 0);
+  TEST_ASSERT(species.phi_wall_up.advance_func != 0);
+
+  gk_species_phi_wall_advance(&app, &species.phi_wall_lo, 1.0);
+  gk_species_phi_wall_advance(&app, &species.phi_wall_up, 1.0);
+
+  struct gkyl_range_iter iter;
+  gkyl_range_iter_init(&iter, &app.local_ext);
+  while (gkyl_range_iter_next(&iter)) {
+    long lidx = gkyl_range_idx(&app.local_ext, iter.idx);
+    const double *phi_c = gkyl_array_cfetch(species.phi_wall_lo.phi, lidx);
+    for (int k=0; k<app.basis.num_basis; ++k)
+      TEST_CHECK(phi_c[k] == 0.0);
+  }
+
+  gk_species_phi_wall_release(&app, &species.phi_wall_lo);
+  gk_species_phi_wall_release(&app, &species.phi_wall_up);
+}
+
+TEST_LIST = {
+  { "test_species_wall_profiles", test_species_wall_profiles },
+  { "test_grounded_and_non_sheath_walls", test_grounded_and_non_sheath_walls },
+  { 0 },
+};