diff --git a/docs/source/tutorials/GenVeg/GenVeg_Dune_Simulation.xlsx b/docs/source/tutorials/GenVeg/GenVeg_Dune_Simulation.xlsx
index d3843f6a84..2d1130b996 100644
Binary files a/docs/source/tutorials/GenVeg/GenVeg_Dune_Simulation.xlsx and b/docs/source/tutorials/GenVeg/GenVeg_Dune_Simulation.xlsx differ
diff --git a/docs/source/user_guide/reference/components.md b/docs/source/user_guide/reference/components.md
index 17fedd3a06..c0dacfa04e 100644
--- a/docs/source/user_guide/reference/components.md
+++ b/docs/source/user_guide/reference/components.md
@@ -60,6 +60,7 @@ following categories of components:
 
 * {mod}`landlab.components.vegetation_dynamics`
 * {mod}`landlab.components.plant_competition_ca`
+* {mod}`landlab.components.genveg`
 
 ## Biota
 
diff --git a/src/landlab/components/genveg/allometry.py b/src/landlab/components/genveg/allometry.py
index fdda6a46c5..bbc3518d72 100644
--- a/src/landlab/components/genveg/allometry.py
+++ b/src/landlab/components/genveg/allometry.py
@@ -79,6 +79,7 @@ def calc_abg_dims(self, abg_biomass, cm=False):
         return (basal_dia, shoot_width, height)
 
     def calc_abg_biomass_from_dim(self, dim, var_name, cm=False):
+        # Use any dimension to estimate abg biomass
         abg_biomass = np.zeros_like(dim)
         filter = np.nonzero(dim > self.morph_params[var_name]["min"])
         if cm is True:
@@ -102,6 +103,7 @@ def __init__(self, params, empirical_coeffs, cm=True):
         super().__init__(params, empirical_coeffs, cm)
 
     def calc_abg_dims(self, abg_biomass, cm=True):
+        # Calculate all abg dimensions from abg biomass
         basal_dia = np.zeros_like(abg_biomass)
         shoot_sys_width = np.zeros_like(abg_biomass)
         height = np.zeros_like(abg_biomass)
@@ -116,6 +118,7 @@ def calc_abg_dims(self, abg_biomass, cm=True):
         return (basal_dia, shoot_sys_width, height)
 
     def _calc_shoot_width_from_basal_dia(self, basal_dia, cm=False):
+        # Apply log-log linear relationship to estimate canopy width from basal d
         canopy_area = np.zeros_like(basal_dia)
         filter = np.nonzero(basal_dia > self.morph_params["basal_dia"]["min"])
         if cm is True:
@@ -127,6 +130,7 @@ def _calc_shoot_width_from_basal_dia(self, basal_dia, cm=False):
         return shoot_width
 
     def _calc_height_from_basal_dia(self, basal_dia, cm=False):
+        # Apply log-log linear relationship to estimate height from basal d
         height = np.zeros_like(basal_dia)
         filter = np.nonzero(basal_dia > self.morph_params["basal_dia"]["min"])
         if cm is True:
@@ -137,6 +141,7 @@ def _calc_height_from_basal_dia(self, basal_dia, cm=False):
         return height
 
     def _calc_basal_dia_from_shoot_width(self, shoot_width, cm=False):
+        # Apply log-log linear relationship to estimate basal d from canopy width
         basal_dia = np.zeros_like(shoot_width)
         filter = np.nonzero(shoot_width > self.morph_params["shoot_sys_width"]["min"])
         canopy_area = 0.25 * np.pi * shoot_width**2
@@ -148,6 +153,8 @@ def _calc_basal_dia_from_shoot_width(self, shoot_width, cm=False):
         return basal_dia
 
 
+"""
+Not implemented in this version but code preserved here for future use
 class Multi_Dimensional(Biomass):
     def __init__(self, params):
         # Not implemented yet but saved here for future versions allowing
@@ -180,3 +187,4 @@ def _apply3_allometry_eq_for_xs(self, ys, zs, predict_var):
         c = self.morph_params["empirical_coeffs"][predict_var]["c"]
         ln_xs = (ln_zs - a - c * ln_ys) / b
         return np.exp(ln_xs)
+"""
diff --git a/src/landlab/components/genveg/dispersal.py b/src/landlab/components/genveg/dispersal.py
index 39fed5698c..7ed1e090c7 100644
--- a/src/landlab/components/genveg/dispersal.py
+++ b/src/landlab/components/genveg/dispersal.py
@@ -7,22 +7,24 @@
 class Repro:
     def __init__(self, params):
         grow_params = params["grow_params"]
+        disperse_params = params["disperse_params"]
         self.min_size = grow_params["growth_biomass"]["min"]
+        self.min_size_for_repro = disperse_params["min_size_dispersal"]
 
 
 class Clonal(Repro):
     def __init__(self, params):
-        disperse_params = params["dispersal_params"]
+        disperse_params = params["disperse_params"]
         super().__init__(params)
         self.unit_cost = disperse_params["unit_cost_dispersal"]
-        self.max_dist_dispersal = disperse_params["max_dist_dispersal"]
+        self.max_dist_dispersal = disperse_params["max_dist_runner"]
 
     def disperse(self, plants):
         """
         This method determines the potential location and cost
-        of clonal dispersal. A separate method in integrator.py
-        determines if dispersal is successful based if the potential location
-        is occupied.
+        of clonal dispersal. A separate method in
+        integrator.py determines if dispersal is successful if the
+        potential location is unoccupied.
         """
         max_runner_length = np.zeros_like(plants["root"])
         available_carb = plants["reproductive"] - 2 * self.min_size
@@ -30,38 +32,100 @@ def disperse(self, plants):
             available_carb[available_carb > 0] / self.unit_cost
         )
         runner_length = rng.uniform(
-            low=0.05,
+            low=0.02,
             high=self.max_dist_dispersal,
             size=plants.size,
         )
         pup_azimuth = np.deg2rad(rng.uniform(low=0.01, high=360, size=plants.size))
-
         filter = np.nonzero(runner_length <= max_runner_length)
-
-        plants["pup_x_loc"][filter] = (
+        # Save to first position of pup subarray
+        plants["dispersal"]["pup_x_loc"][filter] = (
             runner_length[filter] * np.cos(pup_azimuth[filter])
             + plants["x_loc"][filter]
         )
-        plants["pup_y_loc"][filter] = (
+        plants["dispersal"]["pup_y_loc"][filter] = (
             runner_length[filter] * np.sin(pup_azimuth)[filter]
             + plants["y_loc"][filter]
         )
-        plants["pup_cost"][filter] = runner_length[filter] * self.unit_cost
-
+        plants["dispersal"]["pup_cost"][filter] = runner_length[filter] * self.unit_cost
         return plants
 
 
 class Seed(Repro):
-    def __init__(self, params):
-        pass
+    def __init__(self, params, dt):
+        super().__init__(params)
+        params = params["disperse_params"]
+        self.biomass_to_seedlings = params["mass_to_seedling_rate"]
+        self.mean_dispersal_distance = params["mean_seed_distance"]
+        self.max_dispersal_distance = params["max_seed_distance"]
+        self.dispersal_shape = params["seed_distribution_shape_param"]
+        self.seed_size = params["seed_mass"]
+        self.seed_efficiency = params["seed_efficiency"]
+        self.dt = dt
 
-    def disperse(self):
-        print("New plant emerges from seed some distance within parent plant")
+    def disperse(self, plants, total_biomass):
+        """
+        This method implements dispersal using an algorithm similar to Vincenot 2016.
+        The log-normal distribution shape can be altered to accomodate many patterns of seed
+        distribution. Seed production is dependent on plant size and is limited by 
+        biomass in the reproductive organ system.
+        """
+        # Tracking max seeds that can be produced based on reproductive mass available
+        max_seeds = np.floor(plants["reproductive"] / (self.seed_size * self.seed_efficiency))
+        # Partial seedlings are banked to accomodate species with low seedling production rates
+        reserve_part, num_seedlings = np.modf(
+            (total_biomass * self.biomass_to_seedlings * self.dt.astype(int)),
+            out=(np.zeros_like(plants["reproductive"]), np.zeros_like(plants["reproductive"])),
+            where=(total_biomass >= self.min_size_for_repro)
+        )
+        num_seedlings[num_seedlings > max_seeds] = max_seeds
+        plants["reproductive"] -= num_seedlings * self.seed_size / self.seed_efficiency
+        plants["dispersal"]["seedling_reserve"] += reserve_part
+        if plants["dispersal"]["seedling_reserve"] >= 1:
+            num_seedlings += np.floor(plants["dispersal"]["seedling_reserve"])
+            plants["dispersal"]["seedling_reserve"] -= np.floor(plants["dispersal"]["seedling_reserve"])
+        # Loop through each plant to assign seedling locations
+        for i in range(plants.size):
+            count = num_seedlings[i].astype(int)
+            if count == 0:
+                break
+            dist_from_plant = rng.lognormal(self.mean_dispersal_distance, self.dispersal_shape, size=count)
+            dist_from_plant[dist_from_plant > self.max_dispersal_distance] = 0.0
+            pup_azimuth = np.deg2rad(rng.uniform(low=0.01, high=360, size=count))
+            plants["dispersal"]["seedling_x_loc"][i][:count] = (
+                dist_from_plant * np.cos(pup_azimuth) + plants["x_loc"][i]
+            )
+            plants["dispersal"]["seedling_y_loc"][i][:count] = (
+                dist_from_plant * np.sin(pup_azimuth) + plants["y_loc"][i]
+            )
+        return plants
 
 
 class Random(Repro):
-    def __init__(self, params):
-        pass
+    def __init__(self, params, cell_area, extent, reference):
+        super().__init__(params)
+        self.daily_col_prob = params["disperse_params"]["daily_colonization_probability"]
+        self.cell_area = cell_area
+        # Note reference is longitude-latitude and extent is y,x
+        self.grid_boundary = (
+            reference[1],
+            reference[1] + extent[1],
+            reference[0],
+            reference[0] + extent[0]
+        )
 
-    def disperse(self):
-        print("New plant randomly appears")
+    def disperse(self, plants):
+        filter = (rng.uniform(0, 1, size=plants["root"].size) > (self.daily_col_prob * self.cell_area))
+        plants["dispersal"]["rand_x_loc"][:, 0] = rng.uniform(
+            low=self.grid_boundary[0],
+            high=self.grid_boundary[1],
+            size=plants["root"].size
+        )
+        plants["dispersal"]["rand_y_loc"][:, 0] = rng.uniform(
+            low=self.grid_boundary[2],
+            high=self.grid_boundary[3],
+            size=plants["root"].size
+        )
+        plants["dispersal"]["rand_x_loc"][filter] = np.nan
+        plants["dispersal"]["rand_y_loc"][filter] = np.nan
+        return plants
diff --git a/src/landlab/components/genveg/duration.py b/src/landlab/components/genveg/duration.py
index 95ccb9265b..d6ba32bcfb 100644
--- a/src/landlab/components/genveg/duration.py
+++ b/src/landlab/components/genveg/duration.py
@@ -208,9 +208,10 @@ def __init__(self, species_grow_params, duration_params, dt, green_parts):
         )
 
     def emerge(self, plants, available_mass, persistent_total_mass):
-        print("I emerge from dormancy and I am a deciduous perennial")
-        # next steps are to clean this up using same approach as dormancy
-
+        """
+        This method allows deciduous perennial plants to emerge from dormancy
+        using non-structural biomass stored in the persistent parts of the plant.
+        """
         total_mass_new_green = np.zeros_like(plants["root_biomass"])
         new_green_biomass = {}
         for part in self.green_parts:
diff --git a/src/landlab/components/genveg/form.py b/src/landlab/components/genveg/form.py
deleted file mode 100644
index 4a1d1d4f9c..0000000000
--- a/src/landlab/components/genveg/form.py
+++ /dev/null
@@ -1,134 +0,0 @@
-import numpy as np
-
-from .dispersal import Clonal
-from .dispersal import Random
-from .dispersal import Seed
-
-rng = np.random.default_rng()
-
-
-class Abg_repro:
-    def __init__(self, params):
-        self.abg_parts = ("leaf", "stem", "reproductive")
-        self.dead_abg_parts = ("dead_leaf", "dead_stem", "dead_reproductive")
-        sums = {}
-        sum_vars = [
-
-            ["max_abg_biomass", "plant_part_max", self.abg_parts],
-            ["min_abg_biomass", "plant_part_min", self.abg_parts],
-        ]
-        for sum_var in sum_vars:
-            sums[sum_var[0]] = 0
-            for part in sum_var[2]:
-                sums[sum_var[0]] += params[
-                    "grow_params"
-                ][sum_var[1]][part]
-        self.max_abg_biomass = sums["max_abg_biomass"]
-        self.min_abg_biomass = sums["min_abg_biomass"]
-
-
-class Blg_repro:
-    def __init__(self, params):
-        self.abg_parts = ("leaf", "stem")
-        self.dead_abg_parts = ("dead_leaf", "dead_stem")
-        sums = {}
-        sum_vars = [
-
-            ["max_abg_biomass", "plant_part_max", self.abg_parts],
-            ["min_abg_biomass", "plant_part_min", self.abg_parts],
-        ]
-        for sum_var in sum_vars:
-            sums[sum_var[0]] = 0
-            for part in sum_var[2]:
-                sums[sum_var[0]] += params[
-                    "grow_params"
-                ][sum_var[1]][part]
-        self.max_abg_biomass = sums["max_abg_biomass"]
-        self.min_abg_biomass = sums["min_abg_biomass"]
-
-
-# Growth form classes and selection method
-class Bunch(Seed):
-    def __init__(self, params):
-        super().__init__(params)
-
-    def branch(self):
-        print("Limited lateral branching due to clumping")
-
-    def disperse(self, plants):
-        return plants
-
-
-class Colonizing(Random):
-    def __init__(self, params):
-        super().__init__(params)
-
-    def branch(self):
-        print("No branching annual")
-
-    def disperse(self, plants):
-        return plants
-
-
-class Multiplestems(Seed):
-    def __init__(self, params):
-        super().__init__(params)
-
-    def branch(self):
-        print("Create two or more main stems at or near soil surface")
-
-    def disperse(self, plants):
-        return plants
-
-
-class Rhizomatous(Clonal):
-    def __init__(self, params):
-        super().__init__(params)
-
-    def set_initial_branches(self, max_branches, arr_size):
-        n_branches = np.ceil(rng.rayleigh(scale=0.26, size=arr_size) * max_branches)
-        return n_branches
-
-    def branch(self):
-        print("Tiller via rhizomes")
-
-
-class Singlecrown(Seed):
-    def __init__(self, params):
-        super().__init__(params)
-
-    def branch(self):
-        print("Herbaceous plant with one persistent base")
-
-    def disperse(self, plants):
-        return plants
-
-
-class Singlestem(Seed):
-    def __init__(self, params):
-        super().__init__(params)
-
-    def branch(self):
-        print("Plant develops one stem like a tree or a corn plant")
-
-    def disperse(self, plants):
-        return plants
-
-
-class Stoloniferous(Clonal):
-    def __init__(self, params):
-        super().__init__(params)
-
-    def branch(self):
-        print("Branches via stolons")
-
-
-class Thicketforming(Seed):
-    def __init__(self, params):
-        super().__init__(params)
-
-    def branch(self):
-        print("Limited lateral branching due to dense thickets")
-
-    def disperse(self, plants):
-        return plants
diff --git a/src/landlab/components/genveg/growth.py b/src/landlab/components/genveg/growth.py
index 85d40b1247..d319513240 100644
--- a/src/landlab/components/genveg/growth.py
+++ b/src/landlab/components/genveg/growth.py
@@ -168,10 +168,8 @@ def __init__(
         """
         # Initialize species object to get correct species parameter list
         self._grid = grid
-        (_, _latitude) = self._grid.xy_of_reference
-        self._lat_rad = np.radians(_latitude)
         self.dt = dt
-        super().__init__(species_params, self._lat_rad, self.dt)
+        super().__init__(species_params, self._grid, self.dt)
         self.species_name = self.species_plant_factors["species"]
         self.time_ind = 1
         event_flags = self.set_event_flags(_current_jday)
@@ -208,9 +206,16 @@ def __init__(
             np.nan,
             np.nan,
             999999,
-            np.nan,
-            np.nan,
-            np.nan,
+            (
+                np.nan,
+                np.nan,
+                np.nan,
+                [np.nan],
+                [np.nan],
+                np.nan,
+                [np.nan],
+                [np.nan],
+            ),
             999999,
         )
         self.dtypes = [
@@ -224,8 +229,8 @@ def __init__(
             (("stem", "stem_biomass"), float),
             (("reproductive", "repro_biomass"), float),
             ("dead_root", float),
-            ("dead_leaf", float),
             ("dead_stem", float),
+            ("dead_leaf", float),
             ("dead_reproductive", float),
             ("dead_root_age", float),
             ("dead_leaf_age", float),
@@ -240,30 +245,41 @@ def __init__(
             ("live_leaf_area", float),
             ("plant_age", float),
             ("n_stems", int),
-            ("pup_x_loc", float),
-            ("pup_y_loc", float),
-            ("pup_cost", float),
+            ("dispersal", [
+                ("pup_x_loc", float),
+                ("pup_y_loc", float),
+                ("pup_cost", float),
+                ("seedling_x_loc", float, (10,)),
+                ("seedling_y_loc", float, (10,)),
+                ("seedling_reserve", float),
+                ("rand_x_loc", float, (10,)),
+                ("rand_y_loc", float, (10,)),
+            ]),
             ("item_id", int),
         ]
         mask_scalar = 1
         empty_list = []
         mask = []
-        for i in range(max_plants[0]):
+        for _ in range(max_plants[0]):
             empty_list.append(self.no_data_scalar)
             mask.append(mask_scalar)
         self.plants = np.ma.array(empty_list, mask=mask, dtype=self.dtypes)
-        self.plants.fill_value = self.no_data_scalar
+        self.plants[:] = self.no_data_scalar
+        species_cover = kwargs.get(
+            "species_cover",
+            np.zeros_like(self._grid["cell"]["vegetation__total_biomass"])
+        )
         try:
-            (init_plants, self.n_plants) = kwargs.get(
-                ("plant_array", "n_plants"),
+            init_plants = kwargs.get(
+                "plant_array",
                 self._init_plants_from_grid(
-                    _in_growing_season, kwargs["species_cover"]
+                    _in_growing_season, species_cover
                 ),
             )
         except KeyError:
             msg = "GenVeg requires a pre-populated plant array or a species cover."
             raise ValueError(msg)
-
+        self.n_plants = init_plants.size
         self.plants[: self.n_plants] = init_plants
 
         self.call = []
@@ -395,11 +411,16 @@ def species_plants(self):
     def species_get_variable(self, var_name):
         return self.species_grow_params
 
-    def update_plants(self, var_names, pids, var_vals):
-        for idx, var_name in enumerate(var_names):
-            self.plants[var_name][np.isin(self.plants["pid"], pids)] = var_vals[idx]
-        return self.plants
-
+    def update_plants(self, var_names, pids, var_vals, subarray=None):
+        if subarray is None:
+            for idx, var_name in enumerate(var_names):
+                self.plants[var_name][np.isin(self.plants["pid"], pids)] = var_vals[idx]
+            return self.plants
+        else:
+            for idx, var_name in enumerate(var_names):
+                self.plants[subarray][var_name][np.isin(self.plants["pid"], pids)] = var_vals[idx]
+            return self.plants
+        
     def add_new_plants(self, new_plants_list, _rel_time):
 
         # Reassess this. We need the INDEX of the last nanmax PID
@@ -619,7 +640,7 @@ def _init_plants_from_grid(self, in_growing_season, species_cover):
                     plantlist = self.set_initial_cover(cover_area, plant, pidval, cell_index, plantlist)
         plant_array = np.array(plantlist, dtype=self.dtypes)
         plant_array = self.set_initial_biomass(plant_array, in_growing_season)
-        return (plant_array, pidval)
+        return plant_array
 
     def set_event_flags(self, _current_jday):
         """
diff --git a/src/landlab/components/genveg/habit.py b/src/landlab/components/genveg/habit.py
index 6d49b7cb27..34267650c2 100644
--- a/src/landlab/components/genveg/habit.py
+++ b/src/landlab/components/genveg/habit.py
@@ -27,6 +27,7 @@ def __init__(self, params, allometry, dt=1, green_parts=(None)):
         )
 
     def _calc_canopy_area_from_shoot_width(self, shoot_sys_width):
+        # instead of returning error, set to zero
         UnitTestChecks().is_negative_present(shoot_sys_width, "shoot_sys_width")
         canopy_area = 0.25 * np.pi * shoot_sys_width**2
         return canopy_area
@@ -76,14 +77,18 @@ def _select_duration_class(
         duration_params,
         dt,
         duration_val,
-        green_parts=(None),
+        green_parts=("leaf"),
     ):
         duration = {
-            "annual": Annual(species_grow_params, duration_params, dt),
+            "annual": Annual(
+                species_grow_params, duration_params, dt
+            ),
             "perennial deciduous": Deciduous(
                 species_grow_params, duration_params, dt, green_parts
             ),
-            "perennial evergreen": Evergreen(species_grow_params, duration_params, dt),
+            "perennial evergreen": Evergreen(
+                species_grow_params, duration_params, dt
+            ),
         }
         return duration[duration_val]
 
@@ -339,6 +344,8 @@ def set_initial_biomass(self, plants, in_growing_season):
         return plants
 
 
+"""
+This will be added later
 class Tree(Habit):
     def __init__(self, params, dt, empirical_coeffs={"root_dia_coeffs": {"a": 0.35, "b": 0.31, }}):
         green_parts = ("leaf")
@@ -355,3 +362,4 @@ def __init__(self, params, dt, empirical_coeffs={"root_dia_coeffs": {"a": 0.35,
         green_parts = ("leaf")
         allometry = self._select_allometry_class(params, empirical_coeffs)
         super().__init__(params, allometry, dt, green_parts)
+"""
diff --git a/src/landlab/components/genveg/integrator.py b/src/landlab/components/genveg/integrator.py
index 8ffe4adde8..ac85c1413e 100644
--- a/src/landlab/components/genveg/integrator.py
+++ b/src/landlab/components/genveg/integrator.py
@@ -46,7 +46,8 @@ class GenVeg(Component, PlantGrowth):
     }
 
     def __init__(
-        self, grid, dt, current_day, vegparams, plant_array=np.empty((0, 30), dtype=[])
+        #self, grid, dt, current_day, vegparams, plant_array=np.empty((0, 28), dtype=[])
+        self, grid, dt, current_day, vegparams, **kwargs
     ):
         # save grid object to class
         super().__init__(grid)
@@ -70,7 +71,10 @@ def __init__(
         _ = self._grid.add_zeros("vegetation__n_plants", at="cell", clobber=True)
         _ = self._grid.add_zeros("vegetation__plant_height", at="cell", clobber=True)
         _ = self._grid.add_zeros("vegetation__lai", at="cell", clobber=True)
-
+        """
+        This is being instantiated incorrectly. We need to check the kwarg for a list containing
+        the plant arrays and then instantiate if it doesn't exist
+        """
         # Instantiate a PlantGrowth object and
         # summarize number of plants and biomass per cell
         if plant_array.size == 0:
@@ -487,17 +491,17 @@ def check_if_loc_unocc(self, plant_loc, plant_width, all_plants, check_type):
         return is_center_unocc
 
     def check_for_dispersal_success(self, all_plants):
-        new_pups = all_plants[~np.isnan(all_plants["pup_x_loc"])]
+        new_pups = all_plants[~np.isnan(all_plants["dispersal"]["pup_x_loc"])]
 
         if new_pups.size != 0:
-            pup_locs = tuple(zip(new_pups["pup_x_loc"], new_pups["pup_y_loc"]))
-            pup_widths = np.zeros_like(new_pups["pup_x_loc"])
+            pup_locs = tuple(zip(new_pups["dispersal"]["pup_x_loc"], new_pups["dispersal"]["pup_y_loc"]))
+            pup_widths = np.zeros_like(new_pups["dispersal"]["pup_x_loc"])
             loc_unoccupied = self.check_if_loc_unocc(
                 pup_locs, pup_widths, all_plants, "below"
             )
             new_pups = new_pups[np.nonzero(loc_unoccupied)]
-            new_pups["x_loc"] = new_pups["pup_x_loc"]
-            new_pups["y_loc"] = new_pups["pup_y_loc"]
+            new_pups["x_loc"] = new_pups["dispersal"]["pup_x_loc"]
+            new_pups["y_loc"] = new_pups["dispersal"]["pup_y_loc"]
 
         for species_obj in self.plant_species:
             species = species_obj.species_name
@@ -522,7 +526,7 @@ def check_for_dispersal_success(self, all_plants):
                     species_new_pups["root"]
                     + species_new_pups["leaf"]
                     + species_new_pups["stem"]
-                    + species_parents["pup_cost"]
+                    + species_parents["dispersal"]["pup_cost"]
                 )
 
                 species_obj.update_plants(
@@ -543,8 +547,8 @@ def check_for_dispersal_success(self, all_plants):
                         np.full_like(species_plants["root"], np.nan),
                     )
                 ),
+                subarray="dispersal"
             )
-
         return self.combine_plant_arrays()
 
     def combine_plant_arrays(self):
diff --git a/src/landlab/components/genveg/species.py b/src/landlab/components/genveg/species.py
index 1b73807fe8..7f24cd8840 100644
--- a/src/landlab/components/genveg/species.py
+++ b/src/landlab/components/genveg/species.py
@@ -12,19 +12,10 @@
 from sympy import symbols
 
 from .check_objects import UnitTestChecks
-from .form import Bunch
-from .form import Colonizing
-from .form import Multiplestems
-from .form import Rhizomatous
-from .form import Singlecrown
-from .form import Singlestem
-from .form import Stoloniferous
-from .form import Thicketforming
+from .dispersal import Clonal, Seed, Random
 from .habit import Forbherb
 from .habit import Graminoid
 from .habit import Shrub
-from .habit import Tree
-from .habit import Vine
 from .photosynthesis import C3
 from .photosynthesis import C4
 from .photosynthesis import Cam
@@ -34,21 +25,25 @@
 
 # Define species class that inherits composite class methods
 class Species:
-    def __init__(self, species_params, latitude, dt=1):
+    def __init__(self, species_params, grid, dt=1):
         self.dt = dt
+        extent = grid.extent
+        cell_area = grid.area_of_cell
+        reference = grid.xy_of_reference
+        _lat_rad = np.radians(reference[1])
         self.validate_plant_factors(species_params["plant_factors"])
         self.validate_duration_params(species_params["duration_params"])
         self.define_plant_parts(species_params)
         species_params = self.calculate_derived_params(species_params)
-        self.form = self.select_form_class(species_params)
         self.habit = self.select_habit_class(species_params)
-        self.photosynthesis = self.select_photosythesis_type(species_params, latitude)
+        self.photosynthesis = self.select_photosythesis_type(species_params, _lat_rad)
+        self.dispersal = self.select_dispersal_type(species_params, cell_area, extent, reference)
         # check these below to see if we need to save the composition dictionary not the original
         self.species_plant_factors = species_params["plant_factors"]
         self.species_duration_params = species_params["duration_params"]
         self.species_grow_params = species_params["grow_params"]
         self.species_photo_params = species_params["photo_params"]
-        self.species_dispersal_params = species_params["dispersal_params"]
+        self.species_dispersal_params = species_params["disperse_params"]
         self.species_mort_params = species_params["mortality_params"]
         self.species_morph_params = self.habit.morph_params
         self.populate_biomass_allocation_array()
@@ -56,20 +51,11 @@ def __init__(self, species_params, latitude, dt=1):
     def validate_plant_factors(self, plant_factors):
         plant_factor_options = {
             "species": [],
-            "growth_habit": ["forb_herb", "graminoid", "shrub", "tree", "vine"],
+            "growth_habit": ["forb_herb", "graminoid", "shrub"],
             "monocot_dicot": ["monocot", "dicot"],
             "angio_gymno": ["angiosperm", "gymnosperm"],
             "duration": ["annual", "perennial deciduous", "perennial evergreen"],
-            "growth_form": [
-                "bunch",
-                "colonizing",
-                "multiple_stems",
-                "rhizomatous",
-                "single_crown",
-                "single_stem",
-                "stoloniferous",
-                "thicket_forming",
-            ],
+            "reproductive_modes": ["clonal", "seed", "random"],
             "storage": ["aboveground", "belowground"],
             "p_type": ["C3", "C4", "Cam"],
         }
@@ -78,9 +64,14 @@ def validate_plant_factors(self, plant_factors):
             try:
                 opt_list = plant_factor_options[key]
                 if opt_list:
-                    if plant_factors[key] not in opt_list:
-                        msg = "Invalid " + str(key) + " option"
-                        raise ValueError(msg)
+                    msg = "Invalid " + str(key) + " option"
+                    if isinstance(plant_factors[key], str):
+                        if plant_factors[key] not in opt_list:
+                            raise ValueError(msg)
+                    elif isinstance(plant_factors[key], list):
+                        for item in plant_factors[key]:
+                            if item not in opt_list:
+                                raise ValueError(msg)
             except KeyError:
                 raise KeyError(
                     "Unexpected variable name in species parameter dictionary. Please check input parameter file"
@@ -275,31 +266,32 @@ def select_photosythesis_type(self, species_params, latitude):
             latitude, photo_params=species_params["photo_params"]
         )
 
+    def select_dispersal_type(self, species_params, cell_area, extent, reference):
+        repro_modes = species_params["plant_factors"]["reproductive_modes"]
+        repro_options = {
+            "clonal": Clonal(species_params),
+            "seed": Seed(species_params, dt=self.dt),
+            "random": Random(species_params, cell_area, extent, reference)
+        }
+        dispersal_classes = []
+
+        for mode in repro_modes:
+            if mode in repro_options:
+                dispersal_classes.append(repro_options[mode])
+            else:
+                print(f"Warning: Invalid reproductive mode '{mode}'. Skipping.")
+
+        return dispersal_classes
+
     def select_habit_class(self, species_params):
         habit_type = species_params["plant_factors"]["growth_habit"]
         habit = {
             "forb_herb": Forbherb,
             "graminoid": Graminoid,
             "shrub": Shrub,
-            "tree": Tree,
-            "vine": Vine,
         }
         return habit[habit_type](species_params, self.dt)
 
-    def select_form_class(self, species_params):
-        form_type = species_params["plant_factors"]["growth_form"]
-        form = {
-            "bunch": Bunch,
-            "colonizing": Colonizing,
-            "multiple_stems": Multiplestems,
-            "rhizomatous": Rhizomatous,
-            "single_crown": Singlecrown,
-            "single_stem": Singlestem,
-            "stoloniferous": Stoloniferous,
-            "thicket_forming": Thicketforming,
-        }
-        return form[form_type](species_params)
-
     def populate_biomass_allocation_array(self):
         # This method precalculates the biomass allocation array based on plant
         # type (angiosperm/gymnosperm, monocot/dicot) or based on user-defined
@@ -528,12 +520,10 @@ def _adjust_biomass_allocation_towards_ideal(self, _new_biomass):
         _new_biomass["stem_biomass"] = _new_stem_mass_frac * _total_biomass
         return _new_biomass
 
-    def branch(self):
-        self.form.branch()
-
-    def disperse(self, plants, jday):
+    def disperse(self, plants):
         # decide how to parameterize reproductive schedule, make repro event
-        plants = self.form.disperse(plants)
+        for disperse_method in self.dispersal:
+            plants = self.dispersal.disperse(plants)
         return plants
 
     def enter_dormancy(
diff --git a/tests/components/genveg/conftest.py b/tests/components/genveg/conftest.py
index deff3dd8e5..bb643fbbef 100644
--- a/tests/components/genveg/conftest.py
+++ b/tests/components/genveg/conftest.py
@@ -2,6 +2,13 @@
 import pytest
 
 from landlab import RasterModelGrid
+from landlab.components.genveg.integrator import GenVeg
+from landlab.components.genveg.growth import PlantGrowth
+from landlab.components.genveg.allometry import Biomass, Dimensional
+from landlab.components.genveg.duration import Annual, Deciduous, Evergreen
+from landlab.components.genveg.habit import Habit, Graminoid
+from landlab.components.genveg.species import Species
+from landlab.components.genveg.dispersal import Clonal, Seed, Random
 
 
 @pytest.fixture(autouse=True)
@@ -15,10 +22,17 @@ def example_input_params():
     param_dict = {
         "BTS": {
             "col_params": {"prob_colonization": 0.01, "time_to_colonization": 365},
-            "dispersal_params": {
-                "max_dist_dispersal": 0.4,
+            "disperse_params": {
+                "max_dist_runner": 0.4,
                 "min_size_dispersal": 0.5,
                 "unit_cost_dispersal": 1.2,
+                "mass_to_seedling_rate": 0.02,
+                "mean_seed_distance": 0.6,
+                "max_seed_distance": 10,
+                "seed_distribution_shape_param": 1.5,
+                "seed_mass": 0.002,
+                "seed_efficiency": (300 / 5000),
+                "daily_colonization_probability": 0.01,
             },
             "duration_params": {
                 "growing_season_end": 305,
@@ -205,10 +219,10 @@ def example_input_params():
             "plant_factors": {
                 "angio_gymno": "angiosperm",
                 "duration": "perennial deciduous",
-                "growth_form": "rhizomatous",
                 "growth_habit": "graminoid",
                 "monocot_dicot": "monocot",
                 "p_type": "C3",
+                "reproductive_modes": ["clonal", "seed", "random"],
                 "species": "BTS",
                 "storage": "belowground"
             },
@@ -247,9 +261,16 @@ def example_plant():
         ("live_leaf_area", float),
         ("plant_age", float),
         ("n_stems", int),
-        ("pup_x_loc", float),
-        ("pup_y_loc", float),
-        ("pup_cost", float),
+        ("dispersal", [
+            ("pup_x_loc", float),
+            ("pup_y_loc", float),
+            ("pup_cost", float),
+            ("seedling_x_loc", float, (10,)),
+            ("seedling_y_loc", float, (10,)),
+            ("seedling_reserve", float),
+            ("random_x_loc", float, (10,)),
+            ("random_y_loc", float, (10,)),
+        ]),
         ("item_id", int),
     ]
     plants = np.empty(1, dtype=dtypes)
@@ -287,9 +308,16 @@ def example_plant():
                 0.0,
                 0.0,
                 0,
-                np.nan,
-                np.nan,
-                np.nan,
+                (
+                    np.nan,
+                    np.nan,
+                    np.nan,
+                    [np.nan],
+                    [np.nan],
+                    np.nan,
+                    [np.nan],
+                    [np.nan],
+                ),
                 i,
             )
         )
@@ -297,6 +325,7 @@ def example_plant():
     plants["root"] = np.array([0.80000000000000004])
     plants["stem"] = np.array([0.29999999999999999])
     plants["leaf"] = np.array([0.50000000000000000])
+    plants["reproductive"] = np.array([0.30])
     plants["live_leaf_area"] = 0.022 * plants["leaf"]
     plants["shoot_sys_width"] = np.array(
         [(4 * (plants["live_leaf_area"] / 0.012) / np.pi) ** 0.5]
@@ -338,9 +367,16 @@ def example_plant_array(set_random_seed):
         ("live_leaf_area", float),
         ("plant_age", float),
         ("n_stems", int),
-        ("pup_x_loc", float),
-        ("pup_y_loc", float),
-        ("pup_cost", float),
+        ("dispersal", [
+            ("pup_x_loc", float),
+            ("pup_y_loc", float),
+            ("pup_cost", float),
+            ("seedling_x_loc", float, (10,)),
+            ("seedling_y_loc", float, (10,)),
+            ("seedling_reserve", float),
+            ("rand_x_loc", float, (10,)),
+            ("rand_y_loc", float, (10,)),
+        ]),
         ("item_id", int),
     ]
     plants = np.empty(100, dtype=dtypes)
@@ -378,9 +414,16 @@ def example_plant_array(set_random_seed):
                 0.0,
                 0.0,
                 0,
-                np.nan,
-                np.nan,
-                np.nan,
+                (
+                    np.nan,
+                    np.nan,
+                    np.nan,
+                    [np.nan],
+                    [np.nan],
+                    np.nan,
+                    [np.nan],
+                    [np.nan],
+                ),
                 i,
             )
         )
@@ -411,9 +454,9 @@ def example_plant_array(set_random_seed):
     plants["live_leaf_area"] = plants["total_leaf_area"]
     plants["plant_age"] = rng.uniform(low=1 / 365, high=5, size=plants.size)
     plants["n_stems"] = rng.integers(1, 6, size=plants.size)
-    plants["pup_x_loc"][plants["reproductive"] > 0.1] = 0.0
-    plants["pup_y_loc"][plants["reproductive"] > 0.1] = 0.0
-    plants["pup_cost"][plants["reproductive"] > 0.1] = 0.75
+    plants["dispersal"]["pup_x_loc"][plants["reproductive"] > 0.1] = 0.1
+    plants["dispersal"]["pup_y_loc"][plants["reproductive"] > 0.1] = 0.2
+    plants["dispersal"]["pup_cost"][plants["reproductive"] > 0.1] = 0.75
     plants["item_id"] = np.array([1, 2, 3, 4, 5, 6, 7, 8])
     return plants
 
@@ -452,6 +495,12 @@ def one_cell_grid():
         np.full(grid.number_of_cells, "Corn"),
         at="cell",
     )
+    _ = grid.add_field(
+        "vegetation__total_biomass",
+        np.zeros(grid.number_of_cells),
+        at="cell",
+        units="g",
+    )
     _ = grid.add_field(
         "soil_water__volume_fraction",
         water_content * np.ones(grid.number_of_cells),
@@ -495,6 +544,12 @@ def two_cell_grid():
         np.full(grid.number_of_cells, "Corn"),
         at="cell",
     )
+    _ = grid.add_field(
+        "vegetation__total_biomass",
+        np.zeros(grid.number_of_cells),
+        at="cell",
+        units="g",
+    )
     _ = grid.add_field(
         "soil_water__volume_fraction",
         water_content * np.ones(grid.number_of_cells),
@@ -502,3 +557,148 @@ def two_cell_grid():
         units="m**3/m**3",
     )
     return grid
+
+
+@pytest.fixture
+def biomass_minmax(example_input_params):
+    params = example_input_params["BTS"]
+    params["grow_params"]["abg_biomass"] = {}
+    params["morph_params"]["canopy_area"] = {}
+    params["grow_params"]["abg_biomass"]["min"] = params["grow_params"]["plant_part_min"]["leaf"] + params["grow_params"]["plant_part_min"]["stem"]
+    params["grow_params"]["abg_biomass"]["max"] = params["grow_params"]["plant_part_max"]["leaf"] + params["grow_params"]["plant_part_max"]["stem"]
+    params["morph_params"]["canopy_area"]["min"] = params["morph_params"]["shoot_sys_width"]["min"]**2 * 0.25 * np.pi
+    params["morph_params"]["canopy_area"]["max"] = params["morph_params"]["shoot_sys_width"]["max"]**2 * 0.25 * np.pi
+    return Biomass(params)
+
+
+@pytest.fixture
+def dimensional_default(example_input_params):
+    params = example_input_params["BTS"]
+    params["morph_params"]["allometry_method"] = "default"
+    params["grow_params"]["abg_biomass"] = {}
+    params["morph_params"]["canopy_area"] = {}
+    params["grow_params"]["abg_biomass"]["min"] = params["grow_params"]["plant_part_min"]["leaf"] + params["grow_params"]["plant_part_min"]["stem"]
+    params["grow_params"]["abg_biomass"]["max"] = params["grow_params"]["plant_part_max"]["leaf"] + params["grow_params"]["plant_part_max"]["stem"]
+    params["morph_params"]["canopy_area"]["min"] = params["morph_params"]["shoot_sys_width"]["min"]**2 * 0.25 * np.pi
+    params["morph_params"]["canopy_area"]["max"] = params["morph_params"]["shoot_sys_width"]["max"]**2 * 0.25 * np.pi
+    empirical_coeffs = {
+        "basal_dia_coeffs": {"a": 0.5093, "b": 0.47},
+        "height_coeffs": {"a": np.log(0.232995), "b": 0.619077},
+        "canopy_area_coeffs": {"a": np.log(0.23702483 * 0.2329925**0.9459644), "b": 0.72682 + (0.619077 * 0.9459644)},
+    }
+    yield Dimensional(params, empirical_coeffs)
+
+
+@pytest.fixture
+def dimensional_user(example_input_params):
+    params = example_input_params["BTS"]
+    params["morph_params"]["allometry_method"] = "user-defined"
+    params["grow_params"]["abg_biomass"] = {}
+    params["morph_params"]["canopy_area"] = {}
+    params["grow_params"]["abg_biomass"]["min"] = params["grow_params"]["plant_part_min"]["leaf"] + params["grow_params"]["plant_part_min"]["stem"]
+    params["grow_params"]["abg_biomass"]["max"] = params["grow_params"]["plant_part_max"]["leaf"] + params["grow_params"]["plant_part_max"]["stem"]
+    params["morph_params"]["canopy_area"]["min"] = params["morph_params"]["shoot_sys_width"]["min"]**2 * 0.25 * np.pi
+    params["morph_params"]["canopy_area"]["max"] = params["morph_params"]["shoot_sys_width"]["max"]**2 * 0.25 * np.pi
+    empirical_coeffs = {
+        "basal_dia_coeffs": {"a": 0.5093, "b": 0.47},
+        "height_coeffs": {"a": np.log(0.232995), "b": 0.619077},
+        "canopy_area_coeffs": {"a": np.log(0.23702483 * 0.2329925**0.9459644), "b": 0.72682 + (0.619077 * 0.9459644)},
+    }
+    yield Dimensional(params, empirical_coeffs)
+
+
+@pytest.fixture
+def habit_object(example_input_params):
+    params = example_input_params["BTS"]
+    params["grow_params"]["abg_biomass"] = {}
+    params["morph_params"]["canopy_area"] = {}
+    params["grow_params"]["abg_biomass"]["min"] = params["grow_params"]["plant_part_min"]["leaf"] + params["grow_params"]["plant_part_min"]["stem"]
+    params["grow_params"]["abg_biomass"]["max"] = params["grow_params"]["plant_part_max"]["leaf"] + params["grow_params"]["plant_part_max"]["stem"]
+    params["morph_params"]["canopy_area"]["min"] = params["morph_params"]["shoot_sys_width"]["min"]**2 * 0.25 * np.pi
+    params["morph_params"]["canopy_area"]["max"] = params["morph_params"]["shoot_sys_width"]["max"]**2 * 0.25 * np.pi
+    allometry = Biomass(params, empirical_coeffs={"root_dia_coeffs": {"a": 0.08, "b": 0.24}})
+    yield Habit(params, allometry, dt=1, green_parts=("leaf", "stem"))
+
+
+@pytest.fixture
+def graminoid_object(example_input_params):
+    params = example_input_params["BTS"]
+    params["grow_params"]["abg_biomass"] = {}
+    params["morph_params"]["canopy_area"] = {}
+    params["grow_params"]["abg_biomass"]["min"] = params["grow_params"]["plant_part_min"]["leaf"] + params["grow_params"]["plant_part_min"]["stem"]
+    params["grow_params"]["abg_biomass"]["max"] = params["grow_params"]["plant_part_max"]["leaf"] + params["grow_params"]["plant_part_max"]["stem"]
+    params["morph_params"]["canopy_area"]["min"] = params["morph_params"]["shoot_sys_width"]["min"]**2 * 0.25 * np.pi
+    params["morph_params"]["canopy_area"]["max"] = params["morph_params"]["shoot_sys_width"]["max"]**2 * 0.25 * np.pi
+    yield Graminoid(params, dt=1)
+
+
+@pytest.fixture
+def species_object(example_input_params, one_cell_grid):
+    dt = np.timedelta64(1, 'D')
+    yield Species(example_input_params["BTS"], one_cell_grid, dt=dt)
+
+
+@pytest.fixture
+def clonal(example_input_params):
+    example_input_params["BTS"]["grow_params"]["growth_biomass"] = {}
+    example_input_params["BTS"]["grow_params"]["growth_biomass"]["min"] = (
+        example_input_params["BTS"]["grow_params"]["plant_part_min"]["root"]
+        + example_input_params["BTS"]["grow_params"]["plant_part_min"]["leaf"]
+        + example_input_params["BTS"]["grow_params"]["plant_part_min"]["stem"]
+    )
+    yield Clonal(example_input_params["BTS"])
+
+
+@pytest.fixture
+def random(example_input_params, one_cell_grid):
+    example_input_params["BTS"]["grow_params"]["growth_biomass"] = {}
+    example_input_params["BTS"]["grow_params"]["growth_biomass"]["min"] = (
+        example_input_params["BTS"]["grow_params"]["plant_part_min"]["root"]
+        + example_input_params["BTS"]["grow_params"]["plant_part_min"]["leaf"]
+        + example_input_params["BTS"]["grow_params"]["plant_part_min"]["stem"]
+    )
+    cell_area = one_cell_grid.area_of_cell
+    extent = one_cell_grid.extent
+    reference = one_cell_grid.xy_of_reference
+    yield Random(example_input_params["BTS"], cell_area, extent, reference)
+
+
+@pytest.fixture
+def seed(example_input_params, one_cell_grid):
+    example_input_params["BTS"]["grow_params"]["growth_biomass"] = {}
+    example_input_params["BTS"]["grow_params"]["growth_biomass"]["min"] = (
+        example_input_params["BTS"]["grow_params"]["plant_part_min"]["root"]
+        + example_input_params["BTS"]["grow_params"]["plant_part_min"]["leaf"]
+        + example_input_params["BTS"]["grow_params"]["plant_part_min"]["stem"]
+    )
+    dt = np.timedelta64(1, 'D')
+    yield Seed(example_input_params["BTS"], dt)
+
+
+@pytest.fixture
+def growth_obj(one_cell_grid, example_input_params, example_plant_array):
+    dt = np.timedelta64(1, 'D')
+    jday = 195
+    rel_time = 194
+    array = example_plant_array.copy()
+    yield PlantGrowth(
+        one_cell_grid,
+        dt,
+        rel_time,
+        jday,
+        species_params=example_input_params["BTS"],
+        plant_array=array,
+    )
+
+
+@pytest.fixture
+def genveg_obj(one_cell_grid, example_input_params, example_plant_array):
+    dt = np.timedelta64(1, 'D')
+    current_day = np.datetime64("2019-07-14", "D")
+    yield GenVeg(
+        one_cell_grid,
+        dt,
+        current_day,
+        example_input_params,
+        plant_array=example_plant_array
+    )
diff --git a/tests/components/genveg/test_allometry.py b/tests/components/genveg/test_allometry.py
index f9a5eb6c67..450a60ac87 100644
--- a/tests/components/genveg/test_allometry.py
+++ b/tests/components/genveg/test_allometry.py
@@ -3,94 +3,39 @@
 from numpy.testing import assert_allclose
 from numpy.testing import assert_array_almost_equal
 
-from landlab.components.genveg.allometry import Biomass, Dimensional
-from landlab.components.genveg.species import Species
 
 dt = np.timedelta64(1, 'D')
 
 
-def create_species_object(example_input_params):
-    return Species(example_input_params["BTS"], dt=dt, latitude=0.9074)
-
-
-def create_biomass_object(example_input_params):
-    params = example_input_params["BTS"]
-    params["grow_params"]["abg_biomass"] = {}
-    params["morph_params"]["canopy_area"] = {}
-    params["grow_params"]["abg_biomass"]["min"] = params["grow_params"]["plant_part_min"]["leaf"] + params["grow_params"]["plant_part_min"]["stem"]
-    params["grow_params"]["abg_biomass"]["max"] = params["grow_params"]["plant_part_max"]["leaf"] + params["grow_params"]["plant_part_max"]["stem"]
-    params["morph_params"]["canopy_area"]["min"] = params["morph_params"]["shoot_sys_width"]["min"]**2 * 0.25 * np.pi
-    params["morph_params"]["canopy_area"]["max"] = params["morph_params"]["shoot_sys_width"]["max"]**2 * 0.25 * np.pi
-    return Biomass(params)
-
-
-def create_default_dimensional_object(example_input_params):
-    params = example_input_params["BTS"]
-    params["morph_params"]["allometry_method"] = "default"
-    params["grow_params"]["abg_biomass"] = {}
-    params["morph_params"]["canopy_area"] = {}
-    params["grow_params"]["abg_biomass"]["min"] = params["grow_params"]["plant_part_min"]["leaf"] + params["grow_params"]["plant_part_min"]["stem"]
-    params["grow_params"]["abg_biomass"]["max"] = params["grow_params"]["plant_part_max"]["leaf"] + params["grow_params"]["plant_part_max"]["stem"]
-    params["morph_params"]["canopy_area"]["min"] = params["morph_params"]["shoot_sys_width"]["min"]**2 * 0.25 * np.pi
-    params["morph_params"]["canopy_area"]["max"] = params["morph_params"]["shoot_sys_width"]["max"]**2 * 0.25 * np.pi
-    empirical_coeffs = {
-        "basal_dia_coeffs": {"a": 0.5093, "b": 0.47},
-        "height_coeffs": {"a": np.log(0.232995), "b": 0.619077},
-        "canopy_area_coeffs": {"a": np.log(0.23702483 * 0.2329925**0.9459644), "b": 0.72682 + (0.619077 * 0.9459644)},
-    }
-    return Dimensional(params, empirical_coeffs)
-
-
-def test__calc2_allometry_coeffs(example_input_params):
-    b = create_biomass_object(example_input_params)
+def test__calc2_allometry_coeffs(biomass_minmax):
     x_min = np.array([1.])
     x_max = np.array([3.])
     y_min = np.array([2.])
     y_max = np.array([6.])
     slope = np.array([1])
     y_intercept = np.array([np.log(2)])
-    (b, m) = b._calc2_allometry_coeffs(x_min, x_max, y_min, y_max)
+    (b, m) = biomass_minmax._calc2_allometry_coeffs(x_min, x_max, y_min, y_max)
     assert_array_almost_equal(y_intercept, b, decimal=5)
     assert_array_almost_equal(slope, m, decimal=5)
 
 
-def test_apply2_allometry_eq_for_xs(example_input_params):
-    b = create_biomass_object(example_input_params)
-    b.morph_params["empirical_coeffs"]["basal_dia_coeffs"] = {"a": np.log(2), "b": 1}
+def test_apply2_allometry_eq_for_xs(biomass_minmax):
+    biomass_minmax.morph_params["empirical_coeffs"]["basal_dia_coeffs"] = {"a": np.log(2), "b": 1}
     abg_biomass = np.array([2.1, 0.2, 15.2])
-    pred_basal_diameters = b._apply2_allometry_eq_for_xs(abg_biomass, "basal_dia_coeffs")
+    pred_basal_diameters = biomass_minmax._apply2_allometry_eq_for_xs(abg_biomass, "basal_dia_coeffs")
     basal_diameters = np.array([1.05, 0.1, 7.6])
     assert_array_almost_equal(pred_basal_diameters, basal_diameters, decimal=5)
 
 
-def test_apply2_allometry_eq_for_ys(example_input_params):
-    b = create_biomass_object(example_input_params)
-    b.morph_params["empirical_coeffs"]["basal_dia_coeffs"] = {"a": np.log(2), "b": 1}
+def test_apply2_allometry_eq_for_ys(biomass_minmax):
+    biomass_minmax.morph_params["empirical_coeffs"]["basal_dia_coeffs"] = {"a": np.log(2), "b": 1}
     basal_diameter = np.array([0.5, 3.5, 12])
-    pred_abg_biomass = b._apply2_allometry_eq_for_ys(basal_diameter, "basal_dia_coeffs")
+    pred_abg_biomass = biomass_minmax._apply2_allometry_eq_for_ys(basal_diameter, "basal_dia_coeffs")
     abg_biomass = np.array([1., 7., 24.])
     assert_array_almost_equal(pred_abg_biomass, abg_biomass, decimal=5)
 
 
-def test_allometric_coeffs_are_user_defined(example_input_params):
-    example_input_params["BTS"]["morph_params"]["allometry_method"] = "user-defined"
-    b = create_biomass_object(example_input_params)
-    assert b.morph_params["empirical_coeffs"] == example_input_params["BTS"]["morph_params"]["empirical_coeffs"]
-
-
-def test_allometric_coeffs_are_default(example_input_params):
-    b = create_default_dimensional_object(example_input_params)
-    empirical_coeffs = {
-        "basal_dia_coeffs": {"a": 0.5093, "b": 0.47},
-        "height_coeffs": {"a": np.log(0.232995), "b": 0.619077},
-        "canopy_area_coeffs": {"a": np.log(0.23702483 * 0.2329925**0.9459644), "b": 0.72682 + (0.619077 * 0.9459644)},
-    }
-    for item in empirical_coeffs:
-        assert b.morph_params["empirical_coeffs"][item] == pytest.approx(empirical_coeffs[item])
-
-
-def test_allometric_coeffs_calculated(example_input_params):
-    b = create_biomass_object(example_input_params)
+def test_allometric_coeffs_calculated(biomass_minmax):
     # Check calculation for min-max
     empirical_coeffs = {
         "basal_dia_coeffs": {"a": 6.26936, "b": 1.74048},
@@ -98,69 +43,81 @@ def test_allometric_coeffs_calculated(example_input_params):
         "canopy_area_coeffs": {"a": 4.2926, "b": 0.766496},
     }
     for item in empirical_coeffs:
-        assert b.morph_params["empirical_coeffs"][item] == pytest.approx(empirical_coeffs[item])
+        assert biomass_minmax.morph_params["empirical_coeffs"][item] == pytest.approx(empirical_coeffs[item])
 
 
-def test_calc_abg_dims(example_input_params):
+def test_calc_abg_dims_minmax(biomass_minmax):
     # Biomass relationships
-    b = create_biomass_object(example_input_params)
     abg_biomass = np.array([2.1, 0.5, 15.2])
     basal_diameter = np.array([0.041761, 0.018309, 0.130217])
     height = np.array([0.5572, 0.3432, 1.0873])
     shoot_sys_width = np.array([0.11132, 0.04365, 0.40486])
-    pred_bd, pred_ssw, pred_h = b.calc_abg_dims(abg_biomass, cm=False)
-    assert_allclose(basal_diameter, pred_bd, rtol=0.001)
-    assert_allclose(height, pred_h, rtol=0.001)
-    assert_allclose(shoot_sys_width, pred_ssw, rtol=0.001)
-    # Dimensional relationships
+    pred_bd, pred_ssw, pred_h = biomass_minmax.calc_abg_dims(abg_biomass, cm=False)
+    print(biomass_minmax.morph_params["allometry_method"])
+    assert_allclose(pred_bd, basal_diameter, rtol=0.001)
+    assert_allclose(pred_h, height, rtol=0.001)
+    assert_allclose(pred_ssw, shoot_sys_width, rtol=0.001)
+
 
-    d = create_default_dimensional_object(example_input_params)
+def test_calc_abg_biomass_from_dim(biomass_minmax):
+    basal_dia = np.array([0.005, 0.035, 0.12])
+    shoot_sys_width = np.array([0.0640, 0.2, 0.48])
+    abg_from_bd = np.array([0, 1.5445, 13.1851])
+    abg_from_ssw = np.array([0.899, 5.156, 19.7325])
+    pred_abg_bd = biomass_minmax.calc_abg_biomass_from_dim(basal_dia, "basal_dia", cm=False)
+    assert_allclose(abg_from_bd, pred_abg_bd, rtol=0.001)
+    pred_abg_ssw = biomass_minmax.calc_abg_biomass_from_dim((0.25 * np.pi * shoot_sys_width**2), "canopy_area")
+    assert_allclose(abg_from_ssw, pred_abg_ssw, rtol=0.001)
+
+
+def test_allometric_coeffs_are_user_defined(dimensional_user, example_input_params):
+    empirical_coeffs = example_input_params["BTS"]["morph_params"]["empirical_coeffs"]
+    assert dimensional_user.morph_params["empirical_coeffs"] == empirical_coeffs
+
+
+def test_allometric_coeffs_are_default(dimensional_default):
+    empirical_coeffs = {
+        "basal_dia_coeffs": {"a": 0.5093, "b": 0.47},
+        "height_coeffs": {"a": np.log(0.232995), "b": 0.619077},
+        "canopy_area_coeffs": {"a": np.log(0.23702483 * 0.2329925**0.9459644), "b": 0.72682 + (0.619077 * 0.9459644)},
+    }
+    for item in empirical_coeffs:
+        assert dimensional_default.morph_params["empirical_coeffs"][item] == pytest.approx(empirical_coeffs[item])
+
+
+def test_calc_abg_dims_default(dimensional_default):
+    # Dimensional relationships
     basal_diameter = np.array([0.009, 0.036, 0.05])
     height = np.array([0.218283, 0.514921, 0.631049])
     shoot_sys_width = np.array([0.257388353, 0.639255015, 0.79304121])
     abg_biomass = np.array([1.58372641, 3.03842384, 3.54570077])
-    pred_bd, pred_ssw, pred_h = d.calc_abg_dims(abg_biomass, cm=True)
+    pred_bd, pred_ssw, pred_h = dimensional_default.calc_abg_dims(abg_biomass, cm=True)
     assert_allclose(pred_bd, basal_diameter, rtol=0.001)
     assert_allclose(pred_ssw, shoot_sys_width, rtol=0.001)
     assert_allclose(pred_h, height, rtol=0.001)
 
 
-def test__calc_shoot_width_from_basal_dia(example_input_params):
-    d = create_default_dimensional_object(example_input_params)
+def test__calc_shoot_width_from_basal_dia(dimensional_default):
     basal_dia = np.array([0.009, 0.036, 0.05])
     shoot_width = np.array([0.257388353, 0.639255015, 0.79304121])
-    pred_shoot_width = d._calc_shoot_width_from_basal_dia(basal_dia, cm=True)
+    pred_shoot_width = dimensional_default._calc_shoot_width_from_basal_dia(basal_dia, cm=True)
     assert_allclose(pred_shoot_width, shoot_width, rtol=0.001)
 
 
-def test__calc_height_from_basal_dia(example_input_params):
-    d = create_default_dimensional_object(example_input_params)
+def test__calc_height_from_basal_dia(dimensional_default):
     basal_dia = np.array([0.009, 0.036, 0.05])
     height = np.array([0.218283, 0.514921, 0.631049])
-    pred_height = d._calc_height_from_basal_dia(basal_dia, cm=True)
+    pred_height = dimensional_default._calc_height_from_basal_dia(basal_dia, cm=True)
     assert_allclose(pred_height, height, rtol=0.001)
 
 
-def test__calc_basal_dia_from_shoot_width(example_input_params):
-    d = create_default_dimensional_object(example_input_params)
+def test__calc_basal_dia_from_shoot_width(dimensional_default):
     basal_dia = np.array([0.009, 0.036, 0.05])
     shoot_width = np.array([0.257388353, 0.639255015, 0.79304121])
-    pred_basal_dia = d._calc_basal_dia_from_shoot_width(shoot_width, cm=True)
+    pred_basal_dia = dimensional_default._calc_basal_dia_from_shoot_width(shoot_width, cm=True)
     assert_allclose(pred_basal_dia, basal_dia, rtol=0.001)
 
 
-def test_calc_abg_biomass_from_dim(example_input_params):
-    b = create_biomass_object(example_input_params)
-    basal_dia = np.array([0.005, 0.035, 0.12])
-    shoot_sys_width = np.array([0.0640, 0.2, 0.48])
-    abg_from_bd = np.array([0, 1.5445, 13.1851])
-    abg_from_ssw = np.array([0.899, 5.156, 19.7325])
-    pred_abg_bd = b.calc_abg_biomass_from_dim(basal_dia, "basal_dia", cm=False)
-    assert_allclose(abg_from_bd, pred_abg_bd, rtol=0.001)
-    pred_abg_ssw = b.calc_abg_biomass_from_dim((0.25 * np.pi * shoot_sys_width**2), "canopy_area")
-    assert_allclose(abg_from_ssw, pred_abg_ssw, rtol=0.001)
-
-
 """
 
 
diff --git a/tests/components/genveg/test_dispersal.py b/tests/components/genveg/test_dispersal.py
new file mode 100644
index 0000000000..c9fce34e65
--- /dev/null
+++ b/tests/components/genveg/test_dispersal.py
@@ -0,0 +1,66 @@
+import numpy as np
+import pytest
+from numpy.testing import assert_allclose, assert_array_less
+from landlab.components.genveg.dispersal import Clonal, Seed, Random
+
+
+def test_clonal_disperse(clonal, example_plant):
+    # Test distance from 0, 0 is equal to pup_cost / unit_cost
+    x_loc = example_plant["x_loc"]
+    y_loc = example_plant["y_loc"]
+    plant_pup = clonal.disperse(example_plant)
+    distance = np.sqrt((plant_pup["dispersal"]["pup_x_loc"] - x_loc)**2 + (plant_pup["dispersal"]["pup_y_loc"] - y_loc)**2)
+    est_distance = plant_pup["dispersal"]["pup_cost"] / clonal.unit_cost
+    assert_allclose(distance, est_distance, rtol=0.0001)
+    # Test pup_cost is less than reproductive
+    assert (plant_pup["dispersal"]["pup_cost"] <= example_plant["reproductive"]) or np.isnan(plant_pup["dispersal"]["pup_cost"])
+
+
+def test_seed_disperse(seed, example_plant):
+    repro_initial = example_plant["reproductive"].copy()
+    x_loc = example_plant["x_loc"]
+    y_loc = example_plant["y_loc"]
+    size = seed.seed_size
+    efficiency = seed.seed_efficiency
+    total_biomass = example_plant["root"] + example_plant["leaf"] + example_plant["stem"]
+    updated_plants = seed.disperse(example_plant, total_biomass)
+    count = np.count_nonzero(~np.isnan(updated_plants["dispersal"][0]["seedling_x_loc"][:10]))
+    # Plant biomass too small to produce seeds
+    assert (count == 0)
+    biomass_seeds = count * size / efficiency
+    assert_allclose(updated_plants["reproductive"], repro_initial, rtol=0.001)
+    # Plant total biomass should produce multiple seeds
+    total_biomass = np.array([350])
+    updated_plants = seed.disperse(example_plant, total_biomass)
+    count = np.count_nonzero(~np.isnan(updated_plants["dispersal"][0]["seedling_x_loc"][:10]))
+    assert (count != 0)
+    biomass_seeds = count * size / efficiency
+    new_repro = repro_initial - biomass_seeds
+    distance = np.sqrt(
+        (updated_plants["dispersal"]["seedling_x_loc"] - x_loc)**2
+        + (updated_plants["dispersal"]["seedling_y_loc"] - y_loc)**2
+    )
+    filter = ~np.isnan(distance)
+    assert_allclose(updated_plants["reproductive"], new_repro, rtol=0.001)
+    assert_array_less(distance[filter], seed.max_dispersal_distance)
+    # Plant can use reserves to produce seedling rarely
+    example_plant["reproductive"] = 0.5
+    example_plant["dispersal"]["seedling_reserve"] = np.array([0.999996])
+    total_biomass = np.array([0.3])
+    updated_plants = seed.disperse(example_plant, total_biomass)
+    count = np.count_nonzero(~np.isnan(updated_plants["dispersal"][0]["seedling_x_loc"][:10]))
+    assert (count > 0)
+
+
+def test_random_disperse(random, one_cell_grid, example_plant_array):
+    # Test we have 
+    extent = one_cell_grid.extent
+    reference = one_cell_grid.xy_of_reference
+    random.daily_col_prob = 0.5
+    update_plant_array = random.disperse(example_plant_array)
+    filter = ~np.isnan(update_plant_array["dispersal"]["rand_x_loc"])
+    x_upper_bound = reference[1] + extent[1]
+    y_upper_bound = reference[0] + extent[0]
+    assert_array_less(update_plant_array["dispersal"]["rand_x_loc"][filter], x_upper_bound)
+    assert_array_less(update_plant_array["dispersal"]["rand_y_loc"][filter], y_upper_bound)
+
diff --git a/tests/components/genveg/test_growth.py b/tests/components/genveg/test_growth.py
new file mode 100644
index 0000000000..7695517e90
--- /dev/null
+++ b/tests/components/genveg/test_growth.py
@@ -0,0 +1,47 @@
+import numpy as np
+from numpy.testing import assert_allclose
+
+
+def test__init__methods_run(growth_obj, example_plant):
+    pass
+
+
+def test_species_plants(growth_obj, example_plant):
+    # Test to make sure this gets plant array
+    pass
+
+
+def test_species_get_variables(growth_obj, example_input_params):
+    pass
+
+
+def test_update_plants(growth_obj, example_plant_array):
+    pass
+
+
+def test_add_new_plants(growth_obj, example_plant_array, example_plant):
+    pass
+
+
+def test_grow(growth_obj, example_plant, one_cell_grid):
+    pass
+
+
+def test__init_plants_from_grid(growth_obj, one_cell_grid):
+    pass
+
+
+def test_set_event_flags(growth_obj):
+    pass
+
+
+def test_kill_small_plants(growth_obj, example_plant):
+    pass
+
+
+def test_remove_plants(growth_obj, example_plant_array):
+    pass
+
+
+def test_save_plant_output(growth_obj):
+    pass
diff --git a/tests/components/genveg/test_habit.py b/tests/components/genveg/test_habit.py
index 5d441be964..14f49c9e05 100644
--- a/tests/components/genveg/test_habit.py
+++ b/tests/components/genveg/test_habit.py
@@ -4,54 +4,25 @@
 from numpy.testing import assert_array_almost_equal
 from numpy.testing import assert_equal
 
-from landlab.components.genveg.allometry import Biomass
-from landlab.components.genveg.habit import Habit, Graminoid
-from landlab.components.genveg.species import Species
+from landlab.components.genveg.duration import Annual, Deciduous, Evergreen
 
-dt = np.timedelta64(1, 'D')
-
-
-def create_species_object(example_input_params):
-    return Species(example_input_params["BTS"], dt=dt, latitude=0.9074)
-
-
-def create_habit_object(example_input_params):
-    params = example_input_params["BTS"]
-    params["grow_params"]["abg_biomass"] = {}
-    params["morph_params"]["canopy_area"] = {}
-    params["grow_params"]["abg_biomass"]["min"] = params["grow_params"]["plant_part_min"]["leaf"] + params["grow_params"]["plant_part_min"]["stem"]
-    params["grow_params"]["abg_biomass"]["max"] = params["grow_params"]["plant_part_max"]["leaf"] + params["grow_params"]["plant_part_max"]["stem"]
-    params["morph_params"]["canopy_area"]["min"] = params["morph_params"]["shoot_sys_width"]["min"]**2 * 0.25 * np.pi
-    params["morph_params"]["canopy_area"]["max"] = params["morph_params"]["shoot_sys_width"]["max"]**2 * 0.25 * np.pi
-    allometry = Biomass(params, empirical_coeffs={"root_dia_coeffs": {"a": 0.08, "b": 0.24}})
-    return Habit(params, allometry, dt=1, green_parts=("leaf", "stem"))
 
-
-def create_graminoid_object(example_input_params):
-    params = example_input_params["BTS"]
-    params["grow_params"]["abg_biomass"] = {}
-    params["morph_params"]["canopy_area"] = {}
-    params["grow_params"]["abg_biomass"]["min"] = params["grow_params"]["plant_part_min"]["leaf"] + params["grow_params"]["plant_part_min"]["stem"]
-    params["grow_params"]["abg_biomass"]["max"] = params["grow_params"]["plant_part_max"]["leaf"] + params["grow_params"]["plant_part_max"]["stem"]
-    params["morph_params"]["canopy_area"]["min"] = params["morph_params"]["shoot_sys_width"]["min"]**2 * 0.25 * np.pi
-    params["morph_params"]["canopy_area"]["max"] = params["morph_params"]["shoot_sys_width"]["max"]**2 * 0.25 * np.pi
-    return Graminoid(params, dt=1)
+dt = np.timedelta64(1, 'D')
 
 
-def test_calc_canopy_area_from_shoot_width(example_input_params):
-    h = create_habit_object(example_input_params)
+def test_calc_canopy_area_from_shoot_width(habit_object):
     # zero array returns zero array
     assert_array_almost_equal(
-        h._calc_canopy_area_from_shoot_width(shoot_sys_width=np.array([0, 0, 0])),
+        habit_object._calc_canopy_area_from_shoot_width(shoot_sys_width=np.array([0, 0, 0])),
         np.array([0, 0, 0]),
     )
     # single input returns correct input
     assert_equal(
-        h._calc_canopy_area_from_shoot_width(shoot_sys_width=0.325), 0.08295768100885548
+        habit_object._calc_canopy_area_from_shoot_width(shoot_sys_width=0.325), 0.08295768100885548
     )
     # an array of values
     assert_allclose(
-        h._calc_canopy_area_from_shoot_width(
+        habit_object._calc_canopy_area_from_shoot_width(
             shoot_sys_width=np.array([0, 0.0004, 0.678, 1.5, 3])
         ),
         np.array(
@@ -66,49 +37,63 @@ def test_calc_canopy_area_from_shoot_width(example_input_params):
     )
 
 
-def test__calc_diameter_from_area(example_input_params):
-    h = create_habit_object(example_input_params)
+def test__calc_diameter_from_area(habit_object):
     canopy_area = np.array([0.1, 1.28, 3.7])
     shoot_width = np.array([0.35682, 1.27662, 2.17048])
-    pred_shoot_width = h._calc_diameter_from_area(canopy_area)
+    pred_shoot_width = habit_object._calc_diameter_from_area(canopy_area)
     assert_array_almost_equal(shoot_width, pred_shoot_width, decimal=5)
 
 
-def test__calc_canopy_volume(example_input_params):
-    h = create_habit_object(example_input_params)
+def test__calc_canopy_volume(habit_object):
     shoot_width = np.array([0.08, 0.35, 0.72])
     height = np.array([0.1, 0.5, 1.7])
     basal_dia = np.array([0.005, 0.02, 0.1])
     volume = np.array([0.00017868, 0.017004, 0.2672])
-    pred_vol = h._calc_canopy_volume(shoot_width, basal_dia, height)
+    pred_vol = habit_object._calc_canopy_volume(shoot_width, basal_dia, height)
     assert_allclose(pred_vol, volume, rtol=0.0001)
 
 
-def test_calc_root_sys_width(example_input_params):
-    h = create_habit_object(example_input_params)
+def test_calc_root_sys_width(habit_object):
     shoot_width = np.array([0.08, 0.35, 0.72])
     height = np.array([0.1, 0.5, 1.7])
     basal_dia = np.array([0.005, 0.02, 0.1])
     root_width = np.array([0.080043, 0.084081, 0.144128])
-    pred_root_width = h.calc_root_sys_width(shoot_width, basal_dia, height)
+    pred_root_width = habit_object.calc_root_sys_width(shoot_width, basal_dia, height)
     assert_allclose(pred_root_width, root_width, rtol=0.0001)
 
 
-def test_estimate_abg_biomass_from_cover(example_input_params, example_plant_array):
-    h = create_habit_object(example_input_params)
+def test_estimate_abg_biomass_from_cover_graminoid(graminoid_object, example_plant_array):
     abg_biomass = example_plant_array["leaf"] + example_plant_array["stem"]
-    example_plant_array["basal_dia"], example_plant_array["shoot_sys_width"], example_plant_array["shoot_sys_height"] = h.calc_abg_dims_from_biomass(abg_biomass)
-    est_abg_biomass = h.estimate_abg_biomass_from_cover(example_plant_array)
+    example_plant_array["basal_dia"], example_plant_array["shoot_sys_width"], example_plant_array["shoot_sys_height"] = graminoid_object.calc_abg_dims_from_biomass(abg_biomass)
+    est_abg_biomass = graminoid_object.estimate_abg_biomass_from_cover(example_plant_array)
     assert_allclose(est_abg_biomass, abg_biomass, rtol=0.0001)
-    g = create_graminoid_object(example_input_params)
+
+
+def test_estimate_abg_biomass_from_cover(habit_object, example_plant_array):
     abg_biomass = example_plant_array["leaf"] + example_plant_array["stem"]
-    example_plant_array["basal_dia"], example_plant_array["shoot_sys_width"], example_plant_array["shoot_sys_height"] = h.calc_abg_dims_from_biomass(abg_biomass)
-    est_abg_biomass = g.estimate_abg_biomass_from_cover(example_plant_array)
+    example_plant_array["basal_dia"], example_plant_array["shoot_sys_width"], example_plant_array["shoot_sys_height"] = habit_object.calc_abg_dims_from_biomass(abg_biomass)
+    est_abg_biomass = habit_object.estimate_abg_biomass_from_cover(example_plant_array)
     assert_allclose(est_abg_biomass, abg_biomass, rtol=0.0001)
 
 
-def test_calc_canopy_area_from_shoot_width_raises_error(example_input_params):
-    h = create_habit_object(example_input_params)
+# change this to ensure zero in case of negative
+def test_calc_canopy_area_from_shoot_width_raises_error(habit_object):
     with pytest.raises(ValueError):
-        h._calc_canopy_area_from_shoot_width(-1.5)
-        h._calc_canopy_area_from_shoot_width(np.array([0, 0.004, -0.678, 1.5, 3]))
+        habit_object._calc_canopy_area_from_shoot_width(-1.5)
+        habit_object._calc_canopy_area_from_shoot_width(np.array([0, 0.004, -0.678, 1.5, 3]))
+
+
+def test_select_duration_class(habit_object, example_input_params):
+    dummy_habit = example_input_params
+    dt = np.timedelta64(1, 'D')
+    for duration_opt, cls in zip(["annual", "perennial deciduous", "perennial evergreen"], [Annual, Deciduous, Evergreen]):
+        dummy_habit["BTS"]["plant_factors"]["duration"] = duration_opt
+        assert isinstance(
+            habit_object._select_duration_class(
+                dummy_habit["BTS"]["grow_params"],
+                dummy_habit["BTS"]["duration_params"],
+                dt,
+                duration_opt, green_parts=("leaf")
+            ),
+            cls,
+        )
diff --git a/tests/components/genveg/test_integrator.py b/tests/components/genveg/test_integrator.py
index 62ff40c490..9ca9197312 100644
--- a/tests/components/genveg/test_integrator.py
+++ b/tests/components/genveg/test_integrator.py
@@ -18,3 +18,23 @@ def test_check_for_dispersal_success(
     # filter = np.nonzero(not np.isnan(example_plant_array["pup_x_loc"]))
     (dispersed_size,) = dispersed_plant_array.shape
     return example_plant_array.shape
+
+
+def test__init__(genveg_obj, one_cell_grid):
+    # test initialization routine - consider breaking this up
+    pass
+
+
+def test_calculate_grid_vars(genveg_obj, one_cell_grid):
+    pass
+
+
+def test_check_for_grid_fields(genveg_obj, one_cell_grid):
+    pass
+
+
+def test_run_one_step(genveg_obj, one_cell_grid, example_plant):
+    pass
+
+
+
diff --git a/tests/components/genveg/test_species.py b/tests/components/genveg/test_species.py
index 85c632908c..5eee07923c 100644
--- a/tests/components/genveg/test_species.py
+++ b/tests/components/genveg/test_species.py
@@ -2,34 +2,16 @@
 import pytest
 from numpy.testing import assert_allclose, assert_almost_equal, assert_array_less
 
-from landlab.components.genveg.form import Bunch
-from landlab.components.genveg.form import Colonizing
-from landlab.components.genveg.form import Multiplestems
-from landlab.components.genveg.form import Rhizomatous
-from landlab.components.genveg.form import Singlecrown
-from landlab.components.genveg.form import Singlestem
-from landlab.components.genveg.form import Stoloniferous
-from landlab.components.genveg.form import Thicketforming
-from landlab.components.genveg.habit import Forbherb
-from landlab.components.genveg.habit import Graminoid
-from landlab.components.genveg.habit import Shrub
-from landlab.components.genveg.habit import Tree
-from landlab.components.genveg.habit import Vine
+from landlab.components.genveg.species import Species
 from landlab.components.genveg.photosynthesis import C3
 from landlab.components.genveg.photosynthesis import C4
 from landlab.components.genveg.photosynthesis import Cam
-from landlab.components.genveg.species import Species
-
-dt = np.timedelta64(1, 'D')
-
+from landlab.components.genveg.dispersal import Clonal, Seed, Random
+from landlab.components.genveg.habit import Forbherb, Graminoid, Shrub
 
-def create_species_object(example_input_params, dt=dt):
-    return Species(example_input_params["BTS"], dt=dt, latitude=0.9074)
 
-
-def test_get_daily_nsc_concentration(example_input_params):
-    species_object = create_species_object(example_input_params)
-    days = example_input_params["BTS"]["duration_params"]
+def test_get_daily_nsc_concentration(species_object):
+    days = species_object.species_duration_params
     parts = ["root", "leaf", "reproductive", "stem"]
     nsc_content_gs_start_actual = {
         "root": 276.9606782 / 1000,
@@ -55,9 +37,8 @@ def test_get_daily_nsc_concentration(example_input_params):
         )
 
 
-def test_calc_area_of_circle(example_input_params):
-    species_object = create_species_object(example_input_params)
-    morph_params = example_input_params["BTS"]["morph_params"]
+def test_calc_area_of_circle(species_object):
+    morph_params = species_object.species_morph_params
     m_params = [
         "shoot_sys_width",
         "root_sys_width",
@@ -77,21 +58,18 @@ def test_calc_area_of_circle(example_input_params):
         idx += 1
 
 
-def test_calculate_dead_age(example_input_params):
+def test_calculate_dead_age(species_object):
     age_t1 = np.array([60, 60, 60])
     mass_t1 = np.array([2, 2, 2])
     mass_t2 = np.array([0, 3, 4])
     age_t2 = np.array([60, 40, 30])
-    calc_age_t2 = create_species_object(example_input_params).calculate_dead_age(age_t1, mass_t1, mass_t2)
+    calc_age_t2 = species_object.calculate_dead_age(age_t1, mass_t1, mass_t2)
     assert_almost_equal(age_t2, calc_age_t2, decimal=6)
 
 
-def test_calculate_shaded_leaf_mortality(example_plant_array, example_plant, example_input_params):
-    species_object = create_species_object(example_input_params)
+def test_calculate_shaded_leaf_mortality(example_plant_array, example_plant, species_object):
     # Check to make sure no leaf mortality occurred for LAI 0.012 < LAI_crit 2
     init_leaf_weight = example_plant_array["leaf"].copy()
-    print(example_plant_array["total_leaf_area"])
-    print(example_plant_array["total_leaf_area"] / (np.pi / 4 * example_plant_array["shoot_sys_width"]**2))
     plant_out = species_object.calculate_shaded_leaf_mortality(example_plant_array)
     assert_almost_equal(example_plant_array["leaf"][0:2], init_leaf_weight[0:2])
     assert_array_less(plant_out["leaf"][2], init_leaf_weight[2])
@@ -105,33 +83,57 @@ def test_calculate_shaded_leaf_mortality(example_plant_array, example_plant, exa
     assert_almost_equal(weight_change, wofost_leaf_weight_change, 6)
 
 
-def test_calculate_whole_plant_mortality(example_plant_array, one_cell_grid, example_input_params):
-    max_temp_dt = np.timedelta64(example_input_params["BTS"]["mortality_params"]["duration"]["1"], 'D')
-    min_temp_dt = np.timedelta64(example_input_params["BTS"]["mortality_params"]["duration"]["1"], 'D')
-    species_object_max = create_species_object(example_input_params, dt=max_temp_dt)
-    species_object_min = create_species_object(example_input_params, dt=min_temp_dt)
-    max_temp = one_cell_grid["cell"]["air__max_temperature_C"][example_plant_array["cell_index"]]
-    min_temp = one_cell_grid["cell"]["air__min_temperature_C"][example_plant_array["cell_index"]]
-    # Make sure plant doesn't die at moderate ambient temp
-    new_biomass = species_object_max.calculate_whole_plant_mortality(example_plant_array, max_temp, "1")
-    assert_allclose(new_biomass["root"], example_plant_array["root"], rtol=0.0001)
-    # Change max temp and make sure plant dies
+def test_calculate_whole_plant_mortality_max_temp(example_plant, one_cell_grid, species_object):
+    max_temp_dt = species_object.species_mort_params["duration"]["1"]
+    max_temp = one_cell_grid["cell"]["air__max_temperature_C"][example_plant["cell_index"]]
+    new_biomass = example_plant
+    n_runs = 0
+    while new_biomass["root"] > np.zeros_like(example_plant["root"]):
+        new_biomass = species_object.calculate_whole_plant_mortality(example_plant, max_temp, "1")
+        n_runs += 1
+        if n_runs >= (1.5 * max_temp_dt):
+            break
+    assert_allclose(n_runs, (max_temp_dt * 1.5), rtol=0.01)
     one_cell_grid["cell"]["air__max_temperature_C"] = 38 * np.ones(one_cell_grid.number_of_cells)
-    max_temp = one_cell_grid["cell"]["air__max_temperature_C"]
-    new_biomass = species_object_max.calculate_whole_plant_mortality(example_plant_array, max_temp, "1")
-    assert_allclose(new_biomass["root"], np.zeros_like(new_biomass["root"]), rtol=0.0001)
+    max_temp = one_cell_grid["cell"]["air__max_temperature_C"][example_plant["cell_index"]]
+    n_runs = 0
+    while new_biomass["root"] > np.zeros_like(example_plant["root"]):
+        new_biomass = species_object.calculate_whole_plant_mortality(example_plant, max_temp, "1")
+        n_runs += 1
+        if n_runs >= (1.5 * max_temp_dt):
+            break
+    assert n_runs < (1.5 * max_temp_dt)
+
+
+def test_calculate_whole_plant_mortality_min_temp(example_plant, one_cell_grid, species_object):
     # Change min temp and make sure plant dies
+    min_temp_dt = species_object.species_mort_params["duration"]["2"]
+    min_temp = one_cell_grid["cell"]["air__min_temperature_C"][example_plant["cell_index"]]
+    new_biomass = example_plant
+    n_runs = 0
+    while new_biomass["root"] > np.zeros_like(example_plant["root"]):
+        new_biomass = species_object.calculate_whole_plant_mortality(example_plant, min_temp, "2")
+        n_runs += 1
+        if n_runs >= (1.5 * min_temp_dt):
+            break
+    assert_allclose(n_runs, (min_temp_dt * 1.5), rtol=0.01)
     one_cell_grid["cell"]["air__min_temperature_C"] = -38 * np.ones(one_cell_grid.number_of_cells)
-    min_temp = one_cell_grid["cell"]["air__min_temperature_C"]
-    new_biomass = species_object_min.calculate_whole_plant_mortality(example_plant_array, min_temp, "2")
-    assert_allclose(new_biomass["root"], np.zeros_like(new_biomass["root"]), rtol=0.0001)
+    min_temp = one_cell_grid["cell"]["air__min_temperature_C"][example_plant["cell_index"]]
+    n_runs = 0
+    while new_biomass["root"] > np.zeros_like(example_plant["root"]):
+        new_biomass = species_object.calculate_whole_plant_mortality(example_plant, min_temp, "2")
+        n_runs += 1
+        if n_runs >= (1.5 * min_temp_dt):
+            break
+    assert n_runs < (1.5 * min_temp_dt)
 
 
-def test_update_morphology(example_input_params, example_plant_array):
+def test_update_morphology(example_input_params, example_plant_array, one_cell_grid):
     # This method should update basal diameter, shoot width, shoot height, root
     # width, live leaf area, dead leaf area
     example_input_params["BTS"]["morph_params"]["allometry_method"] = "default"
-    sp_obj = create_species_object(example_input_params)
+    dt = np.timedelta64(1, 'D')
+    sp_obj = Species(example_input_params["BTS"], one_cell_grid, dt=dt)
     # change to default
     pred_example_plant_array = sp_obj.update_morphology(example_plant_array)
     known_plants = {}
@@ -148,14 +150,12 @@ def test_update_morphology(example_input_params, example_plant_array):
         "live_leaf_area",
     ]
     for var in update_vars:
-        print(var)
         assert_allclose(known_plants[var], pred_example_plant_array[var], rtol=0.001)
 
 
-def test_populate_biomass_allocation_array(example_input_params):
-    s = create_species_object(example_input_params)
+def test_populate_biomass_allocation_array(species_object):
     biomass_array_from_excel = np.array([
-        [0.01, 0.11,0.21,0.31,0.41,0.51,0.61,0.71,0.81,0.91,1.01,1.11,1.21,1.31,1.41,1.51,1.61,1.71,1.81,1.91,2.01,2.11,2.21,2.31,2.41,2.51,2.61,2.71,2.81,2.91,3.01,3.11,3.21,3.31,3.41,3.51,3.61,3.71,3.81,3.91,4.01,4.11,4.21,4.31],
+        [0.01, 0.11, 0.21,0.31,0.41,0.51,0.61,0.71,0.81,0.91,1.01,1.11,1.21,1.31,1.41,1.51,1.61,1.71,1.81,1.91,2.01,2.11,2.21,2.31,2.41,2.51,2.61,2.71,2.81,2.91,3.01,3.11,3.21,3.31,3.41,3.51,3.61,3.71,3.81,3.91,4.01,4.11,4.21,4.31],
         [0.029531558, 0.314126222,0.597580079,0.881430388,1.165852132,1.450873064,1.736484358,2.022666426,2.309396945,2.596653692,2.884415563,3.1726629,3.461377538,3.750542729,4.040143031,4.330164189,4.620593014,4.911417273,5.202625596,5.494207382,5.786152728,6.078452358,6.371097567,6.664080162,6.957392425,7.251027061,7.544977173,7.839236221,8.133797995,8.428656593,8.723806393,9.019242035,9.314958401,9.610950598,9.907213944,10.20374395,10.50053631,10.7975869,11.09489174,11.39244701,11.69024903,11.98829426,12.28657927,12.58510076],
         [1.283164616, 1.143451734,1.118287452,1.091405185,1.073974314,1.061083353,1.050870103,1.042423648,1.035229826,1.028969927,1.023432791,1.018471342,1.013979073,1.009876397,1.006102299,1.002608985,0.999358326,0.996319423,0.993466895,0.990779636,0.988239915,0.98583269,0.983545095,0.981366039,0.979285895,0.977296253,0.975389729,0.973559797,0.971800668,0.970107179,0.968474711,0.96689911,0.965376631,0.963903883,0.962477789,0.961095547,0.959754597,0.958452596,0.957187394,0.955957012,0.954759627,0.953593552,0.952457226,0.951349204],
         [0.613402864, 0.691408432,0.716251113,0.747097905,0.770243126,0.789125974,0.805242837,0.819397029,0.832075367,0.84359754,0.854185916,0.864002036,0.873167308,0.881775509,0.889900723,0.897602596,0.90492992,0.911923172,0.918616329,0.925038223,0.931213544,0.937163614,0.942906988,0.948459918,0.953836728,0.959050116,0.96411139,0.969030676,0.973817075,0.978478798,0.983023289,0.98745731,0.991787029,0.99601809,1.000155667,1.004204522,1.008169042,1.012053284,1.015861004,1.01959569,1.023260584,1.026858708,1.030392882,1.033865742],
@@ -176,24 +176,23 @@ def test_populate_biomass_allocation_array(example_input_params):
     idx = 0
     for var in allocation_array_cols:
         # Using larger relative tolerance because knowns are approximations not direct solve
-        assert_allclose(s.biomass_allocation_array[var], biomass_array_from_excel[idx], rtol=0.05)
+        assert_allclose(species_object.biomass_allocation_array[var], biomass_array_from_excel[idx], rtol=0.05)
         idx += 1
 
 
-def test_set_initial_cover(example_input_params):
-    s = create_species_object(example_input_params)
+def test_set_initial_cover(species_object):
     species_name = "BTS"
     pidval = 0
     cell_index = 1
     plantlist = []
     # Cover area is too small for any plants
     cover_area = 0.00001
-    plantlist = s.set_initial_cover(cover_area, species_name, pidval, cell_index, plantlist)
+    plantlist = species_object.set_initial_cover(cover_area, species_name, pidval, cell_index, plantlist)
     assert len(plantlist) == 0
     assert pidval == 0
     # Correct dimension is populated
     cover_area = 0.2
-    plantlist = s.set_initial_cover(cover_area, species_name, pidval, cell_index, plantlist)
+    plantlist = species_object.set_initial_cover(cover_area, species_name, pidval, cell_index, plantlist)
     basal_dias = plantlist[0][18]
     assert basal_dias > 0
     # Area being used properly
@@ -201,26 +200,25 @@ def test_set_initial_cover(example_input_params):
     for i in range(len(plantlist)):
         inc_area = 0.25 * np.pi * plantlist[i][18]**2
         sum_area += inc_area
-    min_cover_area = 0.25 * np.pi * s.species_morph_params["basal_dia"]["min"]**2
+    min_cover_area = 1.2 * 0.25 * np.pi * species_object.species_morph_params["basal_dia"]["min"]**2
     assert sum_area < cover_area
     assert (cover_area - sum_area) < min_cover_area
 
 
-def test_set_initial_biomass(example_input_params, example_plant):
+def test_set_initial_biomass(species_object, example_plant):
     # Testing for Deciduous set_initial_biomass, testing for other classes handled under test_emerge
-    s = create_species_object(example_input_params)
     in_growing_season = True
     abg_biomass = example_plant["leaf"] + example_plant["stem"]
     root = 0.433864
     leaf = 0.485539
     stem = 0.314443
-    min_repro = example_input_params["BTS"]["grow_params"]["plant_part_min"]["reproductive"]
-    max_repro = example_input_params["BTS"]["grow_params"]["plant_part_max"]["reproductive"]
-    example_plant["basal_dia"], shoot_sys_width, height = s.habit.calc_abg_dims_from_biomass(abg_biomass)
+    min_repro = species_object.species_grow_params["plant_part_min"]["reproductive"]
+    max_repro = species_object.species_grow_params["plant_part_max"]["reproductive"]
+    example_plant["basal_dia"], shoot_sys_width, height = species_object.habit.calc_abg_dims_from_biomass(abg_biomass)
     example_plant["root"] = np.zeros_like(example_plant["root"])
     example_plant["leaf"] = np.zeros_like(example_plant["leaf"])
     example_plant["stem"] = np.zeros_like(example_plant["stem"])
-    emerged_plant = s.set_initial_biomass(example_plant, in_growing_season)
+    emerged_plant = species_object.set_initial_biomass(example_plant, in_growing_season)
     assert_allclose(emerged_plant["root"], root, rtol=0.001)
     assert_allclose(emerged_plant["leaf"], leaf, rtol=0.001)
     assert_allclose(emerged_plant["stem"], stem, rtol=0.001)
@@ -229,7 +227,7 @@ def test_set_initial_biomass(example_input_params, example_plant):
     assert_allclose(emerged_plant["shoot_sys_width"], shoot_sys_width, rtol=0.001)
     assert_allclose(emerged_plant["shoot_sys_height"], height, rtol=0.001)
     in_growing_season = False
-    emerged_plant = s.set_initial_biomass(example_plant, in_growing_season)
+    emerged_plant = species_object.set_initial_biomass(example_plant, in_growing_season)
     zeros = np.zeros_like(example_plant["root"])
     assert_allclose(emerged_plant["root"], root, rtol=0.001)
     assert_allclose(emerged_plant["leaf"], zeros, rtol=0.001)
@@ -238,8 +236,7 @@ def test_set_initial_biomass(example_input_params, example_plant):
     assert_allclose(emerged_plant["shoot_sys_height"], zeros, rtol=0.001)
 
 
-def test_sum_plant_parts(example_plant, example_input_params):
-    s = create_species_object(example_input_params)
+def test_sum_plant_parts(example_plant, species_object):
     example_plant["reproductive"] = np.array([0.32])
     example_plant["dead_root"] = np.array([0.15])
     example_plant["dead_leaf"] = np.array([0.12])
@@ -255,12 +252,11 @@ def test_sum_plant_parts(example_plant, example_input_params):
         "dead_aboveground": 0.21,
     }
     for sum_type in sums:
-        sum_calc = s.sum_plant_parts(example_plant, sum_type)
+        sum_calc = species_object.sum_plant_parts(example_plant, sum_type)
         assert_allclose(sum_calc, sums[sum_type], rtol=0.001)
 
 
-def test_emerge(example_plant, example_input_params):
-    species_object = create_species_object(example_input_params)
+def test_emerge(example_plant, species_object):
     example_plant["leaf"] = np.zeros_like(example_plant["leaf"])
     example_plant["stem"] = np.zeros_like(example_plant["stem"])
     jday = 195
@@ -303,12 +299,11 @@ def test_emerge(example_plant, example_input_params):
         assert_allclose(check_plant[var], emerged_plant[var], rtol=0.001)
 
 
-def test__adjust_biomass_allocation_towards_ideal(example_input_params, example_plant):
-    s = create_species_object(example_input_params)
+def test__adjust_biomass_allocation_towards_ideal(species_object, example_plant):
     new_root = np.array([0.783358])
     new_leaf = np.array([0.508459])
     new_stem = np.array([0.308183])
-    adjust_plant = s._adjust_biomass_allocation_towards_ideal(example_plant)
+    adjust_plant = species_object._adjust_biomass_allocation_towards_ideal(example_plant)
     assert_allclose(adjust_plant["root"], new_root, rtol=0.001)
     assert_allclose(adjust_plant["leaf"], new_leaf, rtol=0.001)
     assert_allclose(adjust_plant["stem"], new_stem, rtol=0.001)
@@ -317,23 +312,22 @@ def test__adjust_biomass_allocation_towards_ideal(example_input_params, example_
     small_root_plant["root"] = np.array([0.005])
     small_root_plant["leaf"] = np.array([0.5])
     small_root_plant["stem"] = np.array([0.2999])
-    s.species_grow_params["plant_part_min"]["root"] = 0.03
+    species_object.species_grow_params["plant_part_min"]["root"] = 0.03
     new_root = np.array([0.030000])
     new_leaf = np.array([0.484195])
     new_stem = np.array([0.290805])
-    adjust_plant = s._adjust_biomass_allocation_towards_ideal(small_root_plant)
+    adjust_plant = species_object._adjust_biomass_allocation_towards_ideal(small_root_plant)
     assert_allclose(adjust_plant["root"], new_root, rtol=0.001)
     assert_allclose(adjust_plant["leaf"], new_leaf, rtol=0.001)
     assert_allclose(adjust_plant["stem"], new_stem, rtol=0.001)
 
 
-def test_allocate_biomass_dynamically(example_input_params, example_plant):
-    s = create_species_object(example_input_params)
+def test_allocate_biomass_dynamically(species_object, example_plant):
     new_root = np.array([0.789284])
     new_leaf = np.array([0.514743])
     new_stem = np.array([0.312879])
-    new_repro = np.array([0.000104917])
-    allocated_plant = s.allocate_biomass_dynamically(example_plant, np.array([0.025]))
+    new_repro = np.array([0.300104917])
+    allocated_plant = species_object.allocate_biomass_dynamically(example_plant, np.array([0.025]))
     assert_allclose(allocated_plant["root"], new_root, rtol=0.001)
     assert_allclose(allocated_plant["leaf"], new_leaf, rtol=0.001)
     assert_allclose(allocated_plant["stem"], new_stem, rtol=0.001)
@@ -342,7 +336,7 @@ def test_allocate_biomass_dynamically(example_input_params, example_plant):
 
 def test_mortality(example_input_params, two_cell_grid, example_plant_array):
     temp_dt = np.timedelta64(example_input_params["BTS"]["mortality_params"]["duration"]["2"], 'D')
-    species_object = create_species_object(example_input_params, dt=temp_dt)
+    species_object = Species(example_input_params["BTS"], two_cell_grid, dt=temp_dt)
     example_plant_array["cell_index"][4:] = np.array([1, 1, 1, 1])
     inital_leaf = example_plant_array["leaf"].copy()
     initial_plants = example_plant_array["root"].copy()
@@ -366,8 +360,7 @@ def test_mortality(example_input_params, two_cell_grid, example_plant_array):
     assert_almost_equal(plant_out["leaf"], np.zeros_like(plant_out["leaf"]), decimal=6)
 
 
-def test_respire(example_plant, one_cell_grid, example_input_params):
-    species_object = create_species_object(example_input_params)
+def test_respire(example_plant, one_cell_grid, species_object):
     # Respire with no water growth limitation from WOFOST in PCSE and crop development stage=1 divided by carb conversion
     wofost_resp_root = np.array([0.0033926])
     wofost_resp_stem = np.array([0.00121421])
@@ -398,8 +391,7 @@ def test_respire(example_plant, one_cell_grid, example_input_params):
     assert_almost_equal(wofost_resp_stem, mass_change_resp_stem_half_sat, 6)
 
 
-def test_sum_vars_in_calculate_derived_params(example_input_params):
-    species_object = create_species_object(example_input_params)
+def test_sum_vars_in_calculate_derived_params(species_object, example_input_params):
     species_param = species_object.calculate_derived_params(example_input_params["BTS"])
     # Checked via excel
     # Max total Biomass
@@ -418,15 +410,15 @@ def test_sum_vars_in_calculate_derived_params(example_input_params):
     assert_almost_equal(species_param["grow_params"]["nsc_biomass"]["min"], 0.03369)
 
 
-def test_senesce(example_input_params, example_plant):
-    species_object = create_species_object(example_input_params)
+def test_senesce(species_object, example_plant):
     jday = 195
     # leaves and stems should move nonstructural carb content to roots at a fixed rate
     # calculated values from Excel at day 195 for one plant
+    example_plant["reproductive"] = np.array([0.0])
     end_root = np.array([0.803921028])
     end_stem = np.array([0.29399902])
     end_leaf = np.array([0.485])
-    end_repro = np.array([0.])
+    end_repro = np.array([0.0])
     plant_out = species_object.senesce(example_plant, jday)
     assert_almost_equal(plant_out["reproductive"], end_repro, decimal=6)
     assert_almost_equal(plant_out["root"], end_root, decimal=6)
@@ -434,8 +426,7 @@ def test_senesce(example_input_params, example_plant):
     assert_almost_equal(plant_out["leaf"], end_leaf, decimal=6)
 
 
-def test_nsc_rate_change_per_season_and_part(example_input_params):
-    species_object = create_species_object(example_input_params)
+def test_nsc_rate_change_per_season_and_part(species_object, example_input_params):
     species_param = species_object.calculate_derived_params(example_input_params["BTS"])
     ncs_rate_change = species_param["duration_params"]["nsc_rate_change"]
     # winter_nsc_rate
@@ -480,51 +471,35 @@ def test_nsc_rate_change_per_season_and_part(example_input_params):
     assert_almost_equal(ncs_rate_change["fall_nsc_rate"]["stem"], 0.015625)
 
 
-def test_select_habit_class(example_input_params):
-    species_object = create_species_object(example_input_params)
+def test_select_habit_class(species_object, example_input_params):
     dummy_species = example_input_params
     for spec, cls in zip(
-        ["forb_herb", "graminoid", "shrub", "tree", "vine"],
-        [Forbherb, Graminoid, Shrub, Tree, Vine],
+        ["forb_herb", "graminoid", "shrub"],
+        [Forbherb, Graminoid, Shrub],
     ):
         dummy_species["BTS"]["plant_factors"]["growth_habit"] = spec
         print(dummy_species["BTS"]["plant_factors"]["growth_habit"])
         assert isinstance(species_object.select_habit_class(dummy_species["BTS"]), cls)
 
 
-def test_select_form_class(example_input_params):
-    species_object = create_species_object(example_input_params)
-    dummy_growth_form = example_input_params
+def test_select_dispersal_type(species_object, example_input_params, one_cell_grid):
+    dummy_dispersal = example_input_params
     for growth, cls in zip(
-        [
-            "bunch",
-            "colonizing",
-            "multiple_stems",
-            "rhizomatous",
-            "single_crown",
-            "single_stem",
-            "stoloniferous",
-            "thicket_forming",
-        ],
-        [
-            Bunch,
-            Colonizing,
-            Multiplestems,
-            Rhizomatous,
-            Singlecrown,
-            Singlestem,
-            Stoloniferous,
-            Thicketforming,
-        ],
+        ["clonal", "seed", "random",],
+        [Clonal, Seed, Random,],
     ):
-        dummy_growth_form["BTS"]["plant_factors"]["growth_form"] = growth
-        assert isinstance(
-            species_object.select_form_class(dummy_growth_form["BTS"]), cls
+        dummy_dispersal["BTS"]["plant_factors"]["reproductive_modes"] = [growth]
+        dispersal_classes = species_object.select_dispersal_type(
+            dummy_dispersal["BTS"],
+            one_cell_grid.area_of_cell,
+            one_cell_grid.extent,
+            one_cell_grid.xy_of_reference,
         )
+        for dispersal_class in dispersal_classes:
+            assert isinstance(dispersal_class, cls)
 
 
-def test_select_photosythesis_type(example_input_params):
-    species_object = create_species_object(example_input_params)
+def test_select_photosythesis_type(species_object, example_input_params):
     dummy_photosythesis = example_input_params
     for photosynthesis_opt, cls in zip(["C3", "C4", "cam"], [C3, C4, Cam]):
         dummy_photosythesis["BTS"]["plant_factors"]["p_type"] = photosynthesis_opt
@@ -536,22 +511,19 @@ def test_select_photosythesis_type(example_input_params):
         )
 
 
-def test_update_dead_biomass(example_input_params, example_plant):
+def test_update_dead_biomass(species_object, example_plant):
     example_plant_new = example_plant.copy()
     example_plant_new["leaf"] *= 0.5
     example_plant_new["root"] *= 0.5
     example_plant_new["stem"] *= 0.5
-    species_object = create_species_object(example_input_params)
     example_plant_new = species_object.update_dead_biomass(example_plant_new, example_plant)
     assert_almost_equal(example_plant_new["dead_leaf"], example_plant_new["leaf"], decimal=6)
     assert_almost_equal(example_plant_new["dead_root"], example_plant_new["root"])
     assert_almost_equal(example_plant_new["dead_stem"], example_plant_new["stem"])
 
 
-def test_validate_plant_factors_raises_errors(example_input_params):
+def test_validate_plant_factors_raises_errors(species_object):
     # create species class (Note this will run test_validate_plant_factors during init
-    species_object = create_species_object(example_input_params)
-
     # create an unexpected variable
     dummy_plant_factor = {"unexpected-variable-name": []}
     with pytest.raises(KeyError) as excinfo:
@@ -571,10 +543,7 @@ def test_validate_plant_factors_raises_errors(example_input_params):
     assert str(excinfo.value) == "Invalid p_type option"
 
 
-def test_validate_duration_params_raises_errors(example_input_params):
-    # create Species class (note this will run validate_duration_params during init)
-    species_object = create_species_object(example_input_params)
-
+def test_validate_duration_params_raises_errors(species_object):
     # ValueError msg for growing_season_start is between 1-365
     start_below_bound = {"growing_season_start": 0}
     start_above_bound = {"growing_season_start": 367}
@@ -602,7 +571,7 @@ def test_validate_duration_params_raises_errors(example_input_params):
     assert str(excinfo.value) == "Start of senescence must be within the growing season"
 
 
-def test_litter_decomp_new_biomass_values(example_input_params, example_plant_array):
+def test_litter_decomp_new_biomass_values(species_object, example_plant_array):
     # expected values obtained from excel
     expected_new_biomass = {
         'dead_root': np.array([0.011654923, 0.011654923, 0.011654923, 0.011654923, 1.631689185, 1.631689185, 1.631689185, 1.631689185]),
@@ -614,14 +583,8 @@ def test_litter_decomp_new_biomass_values(example_input_params, example_plant_ar
         'dead_leaf_age': np.array([1330.561664, 1353.338509, 1239.787392, 1434.019734, 1349.797247, 736.8478796, 1178.379645, 235.6844825]),
         'dead_reproductive_age': np.array([508.8777074, 1122.713079, 205.4386012, 1074.25249, 661.3439315, 427.4383309, 1262.53064, 1086.549562])
     }
-
-    # initialize class
-    species_object = create_species_object(example_input_params)
-    # set dt
-    species_object.dt = np.timedelta64(1, "D")
     # call litter_decomp
     new_biomass = species_object.litter_decomp(example_plant_array)
-
     # tests
     for k, expected_value in expected_new_biomass.items():
         assert_allclose(
@@ -629,11 +592,7 @@ def test_litter_decomp_new_biomass_values(example_input_params, example_plant_ar
         )
 
 
-def test_litter_decomp_bad_values_replaced(example_input_params, example_plant_array):
-    # initialize class
-    species_object = create_species_object(example_input_params)
-    # set dt
-    species_object.dt = np.timedelta64(1, "D")
+def test_litter_decomp_bad_values_replaced(species_object, example_plant_array):
     # make a nan value, negative value, and inf value
     bad_epa_values = example_plant_array
     bad_epa_values["dead_root"][0] = -0.0125