diff --git a/main.cpp b/main.cpp
index b2f12ac..55ce06c 100644
--- a/main.cpp
+++ b/main.cpp
@@ -262,6 +262,7 @@ struct PoPSOptions
     struct Option *probability, *probability_series;
     struct Option* spread_rate_output;
     struct Option *quarantine, *quarantine_output, *quarantine_directions;
+    struct Option *dispersers_output, *established_dispersers_output;
     struct Option *output_frequency, *output_frequency_n;
 };
 
@@ -739,6 +740,20 @@ int main(int argc, char* argv[])
     opt.percent_natural_dispersal->options = "0-1";
     opt.percent_natural_dispersal->guisection = _("Dispersal");
 
+    opt.dispersers_output = G_define_standard_option(G_OPT_R_OUTPUT);
+    opt.dispersers_output->key = "dispersers_output";
+    opt.dispersers_output->label = _("Output raster of disperses");
+    opt.dispersers_output->description = _("Dispersers are accumulated over all steps and stochastic runs");
+    opt.dispersers_output->required = NO;
+    opt.dispersers_output->guisection = _("Output");
+
+    opt.established_dispersers_output = G_define_standard_option(G_OPT_R_OUTPUT);
+    opt.established_dispersers_output->key = "established_dispersers_output";
+    opt.established_dispersers_output->label = _("Output raster of established disperses");
+    opt.established_dispersers_output->description = _("Dispersers are accumulated over all steps and stochastic runs");
+    opt.established_dispersers_output->required = NO;
+    opt.established_dispersers_output->guisection = _("Output");
+
     opt.infected_to_dead_rate = G_define_option();
     opt.infected_to_dead_rate->type = TYPE_DOUBLE;
     opt.infected_to_dead_rate->key = "mortality_rate";
@@ -1274,6 +1289,9 @@ int main(int argc, char* argv[])
         established_dispersers.emplace_back(
             I_species_rast.rows(), I_species_rast.cols());
     }
+    std::vector<Img> dispersers_rasts(num_runs, Img(S_species_rast, 0));
+    std::vector<Img> established_dispersers_rasts(num_runs, Img(S_species_rast, 0));
+
     std::vector<std::vector<Img>> soil_reservoirs(
         use_soils ? num_runs : 0,
         std::vector<Img>(
@@ -1401,6 +1419,10 @@ int main(int argc, char* argv[])
                         Network<Img::IndexType>::null_network(),
                         suitable_cells);
                     ++weather_step;
+                    if (opt.dispersers_output->answer)
+                        dispersers_rasts[run] += dispersers[run];
+                    if (opt.established_dispersers_output->answer)
+                        established_dispersers_rasts[run] += established_dispersers[run];
                 }
             }
 
@@ -1604,6 +1626,31 @@ int main(int argc, char* argv[])
         fprintf(fp, "%s", output.c_str());
         G_close_option_file(fp);
     }
-
+    if (opt.dispersers_output->answer) {
+        Img dispersers_rasts_sum(I_species_rast.rows(), I_species_rast.cols(), 0);
+        for (unsigned i = 0; i < num_runs; i++)
+            dispersers_rasts_sum += dispersers_rasts[i];
+        if (opt.dispersers_output->answer) {
+            // write final result
+            raster_to_grass(
+                        dispersers_rasts_sum,
+                        opt.dispersers_output->answer,
+                        "Sum of all dispersers",
+                        interval.end_date());
+        }
+    }
+    if (opt.established_dispersers_output->answer) {
+        Img established_dispersers_rasts_sum(I_species_rast.rows(), I_species_rast.cols(), 0);
+        for (unsigned i = 0; i < num_runs; i++)
+            established_dispersers_rasts_sum += established_dispersers_rasts[i];
+        if (opt.established_dispersers_output->answer) {
+            // write final result
+            raster_to_grass(
+                        established_dispersers_rasts_sum,
+                        opt.established_dispersers_output->answer,
+                        "Sum of all established dispersers",
+                        interval.end_date());
+        }
+    }
     return 0;
 }
diff --git a/testsuite/test_r_pops_spread.py b/testsuite/test_r_pops_spread.py
index 08fcd06..07598ff 100755
--- a/testsuite/test_r_pops_spread.py
+++ b/testsuite/test_r_pops_spread.py
@@ -191,7 +191,7 @@ def tearDown(self):
             "g.remove",
             flags="f",
             type="raster",
-            pattern="average*,single*,stddev*,probability*,dead*",
+            pattern="average*,single*,stddev*,probability*,dead*,*dispersers",
         )
 
     def test_outputs(self):
@@ -1406,6 +1406,54 @@ def test_survival_rate_0_percent(self):
         # Even with multiple runs, stddev should be still zero.
         self.assertRasterFitsUnivar(raster="stddev", reference=values, precision=0.001)
 
+    def test_outputs_dispersers(self):
+        """Check dead output of mortality"""
+        start = "2019-01-01"
+        end = "2022-12-31"
+        self.assertModule(
+            "r.pops.spread",
+            host="host",
+            total_plants="max_host",
+            infected="infection",
+            average="average",
+            average_series="average",
+            single_series="single",
+            stddev="stddev",
+            stddev_series="stddev",
+            probability="probability",
+            probability_series="probability",
+            start_date=start,
+            end_date=end,
+            seasonality=[1, 12],
+            step_unit="week",
+            step_num_units=1,
+            reproductive_rate=1,
+            natural_dispersal_kernel="exponential",
+            natural_distance=50,
+            natural_direction="W",
+            natural_direction_strength=3,
+            anthropogenic_dispersal_kernel="cauchy",
+            anthropogenic_distance=1000,
+            anthropogenic_direction_strength=0,
+            percent_natural_dispersal=0.95,
+            random_seed=1,
+            runs=5,
+            nprocs=5,
+            dispersers_output="dispersers",
+            established_dispersers_output="established_dispersers",
+        )
+        self.assertRasterExists("dispersers")
+        self.assertRasterExists("established_dispersers")
+
+        values = dict(null_cells=0, min=0, max=15530, mean=522.275)
+        self.assertRasterFitsUnivar(
+            raster="dispersers", reference=values, precision=0.001
+        )
+        values = dict(null_cells=0, min=0, max=129, mean=8.833)
+        self.assertRasterFitsUnivar(
+            raster="established_dispersers", reference=values, precision=0.001
+        )
+
     def test_with_and_without_anthropogenic_dispersal_multiple_seeds(self):
         """Check that multiple seeds keep anthropogenic dispersal separate