diff --git a/R/visualization.R b/R/visualization.R
index 55cfef162..c175b3acb 100644
--- a/R/visualization.R
+++ b/R/visualization.R
@@ -264,6 +264,8 @@ sort_splits <- function(model) {
 #'   vector of population names. If \code{NULL} (the default), the default plotting
 #'   algorithm will be used, ordering populations from the most ancestral to the
 #'   most recent using an in-order tree traversal.
+#' @param file Output file for a figure saved via \code{ggsave}
+#' @param ... Optional argument which will be passed to \code{ggsave}
 #'
 #' @return A ggplot2 object with the visualized slendr model
 #'
@@ -282,7 +284,7 @@ sort_splits <- function(model) {
 #'   labs geom_segment arrow
 #' @export
 plot_model <- function(model, sizes = TRUE, proportions = FALSE, gene_flow = TRUE, log = FALSE,
-                       order = NULL) {
+                       order = NULL, file = NULL, ...) {
   populations <- model$populations
 
   log10_ydelta <- 0.001
@@ -519,7 +521,10 @@ plot_model <- function(model, sizes = TRUE, proportions = FALSE, gene_flow = TRU
 
   p <- p + scale_y_continuous(trans = trans)
 
-  p
+  if (!is.null(file))
+    ggplot2::ggsave(file, plot = p, ...)
+  else
+    p
 }
 
 #' Animate the simulated population dynamics
diff --git a/docs/reference/plot_model.html b/docs/reference/plot_model.html
index 2f45ad445..ab06b00dc 100644
--- a/docs/reference/plot_model.html
+++ b/docs/reference/plot_model.html
@@ -92,7 +92,9 @@ <h2 id="ref-usage">Usage<a class="anchor" aria-label="anchor" href="#ref-usage">
 <span>  proportions <span class="op">=</span> <span class="cn">FALSE</span>,</span>
 <span>  gene_flow <span class="op">=</span> <span class="cn">TRUE</span>,</span>
 <span>  log <span class="op">=</span> <span class="cn">FALSE</span>,</span>
-<span>  order <span class="op">=</span> <span class="cn">NULL</span></span>
+<span>  order <span class="op">=</span> <span class="cn">NULL</span>,</span>
+<span>  file <span class="op">=</span> <span class="cn">NULL</span>,</span>
+<span>  <span class="va">...</span></span>
 <span><span class="op">)</span></span></code></pre></div>
     </div>
 
@@ -126,6 +128,14 @@ <h2 id="arguments">Arguments<a class="anchor" aria-label="anchor" href="#argumen
 algorithm will be used, ordering populations from the most ancestral to the
 most recent using an in-order tree traversal.</p></dd>
 
+
+<dt>file</dt>
+<dd><p>Output file for a figure saved via <code>ggsave</code></p></dd>
+
+
+<dt>...</dt>
+<dd><p>Optional argument which will be passed to <code>ggsave</code></p></dd>
+
 </dl></div>
     <div class="section level2">
     <h2 id="value">Value<a class="anchor" aria-label="anchor" href="#value"></a></h2>
diff --git a/man/plot_model.Rd b/man/plot_model.Rd
index dbdd62ea9..b9a47f0fb 100644
--- a/man/plot_model.Rd
+++ b/man/plot_model.Rd
@@ -10,7 +10,9 @@ plot_model(
   proportions = FALSE,
   gene_flow = TRUE,
   log = FALSE,
-  order = NULL
+  order = NULL,
+  file = NULL,
+  ...
 )
 }
 \arguments{
@@ -30,6 +32,10 @@ over very long time-scales.}
 vector of population names. If \code{NULL} (the default), the default plotting
 algorithm will be used, ordering populations from the most ancestral to the
 most recent using an in-order tree traversal.}
+
+\item{file}{Output file for a figure saved via \code{ggsave}}
+
+\item{...}{Optional argument which will be passed to \code{ggsave}}
 }
 \value{
 A ggplot2 object with the visualized slendr model