diff --git a/NEWS.md b/NEWS.md
index 7affeaf9e..dc97604c8 100644
--- a/NEWS.md
+++ b/NEWS.md
@@ -134,6 +134,14 @@ individual variability in growth to be modelled.
 - `plotSpectra()` now accepts `log_x`, `log_y`, and `log` arguments for
   controlling axis scaling, matching the mizer array `plot()` methods.
 
+- Size-based plots now accept `size_axis = "l"` to show length in cm on the
+  size axis instead of weight in grams, using the species' allometric
+  weight-length relationship.
+
+- Size-based plots with a `size_axis` argument now accept `llim`, the
+  length-axis equivalent of `wlim`, for filtering and limiting plots when
+  `size_axis = "l"`.
+
 - The `plot()` and `summary()` methods for `MizerParams`, `MizerSim`, and the
   mizer array classes are now registered as S3 methods rather than S4 methods,
   so `plot()` and `summary()` remain plain S3 generics when mizer is loaded,
diff --git a/R/ArraySpeciesBySize-class.R b/R/ArraySpeciesBySize-class.R
index 8a0d25ed7..560485979 100644
--- a/R/ArraySpeciesBySize-class.R
+++ b/R/ArraySpeciesBySize-class.R
@@ -159,6 +159,11 @@ print.summary.ArraySpeciesBySize <- function(x, ...) {
 #' @param wlim A numeric vector of length two providing lower and upper limits
 #'   for the weight (x) axis. Use `NA` to refer to the existing minimum or
 #'   maximum. Only applies to `ArraySpeciesBySize`.
+#' @param llim A numeric vector of length two providing lower and upper limits
+#'   for the length (x) axis when `size_axis = "l"`. Use `NA` to refer to the
+#'   existing minimum or maximum. Only applies to `ArraySpeciesBySize`.
+#' @param size_axis Whether to plot size as weight (`"w"`, default) or length
+#'   (`"l"`), using the allometric weight-length relationship.
 #' @param ylim A numeric vector of length two providing lower and upper limits
 #'   for the value (y) axis. Use `NA` to refer to the existing minimum or
 #'   maximum.
@@ -180,24 +185,36 @@ print.summary.ArraySpeciesBySize <- function(x, ...) {
 #' \donttest{
 #' plot(getEncounter(NS_params))
 #' plot(getFeedingLevel(NS_params), species = c("Cod", "Herring"))
+#' plot(getPredMort(NS_params), species = c("Cod", "Herring"),
+#'      size_axis = "l")
 #' }
 plot.ArraySpeciesBySize <- function(x, species = NULL,
                             all.sizes = FALSE, highlight = NULL,
                             return_data = FALSE, log_x = TRUE, log_y = FALSE,
                             log = NULL,
-                            wlim = c(NA, NA), ylim = c(NA, NA),
+                            wlim = c(NA, NA), llim = c(NA, NA),
+                            ylim = c(NA, NA),
+                            size_axis = c("w", "l"),
                             total = FALSE, background = TRUE,
                             y_ticks = 6, ...) {
+    size_axis <- plot_size_axis(size_axis)
     log_axes <- parsePlotLog(log, log_x = log_x, log_y = log_y)
     log_x <- log_axes$log_x
     log_y <- log_axes$log_y
 
+    assert_that(length(wlim) == 2,
+                length(llim) == 2,
+                length(ylim) == 2)
     value_name <- attr(x, "value_name") %||% "Rate"
     units_str <- attr(x, "units")
     params <- attr(x, "params")
     plot_dat <- prepare_ArraySpeciesBySize_plot_data(
         x, species = species, all.sizes = all.sizes, wlim = wlim,
         total = total, background = background)
+    plot_dat <- convert_plot_size_axis(plot_dat, params, size_axis)
+    if (identical(size_axis, "l")) {
+        plot_dat <- filter_plot_length_limits(plot_dat, llim)
+    }
 
     if (return_data) return(plot_dat)
 
@@ -206,11 +223,13 @@ plot.ArraySpeciesBySize <- function(x, species = NULL,
         y_label <- paste0(value_name, " [", units_str, "]")
     }
 
-    plotDataFrame(plot_dat, params, xlab = "Size [g]", ylab = y_label,
+    plotDataFrame(plot_dat, params, xlab = plot_size_xlab(size_axis),
+                  ylab = y_label,
                   xtrans = if (log_x) "log10" else "identity",
                   ytrans = if (log_y) "log10" else "identity",
-                  xlim = wlim, ylim = ylim, highlight = highlight,
-                  y_ticks = y_ticks, legend_var = "Legend")
+                  xlim = plot_size_xlim(wlim, size_axis, llim), ylim = ylim,
+                  highlight = highlight, y_ticks = y_ticks,
+                  legend_var = "Legend")
 }
 
 parsePlotLog <- function(log, log_x = FALSE, log_y = FALSE) {
@@ -257,14 +276,20 @@ plot2 <- function(x, y, ...) {
 plot2.ArraySpeciesBySize <- function(x, y, name1 = "First", name2 = "Second",
                                      species = NULL, all.sizes = FALSE,
                                      log_x = TRUE, log_y = FALSE, log = NULL,
-                                     wlim = c(NA, NA), ylim = c(NA, NA),
+                                     wlim = c(NA, NA), llim = c(NA, NA),
+                                     ylim = c(NA, NA),
+                                     size_axis = c("w", "l"),
                                      total = FALSE, background = TRUE,
                                      y_ticks = 6, ...) {
     check_plot2_compatible(x, y, "ArraySpeciesBySize")
     compare_array_metadata(x, y)
+    size_axis <- plot_size_axis(size_axis)
     log_axes <- parsePlotLog(log, log_x = log_x, log_y = log_y)
     log_x <- log_axes$log_x
     log_y <- log_axes$log_y
+    assert_that(length(wlim) == 2,
+                length(llim) == 2,
+                length(ylim) == 2)
 
     params <- attr(x, "params")
     y_label <- array_y_label(x, default = "Rate")
@@ -277,11 +302,13 @@ plot2.ArraySpeciesBySize <- function(x, y, name1 = "First", name2 = "Second",
 
     plotComparisonDataFrame(plot_dat1, plot_dat2, params,
                             name1 = name1, name2 = name2,
-                            xlab = "Size [g]", ylab = y_label,
+                            xlab = plot_size_xlab(size_axis), ylab = y_label,
                             xtrans = if (log_x) "log10" else "identity",
                             ytrans = if (log_y) "log10" else "identity",
-                            xlim = wlim, ylim = ylim,
-                            y_ticks = y_ticks, legend_var = "Legend")
+                            xlim = plot_size_xlim(wlim, size_axis, llim),
+                            ylim = ylim,
+                            y_ticks = y_ticks, legend_var = "Legend",
+                            size_axis = size_axis)
 }
 
 #' Plot the relative difference between two mizer array objects
@@ -315,11 +342,17 @@ plotRelative.ArraySpeciesBySize <- function(x, y, species = NULL,
                                             all.sizes = FALSE,
                                             log_x = TRUE,
                                             wlim = c(NA, NA),
+                                            llim = c(NA, NA),
                                             ylim = c(NA, NA),
+                                            size_axis = c("w", "l"),
                                             total = FALSE,
                                             background = TRUE, ...) {
     check_plot2_compatible(x, y, "ArraySpeciesBySize")
     compare_array_metadata(x, y)
+    size_axis <- plot_size_axis(size_axis)
+    assert_that(length(wlim) == 2,
+                length(llim) == 2,
+                length(ylim) == 2)
     params <- attr(x, "params")
     plot_dat1 <- prepare_ArraySpeciesBySize_plot_data(
         x, species = species, all.sizes = all.sizes, wlim = wlim,
@@ -329,9 +362,11 @@ plotRelative.ArraySpeciesBySize <- function(x, y, species = NULL,
         total = total, background = background)
 
     plotRelativeDataFrame(plot_dat1, plot_dat2, params,
-                          xlab = "Size [g]",
+                          xlab = plot_size_xlab(size_axis),
                           xtrans = if (log_x) "log10" else "identity",
-                          xlim = wlim, ylim = ylim, legend_var = "Legend")
+                          xlim = plot_size_xlim(wlim, size_axis, llim),
+                          ylim = ylim,
+                          legend_var = "Legend", size_axis = size_axis)
 }
 
 check_plot2_compatible <- function(x, y, class) {
@@ -405,6 +440,8 @@ addPlot <- function(plot, x, ...) {
 addPlot.ArraySpeciesBySize <- function(plot, x, species = NULL,
                                        all.sizes = FALSE,
                                        wlim = c(NA, NA),
+                                       llim = c(NA, NA),
+                                       size_axis = c("w", "l"),
                                        total = FALSE,
                                        background = TRUE,
                                        colour = NULL,
@@ -419,15 +456,24 @@ addPlot.ArraySpeciesBySize <- function(plot, x, species = NULL,
                 is.number(alpha),
                 alpha >= 0,
                 alpha <= 1)
+    size_axis <- plot_size_axis(size_axis)
+    assert_that(length(wlim) == 2,
+                length(llim) == 2)
 
     plot <- deep_copy(plot)
     plot_dat <- prepare_ArraySpeciesBySize_plot_data(
         x, species = species, all.sizes = all.sizes, wlim = wlim,
         total = total, background = background)
-    check_addPlot_compatible(plot, x_var = "w", y_var = "value",
+    params <- attr(x, "params")
+    plot_dat <- convert_plot_size_axis(plot_dat, params, size_axis)
+    if (identical(size_axis, "l")) {
+        plot_dat <- filter_plot_length_limits(plot_dat, llim)
+    }
+    x_var <- plot_size_x_var(size_axis)
+    check_addPlot_compatible(plot, x_var = x_var, y_var = "value",
                              units = attr(x, "units"))
 
-    mapping <- aes(x = .data[["w"]], y = .data[["value"]],
+    mapping <- aes(x = .data[[x_var]], y = .data[["value"]],
                    group = .data[["Species"]])
     if (is.null(colour)) {
         mapping$colour <- rlang::quo(.data[["Legend"]])
diff --git a/R/ArrayTimeBySpeciesBySize-class.R b/R/ArrayTimeBySpeciesBySize-class.R
index 84e743e84..9fffcd315 100644
--- a/R/ArrayTimeBySpeciesBySize-class.R
+++ b/R/ArrayTimeBySpeciesBySize-class.R
@@ -143,7 +143,9 @@ plot.ArrayTimeBySpeciesBySize <- function(x, species = NULL, time = NULL,
                                           all.sizes = FALSE, highlight = NULL,
                                           return_data = FALSE, log_x = TRUE,
                                           log_y = FALSE, log = NULL,
-                                          wlim = c(NA, NA), ylim = c(NA, NA),
+                                          wlim = c(NA, NA), llim = c(NA, NA),
+                                          ylim = c(NA, NA),
+                                          size_axis = c("w", "l"),
                                           total = FALSE, background = TRUE,
                                           y_ticks = 6, ...) {
     params <- attr(x, "params")
@@ -167,8 +169,10 @@ plot.ArrayTimeBySpeciesBySize <- function(x, species = NULL, time = NULL,
     plot.ArraySpeciesBySize(slice, species = species, all.sizes = all.sizes,
                             highlight = highlight, return_data = return_data,
                             log_x = log_x, log_y = log_y, log = log,
-                            wlim = wlim, ylim = ylim, total = total,
-                            background = background, y_ticks = y_ticks, ...)
+                            wlim = wlim, ylim = ylim, llim = llim,
+                            size_axis = size_axis,
+                            total = total, background = background,
+                            y_ticks = y_ticks, ...)
 }
 
 #' @rdname plot2
@@ -183,7 +187,9 @@ plot2.ArrayTimeBySpeciesBySize <- function(x, y, name1 = "First",
                                            log_x = TRUE, log_y = FALSE,
                                            log = NULL,
                                            wlim = c(NA, NA),
+                                           llim = c(NA, NA),
                                            ylim = c(NA, NA),
+                                           size_axis = c("w", "l"),
                                            total = FALSE,
                                            background = TRUE,
                                            y_ticks = 6, ...) {
@@ -194,7 +200,8 @@ plot2.ArrayTimeBySpeciesBySize <- function(x, y, name1 = "First",
     plot2.ArraySpeciesBySize(slice1, slice2, name1 = name1, name2 = name2,
                              species = species, all.sizes = all.sizes,
                              log_x = log_x, log_y = log_y, log = log,
-                             wlim = wlim, ylim = ylim, total = total,
+                             wlim = wlim, ylim = ylim, llim = llim,
+                             size_axis = size_axis, total = total,
                              background = background, y_ticks = y_ticks, ...)
 }
 
@@ -208,7 +215,9 @@ plotRelative.ArrayTimeBySpeciesBySize <- function(x, y, species = NULL,
                                                   all.sizes = FALSE,
                                                   log_x = TRUE,
                                                   wlim = c(NA, NA),
+                                                  llim = c(NA, NA),
                                                   ylim = c(NA, NA),
+                                                  size_axis = c("w", "l"),
                                                   total = FALSE,
                                                   background = TRUE, ...) {
     check_plot2_compatible(x, y, "ArrayTimeBySpeciesBySize")
@@ -217,7 +226,8 @@ plotRelative.ArrayTimeBySpeciesBySize <- function(x, y, species = NULL,
 
     plotRelative.ArraySpeciesBySize(slice1, slice2, species = species,
                                     all.sizes = all.sizes, log_x = log_x,
-                                    wlim = wlim, ylim = ylim,
+                                    wlim = wlim, ylim = ylim, llim = llim,
+                                    size_axis = size_axis,
                                     total = total, background = background,
                                     ...)
 }
@@ -267,8 +277,10 @@ animate.ArrayTimeBySpeciesBySize <- function(x, species = NULL,
                                              time_range = NULL,
                                              log_x = TRUE,
                                              wlim = c(NA, NA),
+                                             llim = c(NA, NA),
                                              ylim = c(NA, NA),
                                              log_y = TRUE,
+                                             size_axis = c("w", "l"),
                                              total = FALSE,
                                              background = TRUE,
                                              frame_duration = 500,
@@ -279,7 +291,8 @@ animate.ArrayTimeBySpeciesBySize <- function(x, species = NULL,
                 is.number(frame_duration), frame_duration >= 0,
                 is.number(transition_duration), transition_duration >= 0,
                 is.string(easing),
-                length(wlim) == 2, length(ylim) == 2)
+                length(wlim) == 2, length(llim) == 2, length(ylim) == 2)
+    size_axis <- plot_size_axis(size_axis)
 
     params <- attr(x, "params")
     value_name <- attr(x, "value_name") %||% "Value"
@@ -349,8 +362,12 @@ animate.ArrayTimeBySpeciesBySize <- function(x, species = NULL,
         y_label <- paste0(value_name, " [", units_str, "]")
     }
 
-    animate_plotly(df, params, log_x, log_y, y_label, wlim, ylim,
-                   frame_duration, transition_duration, easing)
+    animate_plotly(df, params, log_x, log_y, y_label, wlim, llim,
+                   ylim,
+                   size_axis = size_axis,
+                   frame_duration = frame_duration,
+                   transition_duration = transition_duration,
+                   easing = easing)
 }
 
 #' @export
diff --git a/R/animateSpectra.R b/R/animateSpectra.R
index 1b393fbee..ea7c92290 100644
--- a/R/animateSpectra.R
+++ b/R/animateSpectra.R
@@ -34,6 +34,8 @@
 #'   range. Default is the entire time range of `x`.
 #' @param log_x If `TRUE` (default), use a log10 x-axis for body size.
 #' @param log_y If `TRUE` (default), use a log10 y-axis.
+#' @param size_axis Whether to plot size as weight (`"w"`, default) or length
+#'   (`"l"`), using the allometric weight-length relationship.
 #' @param total A boolean value that determines whether the total over all
 #'   selected species is plotted as an additional trace called `"Total"`.
 #'   Default is `FALSE`.
@@ -43,6 +45,9 @@
 #' @param wlim A numeric vector of length two providing lower and upper limits
 #'   for the body-size (x) axis. Use `NA` to refer to the existing minimum or
 #'   maximum.
+#' @param llim A numeric vector of length two providing lower and upper limits
+#'   for the length (x) axis when `size_axis = "l"`. Use `NA` to refer to the
+#'   existing minimum or maximum.
 #' @param ylim A numeric vector of length two providing lower and upper limits
 #'   for the value (y) axis. Use `NA` to refer to the existing minimum or
 #'   maximum. Limits are applied as Plotly axis ranges, so points outside the
@@ -88,19 +93,23 @@ animate <- function(x, ...) UseMethod("animate")
 #' @export
 animate.MizerSim <- function(x, species = NULL, time_range = NULL,
                               log_x = TRUE, log_y = TRUE,
-                              wlim = c(NA, NA), ylim = c(NA, NA),
+                              wlim = c(NA, NA), llim = c(NA, NA),
+                              ylim = c(NA, NA),
+                              size_axis = c("w", "l"),
                               power = 1, total = FALSE, resource = TRUE,
                               background = TRUE,
                               frame_duration = 500,
                               transition_duration = frame_duration,
                               easing = "linear", ...) {
     sim <- x
+    size_axis <- plot_size_axis(size_axis)
     assert_that(is.flag(total), is.flag(resource), is.flag(background),
                 is.number(power),
                 is.number(frame_duration), frame_duration >= 0,
                 is.number(transition_duration), transition_duration >= 0,
                 is.string(easing),
                 length(wlim) == 2,
+                length(llim) == 2,
                 length(ylim) == 2)
 
     species <- valid_species_arg(sim, species)
@@ -155,8 +164,12 @@ animate.MizerSim <- function(x, species = NULL, time_range = NULL,
     }
     nf <- mutate(nf, value = value * w^power)
 
-    animate_plotly(nf, sim@params, log_x, log_y, y_label, wlim, ylim,
-                   frame_duration, transition_duration, easing)
+    animate_plotly(nf, sim@params, log_x, log_y, y_label, wlim, llim,
+                   ylim,
+                   size_axis = size_axis,
+                   frame_duration = frame_duration,
+                   transition_duration = transition_duration,
+                   easing = easing)
 }
 
 # Build a plotly animation from a prepared long-format data frame.
@@ -165,9 +178,14 @@ animate.MizerSim <- function(x, species = NULL, time_range = NULL,
 # by individual Species within each legend group — so background species always
 # appear together and share a single legend entry.
 animate_plotly <- function(df, params, log_x, log_y, y_label,
-                           wlim = c(NA, NA), ylim = c(NA, NA),
+                           wlim = c(NA, NA), llim = c(NA, NA),
+                           ylim = c(NA, NA),
+                           size_axis = "w",
                            frame_duration = 500, transition_duration = 500,
                            easing = "linear") {
+    size_axis <- plot_size_axis(size_axis)
+    df <- convert_plot_size_axis(df, params, size_axis)
+    x_var <- plot_size_x_var(size_axis)
     legend_name_order <- intersect(names(params@linecolour),
                                    unique(df$legend_name))
     sp_order <- unlist(lapply(legend_name_order, function(ln) {
@@ -185,7 +203,7 @@ animate_plotly <- function(df, params, log_x, log_y, y_label,
         p <- plotly::add_lines(
             p,
             data = df_sp,
-            x = ~w, y = ~value,
+            x = stats::as.formula(paste0("~", x_var)), y = ~value,
             frame = ~time,
             name = ln,
             legendgroup = ln,
@@ -194,7 +212,10 @@ animate_plotly <- function(df, params, log_x, log_y, y_label,
         )
     }
     p <- plotly::layout(p,
-                        xaxis = plotly_axis(df$w, wlim, log_x, "Size [g]"),
+                        xaxis = plotly_axis(
+                            df[[x_var]],
+                            plot_size_xlim(wlim, size_axis, llim),
+                            log_x, plot_size_xlab(size_axis)),
                         yaxis = plotly_axis(df$value, ylim, log_y, y_label),
                         legend = list(traceorder = "normal"))
     plotly::animation_opts(p, frame = frame_duration,
diff --git a/R/plots.R b/R/plots.R
index 7d76b0458..d65c3260f 100644
--- a/R/plots.R
+++ b/R/plots.R
@@ -131,7 +131,7 @@
 NULL
 
 # Hackiness to get past the 'no visible binding ... ' warning when running check
-utils::globalVariables(c("time", "value", "Species", "w", "gear", "Age",
+utils::globalVariables(c("time", "value", "Species", "w", "l", "gear", "Age",
                          "x", "y", "Year", "Yield", "Biomass", "Size",
                          "Proportion", "Prey", "Legend", "Type", "Gear",
                          "Predator", "weight", "a", "b", "age", "w_max",
@@ -289,7 +289,8 @@ plotComparisonDataFrame <- function(frame1, frame2, params,
                                     xlab = waiver(), ylab = waiver(),
                                     xtrans = "identity", ytrans = "identity",
                                     xlim = c(NA, NA), ylim = c(NA, NA),
-                                    y_ticks = 6, legend_var = "Legend") {
+                                    y_ticks = 6, legend_var = "Legend",
+                                    size_axis = NULL) {
     assert_that(is.data.frame(frame1),
                 is.data.frame(frame2),
                 is(params, "MizerParams"))
@@ -309,6 +310,12 @@ plotComparisonDataFrame <- function(frame1, frame2, params,
         stop("The `legend_var` argument must be the name of a variable ",
              "in the data frame.")
     }
+    if (!is.null(size_axis)) {
+        size_axis <- plot_size_axis(size_axis)
+        frame <- convert_plot_size_axis(frame, params, size_axis,
+                                        species_col = group_var)
+        x_var <- plot_size_x_var(size_axis)
+    }
 
     legend_levels <- intersect(names(params@linecolour), frame[[legend_var]])
     frame[[legend_var]] <- factor(frame[[legend_var]], levels = legend_levels)
@@ -343,7 +350,8 @@ plotRelativeDataFrame <- function(frame1, frame2, params,
                                   xtrans = "identity",
                                   xlim = c(NA, NA),
                                   ylim = c(NA, NA),
-                                  legend_var = "Legend") {
+                                  legend_var = "Legend",
+                                  size_axis = NULL) {
     assert_that(is.data.frame(frame1),
                 is.data.frame(frame2),
                 is(params, "MizerParams"))
@@ -365,6 +373,12 @@ plotRelativeDataFrame <- function(frame1, frame2, params,
     frame$rel_diff <- relative_difference(frame[[paste0(y_var, ".x")]],
                                           frame[[paste0(y_var, ".y")]])
     frame <- frame[is.finite(frame$rel_diff), ]
+    if (!is.null(size_axis)) {
+        size_axis <- plot_size_axis(size_axis)
+        frame <- convert_plot_size_axis(frame, params, size_axis,
+                                        species_col = group_var)
+        x_var <- plot_size_x_var(size_axis)
+    }
 
     legend_levels <- intersect(names(params@linecolour), frame[[legend_var]])
     frame[[legend_var]] <- factor(frame[[legend_var]], levels = legend_levels)
@@ -416,6 +430,76 @@ log_breaks <- function(n = 6) {
     }
 }
 
+plot_size_axis <- function(size_axis = "w") {
+    match.arg(size_axis, c("w", "l"))
+}
+
+plot_size_x_var <- function(size_axis) {
+    if (identical(plot_size_axis(size_axis), "l")) "l" else "w"
+}
+
+plot_size_xlab <- function(size_axis) {
+    if (identical(plot_size_axis(size_axis), "l")) "Length [cm]" else "Size [g]"
+}
+
+plot_size_xlim <- function(wlim, size_axis, llim = c(NA, NA)) {
+    if (identical(plot_size_axis(size_axis), "l")) llim else wlim
+}
+
+filter_plot_length_limits <- function(plot_dat, llim) {
+    if (!"l" %in% names(plot_dat)) {
+        return(plot_dat)
+    }
+    if (!is.na(llim[1])) plot_dat <- plot_dat[plot_dat$l >= llim[1], ]
+    if (!is.na(llim[2])) plot_dat <- plot_dat[plot_dat$l <= llim[2], ]
+    plot_dat
+}
+
+plot_size_tooltip <- function(size_axis, before = NULL, after = NULL) {
+    c(before, plot_size_x_var(size_axis), after)
+}
+
+convert_plot_size_axis <- function(plot_dat, params, size_axis,
+                                   species_col = "Species",
+                                   drop_w = TRUE) {
+    size_axis <- plot_size_axis(size_axis)
+    if (identical(size_axis, "w")) {
+        return(plot_dat)
+    }
+    if (!is(params, "MizerParams")) {
+        stop("Length-axis plots need a MizerParams object.")
+    }
+    if (!"w" %in% names(plot_dat)) {
+        stop("Length-axis plots need plotting data with a `w` column.")
+    }
+    if (!species_col %in% names(plot_dat)) {
+        stop("Length-axis plots need plotting data with a `", species_col,
+             "` column.")
+    }
+
+    species <- as.character(plot_dat[[species_col]])
+    species_idx <- match(species, as.character(params@species_params$species))
+    plot_dat <- plot_dat[!is.na(species_idx), , drop = FALSE]
+    species_idx <- species_idx[!is.na(species_idx)]
+    if (nrow(plot_dat) == 0) {
+        if (!drop_w) {
+            plot_dat$l <- numeric(0)
+            return(plot_dat[, c("l", "w",
+                                setdiff(names(plot_dat), c("l", "w"))),
+                            drop = FALSE])
+        }
+        return(plot_dat[, setdiff(names(plot_dat), "w"), drop = FALSE])
+    }
+    sp <- params@species_params[species_idx, , drop = FALSE]
+    plot_dat$l <- w2l(plot_dat$w, sp)
+    if (drop_w) {
+        return(plot_dat[, c("l", setdiff(names(plot_dat), c("l", "w"))),
+                        drop = FALSE])
+    }
+    plot_dat[, c("l", "w", setdiff(names(plot_dat), c("l", "w"))),
+             drop = FALSE]
+}
+
 
 #' Plot the biomass of species through time
 #'
@@ -748,6 +832,10 @@ plotlyYieldGear <- function(sim, species = NULL,
 #'   the fish grid) or `min(params@w)` when `resource = FALSE`; the upper
 #'   default is `max(params@w_full)`. Data is filtered to this range and the
 #'   axis limits are set accordingly.
+#' @param llim A numeric vector of length two providing lower and upper limits
+#'   for the length axis when `size_axis = "l"`. Use `NA` to auto-scale to the
+#'   data range. Data is filtered to this range and the axis limits are set
+#'   accordingly.
 #' @param ylim A numeric vector of length two providing lower and upper limits
 #'   for the y axis. Use NA to auto-scale to the data range. Values below 1e-20
 #'   are always filtered out from the data regardless of `ylim[1]`. Data above
@@ -775,13 +863,15 @@ plotlyYieldGear <- function(sim, species = NULL,
 #' @param log Character string specifying which axes should use log10 scales,
 #'   in the same form as the base [plot()] argument. For example, `"x"`,
 #'   `"y"`, `"xy"` or `""`. If supplied, this overrides `log_x` and `log_y`.
+#' @param size_axis Whether to plot size as weight (`"w"`, default) or length
+#'   (`"l"`), using the allometric weight-length relationship.
 #' @param return_data A boolean value that determines whether the formatted data
 #' used for the plot is returned instead of the plot itself. Default value is FALSE
 #' @param ... Other arguments (currently unused)
 #'
 #' @return A ggplot2 object, unless `return_data = TRUE`, in which case a data
-#'   frame with the four variables 'w', 'value', 'Species', 'Legend' is
-#'   returned.
+#'   frame with the four variables 'w' (or 'l' if `size_axis = "l"`), 'value',
+#'   'Species', 'Legend' is returned.
 #' @export
 #' @family plotting functions
 #' @seealso [plotting_functions]
@@ -794,6 +884,7 @@ plotlyYieldGear <- function(sim, species = NULL,
 #' plotSpectra(sim, time_range = 10:20)
 #' plotSpectra(sim, time_range = 10:20, power = 0)
 #' plotSpectra(sim, species = c("Cod", "Herring"), power = 1)
+#' plotSpectra(sim, species = c("Cod", "Herring"), size_axis = "l")
 #'
 #' # Returning the data frame
 #' fr <- plotSpectra(sim, return_data = TRUE)
@@ -810,16 +901,19 @@ plotSpectra <- function(object, ...) {
 plotSpectra.MizerSim <- function(object, species = NULL,
                         time_range,
                         geometric_mean = FALSE,
-                        wlim = c(NA, NA), ylim = c(NA, NA),
+                        wlim = c(NA, NA), llim = c(NA, NA),
+                        ylim = c(NA, NA),
                         power = 1, biomass = TRUE,
                         total = FALSE, resource = TRUE,
                         background = TRUE,
                         highlight = NULL, log_x = TRUE, log_y = TRUE,
-                        log = NULL, return_data = FALSE, ...) {
+                        log = NULL, size_axis = c("w", "l"),
+                        return_data = FALSE, ...) {
     # to deal with old-type biomass argument
     if (missing(power)) {
         power <- as.numeric(biomass)
     }
+    size_axis <- plot_size_axis(size_axis)
     log_axes <- parsePlotLog(log, log_x = log_x, log_y = log_y)
     log_x <- log_axes$log_x
     log_y <- log_axes$log_y
@@ -828,6 +922,7 @@ plotSpectra.MizerSim <- function(object, species = NULL,
                 is.flag(background),
                 is.number(power),
                 length(wlim) == 2,
+                length(llim) == 2,
                 length(ylim) == 2)
     species <- valid_species_arg(object, species)
     if (length(species) == 0 && !total && !resource) {
@@ -847,25 +942,30 @@ plotSpectra.MizerSim <- function(object, species = NULL,
     n_pp <- apply(object@n_pp[time_elements, , drop = FALSE], 2, mean_fn)
     plot_spectra(object@params, n = n, n_pp = n_pp,
                  species = species, wlim = wlim, ylim = ylim,
+                 llim = llim,
                  power = power, total = total, resource = resource,
                  background = background, highlight = highlight,
                  log_x = log_x, log_y = log_y,
+                 size_axis = size_axis,
                  return_data = return_data)
 }
 
 #' @rdname plotSpectra
 #' @export
 plotSpectra.MizerParams <- function(object, species = NULL,
-                        wlim = c(NA, NA), ylim = c(NA, NA),
+                        wlim = c(NA, NA), llim = c(NA, NA),
+                        ylim = c(NA, NA),
                         power = 1, biomass = TRUE,
                         total = FALSE, resource = TRUE,
                         background = TRUE,
                         highlight = NULL, log_x = TRUE, log_y = TRUE,
-                        log = NULL, return_data = FALSE, ...) {
+                        log = NULL, size_axis = c("w", "l"),
+                        return_data = FALSE, ...) {
     # to deal with old-type biomass argument
     if (missing(power)) {
         power <- as.numeric(biomass)
     }
+    size_axis <- plot_size_axis(size_axis)
     log_axes <- parsePlotLog(log, log_x = log_x, log_y = log_y)
     log_x <- log_axes$log_x
     log_y <- log_axes$log_y
@@ -874,6 +974,7 @@ plotSpectra.MizerParams <- function(object, species = NULL,
                 is.flag(background),
                 is.number(power),
                 length(wlim) == 2,
+                length(llim) == 2,
                 length(ylim) == 2)
     species <- valid_species_arg(object, species)
     if (length(species) == 0 && !total && !resource) {
@@ -882,18 +983,21 @@ plotSpectra.MizerParams <- function(object, species = NULL,
     plot_spectra(object, n = object@initial_n,
                  n_pp = object@initial_n_pp,
                  species = species, wlim = wlim, ylim = ylim,
+                 llim = llim,
                  power = power, total = total, resource = resource,
                  background = background, highlight = highlight,
                  log_x = log_x, log_y = log_y,
+                 size_axis = size_axis,
                  return_data = return_data)
 }
 
 
 plot_spectra <- function(params, n, n_pp,
-                         species, wlim, ylim, power,
+                         species, wlim, llim, ylim, power,
                          total, resource, background,
-                         highlight, log_x, log_y, return_data) {
+                         highlight, log_x, log_y, size_axis, return_data) {
     params <- validParams(params)
+    size_axis <- plot_size_axis(size_axis)
     if (is.na(wlim[1])) {
         wlim[1] <- if (resource) min(params@w) / 100 else min(params@w)
     }
@@ -964,13 +1068,18 @@ plot_spectra <- function(params, n, n_pp,
     }
     filter_min <- if (is.na(ylim[1])) 1e-20 else ylim[1]
     plot_dat <- plot_dat[plot_dat$value > filter_min, ]
+    plot_dat <- convert_plot_size_axis(plot_dat, params, size_axis)
+    if (identical(size_axis, "l")) {
+        plot_dat <- filter_plot_length_limits(plot_dat, llim)
+    }
 
     if (return_data) return(plot_dat)
 
-    plotDataFrame(plot_dat, params, xlab = "Size [g]", ylab = y_label,
+    plotDataFrame(plot_dat, params, xlab = plot_size_xlab(size_axis),
+                  ylab = y_label,
                   xtrans = if (log_x) "log10" else "identity",
                   ytrans = if (log_y) "log10" else "identity",
-                  xlim = wlim, ylim = ylim,
+                  xlim = plot_size_xlim(wlim, size_axis, llim), ylim = ylim,
                   highlight = highlight, legend_var = "Legend")
 }
 
@@ -994,8 +1103,8 @@ plot_spectra <- function(params, n, n_pp,
 #'   only `"x"` and `""` are supported. If supplied, this overrides `log_x`.
 #'
 #' @return A ggplot2 object, unless `return_data = TRUE`, in which case a data
-#'   frame with the four variables 'w', 'value', 'Species', 'Legend' is
-#'   returned.
+#'   frame with the four variables 'w' (or 'l' if `size_axis = "l"`), 'value',
+#'   'Species', 'Legend' is returned.
 #' @export
 #' @family plotting functions
 #' @seealso [plotSpectra()]
@@ -1013,13 +1122,16 @@ plotCDF <- function(object, ...) {
 plotCDF.MizerSim <- function(object, species = NULL,
                              time_range,
                              geometric_mean = FALSE,
-                             wlim = c(NA, NA), ylim = c(NA, NA),
+                             wlim = c(NA, NA), llim = c(NA, NA),
+                             ylim = c(NA, NA),
                              power = 1, biomass = TRUE,
                              total = FALSE, resource = FALSE,
                              background = TRUE,
                              highlight = NULL, normalise = TRUE,
                              log_x = TRUE, log = NULL,
+                             size_axis = c("w", "l"),
                              return_data = FALSE, ...) {
+    size_axis <- plot_size_axis(size_axis)
     if (missing(power)) {
         power <- as.numeric(biomass)
     }
@@ -1028,13 +1140,15 @@ plotCDF.MizerSim <- function(object, species = NULL,
                 is.flag(background), is.flag(normalise),
                 is.number(power),
                 length(wlim) == 2,
+                length(llim) == 2,
                 length(ylim) == 2)
 
     args <- list(object = object, species = species,
                  geometric_mean = geometric_mean,
-                 wlim = wlim, ylim = c(NA, NA),
+                 wlim = wlim, llim = llim, ylim = c(NA, NA),
                  power = power, total = total,
                  resource = resource, background = background,
+                 size_axis = "w",
                  return_data = TRUE)
     if (!missing(time_range)) {
         args$time_range <- time_range
@@ -1042,19 +1156,23 @@ plotCDF.MizerSim <- function(object, species = NULL,
     plot_dat <- do.call(plotSpectra, args)
     plot_cdf(plot_dat, object@params, power = power, normalise = normalise,
              log_x = log_x, wlim = wlim, ylim = ylim,
-             highlight = highlight, return_data = return_data)
+             llim = llim, highlight = highlight, size_axis = size_axis,
+             return_data = return_data)
 }
 
 #' @rdname plotCDF
 #' @export
 plotCDF.MizerParams <- function(object, species = NULL,
-                                wlim = c(NA, NA), ylim = c(NA, NA),
+                                wlim = c(NA, NA), llim = c(NA, NA),
+                                ylim = c(NA, NA),
                                 power = 1, biomass = TRUE,
                                 total = FALSE, resource = FALSE,
                                 background = TRUE,
                                 highlight = NULL, normalise = TRUE,
                                 log_x = TRUE, log = NULL,
+                                size_axis = c("w", "l"),
                                 return_data = FALSE, ...) {
+    size_axis <- plot_size_axis(size_axis)
     if (missing(power)) {
         power <- as.numeric(biomass)
     }
@@ -1063,29 +1181,43 @@ plotCDF.MizerParams <- function(object, species = NULL,
                 is.flag(background), is.flag(normalise),
                 is.number(power),
                 length(wlim) == 2,
+                length(llim) == 2,
                 length(ylim) == 2)
 
     plot_dat <- plotSpectra(object, species = species,
                             wlim = wlim, ylim = c(NA, NA),
+                            llim = llim,
                             power = power, total = total,
                             resource = resource, background = background,
+                            size_axis = "w",
                             return_data = TRUE)
     plot_cdf(plot_dat, object, power = power, normalise = normalise,
              log_x = log_x, wlim = wlim, ylim = ylim,
-             highlight = highlight, return_data = return_data)
+             llim = llim, highlight = highlight, size_axis = size_axis,
+             return_data = return_data)
 }
 
-plot_cdf <- function(plot_dat, params, power, normalise, log_x, wlim, ylim,
-                     highlight, return_data) {
+plot_cdf <- function(plot_dat, params, power, normalise, log_x, wlim, llim,
+                     ylim, highlight, size_axis, return_data) {
+    size_axis <- plot_size_axis(size_axis)
+    if (identical(size_axis, "l")) {
+        plot_dat_l <- convert_plot_size_axis(plot_dat, params, size_axis,
+                                             drop_w = FALSE)
+        plot_dat_l <- filter_plot_length_limits(plot_dat_l, llim)
+        plot_dat <- plot_dat_l[, setdiff(names(plot_dat_l), "l"),
+                               drop = FALSE]
+    }
     cdf_dat <- prepare_spectra_cdf_data(plot_dat, params,
                                         normalise = normalise)
+    cdf_dat <- convert_plot_size_axis(cdf_dat, params, size_axis)
     if (return_data) return(cdf_dat)
 
     plotDataFrame(cdf_dat, validParams(params),
-                  xlab = "Size [g]", ylab = cdf_y_label(power, normalise),
+                  xlab = plot_size_xlab(size_axis),
+                  ylab = cdf_y_label(power, normalise),
                   xtrans = if (log_x) "log10" else "identity",
                   ytrans = "identity",
-                  xlim = wlim, ylim = ylim,
+                  xlim = plot_size_xlim(wlim, size_axis, llim), ylim = ylim,
                   highlight = highlight, legend_var = "Legend")
 }
 
@@ -1170,28 +1302,33 @@ parsePlotCDFLog <- function(log, log_x) {
 #' }
 plotCDF2 <- function(object1, object2, name1 = "First", name2 = "Second",
                      power = 1, normalise = TRUE, log_x = TRUE, log = NULL,
-                     resource = FALSE, ...) {
+                     resource = FALSE, llim = c(NA, NA),
+                     size_axis = c("w", "l"), ...) {
+    size_axis <- plot_size_axis(size_axis)
     log_x <- parsePlotCDFLog(log, log_x)
-    assert_that(is.number(power), is.flag(normalise))
+    assert_that(is.number(power), is.flag(normalise), length(llim) == 2)
 
     args <- list(...)
     wlim <- args$wlim %||% c(NA, NA)
     ylim <- args$ylim %||% c(NA, NA)
 
     cf1 <- plotCDF(object1, power = power, normalise = normalise,
+                   size_axis = "w",
                    return_data = TRUE, ...)
     cf2 <- plotCDF(object2, power = power, normalise = normalise,
+                   size_axis = "w",
                    return_data = TRUE, ...)
     params <- if (is(object1, "MizerSim")) object1@params else object1
 
     plotComparisonDataFrame(cf1, cf2, validParams(params),
                             name1 = name1, name2 = name2,
-                            xlab = "Size [g]",
+                            xlab = plot_size_xlab(size_axis),
                             ylab = cdf_y_label(power, normalise),
                             xtrans = if (log_x) "log10" else "identity",
                             ytrans = "identity",
-                            xlim = wlim, ylim = ylim,
-                            legend_var = "Legend")
+                            xlim = plot_size_xlim(wlim, size_axis, llim),
+                            ylim = ylim, legend_var = "Legend",
+                            size_axis = size_axis)
 }
 
 #' Compare two size spectra in the same plot
@@ -1225,27 +1362,33 @@ plotCDF2 <- function(object1, object2, name1 = "First", name2 = "Second",
 #' }
 plotSpectra2 <- function(object1, object2, name1 = "First", name2 = "Second",
                          power = 1, log_x = TRUE, log_y = TRUE,
-                         log = NULL, ...) {
+                         log = NULL, llim = c(NA, NA),
+                         size_axis = c("w", "l"), ...) {
+    size_axis <- plot_size_axis(size_axis)
     log_axes <- parsePlotLog(log, log_x = log_x, log_y = log_y)
     log_x <- log_axes$log_x
     log_y <- log_axes$log_y
+    assert_that(length(llim) == 2)
 
     args <- list(...)
     wlim <- args$wlim %||% c(NA, NA)
     ylim <- args$ylim %||% c(NA, NA)
 
-    sf1 <- plotSpectra(object1, power = power, return_data = TRUE, ...)
-    sf2 <- plotSpectra(object2, power = power, return_data = TRUE, ...)
+    sf1 <- plotSpectra(object1, power = power, size_axis = "w",
+                       return_data = TRUE, ...)
+    sf2 <- plotSpectra(object2, power = power, size_axis = "w",
+                       return_data = TRUE, ...)
     params <- if (is(object1, "MizerSim")) object1@params else object1
 
     plotComparisonDataFrame(sf1, sf2, validParams(params),
                             name1 = name1, name2 = name2,
-                            xlab = "Size [g]",
+                            xlab = plot_size_xlab(size_axis),
                             ylab = spectra_y_label(power),
                             xtrans = if (log_x) "log10" else "identity",
                             ytrans = if (log_y) "log10" else "identity",
-                            xlim = wlim, ylim = ylim,
-                            legend_var = "Legend")
+                            xlim = plot_size_xlim(wlim, size_axis, llim),
+                            ylim = ylim, legend_var = "Legend",
+                            size_axis = size_axis)
 }
 
 spectra_y_label <- function(power) {
@@ -1261,11 +1404,15 @@ spectra_y_label <- function(power) {
 #' @export
 plotlySpectra2 <- function(object1, object2, name1 = "First",
                            name2 = "Second", power = 1,
-                           log_x = TRUE, log_y = TRUE, log = NULL, ...) {
+                           log_x = TRUE, log_y = TRUE, log = NULL,
+                           llim = c(NA, NA),
+                           size_axis = c("w", "l"), ...) {
+    size_axis <- plot_size_axis(size_axis)
     ggplotly(plotSpectra2(object1, object2, name1 = name1, name2 = name2,
                           power = power, log_x = log_x, log_y = log_y,
-                          log = log, ...),
-             tooltip = c("Species", "w", "value", "Model"))
+                          log = log, llim = llim, size_axis = size_axis, ...),
+             tooltip = plot_size_tooltip(size_axis, before = "Species",
+                                         after = c("value", "Model")))
 }
 
 #' Plot the relative difference between two spectra
@@ -1292,6 +1439,8 @@ plotlySpectra2 <- function(object1, object2, name1 = "First",
 #' @param ylim A numeric vector of length two providing lower and upper limits
 #'   for the relative difference (y) axis. Use `NA` to refer to the existing
 #'   minimum or maximum.
+#' @param size_axis Whether to plot size as weight (`"w"`, default) or length
+#'   (`"l"`), using the allometric weight-length relationship.
 #' @param ... Arguments passed to [plotSpectra()] for preparing the spectra
 #'   data, for example `species`, `time_range`, `wlim`, `resource`,
 #'   `background` or `total`.
@@ -1307,28 +1456,38 @@ plotlySpectra2 <- function(object1, object2, name1 = "First",
 #' plotSpectraRelative(sim1, sim2)
 #' }
 plotSpectraRelative <- function(object1, object2, log_x = TRUE,
-                                ylim = c(NA, NA), ...) {
+                                ylim = c(NA, NA),
+                                llim = c(NA, NA),
+                                size_axis = c("w", "l"), ...) {
+    size_axis <- plot_size_axis(size_axis)
+    assert_that(length(llim) == 2)
     args <- list(...)
     wlim <- args$wlim %||% c(NA, NA)
 
-    sf1 <- plotSpectra(object1, return_data = TRUE, ...)
-    sf2 <- plotSpectra(object2, return_data = TRUE, ...)
+    sf1 <- plotSpectra(object1, size_axis = "w", return_data = TRUE, ...)
+    sf2 <- plotSpectra(object2, size_axis = "w", return_data = TRUE, ...)
     params <- if (is(object1, "MizerSim")) object1@params else object1
 
     plotRelativeDataFrame(sf1, sf2, validParams(params),
-                          xlab = "Size [g]",
+                          xlab = plot_size_xlab(size_axis),
                           xtrans = if (log_x) "log10" else "identity",
-                          xlim = wlim, ylim = ylim,
-                          legend_var = "Legend")
+                          xlim = plot_size_xlim(wlim, size_axis, llim),
+                          ylim = ylim,
+                          legend_var = "Legend", size_axis = size_axis)
 }
 
 #' @rdname plotSpectraRelative
 #' @export
 plotlySpectraRelative <- function(object1, object2, log_x = TRUE,
-                                  ylim = c(NA, NA), ...) {
+                                  ylim = c(NA, NA),
+                                  llim = c(NA, NA),
+                                  size_axis = c("w", "l"), ...) {
+    size_axis <- plot_size_axis(size_axis)
     ggplotly(plotSpectraRelative(object1, object2, log_x = log_x,
-                                  ylim = ylim, ...),
-             tooltip = c("Legend", "w", "rel_diff"))
+                                  ylim = ylim, llim = llim,
+                                  size_axis = size_axis, ...),
+             tooltip = plot_size_tooltip(size_axis, before = "Legend",
+                                         after = "rel_diff"))
 }
 
 #' @rdname plotCDF
@@ -1336,19 +1495,23 @@ plotlySpectraRelative <- function(object1, object2, log_x = TRUE,
 #' @export
 plotlyCDF <- function(object, species = NULL,
                       time_range, geometric_mean = FALSE,
-                      wlim = c(NA, NA), ylim = c(NA, NA),
+                      wlim = c(NA, NA), llim = c(NA, NA),
+                      ylim = c(NA, NA),
                       power = 1, biomass = TRUE,
                       total = FALSE, resource = FALSE,
                       background = TRUE,
                       highlight = NULL, normalise = TRUE,
-                      log_x = TRUE, log = NULL, ...) {
+                      log_x = TRUE, log = NULL,
+                      size_axis = c("w", "l"), ...) {
+    size_axis <- plot_size_axis(size_axis)
     args <- list(object = object, species = species,
                  geometric_mean = geometric_mean,
                  wlim = wlim, ylim = ylim,
                  biomass = biomass, total = total,
                  resource = resource, background = background,
                  highlight = highlight, normalise = normalise,
-                 log_x = log_x, log = log, ...)
+                 log_x = log_x, log = log, llim = llim,
+                 size_axis = size_axis, ...)
     if (!missing(time_range)) {
         args$time_range <- time_range
     }
@@ -1356,7 +1519,8 @@ plotlyCDF <- function(object, species = NULL,
         args$power <- power
     }
     ggplotly(do.call("plotCDF", args),
-             tooltip = c("Species", "w", "value"))
+             tooltip = plot_size_tooltip(size_axis, before = "Species",
+                                         after = "value"))
 }
 
 #' @rdname plotCDF2
@@ -1364,30 +1528,39 @@ plotlyCDF <- function(object, species = NULL,
 #' @export
 plotlyCDF2 <- function(object1, object2, name1 = "First", name2 = "Second",
                        power = 1, normalise = TRUE,
-                       log_x = TRUE, log = NULL, resource = FALSE, ...) {
+                       log_x = TRUE, log = NULL, resource = FALSE,
+                       llim = c(NA, NA),
+                       size_axis = c("w", "l"), ...) {
+    size_axis <- plot_size_axis(size_axis)
     args <- list(object1 = object1, object2 = object2,
                  name1 = name1, name2 = name2,
-                 normalise = normalise, log_x = log_x, log = log, ...)
+                 normalise = normalise, log_x = log_x, log = log,
+                 llim = llim, size_axis = size_axis, ...)
     if (!missing(power)) {
         args$power <- power
     }
     ggplotly(do.call("plotCDF2", args),
-             tooltip = c("Species", "w", "value", "Model"))
+             tooltip = plot_size_tooltip(size_axis, before = "Species",
+                                         after = c("value", "Model")))
 }
 
 #' @rdname plotSpectra
 #' @export
 plotlySpectra <- function(object, species = NULL,
                         time_range, geometric_mean = FALSE,
-                        wlim = c(NA, NA), ylim = c(NA, NA),
+                        wlim = c(NA, NA), llim = c(NA, NA),
+                        ylim = c(NA, NA),
                         power = 1, biomass = TRUE,
                         total = FALSE, resource = TRUE,
                         background = TRUE,
                         highlight = NULL, log_x = TRUE, log_y = TRUE,
-                        log = NULL, ...) {
+                        log = NULL,
+                        size_axis = c("w", "l"), ...) {
+    size_axis <- plot_size_axis(size_axis)
     argg <- as.list(environment())
     ggplotly(do.call("plotSpectra", argg),
-             tooltip = c("Species", "w", "value"))
+             tooltip = plot_size_tooltip(size_axis, before = "Species",
+                                         after = "value"))
 }
 
 #' Plot the feeding level of species by size
@@ -1414,9 +1587,10 @@ plotlySpectra <- function(object, species = NULL,
 #'   plotted. Default FALSE.
 #'
 #' @return A ggplot2 object, unless `return_data = TRUE`, in which case a data
-#'   frame with the variables 'w', 'value' and 'Species' is returned. If also
-#'   `include_critical = TRUE` then the data frame contains a fourth variable
-#'   'Type' that distinguishes between 'actual' and 'critical' feeding level.
+#'   frame with the variables 'w' (or 'l' if `size_axis = "l"`), 'value' and
+#'   'Species' is returned. If also `include_critical = TRUE` then the data
+#'   frame contains a fourth variable 'Type' that distinguishes between
+#'   'actual' and 'critical' feeding level.
 #' @export
 #' @family plotting functions
 #' @seealso [plotting_functions], [getFeedingLevel()]
@@ -1443,9 +1617,14 @@ plotFeedingLevel <- function(object, ...) {
 plotFeedingLevel.MizerSim <- function(object, species = NULL,
             time_range, highlight = NULL,
             all.sizes = FALSE, include_critical = FALSE,
+            wlim = c(NA, NA), llim = c(NA, NA),
+            size_axis = c("w", "l"),
             return_data = FALSE, ...) {
+    size_axis <- plot_size_axis(size_axis)
     assert_that(is.flag(all.sizes),
                 is.flag(include_critical),
+                length(wlim) == 2,
+                length(llim) == 2,
                 is.flag(return_data))
     if (missing(time_range)) {
         time_range  <- max(as.numeric(dimnames(object@n)$time))
@@ -1459,6 +1638,8 @@ plotFeedingLevel.MizerSim <- function(object, species = NULL,
     plot_feeding_level(params, feed, species = species,
                        highlight = highlight, all.sizes = all.sizes,
                        include_critical = include_critical,
+                       wlim = wlim, llim = llim,
+                       size_axis = size_axis,
                        return_data = return_data)
 }
 
@@ -1467,21 +1648,29 @@ plotFeedingLevel.MizerSim <- function(object, species = NULL,
 plotFeedingLevel.MizerParams <- function(object, species = NULL,
             highlight = NULL,
             all.sizes = FALSE, include_critical = FALSE,
+            wlim = c(NA, NA), llim = c(NA, NA),
+            size_axis = c("w", "l"),
             return_data = FALSE, ...) {
+    size_axis <- plot_size_axis(size_axis)
     assert_that(is.flag(all.sizes),
                 is.flag(include_critical),
+                length(wlim) == 2,
+                length(llim) == 2,
                 is.flag(return_data))
     params <- validParams(object)
     feed <- getFeedingLevel(params, drop = FALSE)
     plot_feeding_level(params, feed, species = species,
                        highlight = highlight, all.sizes = all.sizes,
                        include_critical = include_critical,
+                       wlim = wlim, llim = llim,
+                       size_axis = size_axis,
                        return_data = return_data)
 }
 
 plot_feeding_level <- function(params, feed, species, highlight,
                                all.sizes, include_critical,
-                               return_data) {
+                               wlim, llim, size_axis, return_data) {
+    size_axis <- plot_size_axis(size_axis)
 
     # selector for desired species
     sel_sp <- valid_species_arg(params, species, return.logical = TRUE,
@@ -1513,8 +1702,15 @@ plot_feeding_level <- function(params, feed, species, highlight,
         }
         plot_dat <- plot_dat[complete.cases(plot_dat), ]
     }
+    if (!is.na(wlim[1])) plot_dat <- plot_dat[plot_dat$w >= wlim[1], ]
+    if (!is.na(wlim[2])) plot_dat <- plot_dat[plot_dat$w <= wlim[2], ]
+    plot_dat <- convert_plot_size_axis(plot_dat, params, size_axis)
+    if (identical(size_axis, "l")) {
+        plot_dat <- filter_plot_length_limits(plot_dat, llim)
+    }
 
     if (return_data) return(plot_dat)
+    x_var <- plot_size_x_var(size_axis)
 
     # Need to keep species in order for legend
     legend_levels <-
@@ -1539,9 +1735,10 @@ plot_feeding_level <- function(params, feed, species, highlight,
     } else {
         p <- ggplot(plot_dat, aes(group = Species))
     }
-    p + geom_line(aes(x = w, y = value,
+    p + geom_line(aes(x = .data[[x_var]], y = value,
                       colour = Legend, linetype = Legend, linewidth = Legend)) +
-        scale_x_continuous(name = "Size [g]", trans = "log10") +
+        scale_x_continuous(name = plot_size_xlab(size_axis), trans = "log10",
+                           limits = plot_size_xlim(wlim, size_axis, llim)) +
         scale_y_continuous(name = "Feeding Level") +
         coord_cartesian(ylim = c(0, 1)) +
         scale_colour_manual(values = params@linecolour[legend_levels]) +
@@ -1555,10 +1752,14 @@ plotlyFeedingLevel <- function(object,
                              species = NULL,
                              time_range,
                              highlight = NULL,
-                             include_critical = FALSE, ...) {
+                             include_critical = FALSE,
+                             wlim = c(NA, NA), llim = c(NA, NA),
+                             size_axis = c("w", "l"), ...) {
+    size_axis <- plot_size_axis(size_axis)
     argg <- as.list(environment())
     p <- ggplotly(do.call("plotFeedingLevel", argg),
-                  tooltip = c("Species", "w", "value"))
+                  tooltip = plot_size_tooltip(size_axis, before = "Species",
+                                              after = "value"))
 
     # When critical feeding level is included, ggplotly creates traces split by the
     # interaction of Species and Type, which produces a very long combined legend.
@@ -1638,7 +1839,8 @@ plotlyFeedingLevel <- function(object,
 #'   outside a species' size range. Default FALSE.
 #'
 #' @return  A ggplot2 object, unless `return_data = TRUE`, in which case a data
-#'   frame with the three variables 'w', 'value', 'Species' is returned.
+#'   frame with the three variables 'w' (or 'l' if `size_axis = "l"`), 'value',
+#'   'Species' is returned.
 #' @export
 #' @family plotting functions
 #' @seealso [plotting_functions],  [getPredMort()]
@@ -1663,8 +1865,12 @@ plotPredMort <- function(object, ...) {
 #' @export
 plotPredMort.MizerSim <- function(object, species = NULL,
                          time_range, all.sizes = FALSE,
-                         highlight = NULL, return_data = FALSE,
+                         highlight = NULL,
+                         wlim = c(NA, NA), llim = c(NA, NA),
+                         size_axis = c("w", "l"),
+                         return_data = FALSE,
                          ...) {
+    size_axis <- plot_size_axis(size_axis)
     if (missing(time_range)) {
         time_range <- max(as.numeric(dimnames(object@n)$time))
     }
@@ -1677,7 +1883,9 @@ plotPredMort.MizerSim <- function(object, species = NULL,
                                     units = "1/year",
                                     params = object@params)
     plot(pred_mort, species = species, all.sizes = all.sizes,
-         highlight = highlight, return_data = return_data,
+         highlight = highlight, wlim = wlim, llim = llim,
+         size_axis = size_axis,
+         return_data = return_data,
          ylim = c(0, NA))
 }
 
@@ -1685,10 +1893,15 @@ plotPredMort.MizerSim <- function(object, species = NULL,
 #' @export
 plotPredMort.MizerParams <- function(object, species = NULL,
                          all.sizes = FALSE,
-                         highlight = NULL, return_data = FALSE,
+                         highlight = NULL,
+                         wlim = c(NA, NA), llim = c(NA, NA),
+                         size_axis = c("w", "l"),
+                         return_data = FALSE,
                          ...) {
+    size_axis <- plot_size_axis(size_axis)
     plot(getPredMort(validParams(object)), species = species,
          all.sizes = all.sizes, highlight = highlight,
+         wlim = wlim, llim = llim, size_axis = size_axis,
          return_data = return_data, ylim = c(0, NA))
 }
 
@@ -1705,10 +1918,14 @@ plotM2 <- plotPredMort
 #' @export
 plotlyPredMort <- function(object, species = NULL,
                            time_range,
-                           highlight = NULL, ...) {
+                           highlight = NULL,
+                           wlim = c(NA, NA), llim = c(NA, NA),
+                           size_axis = c("w", "l"), ...) {
+    size_axis <- plot_size_axis(size_axis)
     argg <- as.list(environment())
     ggplotly(do.call("plotPredMort", argg),
-             tooltip = c("Species", "w", "value"))
+             tooltip = plot_size_tooltip(size_axis, before = "Species",
+                                         after = "value"))
 }
 
 #' Plot total fishing mortality of each species by size
@@ -1722,7 +1939,8 @@ plotlyPredMort <- function(object, species = NULL,
 #' @param all.sizes If TRUE, then fishing mortality is plotted also for sizes
 #'   outside a species' size range. Default FALSE.
 #' @return A ggplot2 object, unless `return_data = TRUE`, in which case a data
-#'   frame with the three variables 'w', 'value', 'Species' is returned.
+#'   frame with the three variables 'w' (or 'l' if `size_axis = "l"`), 'value',
+#'   'Species' is returned.
 #' @export
 #' @family plotting functions
 #' @seealso [plotting_functions], [getFMort()]
@@ -1747,8 +1965,12 @@ plotFMort <- function(object, ...) {
 #' @export
 plotFMort.MizerSim <- function(object, species = NULL,
                       time_range, all.sizes = FALSE,
-                      highlight = NULL, return_data = FALSE,
+                      highlight = NULL,
+                      wlim = c(NA, NA), llim = c(NA, NA),
+                      size_axis = c("w", "l"),
+                      return_data = FALSE,
                       ...) {
+    size_axis <- plot_size_axis(size_axis)
     if (missing(time_range)) {
         time_range <- max(as.numeric(dimnames(object@n)$time))
     }
@@ -1759,17 +1981,24 @@ plotFMort.MizerSim <- function(object, species = NULL,
     f <- ArraySpeciesBySize(f, value_name = "Fishing mortality",
                             units = "1/year", params = object@params)
     plot(f, species = species, all.sizes = all.sizes,
-         highlight = highlight, return_data = return_data)
+         highlight = highlight, wlim = wlim, llim = llim,
+         size_axis = size_axis,
+         return_data = return_data)
 }
 
 #' @rdname plotFMort
 #' @export
 plotFMort.MizerParams <- function(object, species = NULL,
                       all.sizes = FALSE,
-                      highlight = NULL, return_data = FALSE,
+                      highlight = NULL,
+                      wlim = c(NA, NA), llim = c(NA, NA),
+                      size_axis = c("w", "l"),
+                      return_data = FALSE,
                       ...) {
+    size_axis <- plot_size_axis(size_axis)
     plot(getFMort(validParams(object)), species = species,
          all.sizes = all.sizes, highlight = highlight,
+         wlim = wlim, llim = llim, size_axis = size_axis,
          return_data = return_data)
 }
 
@@ -1777,10 +2006,14 @@ plotFMort.MizerParams <- function(object, species = NULL,
 #' @export
 plotlyFMort <- function(object, species = NULL,
                         time_range,
-                        highlight = NULL, ...) {
+                        highlight = NULL,
+                        wlim = c(NA, NA), llim = c(NA, NA),
+                        size_axis = c("w", "l"), ...) {
+    size_axis <- plot_size_axis(size_axis)
     argg <- as.list(environment())
     ggplotly(do.call("plotFMort", argg),
-             tooltip = c("Species", "w", "value"))
+             tooltip = plot_size_tooltip(size_axis, before = "Species",
+                                         after = "value"))
 }
 
 
@@ -2052,8 +2285,8 @@ plotlyGrowthCurves <- function(object, species = NULL,
 #' @inheritParams plotSpectra
 #'
 #' @return A ggplot2 object, unless `return_data = TRUE`, in which case a data
-#'   frame with the four variables 'Predator', 'w', 'Proportion', 'Prey' is
-#'   returned.
+#'   frame with the four variables 'Predator', 'w' (or 'l' if
+#'   `size_axis = "l"`), 'Proportion', 'Prey' is returned.
 #' @export
 #' @seealso [getDiet()]
 #' @family plotting functions
@@ -2075,21 +2308,35 @@ plotDiet <- function(object, ...) {
 #' @rdname plotDiet
 #' @export
 plotDiet.MizerSim <- function(object, species = NULL,
+                              wlim = c(NA, NA), llim = c(NA, NA),
+                              size_axis = c("w", "l"),
                               return_data = FALSE, ...) {
-    plotDiet(object@params, species = species, return_data = return_data, ...)
+    size_axis <- plot_size_axis(size_axis)
+    plotDiet(object@params, species = species, wlim = wlim, llim = llim,
+             size_axis = size_axis,
+             return_data = return_data, ...)
 }
 
 #' @rdname plotDiet
 #' @export
-plotDiet.MizerParams <- function(object, species = NULL, return_data = FALSE, ...) {
-    assert_that(is.flag(return_data))
+plotDiet.MizerParams <- function(object, species = NULL,
+                                 wlim = c(NA, NA), llim = c(NA, NA),
+                                 size_axis = c("w", "l"),
+                                 return_data = FALSE, ...) {
+    size_axis <- plot_size_axis(size_axis)
+    assert_that(is.flag(return_data),
+                length(wlim) == 2,
+                length(llim) == 2)
     params <- validParams(object)
     diet <- getDiet(params)
     plot_diet(params, n = params@initial_n, diet = diet, species = species,
+              wlim = wlim, llim = llim, size_axis = size_axis,
               return_data = return_data)
 }
 
-plot_diet <- function(params, n, diet, species, return_data) {
+plot_diet <- function(params, n, diet, species, wlim, llim,
+                      size_axis, return_data) {
+    size_axis <- plot_size_axis(size_axis)
     species <- valid_species_arg(params, species, return.logical = TRUE)
     diet <- diet[species, , , drop = FALSE]
     names(dimnames(diet)) <- c("Predator", "w", "Prey")
@@ -2099,6 +2346,8 @@ plot_diet <- function(params, n, diet, species, return_data) {
     plot_dat$Prey <- factor(plot_dat$Prey, levels = rev(prey))
 
     plot_dat <- plot_dat[plot_dat$Proportion > 0.001, ]
+    if (!is.na(wlim[1])) plot_dat <- plot_dat[plot_dat$w >= wlim[1], ]
+    if (!is.na(wlim[2])) plot_dat <- plot_dat[plot_dat$w <= wlim[2], ]
 
     # Restrict plot to relevant size ranges where abundance is meaningful
     # For each predator species, find the maximum size where density is meaningful
@@ -2120,14 +2369,20 @@ plot_diet <- function(params, n, diet, species, return_data) {
         }
     }
 
+    plot_dat <- convert_plot_size_axis(plot_dat, params, size_axis,
+                                       species_col = "Predator")
+    if (identical(size_axis, "l")) {
+        plot_dat <- filter_plot_length_limits(plot_dat, llim)
+    }
     if (return_data) return(plot_dat)
+    x_var <- plot_size_x_var(size_axis)
 
     legend_levels <-
         intersect(names(params@linecolour), plot_dat$Prey)
     p <- ggplot(plot_dat) +
-        geom_area(aes(x = w, y = Proportion, fill = Prey)) +
-        scale_x_log10() +
-        labs(x = "Size [g]", y = "Proportion") +
+        geom_area(aes(x = .data[[x_var]], y = Proportion, fill = Prey)) +
+        scale_x_log10(limits = plot_size_xlim(wlim, size_axis, llim)) +
+        labs(x = plot_size_xlab(size_axis), y = "Proportion") +
         scale_fill_manual(values = params@linecolour[legend_levels],
                           limits = legend_levels)
     if (sum(species) > 1) {
@@ -2139,9 +2394,14 @@ plot_diet <- function(params, n, diet, species, return_data) {
 #' @rdname plotDiet
 #' @return `plotlyDiet()` returns a plotly object.
 #' @export
-plotlyDiet <- function(object, species = NULL, ...) {
-    ggplotly(plotDiet(object, species = species, ...),
-             tooltip = c("Predator", "w", "Proportion", "Prey"))
+plotlyDiet <- function(object, species = NULL,
+                       wlim = c(NA, NA), llim = c(NA, NA),
+                       size_axis = c("w", "l"), ...) {
+    size_axis <- plot_size_axis(size_axis)
+    ggplotly(plotDiet(object, species = species, wlim = wlim, llim = llim,
+                      size_axis = size_axis, ...),
+             tooltip = plot_size_tooltip(size_axis, before = "Predator",
+                                         after = c("Proportion", "Prey")))
 }
 
 
diff --git a/man/addPlot.Rd b/man/addPlot.Rd
index 602a34132..666b7a4e8 100644
--- a/man/addPlot.Rd
+++ b/man/addPlot.Rd
@@ -15,6 +15,8 @@ addPlot(plot, x, ...)
   species = NULL,
   all.sizes = FALSE,
   wlim = c(NA, NA),
+  llim = c(NA, NA),
+  size_axis = c("w", "l"),
   total = FALSE,
   background = TRUE,
   colour = NULL,
@@ -57,6 +59,13 @@ all species.}
 for the weight (x) axis. Use \code{NA} to refer to the existing minimum or
 maximum. Only applies to \code{ArraySpeciesBySize}.}
 
+\item{llim}{A numeric vector of length two providing lower and upper limits
+for the length (x) axis when \code{size_axis = "l"}. Use \code{NA} to refer to the
+existing minimum or maximum. Only applies to \code{ArraySpeciesBySize}.}
+
+\item{size_axis}{Whether to plot size as weight (\code{"w"}, default) or length
+(\code{"l"}), using the allometric weight-length relationship.}
+
 \item{total}{A boolean value that determines whether the total over all
 selected species is plotted as well. Default is \code{FALSE}.}
 
diff --git a/man/animate.Rd b/man/animate.Rd
index 7e547049e..bb5b8fe8c 100644
--- a/man/animate.Rd
+++ b/man/animate.Rd
@@ -14,8 +14,10 @@
   time_range = NULL,
   log_x = TRUE,
   wlim = c(NA, NA),
+  llim = c(NA, NA),
   ylim = c(NA, NA),
   log_y = TRUE,
+  size_axis = c("w", "l"),
   total = FALSE,
   background = TRUE,
   frame_duration = 500,
@@ -33,7 +35,9 @@ animate(x, ...)
   log_x = TRUE,
   log_y = TRUE,
   wlim = c(NA, NA),
+  llim = c(NA, NA),
   ylim = c(NA, NA),
+  size_axis = c("w", "l"),
   power = 1,
   total = FALSE,
   resource = TRUE,
@@ -62,6 +66,10 @@ range. Default is the entire time range of \code{x}.}
 for the body-size (x) axis. Use \code{NA} to refer to the existing minimum or
 maximum.}
 
+\item{llim}{A numeric vector of length two providing lower and upper limits
+for the length (x) axis when \code{size_axis = "l"}. Use \code{NA} to refer to the
+existing minimum or maximum.}
+
 \item{ylim}{A numeric vector of length two providing lower and upper limits
 for the value (y) axis. Use \code{NA} to refer to the existing minimum or
 maximum. Limits are applied as Plotly axis ranges, so points outside the
@@ -70,6 +78,9 @@ frames.}
 
 \item{log_y}{If \code{TRUE} (default), use a log10 y-axis.}
 
+\item{size_axis}{Whether to plot size as weight (\code{"w"}, default) or length
+(\code{"l"}), using the allometric weight-length relationship.}
+
 \item{total}{A boolean value that determines whether the total over all
 selected species is plotted as an additional trace called \code{"Total"}.
 Default is \code{FALSE}.}
diff --git a/man/plot.Rd b/man/plot.Rd
index 834ebffd6..13f10773a 100644
--- a/man/plot.Rd
+++ b/man/plot.Rd
@@ -21,7 +21,9 @@
   log_y = FALSE,
   log = NULL,
   wlim = c(NA, NA),
+  llim = c(NA, NA),
   ylim = c(NA, NA),
+  size_axis = c("w", "l"),
   total = FALSE,
   background = TRUE,
   y_ticks = 6,
@@ -80,7 +82,9 @@
   log_y = FALSE,
   log = NULL,
   wlim = c(NA, NA),
+  llim = c(NA, NA),
   ylim = c(NA, NA),
+  size_axis = c("w", "l"),
   total = FALSE,
   background = TRUE,
   y_ticks = 6,
@@ -126,10 +130,17 @@ in the same form as the base \code{\link[=plot]{plot()}} argument. For example,
 for the weight (x) axis. Use \code{NA} to refer to the existing minimum or
 maximum. Only applies to \code{ArraySpeciesBySize}.}
 
+\item{llim}{A numeric vector of length two providing lower and upper limits
+for the length (x) axis when \code{size_axis = "l"}. Use \code{NA} to refer to the
+existing minimum or maximum. Only applies to \code{ArraySpeciesBySize}.}
+
 \item{ylim}{A numeric vector of length two providing lower and upper limits
 for the value (y) axis. Use \code{NA} to refer to the existing minimum or
 maximum.}
 
+\item{size_axis}{Whether to plot size as weight (\code{"w"}, default) or length
+(\code{"l"}), using the allometric weight-length relationship.}
+
 \item{total}{A boolean value that determines whether the total over all
 selected species is plotted as well. Default is \code{FALSE}.}
 
@@ -180,6 +191,8 @@ coloured appropriately.
 \donttest{
 plot(getEncounter(NS_params))
 plot(getFeedingLevel(NS_params), species = c("Cod", "Herring"))
+plot(getPredMort(NS_params), species = c("Cod", "Herring"),
+     size_axis = "l")
 }
 \donttest{
 ggplotly(getEncounter(NS_params))
diff --git a/man/plot2.Rd b/man/plot2.Rd
index bd99b0525..7b3e9d1d6 100644
--- a/man/plot2.Rd
+++ b/man/plot2.Rd
@@ -21,7 +21,9 @@ plot2(x, y, ...)
   log_y = FALSE,
   log = NULL,
   wlim = c(NA, NA),
+  llim = c(NA, NA),
   ylim = c(NA, NA),
+  size_axis = c("w", "l"),
   total = FALSE,
   background = TRUE,
   y_ticks = 6,
@@ -58,7 +60,9 @@ plot2(x, y, ...)
   log_y = FALSE,
   log = NULL,
   wlim = c(NA, NA),
+  llim = c(NA, NA),
   ylim = c(NA, NA),
+  size_axis = c("w", "l"),
   total = FALSE,
   background = TRUE,
   y_ticks = 6,
@@ -92,10 +96,17 @@ in the same form as the base \code{\link[=plot]{plot()}} argument. For example,
 for the weight (x) axis. Use \code{NA} to refer to the existing minimum or
 maximum. Only applies to \code{ArraySpeciesBySize}.}
 
+\item{llim}{A numeric vector of length two providing lower and upper limits
+for the length (x) axis when \code{size_axis = "l"}. Use \code{NA} to refer to the
+existing minimum or maximum. Only applies to \code{ArraySpeciesBySize}.}
+
 \item{ylim}{A numeric vector of length two providing lower and upper limits
 for the value (y) axis. Use \code{NA} to refer to the existing minimum or
 maximum.}
 
+\item{size_axis}{Whether to plot size as weight (\code{"w"}, default) or length
+(\code{"l"}), using the allometric weight-length relationship.}
+
 \item{total}{A boolean value that determines whether the total over all
 selected species is plotted as well. Default is \code{FALSE}.}
 
diff --git a/man/plotCDF.Rd b/man/plotCDF.Rd
index 043b2ae46..6b8dbec50 100644
--- a/man/plotCDF.Rd
+++ b/man/plotCDF.Rd
@@ -15,6 +15,7 @@ plotCDF(object, ...)
   time_range,
   geometric_mean = FALSE,
   wlim = c(NA, NA),
+  llim = c(NA, NA),
   ylim = c(NA, NA),
   power = 1,
   biomass = TRUE,
@@ -25,6 +26,7 @@ plotCDF(object, ...)
   normalise = TRUE,
   log_x = TRUE,
   log = NULL,
+  size_axis = c("w", "l"),
   return_data = FALSE,
   ...
 )
@@ -33,6 +35,7 @@ plotCDF(object, ...)
   object,
   species = NULL,
   wlim = c(NA, NA),
+  llim = c(NA, NA),
   ylim = c(NA, NA),
   power = 1,
   biomass = TRUE,
@@ -43,6 +46,7 @@ plotCDF(object, ...)
   normalise = TRUE,
   log_x = TRUE,
   log = NULL,
+  size_axis = c("w", "l"),
   return_data = FALSE,
   ...
 )
@@ -53,6 +57,7 @@ plotlyCDF(
   time_range,
   geometric_mean = FALSE,
   wlim = c(NA, NA),
+  llim = c(NA, NA),
   ylim = c(NA, NA),
   power = 1,
   biomass = TRUE,
@@ -63,6 +68,7 @@ plotlyCDF(
   normalise = TRUE,
   log_x = TRUE,
   log = NULL,
+  size_axis = c("w", "l"),
   ...
 )
 }
@@ -93,6 +99,11 @@ the fish grid) or \code{min(params@w)} when \code{resource = FALSE}; the upper
 default is \code{max(params@w_full)}. Data is filtered to this range and the
 axis limits are set accordingly.}
 
+\item{llim}{A numeric vector of length two providing lower and upper limits
+for the length axis when \code{size_axis = "l"}. Use \code{NA} to auto-scale to the
+data range. Data is filtered to this range and the axis limits are set
+accordingly.}
+
 \item{ylim}{A numeric vector of length two providing lower and upper limits
 for the y axis. Use NA to auto-scale to the data range. Values below 1e-20
 are always filtered out from the data regardless of \code{ylim[1]}. Data above
@@ -132,13 +143,16 @@ integral.}
 scale, in the same form as the base \code{\link[=plot]{plot()}} argument. For \code{plotCDF()},
 only \code{"x"} and \code{""} are supported. If supplied, this overrides \code{log_x}.}
 
+\item{size_axis}{Whether to plot size as weight (\code{"w"}, default) or length
+(\code{"l"}), using the allometric weight-length relationship.}
+
 \item{return_data}{A boolean value that determines whether the formatted data
 used for the plot is returned instead of the plot itself. Default value is FALSE}
 }
 \value{
 A ggplot2 object, unless \code{return_data = TRUE}, in which case a data
-frame with the four variables 'w', 'value', 'Species', 'Legend' is
-returned.
+frame with the four variables 'w' (or 'l' if \code{size_axis = "l"}), 'value',
+'Species', 'Legend' is returned.
 
 \code{plotlyCDF()} returns a plotly object.
 }
diff --git a/man/plotCDF2.Rd b/man/plotCDF2.Rd
index 5224e814a..8d8989e45 100644
--- a/man/plotCDF2.Rd
+++ b/man/plotCDF2.Rd
@@ -15,6 +15,8 @@ plotCDF2(
   log_x = TRUE,
   log = NULL,
   resource = FALSE,
+  llim = c(NA, NA),
+  size_axis = c("w", "l"),
   ...
 )
 
@@ -28,6 +30,8 @@ plotlyCDF2(
   log_x = TRUE,
   log = NULL,
   resource = FALSE,
+  llim = c(NA, NA),
+  size_axis = c("w", "l"),
   ...
 )
 }
@@ -56,6 +60,14 @@ only \code{"x"} and \code{""} are supported. If supplied, this overrides \code{l
 \item{resource}{A boolean value that determines whether resource is included.
 Default is FALSE.}
 
+\item{llim}{A numeric vector of length two providing lower and upper limits
+for the length axis when \code{size_axis = "l"}. Use \code{NA} to auto-scale to the
+data range. Data is filtered to this range and the axis limits are set
+accordingly.}
+
+\item{size_axis}{Whether to plot size as weight (\code{"w"}, default) or length
+(\code{"l"}), using the allometric weight-length relationship.}
+
 \item{...}{Arguments passed to \code{\link[=plotCDF]{plotCDF()}} for preparing the cumulative
 distribution data, for example \code{species}, \code{time_range}, \code{wlim},
 \code{resource}, \code{background} or \code{total}.}
diff --git a/man/plotDiet.Rd b/man/plotDiet.Rd
index c6aeb3d55..3cf9a87bc 100644
--- a/man/plotDiet.Rd
+++ b/man/plotDiet.Rd
@@ -9,11 +9,34 @@
 \usage{
 plotDiet(object, ...)
 
-\method{plotDiet}{MizerSim}(object, species = NULL, return_data = FALSE, ...)
+\method{plotDiet}{MizerSim}(
+  object,
+  species = NULL,
+  wlim = c(NA, NA),
+  llim = c(NA, NA),
+  size_axis = c("w", "l"),
+  return_data = FALSE,
+  ...
+)
 
-\method{plotDiet}{MizerParams}(object, species = NULL, return_data = FALSE, ...)
+\method{plotDiet}{MizerParams}(
+  object,
+  species = NULL,
+  wlim = c(NA, NA),
+  llim = c(NA, NA),
+  size_axis = c("w", "l"),
+  return_data = FALSE,
+  ...
+)
 
-plotlyDiet(object, species = NULL, ...)
+plotlyDiet(
+  object,
+  species = NULL,
+  wlim = c(NA, NA),
+  llim = c(NA, NA),
+  size_axis = c("w", "l"),
+  ...
+)
 }
 \arguments{
 \item{object}{An object of class \linkS4class{MizerSim} or
@@ -26,13 +49,28 @@ species are selected. A vector of species names, or a
 numeric vector with the species indices, or a logical vector indicating for
 each species whether it is to be selected (TRUE) or not.}
 
+\item{wlim}{A numeric vector of length two providing lower and upper limits
+for the w axis. Use NA for the default: the lower default is
+\code{min(params@w) / 100} when \code{resource = TRUE} (to show some resource below
+the fish grid) or \code{min(params@w)} when \code{resource = FALSE}; the upper
+default is \code{max(params@w_full)}. Data is filtered to this range and the
+axis limits are set accordingly.}
+
+\item{llim}{A numeric vector of length two providing lower and upper limits
+for the length axis when \code{size_axis = "l"}. Use \code{NA} to auto-scale to the
+data range. Data is filtered to this range and the axis limits are set
+accordingly.}
+
+\item{size_axis}{Whether to plot size as weight (\code{"w"}, default) or length
+(\code{"l"}), using the allometric weight-length relationship.}
+
 \item{return_data}{A boolean value that determines whether the formatted data
 used for the plot is returned instead of the plot itself. Default value is FALSE}
 }
 \value{
 A ggplot2 object, unless \code{return_data = TRUE}, in which case a data
-frame with the four variables 'Predator', 'w', 'Proportion', 'Prey' is
-returned.
+frame with the four variables 'Predator', 'w' (or 'l' if
+\code{size_axis = "l"}), 'Proportion', 'Prey' is returned.
 
 \code{plotlyDiet()} returns a plotly object.
 }
diff --git a/man/plotFMort.Rd b/man/plotFMort.Rd
index af1c0f968..0a1fd1de0 100644
--- a/man/plotFMort.Rd
+++ b/man/plotFMort.Rd
@@ -15,6 +15,9 @@ plotFMort(object, ...)
   time_range,
   all.sizes = FALSE,
   highlight = NULL,
+  wlim = c(NA, NA),
+  llim = c(NA, NA),
+  size_axis = c("w", "l"),
   return_data = FALSE,
   ...
 )
@@ -24,11 +27,23 @@ plotFMort(object, ...)
   species = NULL,
   all.sizes = FALSE,
   highlight = NULL,
+  wlim = c(NA, NA),
+  llim = c(NA, NA),
+  size_axis = c("w", "l"),
   return_data = FALSE,
   ...
 )
 
-plotlyFMort(object, species = NULL, time_range, highlight = NULL, ...)
+plotlyFMort(
+  object,
+  species = NULL,
+  time_range,
+  highlight = NULL,
+  wlim = c(NA, NA),
+  llim = c(NA, NA),
+  size_axis = c("w", "l"),
+  ...
+)
 }
 \arguments{
 \item{object}{An object of class \linkS4class{MizerSim} or
@@ -51,12 +66,28 @@ outside a species' size range. Default FALSE.}
 
 \item{highlight}{Name or vector of names of the species to be highlighted.}
 
+\item{wlim}{A numeric vector of length two providing lower and upper limits
+for the w axis. Use NA for the default: the lower default is
+\code{min(params@w) / 100} when \code{resource = TRUE} (to show some resource below
+the fish grid) or \code{min(params@w)} when \code{resource = FALSE}; the upper
+default is \code{max(params@w_full)}. Data is filtered to this range and the
+axis limits are set accordingly.}
+
+\item{llim}{A numeric vector of length two providing lower and upper limits
+for the length axis when \code{size_axis = "l"}. Use \code{NA} to auto-scale to the
+data range. Data is filtered to this range and the axis limits are set
+accordingly.}
+
+\item{size_axis}{Whether to plot size as weight (\code{"w"}, default) or length
+(\code{"l"}), using the allometric weight-length relationship.}
+
 \item{return_data}{A boolean value that determines whether the formatted data
 used for the plot is returned instead of the plot itself. Default value is FALSE}
 }
 \value{
 A ggplot2 object, unless \code{return_data = TRUE}, in which case a data
-frame with the three variables 'w', 'value', 'Species' is returned.
+frame with the three variables 'w' (or 'l' if \code{size_axis = "l"}), 'value',
+'Species' is returned.
 }
 \description{
 After running a projection, plot the total fishing mortality of each species
diff --git a/man/plotFeedingLevel.Rd b/man/plotFeedingLevel.Rd
index fecafb388..25121cc3f 100644
--- a/man/plotFeedingLevel.Rd
+++ b/man/plotFeedingLevel.Rd
@@ -16,6 +16,9 @@ plotFeedingLevel(object, ...)
   highlight = NULL,
   all.sizes = FALSE,
   include_critical = FALSE,
+  wlim = c(NA, NA),
+  llim = c(NA, NA),
+  size_axis = c("w", "l"),
   return_data = FALSE,
   ...
 )
@@ -26,6 +29,9 @@ plotFeedingLevel(object, ...)
   highlight = NULL,
   all.sizes = FALSE,
   include_critical = FALSE,
+  wlim = c(NA, NA),
+  llim = c(NA, NA),
+  size_axis = c("w", "l"),
   return_data = FALSE,
   ...
 )
@@ -36,6 +42,9 @@ plotlyFeedingLevel(
   time_range,
   highlight = NULL,
   include_critical = FALSE,
+  wlim = c(NA, NA),
+  llim = c(NA, NA),
+  size_axis = c("w", "l"),
   ...
 )
 }
@@ -63,14 +72,30 @@ outside a species' size range. Default FALSE.}
 \item{include_critical}{If TRUE, then the critical feeding level is also
 plotted. Default FALSE.}
 
+\item{wlim}{A numeric vector of length two providing lower and upper limits
+for the w axis. Use NA for the default: the lower default is
+\code{min(params@w) / 100} when \code{resource = TRUE} (to show some resource below
+the fish grid) or \code{min(params@w)} when \code{resource = FALSE}; the upper
+default is \code{max(params@w_full)}. Data is filtered to this range and the
+axis limits are set accordingly.}
+
+\item{llim}{A numeric vector of length two providing lower and upper limits
+for the length axis when \code{size_axis = "l"}. Use \code{NA} to auto-scale to the
+data range. Data is filtered to this range and the axis limits are set
+accordingly.}
+
+\item{size_axis}{Whether to plot size as weight (\code{"w"}, default) or length
+(\code{"l"}), using the allometric weight-length relationship.}
+
 \item{return_data}{A boolean value that determines whether the formatted data
 used for the plot is returned instead of the plot itself. Default value is FALSE}
 }
 \value{
 A ggplot2 object, unless \code{return_data = TRUE}, in which case a data
-frame with the variables 'w', 'value' and 'Species' is returned. If also
-\code{include_critical = TRUE} then the data frame contains a fourth variable
-'Type' that distinguishes between 'actual' and 'critical' feeding level.
+frame with the variables 'w' (or 'l' if \code{size_axis = "l"}), 'value' and
+'Species' is returned. If also \code{include_critical = TRUE} then the data
+frame contains a fourth variable 'Type' that distinguishes between
+'actual' and 'critical' feeding level.
 }
 \description{
 After running a projection, plot the feeding level of each species by size.
diff --git a/man/plotM2.Rd b/man/plotM2.Rd
index a59aeb884..64ead598e 100644
--- a/man/plotM2.Rd
+++ b/man/plotM2.Rd
@@ -14,7 +14,8 @@ plotM2(object, ...)
 }
 \value{
 A ggplot2 object, unless \code{return_data = TRUE}, in which case a data
-frame with the three variables 'w', 'value', 'Species' is returned.
+frame with the three variables 'w' (or 'l' if \code{size_axis = "l"}), 'value',
+'Species' is returned.
 }
 \description{
 \ifelse{html}{\href{https://lifecycle.r-lib.org/articles/stages.html#deprecated}{\figure{lifecycle-deprecated.svg}{options: alt='[Deprecated]'}}}{\strong{[Deprecated]}}
diff --git a/man/plotPredMort.Rd b/man/plotPredMort.Rd
index 5f8bd1050..ec928a5e6 100644
--- a/man/plotPredMort.Rd
+++ b/man/plotPredMort.Rd
@@ -15,6 +15,9 @@ plotPredMort(object, ...)
   time_range,
   all.sizes = FALSE,
   highlight = NULL,
+  wlim = c(NA, NA),
+  llim = c(NA, NA),
+  size_axis = c("w", "l"),
   return_data = FALSE,
   ...
 )
@@ -24,11 +27,23 @@ plotPredMort(object, ...)
   species = NULL,
   all.sizes = FALSE,
   highlight = NULL,
+  wlim = c(NA, NA),
+  llim = c(NA, NA),
+  size_axis = c("w", "l"),
   return_data = FALSE,
   ...
 )
 
-plotlyPredMort(object, species = NULL, time_range, highlight = NULL, ...)
+plotlyPredMort(
+  object,
+  species = NULL,
+  time_range,
+  highlight = NULL,
+  wlim = c(NA, NA),
+  llim = c(NA, NA),
+  size_axis = c("w", "l"),
+  ...
+)
 }
 \arguments{
 \item{object}{An object of class \linkS4class{MizerSim} or
@@ -51,12 +66,28 @@ outside a species' size range. Default FALSE.}
 
 \item{highlight}{Name or vector of names of the species to be highlighted.}
 
+\item{wlim}{A numeric vector of length two providing lower and upper limits
+for the w axis. Use NA for the default: the lower default is
+\code{min(params@w) / 100} when \code{resource = TRUE} (to show some resource below
+the fish grid) or \code{min(params@w)} when \code{resource = FALSE}; the upper
+default is \code{max(params@w_full)}. Data is filtered to this range and the
+axis limits are set accordingly.}
+
+\item{llim}{A numeric vector of length two providing lower and upper limits
+for the length axis when \code{size_axis = "l"}. Use \code{NA} to auto-scale to the
+data range. Data is filtered to this range and the axis limits are set
+accordingly.}
+
+\item{size_axis}{Whether to plot size as weight (\code{"w"}, default) or length
+(\code{"l"}), using the allometric weight-length relationship.}
+
 \item{return_data}{A boolean value that determines whether the formatted data
 used for the plot is returned instead of the plot itself. Default value is FALSE}
 }
 \value{
 A ggplot2 object, unless \code{return_data = TRUE}, in which case a data
-frame with the three variables 'w', 'value', 'Species' is returned.
+frame with the three variables 'w' (or 'l' if \code{size_axis = "l"}), 'value',
+'Species' is returned.
 }
 \description{
 After running a projection, plot the predation mortality rate of each species
diff --git a/man/plotRelative.Rd b/man/plotRelative.Rd
index 75de88b66..e79390d86 100644
--- a/man/plotRelative.Rd
+++ b/man/plotRelative.Rd
@@ -17,7 +17,9 @@ plotRelative(x, y, ...)
   all.sizes = FALSE,
   log_x = TRUE,
   wlim = c(NA, NA),
+  llim = c(NA, NA),
   ylim = c(NA, NA),
+  size_axis = c("w", "l"),
   total = FALSE,
   background = TRUE,
   ...
@@ -44,7 +46,9 @@ plotRelative(x, y, ...)
   all.sizes = FALSE,
   log_x = TRUE,
   wlim = c(NA, NA),
+  llim = c(NA, NA),
   ylim = c(NA, NA),
+  size_axis = c("w", "l"),
   total = FALSE,
   background = TRUE,
   ...
@@ -68,10 +72,17 @@ spectra and \code{FALSE} for time series.}
 for the weight (x) axis. Use \code{NA} to refer to the existing minimum or
 maximum. Only applies to \code{ArraySpeciesBySize}.}
 
+\item{llim}{A numeric vector of length two providing lower and upper limits
+for the length (x) axis when \code{size_axis = "l"}. Use \code{NA} to refer to the
+existing minimum or maximum. Only applies to \code{ArraySpeciesBySize}.}
+
 \item{ylim}{A numeric vector of length two providing lower and upper limits
 for the value (y) axis. Use \code{NA} to refer to the existing minimum or
 maximum.}
 
+\item{size_axis}{Whether to plot size as weight (\code{"w"}, default) or length
+(\code{"l"}), using the allometric weight-length relationship.}
+
 \item{total}{A boolean value that determines whether the total over all
 selected species is plotted as well. Default is \code{FALSE}.}
 
diff --git a/man/plotSpectra.Rd b/man/plotSpectra.Rd
index e296df713..cd0546010 100644
--- a/man/plotSpectra.Rd
+++ b/man/plotSpectra.Rd
@@ -15,6 +15,7 @@ plotSpectra(object, ...)
   time_range,
   geometric_mean = FALSE,
   wlim = c(NA, NA),
+  llim = c(NA, NA),
   ylim = c(NA, NA),
   power = 1,
   biomass = TRUE,
@@ -25,6 +26,7 @@ plotSpectra(object, ...)
   log_x = TRUE,
   log_y = TRUE,
   log = NULL,
+  size_axis = c("w", "l"),
   return_data = FALSE,
   ...
 )
@@ -33,6 +35,7 @@ plotSpectra(object, ...)
   object,
   species = NULL,
   wlim = c(NA, NA),
+  llim = c(NA, NA),
   ylim = c(NA, NA),
   power = 1,
   biomass = TRUE,
@@ -43,6 +46,7 @@ plotSpectra(object, ...)
   log_x = TRUE,
   log_y = TRUE,
   log = NULL,
+  size_axis = c("w", "l"),
   return_data = FALSE,
   ...
 )
@@ -53,6 +57,7 @@ plotlySpectra(
   time_range,
   geometric_mean = FALSE,
   wlim = c(NA, NA),
+  llim = c(NA, NA),
   ylim = c(NA, NA),
   power = 1,
   biomass = TRUE,
@@ -63,6 +68,7 @@ plotlySpectra(
   log_x = TRUE,
   log_y = TRUE,
   log = NULL,
+  size_axis = c("w", "l"),
   ...
 )
 }
@@ -93,6 +99,11 @@ the fish grid) or \code{min(params@w)} when \code{resource = FALSE}; the upper
 default is \code{max(params@w_full)}. Data is filtered to this range and the
 axis limits are set accordingly.}
 
+\item{llim}{A numeric vector of length two providing lower and upper limits
+for the length axis when \code{size_axis = "l"}. Use \code{NA} to auto-scale to the
+data range. Data is filtered to this range and the axis limits are set
+accordingly.}
+
 \item{ylim}{A numeric vector of length two providing lower and upper limits
 for the y axis. Use NA to auto-scale to the data range. Values below 1e-20
 are always filtered out from the data regardless of \code{ylim[1]}. Data above
@@ -130,13 +141,16 @@ Default is TRUE.}
 in the same form as the base \code{\link[=plot]{plot()}} argument. For example, \code{"x"},
 \code{"y"}, \code{"xy"} or \code{""}. If supplied, this overrides \code{log_x} and \code{log_y}.}
 
+\item{size_axis}{Whether to plot size as weight (\code{"w"}, default) or length
+(\code{"l"}), using the allometric weight-length relationship.}
+
 \item{return_data}{A boolean value that determines whether the formatted data
 used for the plot is returned instead of the plot itself. Default value is FALSE}
 }
 \value{
 A ggplot2 object, unless \code{return_data = TRUE}, in which case a data
-frame with the four variables 'w', 'value', 'Species', 'Legend' is
-returned.
+frame with the four variables 'w' (or 'l' if \code{size_axis = "l"}), 'value',
+'Species', 'Legend' is returned.
 }
 \description{
 Plots the number density multiplied by a power of the weight, with the power
@@ -157,6 +171,7 @@ plotSpectra(sim, wlim = c(1e-6, NA))
 plotSpectra(sim, time_range = 10:20)
 plotSpectra(sim, time_range = 10:20, power = 0)
 plotSpectra(sim, species = c("Cod", "Herring"), power = 1)
+plotSpectra(sim, species = c("Cod", "Herring"), size_axis = "l")
 
 # Returning the data frame
 fr <- plotSpectra(sim, return_data = TRUE)
diff --git a/man/plotSpectra2.Rd b/man/plotSpectra2.Rd
index 4d52bed1a..ca9347f7d 100644
--- a/man/plotSpectra2.Rd
+++ b/man/plotSpectra2.Rd
@@ -14,6 +14,8 @@ plotSpectra2(
   log_x = TRUE,
   log_y = TRUE,
   log = NULL,
+  llim = c(NA, NA),
+  size_axis = c("w", "l"),
   ...
 )
 
@@ -26,6 +28,8 @@ plotlySpectra2(
   log_x = TRUE,
   log_y = TRUE,
   log = NULL,
+  llim = c(NA, NA),
+  size_axis = c("w", "l"),
   ...
 )
 }
@@ -49,6 +53,14 @@ to logarithmic size bins.}
 in the same form as the base \code{\link[=plot]{plot()}} argument. For example, \code{"x"},
 \code{"y"}, \code{"xy"} or \code{""}. If supplied, this overrides \code{log_x} and \code{log_y}.}
 
+\item{llim}{A numeric vector of length two providing lower and upper limits
+for the length axis when \code{size_axis = "l"}. Use \code{NA} to auto-scale to the
+data range. Data is filtered to this range and the axis limits are set
+accordingly.}
+
+\item{size_axis}{Whether to plot size as weight (\code{"w"}, default) or length
+(\code{"l"}), using the allometric weight-length relationship.}
+
 \item{...}{Arguments passed to \code{\link[=plotSpectra]{plotSpectra()}} for preparing the spectra
 data, for example \code{species}, \code{time_range}, \code{wlim}, \code{ylim}, \code{resource},
 \code{background} or \code{total}.}
diff --git a/man/plotSpectraRelative.Rd b/man/plotSpectraRelative.Rd
index 4325841ce..cd19af54f 100644
--- a/man/plotSpectraRelative.Rd
+++ b/man/plotSpectraRelative.Rd
@@ -5,9 +5,25 @@
 \alias{plotlySpectraRelative}
 \title{Plot the relative difference between two spectra}
 \usage{
-plotSpectraRelative(object1, object2, log_x = TRUE, ylim = c(NA, NA), ...)
+plotSpectraRelative(
+  object1,
+  object2,
+  log_x = TRUE,
+  ylim = c(NA, NA),
+  llim = c(NA, NA),
+  size_axis = c("w", "l"),
+  ...
+)
 
-plotlySpectraRelative(object1, object2, log_x = TRUE, ylim = c(NA, NA), ...)
+plotlySpectraRelative(
+  object1,
+  object2,
+  log_x = TRUE,
+  ylim = c(NA, NA),
+  llim = c(NA, NA),
+  size_axis = c("w", "l"),
+  ...
+)
 }
 \arguments{
 \item{object1}{First \code{MizerParams} or \code{MizerSim} object.}
@@ -20,6 +36,9 @@ plotlySpectraRelative(object1, object2, log_x = TRUE, ylim = c(NA, NA), ...)
 for the relative difference (y) axis. Use \code{NA} to refer to the existing
 minimum or maximum.}
 
+\item{size_axis}{Whether to plot size as weight (\code{"w"}, default) or length
+(\code{"l"}), using the allometric weight-length relationship.}
+
 \item{...}{Arguments passed to \code{\link[=plotSpectra]{plotSpectra()}} for preparing the spectra
 data, for example \code{species}, \code{time_range}, \code{wlim}, \code{resource},
 \code{background} or \code{total}.}
diff --git a/tests/testthat/test-plots.R b/tests/testthat/test-plots.R
index d1cfb6597..f662e82df 100644
--- a/tests/testthat/test-plots.R
+++ b/tests/testthat/test-plots.R
@@ -173,6 +173,19 @@ test_that("plotSpectra2 supports base plot log argument", {
                  "`log` must be a character string")
 })
 
+test_that("comparison helpers preserve non-size x variables", {
+    p <- plot2(getBiomass(NS_sim), getBiomass(NS_sim), species = "Cod")
+    expect_s3_class(p, "ggplot")
+    expect_true("Year" %in% names(p$data))
+    expect_error(ggplot2::ggplot_build(p), NA)
+
+    p_rel <- plotRelative(getBiomass(NS_sim), getBiomass(NS_sim),
+                          species = "Cod")
+    expect_s3_class(p_rel, "ggplot")
+    expect_true("Year" %in% names(p_rel$data))
+    expect_error(ggplot2::ggplot_build(p_rel), NA)
+})
+
 test_that("plotSpectraRelative plots symmetric relative difference", {
     params2 <- params
     params2@initial_n[] <- params@initial_n * 2
@@ -277,6 +290,94 @@ test_that("plotCDF2 compares cumulative distributions", {
                  "only supports log scaling on the x axis")
 })
 
+test_that("size-based plots support length axes", {
+    params_len <- params
+    params_len@species_params$a <- 0.01
+    params_len@species_params$b <- 3
+    sim_len <- sim
+    sim_len@params <- params_len
+
+    spectra_w <- plotSpectra(params_len, species = species, resource = FALSE,
+                             total = FALSE, return_data = TRUE)
+    spectra_l <- plotSpectra(params_len, species = species, resource = FALSE,
+                             total = FALSE, size_axis = "l",
+                             return_data = TRUE)
+    expect_true("l" %in% names(spectra_l))
+    expect_false("w" %in% names(spectra_l))
+    sp_idx <- match(as.character(spectra_w$Species),
+                    as.character(params_len@species_params$species))
+    expected_l <- w2l(spectra_w$w,
+                      params_len@species_params[sp_idx, , drop = FALSE])
+    expect_equal(spectra_l$l, expected_l)
+
+    llim <- stats::quantile(spectra_l$l, c(0.25, 0.75), names = FALSE)
+    spectra_l_limited <- plotSpectra(params_len, species = species,
+                                     resource = FALSE, total = FALSE,
+                                     size_axis = "l", llim = llim,
+                                     return_data = TRUE)
+    expect_true(all(spectra_l_limited$l >= llim[1]))
+    expect_true(all(spectra_l_limited$l <= llim[2]))
+
+    spectra_hidden <- plotSpectra(params_len, species = species,
+                                  resource = TRUE, total = TRUE,
+                                  size_axis = "l", return_data = TRUE)
+    expect_false(any(spectra_hidden$Legend %in% c("Resource", "Total")))
+
+    p <- plotSpectra(params_len, species = species, resource = FALSE,
+                     size_axis = "l")
+    expect_identical(p$scales$get_scales("x")$name, "Length [cm]")
+    expect_identical(p$scales$get_scales("x")$trans$name, "log-10")
+
+    expect_true("l" %in% names(plotCDF(params_len, species = species,
+                                       resource = FALSE, size_axis = "l",
+                                       return_data = TRUE)))
+    cdf_l_limited <- plotCDF(params_len, species = species, resource = FALSE,
+                             size_axis = "l", llim = llim,
+                             return_data = TRUE)
+    expect_true(all(cdf_l_limited$l >= llim[1]))
+    expect_true(all(cdf_l_limited$l <= llim[2]))
+    expect_equal(stats::aggregate(value ~ Species, cdf_l_limited, max)$value,
+                 rep(1, length(unique(cdf_l_limited$Species))))
+    expect_true("l" %in% names(plotSpectra2(params_len, params_len,
+                                            species = species,
+                                            resource = FALSE,
+                                            size_axis = "l")$data))
+    expect_true("l" %in% names(plotCDF2(params_len, params_len,
+                                        species = species,
+                                        resource = FALSE,
+                                        size_axis = "l")$data))
+    expect_true("l" %in% names(plotSpectraRelative(params_len, params_len,
+                                                   species = species,
+                                                   resource = FALSE,
+                                                   size_axis = "l")$data))
+
+    expect_true("l" %in% names(plotFeedingLevel(params_len, species = species,
+                                                size_axis = "l",
+                                                return_data = TRUE)))
+    expect_true("l" %in% names(plotPredMort(params_len, species = species,
+                                            size_axis = "l",
+                                            return_data = TRUE)))
+    expect_true("l" %in% names(plotFMort(params_len, species = species,
+                                         size_axis = "l",
+                                         return_data = TRUE)))
+    expect_true("l" %in% names(plotDiet(params_len, species = species[[1]],
+                                        size_axis = "l",
+                                        return_data = TRUE)))
+    expect_true(all(plotDiet(params_len, species = species[[1]],
+                             size_axis = "l", llim = llim,
+                             return_data = TRUE)$l >= llim[1]))
+    expect_true("l" %in% names(plot(getPredMort(params_len),
+                                    species = species, size_axis = "l",
+                                    return_data = TRUE)))
+    rate_l_limited <- plot(getPredMort(params_len), species = species,
+                           size_axis = "l", llim = llim,
+                           return_data = TRUE)
+    expect_true(all(rate_l_limited$l >= llim[1]))
+    expect_true(all(rate_l_limited$l <= llim[2]))
+    expect_true("l" %in% names(plot(getFMort(sim_len), species = species,
+                                    size_axis = "l", return_data = TRUE)))
+})
+
 test_that("yield plotting helpers validate comparison and gear selection", {
     sim_shifted <- sim
     dimnames(sim_shifted@n)$time <- as.character(10:13)
@@ -471,6 +572,10 @@ test_that("axis limits are set correctly", {
     expect_true(is.na(p$scales$scales[[1]]$limits[1]))
     expect_equal(p$scales$scales[[1]]$limits[2], 8)
 
+    p <- plotSpectra(sim, species = species, resource = FALSE,
+                     size_axis = "l", llim = c(10, 100))
+    expect_equal(p$scales$scales[[2]]$limits, c(1, 2))
+
     # Default wlim lower depends on resource argument
     p_res <- plotSpectra(sim, species = species, return_data = TRUE)
     p_nores <- plotSpectra(sim, species = species, resource = FALSE, return_data = TRUE)