diff --git a/NEWS.md b/NEWS.md
index 8561e01e2..b5f3fd03c 100644
--- a/NEWS.md
+++ b/NEWS.md
@@ -159,6 +159,36 @@ stability of steady states.
   copepod really is shorter than a 1 mg fish larva — but it does mean the
   resource and the species are measured differently.
 
+- The total line is now shown on length-based plots, where it used to be
+  dropped. A total can only be formed once every line sits on the same
+  coordinate, and on a length axis the lines do not: each species, and the
+  resource, converts weight to length with its own allometric relationship, so
+  at a given length they sit at different weights. The total is therefore
+  summed *after* the conversion — at equal length rather than at equal weight —
+  interpolating each series onto the union of all the size coordinates,
+  logarithmically in size, with a series contributing nothing outside its own
+  range. Where the series already share a grid, which is always the case on a
+  weight axis, the union is that grid and nothing is approximated: the
+  weight-axis total is unchanged.
+
+  `plotSpectra2()` and `plotSpectraRelative()` keep their total on a length
+  axis too. They used to convert the axis after assembling the two spectra, so
+  the total they had been given — summed at equal weight — arrived at the
+  conversion with no species to convert it by and was dropped. They now let
+  `plotSpectra()` convert, so the total they receive is already the total on
+  the axis being plotted. As a result `plotSpectra2()` also applies `ylim` the
+  way `plotSpectra()` does on a length axis: it could not before, because the
+  values it was filtering were a Jacobian away from the ones the limits
+  described, so with `return_data = TRUE` it returned values outside the limits
+  that a single spectrum plot would have dropped.
+
+  `total = TRUE` also now means the same thing everywhere: the total of
+  everything the object holds. For `plotSpectra()` that was already so — the
+  resource and every species, whatever is drawn — and it stays so. The array
+  plots have been brought into line: `plot(<array>, total = TRUE)` used to sum
+  only the species that were selected for display, and now sums the whole
+  array, so a plot of two species can be read against the community total.
+
 - `plotYieldObservedVsModel()` gains a `gear` argument that restricts the
   comparison to the catch of the selected gears. Both the model yield and the
   observed yield are then taken from those gears only, so in a model where
diff --git a/R/ArraySpeciesBySize-class.R b/R/ArraySpeciesBySize-class.R
index 65c3f9065..8d3264c0f 100644
--- a/R/ArraySpeciesBySize-class.R
+++ b/R/ArraySpeciesBySize-class.R
@@ -364,6 +364,10 @@ plot.ArraySpeciesBySize <- function(x, species = NULL,
     plot_dat <- convert_plot_density_axis(plot_dat, params, size_axis,
                                           density_wrt = array_density_wrt(x),
                                           per_log_size = per_log_size)
+    if (total) {
+        plot_dat <- append_total_line(plot_dat, total_contributors(x, wlim),
+                                      params, size_axis, x, per_log_size)
+    }
     if (identical(size_axis, "l")) {
         plot_dat <- filter_plot_length_limits(plot_dat, llim)
     }
@@ -557,7 +561,13 @@ plot2.ArraySpeciesBySize <- function(x, y, name1 = "First", name2 = "Second",
                             y_ticks = y_ticks, legend_var = "Legend",
                             size_axis = size_axis,
                             density_wrt = array_density_wrt(x),
-                            per_log_size = per_log_size)
+                            per_log_size = per_log_size,
+                            total_dat = if (total) {
+                                rbind(cbind(total_contributors(x, wlim),
+                                            Model = name1),
+                                      cbind(total_contributors(y, wlim),
+                                            Model = name2))
+                            })
 }
 
 #' Plot relative difference between two mizer arrays
@@ -675,7 +685,9 @@ plotRelative.ArraySpeciesBySize <- function(x, y, species = NULL,
                           xtrans = if (log_x) "log10" else "identity",
                           xlim = plot_size_xlim(wlim, size_axis, llim),
                           ylim = ylim,
-                          legend_var = "Legend", size_axis = size_axis)
+                          legend_var = "Legend", size_axis = size_axis,
+                          total_dat1 = if (total) total_contributors(x, wlim),
+                          total_dat2 = if (total) total_contributors(y, wlim))
 }
 
 check_plot2_compatible <- function(x, y, class) {
@@ -973,6 +985,10 @@ addPlot.ArraySpeciesBySize <- function(plot, x, species = NULL,
     plot_dat <- convert_plot_density_axis(plot_dat, params, size_axis,
                                           density_wrt = array_density_wrt(x),
                                           per_log_size = per_log_size)
+    if (total) {
+        plot_dat <- append_total_line(plot_dat, total_contributors(x, wlim),
+                                      params, size_axis, x, per_log_size)
+    }
     if (identical(size_axis, "l")) {
         plot_dat <- filter_plot_length_limits(plot_dat, llim)
     }
@@ -1127,11 +1143,6 @@ prepare_ArraySpeciesBySize_plot_data <- function(x, species = NULL,
     sel <- all_species %in% species
     mat <- unclass(x)[sel, , drop = FALSE]
 
-    # Compute total across all selected species before size-range trimming
-    if (total) {
-        total_row <- colSums(mat, na.rm = TRUE)
-    }
-
     plot_dat <- data.frame(
         w = rep(w, each = sum(sel)),
         value = c(mat),
@@ -1165,23 +1176,32 @@ prepare_ArraySpeciesBySize_plot_data <- function(x, species = NULL,
         }
     }
 
-    # Add total line
-    if (total) {
-        total_dat <- data.frame(
-            w = w,
-            value = total_row,
-            Species = "Total",
-            Legend = "Total"
-        )
-        total_dat <- apply_wlim(total_dat, wlim)
-        plot_dat <- rbind(plot_dat, total_dat)
-    }
-
     names(plot_dat)[2] <- value_name
 
     plot_dat
 }
 
+#' Assemble the contributors to the total of a species-by-size array
+#'
+#' The total is the total of everything the array holds: every species, whether
+#' or not it was selected for display, and every size, whether or not it falls
+#' in a species' own size range. It is a property of the array rather than of
+#' the plot, so that a plot of two species can still be read against the
+#' community total.
+#'
+#' The rows are returned unsummed, because the sum has to be taken after the
+#' size coordinate has been converted — on a length axis the species no longer
+#' share a grid; see [add_total_line()].
+#'
+#' @param x An `ArraySpeciesBySize` object.
+#' @param wlim Numeric vector of length two giving the weight limits.
+#' @return A data frame of plotting data holding every value in the array.
+#' @keywords internal
+total_contributors <- function(x, wlim = c(NA, NA)) {
+    prepare_ArraySpeciesBySize_plot_data(x, species = NULL, all.sizes = TRUE,
+                                         wlim = wlim, background = TRUE)
+}
+
 #' @rdname plotHover
 #' @usage NULL
 #' @examples
diff --git a/R/ArrayTimeBySpeciesBySize-class.R b/R/ArrayTimeBySpeciesBySize-class.R
index ade6d64b4..50b124780 100644
--- a/R/ArrayTimeBySpeciesBySize-class.R
+++ b/R/ArrayTimeBySpeciesBySize-class.R
@@ -434,12 +434,15 @@ animate.ArrayTimeBySpeciesBySize <- function(x, species = NULL,
                       stringsAsFactors = FALSE)
     df$value <- c(sub)
 
-    # Compute total across ALL selected species (including background) before
-    # any background filtering, matching the behaviour of plot.ArraySpeciesBySize
+    # The total is the total of everything the array holds, so its contributors
+    # are every species rather than the selected ones, matching
+    # `plot.ArraySpeciesBySize()`. They are summed only after the size axis has
+    # been converted, inside `animate_plotly()`.
+    total_dat <- NULL
     if (total) {
-        total_sums <- stats::aggregate(value ~ time + w, data = df, FUN = sum,
-                                       na.rm = TRUE)
-        total_sums$Species <- "Total"
+        total_dat <- expand.grid(time = times, Species = all_species, w = w,
+                                 stringsAsFactors = FALSE)
+        total_dat$value <- c(arr)
     }
 
     # Now handle background: exclude rows or group under "Background" legend
@@ -452,10 +455,7 @@ animate.ArrayTimeBySpeciesBySize <- function(x, species = NULL,
         }
     }
 
-    if (total) {
-        total_sums$legend_name <- "Total"
-        df <- rbind(df, total_sums[, names(df)])
-    }
+    if (!is.null(total_dat)) total_dat$legend_name <- "Total"
 
     y_label <- array_y_label(x, default = "Value", size_axis = size_axis,
                              per_log_size = per_log_size)
@@ -465,6 +465,7 @@ animate.ArrayTimeBySpeciesBySize <- function(x, species = NULL,
                    size_axis = size_axis,
                    density_wrt = array_density_wrt(x),
                    per_log_size = per_log_size,
+                   total_dat = total_dat,
                    frame_duration = frame_duration,
                    transition_duration = transition_duration,
                    easing = easing)
diff --git a/R/animateSpectra.R b/R/animateSpectra.R
index 8b0517f66..c45e1ceaf 100644
--- a/R/animateSpectra.R
+++ b/R/animateSpectra.R
@@ -197,17 +197,19 @@ animate.MizerSim <- function(x, species = NULL,
         nf <- rbind(nf, nf_back)
     }
     # Add total ----
+    # The contributors are assembled here but summed only after the size axis
+    # has been converted, inside `animate_plotly()`. The total is the total of
+    # everything the model holds: every species, whether or not it was selected
+    # for display, and the resource, whether or not it is drawn — matching
+    # `plotSpectra()`.
+    total_dat <- NULL
     if (total) {
-        # Calculate total community abundance
-        fish_idx <- (length(sim@params@w_full) -
-                         length(sim@params@w) + 1):length(sim@params@w_full)
-        total_n <- sim@n_pp
-        total_n[, fish_idx] <- total_n[, fish_idx] +
-            rowSums(aperm(sim@n, c(1, 3, 2)), dims = 2)
-        nf_total <- melt(total_n[time_elements, , drop = FALSE])
-        nf_total$Species <- "Total"
-        nf_total$legend_name <- "Total"
-        nf <- rbind(nf, nf_total)
+        total_dat <- melt(sim@n[time_elements, , , drop = FALSE])
+        names(total_dat)[names(total_dat) == "sp"] <- "Species"
+        nf_pp_total <- melt(sim@n_pp[time_elements, , drop = FALSE])
+        nf_pp_total$Species <- "Resource"
+        total_dat <- rbind(total_dat, nf_pp_total)
+        total_dat$legend_name <- "Total"
     }
 
     y_label <- spectra_y_label(power, size_axis,
@@ -220,13 +222,18 @@ animate.MizerSim <- function(x, species = NULL,
     if (isTRUE(sim@params@second_order_w[["bin_average"]])) {
         beta <- sim@params@w_full[2] / sim@params@w_full[1]
         nf$w <- nf$w * sqrt(beta)
+        if (!is.null(total_dat)) total_dat$w <- total_dat$w * sqrt(beta)
     }
     nf <- mutate(nf, value = value * w^power)
+    if (!is.null(total_dat)) {
+        total_dat <- mutate(total_dat, value = value * w^power)
+    }
 
     animate_plotly(nf, sim@params, log_x, log_y, y_label, wlim, llim,
                    ylim,
                    size_axis = size_axis,
                    density_wrt = spectrum_density_wrt(spectrum$per_log_size),
+                   total_dat = total_dat,
                    frame_duration = frame_duration,
                    transition_duration = transition_duration,
                    easing = easing)
@@ -243,14 +250,27 @@ animate_plotly <- function(df, params, log_x, log_y, y_label,
                            size_axis = "w",
                            density_wrt = NA_character_,
                            per_log_size = NULL,
+                           total_dat = NULL,
                            frame_duration = 500, transition_duration = 500,
                            easing = "linear") {
     size_axis <- plot_size_axis(size_axis)
-    df <- convert_plot_density_axis(df, params, size_axis,
-                                    density_wrt = density_wrt,
-                                    per_log_size = per_log_size,
-                                    value_col = "value")
+    convert <- function(d) {
+        convert_plot_density_axis(d, params, size_axis,
+                                  density_wrt = density_wrt,
+                                  per_log_size = per_log_size,
+                                  value_col = "value")
+    }
+    df <- convert(df)
     x_var <- plot_size_x_var(size_axis)
+    # The total is summed over the converted series, so that on a length axis
+    # it is a sum at equal length rather than at equal weight.
+    if (!is.null(total_dat)) {
+        total_dat <- add_total_line(convert(total_dat), x_var, "value",
+                                    by = "time")
+        total_dat <- total_dat[total_dat$Species == "Total", ]
+        total_dat$legend_name <- "Total"
+        df <- rbind(df, total_dat[, names(df), drop = FALSE])
+    }
     legend_name_order <- intersect(names(params@linecolour),
                                    unique(df$legend_name))
     sp_order <- unlist(lapply(legend_name_order, function(ln) {
diff --git a/R/plots.R b/R/plots.R
index 43ad512e7..210b6d163 100644
--- a/R/plots.R
+++ b/R/plots.R
@@ -408,6 +408,10 @@ plotHover.mizer_plot <- function(x = ggplot2::last_plot(), ...,
 #'   and are left alone when the size axis changes.
 #' @param per_log_size Whether to express a density per logarithmic size.
 #'   `NULL` (the default) keeps whichever the values already are.
+#' @param total_dat Optional data frame of the contributors to a total, with a
+#'   `Model` column identifying which of the two they belong to. The total is
+#'   summed from them after the size axis has been converted, see
+#'   [add_total_line()].
 #' @return A `mizer_plot` (ggplot2) object.
 #' @keywords internal
 plotComparisonDataFrame <- function(frame1, frame2, params,
@@ -419,7 +423,8 @@ plotComparisonDataFrame <- function(frame1, frame2, params,
                                     legend_var = "Legend",
                                     size_axis = NULL,
                                     density_wrt = NA_character_,
-                                    per_log_size = NULL) {
+                                    per_log_size = NULL,
+                                    total_dat = NULL) {
     assert_that(is.data.frame(frame1),
                 is.data.frame(frame2),
                 is(params, "MizerParams"))
@@ -448,6 +453,19 @@ plotComparisonDataFrame <- function(frame1, frame2, params,
                                            value_col = y_var)
         x_var <- plot_size_x_var(size_axis)
     }
+    if (!is.null(total_dat)) {
+        # The contributors arrive unconverted and are put on the same axis as
+        # the data they join before being summed there.
+        total_dat <- convert_plot_density_axis(total_dat, params, size_axis,
+                                               density_wrt = density_wrt,
+                                               per_log_size = per_log_size,
+                                               species_col = group_var,
+                                               value_col = y_var)
+        total_dat <- add_total_line(total_dat, x_var, y_var, by = "Model")
+        total_dat <- total_dat[total_dat[[group_var]] == "Total", ]
+        frame <- rbind(frame, total_dat[, names(frame), drop = FALSE])
+        frame$Model <- factor(frame$Model, levels = c(name1, name2))
+    }
 
     legend_levels <- intersect(names(params@linecolour), frame[[legend_var]])
     frame[[legend_var]] <- factor(frame[[legend_var]], levels = legend_levels)
@@ -508,7 +526,8 @@ plotRelativeDataFrame <- function(frame1, frame2, params,
                                   ylim = c(NA, NA),
                                   highlight = NULL,
                                   legend_var = "Legend",
-                                  size_axis = NULL) {
+                                  size_axis = NULL,
+                                  total_dat1 = NULL, total_dat2 = NULL) {
     assert_that(is.data.frame(frame1),
                 is.data.frame(frame2),
                 is(params, "MizerParams"))
@@ -524,18 +543,38 @@ plotRelativeDataFrame <- function(frame1, frame2, params,
              "in the data frame.")
     }
 
+    # The size axis is converted before the two frames are joined, so that a
+    # total can be formed on the axis it is plotted against. The relative
+    # difference itself is unaffected: a common Jacobian cancels.
+    if (!is.null(size_axis)) {
+        size_axis <- plot_size_axis(size_axis)
+        frame1 <- convert_plot_size_axis(frame1, params, size_axis,
+                                         species_col = group_var)
+        frame2 <- convert_plot_size_axis(frame2, params, size_axis,
+                                         species_col = group_var)
+        x_var <- plot_size_x_var(size_axis)
+    }
+    if (!is.null(total_dat1) && !is.null(total_dat2)) {
+        add_total <- function(frame, total_dat) {
+            if (!is.null(size_axis)) {
+                total_dat <- convert_plot_size_axis(total_dat, params,
+                                                    size_axis,
+                                                    species_col = group_var)
+            }
+            total_dat <- add_total_line(total_dat, x_var, y_var)
+            total_dat <- total_dat[total_dat[[group_var]] == "Total", ]
+            rbind(frame, total_dat[, names(frame), drop = FALSE])
+        }
+        frame1 <- add_total(frame1, total_dat1)
+        frame2 <- add_total(frame2, total_dat2)
+    }
+
     by_vars <- c(x_var, group_var, legend_var)
     frame <- dplyr::inner_join(frame1, frame2, by = by_vars,
                                suffix = c(".x", ".y"))
     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)
@@ -726,6 +765,97 @@ array_log_y <- function(x, log_y, log, given) {
     log_y
 }
 
+#' Add a total line to plotting data by summing over its series
+#'
+#' The total has to be formed *after* the size coordinate has been converted,
+#' not before. On a weight axis every series shares the model's weight grid, so
+#' summing at equal weight and summing at equal position are the same thing. On
+#' a length axis they are not: each species, and the resource, converts weight
+#' to length with its own allometric relationship, so at a given length the
+#' series sit at different weights and their grids no longer coincide. The sum
+#' that means something there is the sum at equal *length* — the number of
+#' organisms per unit length, whatever they are — which is what this computes.
+#'
+#' Each series is interpolated onto the sorted union of all the size
+#' coordinates, linearly in the logarithm of size, since the grid is
+#' logarithmic. A series contributes nothing outside its own range. When the
+#' series already share a grid — always on a weight axis, and on a length axis
+#' whenever the weight-length parameters agree — the union is that grid and the
+#' interpolation reproduces the values exactly, so nothing is approximated in
+#' the cases where nothing needs to be.
+#'
+#' @param plot_dat A data frame of plotting data with a size column, a value
+#'   column and a `Species` column.
+#' @param x_var Name of the size column. Defaults to the first column.
+#' @param value_col Name or index of the value column. Defaults to the second.
+#' @param by Names of further columns identifying separate plots, such as the
+#'   time of an animation frame or the model of a comparison. A total is formed
+#'   within each of their combinations.
+#' @return `plot_dat` with the total appended as a series named `"Total"`.
+#' @keywords internal
+add_total_line <- function(plot_dat, x_var = names(plot_dat)[[1]],
+                           value_col = 2, by = NULL) {
+    if (nrow(plot_dat) == 0) return(plot_dat)
+    if (is.numeric(value_col)) value_col <- names(plot_dat)[[value_col]]
+    if (length(by) > 0) {
+        groups <- split(plot_dat, plot_dat[, by, drop = FALSE], drop = TRUE)
+        out <- lapply(groups, add_total_line, x_var = x_var,
+                      value_col = value_col)
+        out <- do.call(rbind, unname(out))
+        rownames(out) <- NULL
+        return(out)
+    }
+    x_all <- sort(unique(plot_dat[[x_var]]))
+    total <- numeric(length(x_all))
+    for (series in split(plot_dat, plot_dat$Species)) {
+        series <- series[order(series[[x_var]]), , drop = FALSE]
+        x <- series[[x_var]]
+        y <- series[[value_col]]
+        if (length(x) < 2) {
+            # `approx()` needs two points; a single one contributes only where
+            # it sits.
+            hit <- match(x, x_all)
+            total[hit] <- total[hit] + y
+            next
+        }
+        contribution <- stats::approx(log(x), y, xout = log(x_all))$y
+        contribution[is.na(contribution)] <- 0
+        total <- total + contribution
+    }
+    total_dat <- plot_dat[rep(1, length(x_all)), , drop = FALSE]
+    total_dat[[x_var]] <- x_all
+    total_dat[[value_col]] <- total
+    total_dat$Species <- "Total"
+    if ("Legend" %in% names(total_dat)) total_dat$Legend <- "Total"
+    rownames(total_dat) <- NULL
+    rbind(plot_dat, total_dat)
+}
+
+#' Convert the contributors to a total and append the total line
+#'
+#' Wraps the two steps that every plot of an array with `total = TRUE` needs:
+#' the contributors are put on the same size axis as the data they will join,
+#' and then summed there, see [add_total_line()].
+#'
+#' @param plot_dat The converted plotting data to append the total to.
+#' @param total_dat The unconverted contributors, see [total_contributors()].
+#' @param params A MizerParams object.
+#' @param size_axis Either `"w"` (weight) or `"l"` (length).
+#' @param x The mizer array the data came from, which says whether the values
+#'   are a density.
+#' @param per_log_size Whether to express a density per logarithmic size.
+#' @return `plot_dat` with the total appended as a series named `"Total"`.
+#' @keywords internal
+append_total_line <- function(plot_dat, total_dat, params, size_axis, x,
+                              per_log_size = NULL) {
+    total_dat <- convert_plot_density_axis(total_dat, params, size_axis,
+                                           density_wrt = array_density_wrt(x),
+                                           per_log_size = per_log_size)
+    total_dat <- add_total_line(total_dat)
+    total_dat <- total_dat[total_dat$Species == "Total", , drop = FALSE]
+    rbind(plot_dat, total_dat[, names(plot_dat), drop = FALSE])
+}
+
 #' The weight-length parameters to plot each row of plotting data with
 #'
 #' A length axis needs an allometric weight-length relationship for every line
@@ -1833,13 +1963,6 @@ plot_spectra <- function(params, n, n_pp,
         wlim[2] <- max(w_full_grid)
     }
 
-    if (total) {
-        # Calculate total community abundance
-        fish_idx <- (length(params@w_full) - length(params@w) + 1):length(params@w_full)
-        total_n <- n_pp
-        total_n[fish_idx] <- total_n[fish_idx] + colSums(n)
-        total_n <- total_n * w_full_grid^power
-    }
     species <- valid_species_arg(params, species)
     y_label <- spectra_y_label(power, size_axis, biomass = biomass,
                                per_log_size = per_log_size)
@@ -1867,14 +1990,6 @@ plot_spectra <- function(params, n, n_pp,
             )
         }
     }
-    if (total) {
-        plot_dat <- rbind(plot_dat,
-                          data.frame(w = w_full_grid,
-                                     value = c(total_n),
-                                     Species = "Total",
-                                     Legend = "Total")
-                          )
-    }
     if (background && any(params@species_params$is_background)) {
         back_n <- n[params@species_params$is_background, , drop = FALSE]
         plot_dat <-
@@ -1892,6 +2007,36 @@ plot_spectra <- function(params, n, n_pp,
                              (plot_dat$w <= wlim[2]), ]
     plot_dat <- convert_plot_spectrum_axis(plot_dat, params, size_axis, power,
                                            per_log_size = per_log_size)
+    if (total) {
+        # The total is the total of everything the model holds: every species,
+        # whether or not it was selected for display, and the resource, whether
+        # or not it is drawn. It is summed over the series as they are plotted,
+        # which on a length axis means at equal length rather than at equal
+        # weight.
+        # `n` has already been multiplied by the power of weight above.
+        resource_dat <- data.frame(w = w_full_grid,
+                                   value = c(n_pp * w_full_grid^power),
+                                   Species = "Resource",
+                                   Legend = "Total")
+        total_dat <- rbind(data.frame(w = rep(w_grid, each = nrow(n)),
+                                      value = c(unclass(n)),
+                                      Species = rep(dimnames(n)[[1]],
+                                                    times = length(w_grid)),
+                                      Legend = "Total"),
+                           resource_dat)
+        # Zeros are kept here, unlike in the plotted data: they are part of
+        # each series and dropping them would let the interpolation bridge a
+        # gap that is really empty.
+        total_dat <- total_dat[total_dat$w >= wlim[1] &
+                                   total_dat$w <= wlim[2], ]
+        total_dat <- convert_plot_spectrum_axis(total_dat, params, size_axis,
+                                                power,
+                                                per_log_size = per_log_size)
+        total_dat <- add_total_line(total_dat)
+        total_dat <- total_dat[total_dat$Species == "Total" &
+                                   total_dat$value > 0, ]
+        plot_dat <- rbind(plot_dat, total_dat[, names(plot_dat), drop = FALSE])
+    }
     if (identical(size_axis, "l")) {
         plot_dat <- filter_plot_length_limits(plot_dat, llim)
     }
@@ -2356,19 +2501,21 @@ plotSpectra2 <- function(object1, object2, name1 = "First", name2 = "Second",
     log_y <- log_axes$log_y
     assert_that(length(wlim) == 2, length(llim) == 2, length(ylim) == 2)
 
-    data_ylim <- if (identical(size_axis, "l")) c(0, NA) else ylim
-    sf1 <- plotSpectra(object1, species = species, wlim = wlim,
-                       ylim = data_ylim,
+    # The size axis is converted by `plotSpectra()` rather than afterwards, so
+    # that the total it forms is the total on the axis being plotted. A total
+    # summed at equal weight cannot be converted to a length axis afterwards:
+    # it is a sum over species that each convert differently. It also means
+    # `ylim` can be handed straight over, because the values it is applied to
+    # are then the values it describes.
+    sf1 <- plotSpectra(object1, species = species, wlim = wlim, ylim = ylim,
                        power = power, total = total, resource = resource,
-                       background = background, size_axis = "w",
+                       background = background, size_axis = size_axis,
                        return_data = TRUE, ...)
-    sf2 <- plotSpectra(object2, species = species, wlim = wlim,
-                       ylim = data_ylim,
+    sf2 <- plotSpectra(object2, species = species, wlim = wlim, ylim = ylim,
                        power = power, total = total, resource = resource,
-                       background = background, size_axis = "w",
+                       background = background, size_axis = size_axis,
                        return_data = TRUE, ...)
     params <- if (is(object1, "MizerSim")) object1@params else object1
-    density_wrt <- spectrum_density_wrt(spectrum$per_log_size)
 
     plotComparisonDataFrame(sf1, sf2, validParams(params),
                             name1 = name1, name2 = name2,
@@ -2381,9 +2528,7 @@ plotSpectra2 <- function(object1, object2, name1 = "First", name2 = "Second",
                             ytrans = if (log_y) "log10" else "identity",
                             xlim = plot_size_xlim(wlim, size_axis, llim),
                             ylim = ylim, highlight = highlight,
-                            legend_var = "Legend",
-                            size_axis = size_axis,
-                            density_wrt = density_wrt)
+                            legend_var = "Legend")
 }
 
 #' Resolve the power of weight multiplying a spectrum
@@ -2581,13 +2726,15 @@ plotSpectraRelative <- function(object1, object2,
     size_axis <- plot_size_axis(size_axis)
     assert_that(length(wlim) == 2, length(llim) == 2, length(ylim) == 2)
 
+    # As in `plotSpectra2()`, the conversion is left to `plotSpectra()` so that
+    # the total is summed on the axis being plotted.
     sf1 <- plotSpectra(object1, species = species, wlim = wlim,
                        power = 1, total = total, resource = resource,
-                       background = background, size_axis = "w",
+                       background = background, size_axis = size_axis,
                        return_data = TRUE, ...)
     sf2 <- plotSpectra(object2, species = species, wlim = wlim,
                        power = 1, total = total, resource = resource,
-                       background = background, size_axis = "w",
+                       background = background, size_axis = size_axis,
                        return_data = TRUE, ...)
     params <- if (is(object1, "MizerSim")) object1@params else object1
 
@@ -2597,7 +2744,7 @@ plotSpectraRelative <- function(object1, object2,
                           xlim = plot_size_xlim(wlim, size_axis, llim),
                           ylim = ylim,
                           highlight = highlight,
-                          legend_var = "Legend", size_axis = size_axis)
+                          legend_var = "Legend")
 }
 
 #' @rdname plotSpectraRelative
diff --git a/inst/skills/analyse-and-plot/SKILL.md b/inst/skills/analyse-and-plot/SKILL.md
index 360273e4f..73bee0d23 100644
--- a/inst/skills/analyse-and-plot/SKILL.md
+++ b/inst/skills/analyse-and-plot/SKILL.md
@@ -224,7 +224,7 @@ dedicated `plot…()` functions below — share these optional arguments:
 | `wlim`/`llim` | numeric vector `c(min, max)` — restrict the size (x) axis |
 | `ylim` | numeric vector `c(min, max)` — restrict the value (y) axis |
 | `highlight` | character vector — draw named species with thicker lines |
-| `total` | logical — add a line for the community total |
+| `total` | logical — add a line for the community total. The total of *everything the object holds*, so it does not change when you select species or hide the resource; on a length axis it is summed at equal length |
 | `log_x`, `log_y` | logical — log-scale the x or y axis |
 | `size_axis` | `"w"` (default) or `"l"` — plot against weight or against length |
 
diff --git a/inst/skills/upgrade-mizer-code/SKILL.md b/inst/skills/upgrade-mizer-code/SKILL.md
index f45f82695..9df6d816e 100644
--- a/inst/skills/upgrade-mizer-code/SKILL.md
+++ b/inst/skills/upgrade-mizer-code/SKILL.md
@@ -52,6 +52,10 @@ plots) are in the changelog and are not repeated here.
 | A custom array plotted with `size_axis = "l"` is not transformed, or is transformed when it should not be | an array declares what it holds with `type` now | Arrays say what kind of value they hold (3.3) |
 | `array_spectrum_power()`, or `spectrum_power =` in an internal plot helper, is no longer found | replaced by the `type` metadata | Arrays say what kind of value they hold (3.3) |
 | New warning that two arrays hold "a value of type" different things in `plot2()` or `plotRelative()` | the two arrays disagree about what they hold | Arrays say what kind of value they hold (3.3) |
+| A `Total` line appears on a plot with `size_axis = "l"` where there used to be none | the total is now summed after the length conversion | The total is summed on the axis it is plotted against (3.3) |
+| `plot(<array>, species = ..., total = TRUE)` gives a bigger total than before | the total is the total of the whole array, not of the selected species | The total is summed on the axis it is plotted against (3.3) |
+| `plotSpectra2()` or `plotSpectraRelative()` with `size_axis = "l"` shows a `Total` line where it used to show none | the total is now formed on the axis being plotted | The total is summed on the axis it is plotted against (3.3) |
+| `plotSpectra2(size_axis = "l", ylim = ..., return_data = TRUE)` returns fewer rows | `ylim` now filters the data there as it does for `plotSpectra()` | The total is summed on the axis it is plotted against (3.3) |
 | New warning that a change to a species or resource parameter "has not taken effect" | the rate it feeds was set by hand and is no longer calculated | A change that cannot take effect now warns (3.3) |
 | That warning appears with `given_species_params<-()` but not with `species_params<-()` | the diagnostics belong to the given species parameter setter | The two species parameter setters divide the diagnostics between them (3.3) |
 | Setting a given species parameter to `NA` now warns that the change has not taken effect | clearing a value counts as a change, and a frozen array blocks it | The two species parameter setters divide the diagnostics between them (3.3) |
@@ -239,6 +243,50 @@ single `density_wrt` argument in place of `spectrum_power` and
 `spectrum_per_log_size`, and that the internal `array_spectrum_power()` is gone.
 The `power`-based interface of `plotSpectra()` and friends is unchanged.
 
+### The total is summed on the axis it is plotted against
+
+A total can only be formed once every line sits on the same coordinate. On a
+weight axis they do: every species shares the model's weight grid. On a length
+axis they do not, because each species — and now the resource — converts weight
+to length with its own allometric relationship, so at a given length the lines
+sit at different weights. That is why the `Total` line used to be dropped from
+length-based plots.
+
+It is now summed *after* the conversion, at equal length rather than at equal
+weight, interpolating each line onto the union of all the size coordinates
+(logarithmically in size, with a line contributing nothing outside its own
+range). Where the lines already share a grid — always on a weight axis, and on
+a length axis whenever the weight-length parameters agree — the union is that
+grid and the interpolation reproduces the values exactly. **The weight-axis
+total is unchanged**, for every power.
+
+`plotSpectra2()` and `plotSpectraRelative()` are fixed by the same change.
+They used to convert the size axis after assembling the two spectra, so the
+total they had been handed — already summed at equal weight — reached the
+conversion with no species to convert it by and was silently dropped. They now
+let `plotSpectra()` do the conversion, so the total they receive is the total
+on the axis being plotted.
+
+That also settles what `ylim` does there. `plotSpectra()` applies it both as
+the axis limits and as a filter on the data, with a hard floor at 1e-20.
+`plotSpectra2()` could not do the same on a length axis, because the values it
+was filtering were a Jacobian away from the ones the limits described, so it
+skipped the filter. Now that it converts first, the filter applies as it does
+everywhere else. The plot is unchanged — the axis limits hid those points
+anyway — but `return_data = TRUE` no longer hands back values outside the
+limits.
+
+One thing does change on the weight axis. `total = TRUE` now means the same
+thing everywhere: **the total of everything the object holds**, whatever is
+drawn. `plotSpectra()` always worked that way — the resource and every species
+count, whether or not the resource is shown and whichever species were
+selected — and it still does. The array plots did not: `plot(<array>,
+total = TRUE)` summed only the species selected for display, and only the sizes
+inside each species' own range. It now sums the whole array, so the total no
+longer moves when you change `species`, `all.sizes` or `background`, and a plot
+of two species can be read against the community total. If you were relying on
+the total of a selection, sum the selection yourself.
+
 ### Length and weight parameters follow the one you gave last
 
 A size can be given either as a weight (`w_mat`, `w_max`, …) or as the length it
diff --git a/man/add_total_line.Rd b/man/add_total_line.Rd
new file mode 100644
index 000000000..b257a4026
--- /dev/null
+++ b/man/add_total_line.Rd
@@ -0,0 +1,48 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/plots.R
+\name{add_total_line}
+\alias{add_total_line}
+\title{Add a total line to plotting data by summing over its series}
+\usage{
+add_total_line(
+  plot_dat,
+  x_var = names(plot_dat)[[1]],
+  value_col = 2,
+  by = NULL
+)
+}
+\arguments{
+\item{plot_dat}{A data frame of plotting data with a size column, a value
+column and a \code{Species} column.}
+
+\item{x_var}{Name of the size column. Defaults to the first column.}
+
+\item{value_col}{Name or index of the value column. Defaults to the second.}
+
+\item{by}{Names of further columns identifying separate plots, such as the
+time of an animation frame or the model of a comparison. A total is formed
+within each of their combinations.}
+}
+\value{
+\code{plot_dat} with the total appended as a series named \code{"Total"}.
+}
+\description{
+The total has to be formed \emph{after} the size coordinate has been converted,
+not before. On a weight axis every series shares the model's weight grid, so
+summing at equal weight and summing at equal position are the same thing. On
+a length axis they are not: each species, and the resource, converts weight
+to length with its own allometric relationship, so at a given length the
+series sit at different weights and their grids no longer coincide. The sum
+that means something there is the sum at equal \emph{length} — the number of
+organisms per unit length, whatever they are — which is what this computes.
+}
+\details{
+Each series is interpolated onto the sorted union of all the size
+coordinates, linearly in the logarithm of size, since the grid is
+logarithmic. A series contributes nothing outside its own range. When the
+series already share a grid — always on a weight axis, and on a length axis
+whenever the weight-length parameters agree — the union is that grid and the
+interpolation reproduces the values exactly, so nothing is approximated in
+the cases where nothing needs to be.
+}
+\keyword{internal}
diff --git a/man/append_total_line.Rd b/man/append_total_line.Rd
new file mode 100644
index 000000000..425a302f2
--- /dev/null
+++ b/man/append_total_line.Rd
@@ -0,0 +1,38 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/plots.R
+\name{append_total_line}
+\alias{append_total_line}
+\title{Convert the contributors to a total and append the total line}
+\usage{
+append_total_line(
+  plot_dat,
+  total_dat,
+  params,
+  size_axis,
+  x,
+  per_log_size = NULL
+)
+}
+\arguments{
+\item{plot_dat}{The converted plotting data to append the total to.}
+
+\item{total_dat}{The unconverted contributors, see \code{\link[=total_contributors]{total_contributors()}}.}
+
+\item{params}{A MizerParams object.}
+
+\item{size_axis}{Either \code{"w"} (weight) or \code{"l"} (length).}
+
+\item{x}{The mizer array the data came from, which says whether the values
+are a density.}
+
+\item{per_log_size}{Whether to express a density per logarithmic size.}
+}
+\value{
+\code{plot_dat} with the total appended as a series named \code{"Total"}.
+}
+\description{
+Wraps the two steps that every plot of an array with \code{total = TRUE} needs:
+the contributors are put on the same size axis as the data they will join,
+and then summed there, see \code{\link[=add_total_line]{add_total_line()}}.
+}
+\keyword{internal}
diff --git a/man/plotComparisonDataFrame.Rd b/man/plotComparisonDataFrame.Rd
index aecaacaa8..06fd35493 100644
--- a/man/plotComparisonDataFrame.Rd
+++ b/man/plotComparisonDataFrame.Rd
@@ -21,7 +21,8 @@ plotComparisonDataFrame(
   legend_var = "Legend",
   size_axis = NULL,
   density_wrt = NA_character_,
-  per_log_size = NULL
+  per_log_size = NULL,
+  total_dat = NULL
 )
 }
 \arguments{
@@ -57,6 +58,11 @@ and are left alone when the size axis changes.}
 
 \item{per_log_size}{Whether to express a density per logarithmic size.
 \code{NULL} (the default) keeps whichever the values already are.}
+
+\item{total_dat}{Optional data frame of the contributors to a total, with a
+\code{Model} column identifying which of the two they belong to. The total is
+summed from them after the size axis has been converted, see
+\code{\link[=add_total_line]{add_total_line()}}.}
 }
 \value{
 A \code{mizer_plot} (ggplot2) object.
diff --git a/man/plotRelativeDataFrame.Rd b/man/plotRelativeDataFrame.Rd
index 1f7ff1414..09d5f0c29 100644
--- a/man/plotRelativeDataFrame.Rd
+++ b/man/plotRelativeDataFrame.Rd
@@ -14,7 +14,9 @@ plotRelativeDataFrame(
   ylim = c(NA, NA),
   highlight = NULL,
   legend_var = "Legend",
-  size_axis = NULL
+  size_axis = NULL,
+  total_dat1 = NULL,
+  total_dat2 = NULL
 )
 }
 \arguments{
diff --git a/man/total_contributors.Rd b/man/total_contributors.Rd
new file mode 100644
index 000000000..8e7d166da
--- /dev/null
+++ b/man/total_contributors.Rd
@@ -0,0 +1,29 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/ArraySpeciesBySize-class.R
+\name{total_contributors}
+\alias{total_contributors}
+\title{Assemble the contributors to the total of a species-by-size array}
+\usage{
+total_contributors(x, wlim = c(NA, NA))
+}
+\arguments{
+\item{x}{An \code{ArraySpeciesBySize} object.}
+
+\item{wlim}{Numeric vector of length two giving the weight limits.}
+}
+\value{
+A data frame of plotting data holding every value in the array.
+}
+\description{
+The total is the total of everything the array holds: every species, whether
+or not it was selected for display, and every size, whether or not it falls
+in a species' own size range. It is a property of the array rather than of
+the plot, so that a plot of two species can still be read against the
+community total.
+}
+\details{
+The rows are returned unsummed, because the sum has to be taken after the
+size coordinate has been converted — on a length axis the species no longer
+share a grid; see \code{\link[=add_total_line]{add_total_line()}}.
+}
+\keyword{internal}
diff --git a/tests/testthat/test-ArraySpeciesBySize-class.R b/tests/testthat/test-ArraySpeciesBySize-class.R
index 94e643e9d..bfcf9650c 100644
--- a/tests/testthat/test-ArraySpeciesBySize-class.R
+++ b/tests/testthat/test-ArraySpeciesBySize-class.R
@@ -276,6 +276,62 @@ test_that("plot() can express a density per logarithmic size", {
     expect_true(any(with_total$Species == "Total"))
 })
 
+test_that("the total is shown on a length axis and matches the lines drawn", {
+    density <- initialN(NS_params_small)
+
+    on_w <- plot(density, total = TRUE, return_data = TRUE)
+    on_l <- plot(density, total = TRUE, size_axis = "l", return_data = TRUE)
+    expect_true("Total" %in% on_w$Species)
+    expect_true("Total" %in% on_l$Species)
+
+    # The species here share one weight-length relationship, so the total is
+    # the plain sum of the converted lines at each length
+    total_l <- on_l[on_l$Species == "Total", ]
+    lines_l <- on_l[on_l$Species != "Total", ]
+    summed <- vapply(total_l$l, function(len) {
+        sum(lines_l[[2]][lines_l$l == len])
+    }, numeric(1))
+    expect_equal(total_l[[2]], summed)
+
+    # The total is the total of everything in the array, so it depends on
+    # neither the species selection nor the size trimming
+    full <- unname(colSums(unclass(density)))
+    expect_equal(on_w[[2]][on_w$Species == "Total"], full)
+    for (args in list(list(all.sizes = TRUE), list(species = "Cod"),
+                      list(background = FALSE))) {
+        d <- do.call(plot, c(list(density, total = TRUE, return_data = TRUE),
+                             args))
+        expect_equal(d[[2]][d$Species == "Total"], full)
+    }
+})
+
+test_that("the total on a length axis interpolates onto a common grid", {
+    # Give two species their own weight-length relationships so that their
+    # length grids no longer coincide
+    params <- NS_params_small
+    species_params(params)$a <- c(0.008, 0.012, 0.01)
+    species_params(params)$b <- c(3.1, 2.9, 3)
+    density <- initialN(params)
+
+    on_w <- plot(density, total = TRUE, return_data = TRUE)
+    on_l <- plot(density, total = TRUE, size_axis = "l", return_data = TRUE)
+    # The weight axis still has one point per size, the length axis needs the
+    # union of the three grids
+    expect_identical(sum(on_w$Species == "Total"), length(params@w))
+    expect_gt(sum(on_l$Species == "Total"), length(params@w))
+
+    # At the largest length only one species reaches, so the total is exactly
+    # that species there. The total covers the whole array, so compare it
+    # against the untrimmed lines.
+    total_l <- on_l[on_l$Species == "Total", ]
+    all_l <- plot(density, all.sizes = TRUE, size_axis = "l",
+                  return_data = TRUE)
+    lines_l <- all_l[all_l$Species != "Total", ]
+    at <- max(total_l$l)
+    expect_identical(sum(lines_l$l == at), 1L)
+    expect_equal(total_l[[2]][total_l$l == at], lines_l[[2]][lines_l$l == at])
+})
+
 test_that("a proportion is plotted against the whole of [0, 1]", {
     feeding_level <- getFeedingLevel(NS_params_small)
     expect_identical(array_type(feeding_level), "proportion")
diff --git a/tests/testthat/test-plots.R b/tests/testthat/test-plots.R
index 00d8e2c82..803c3720b 100644
--- a/tests/testthat/test-plots.R
+++ b/tests/testthat/test-plots.R
@@ -508,14 +508,14 @@ test_that("size-based plots support length axes", {
     expect_true(all(spectra_l_limited$l >= llim[1]))
     expect_true(all(spectra_l_limited$l <= llim[2]))
 
-    # The resource has its own weight-length relationship and so appears on a
-    # length axis; the total still does not, because it is a sum over lines
-    # that each have their own.
+    # The resource and the total both appear on a length axis: the resource
+    # has its own weight-length relationship, and the total is summed after
+    # the conversion, at equal length.
     spectra_all <- plotSpectra(params_len, species = species,
                                resource = TRUE, total = TRUE,
                                size_axis = "l", return_data = TRUE)
     expect_true("Resource" %in% spectra_all$Legend)
-    expect_false("Total" %in% spectra_all$Legend)
+    expect_true("Total" %in% spectra_all$Legend)
 
     p <- plotSpectra(params_len, species = species, resource = FALSE,
                      size_axis = "l")
@@ -616,12 +616,20 @@ test_that("length-axis spectra transform their densities", {
         compared_data$l
     expect_equal(compared_data[[2]], unname(by_weight[[2]] * jacobian))
 
+    # `ylim` is applied to the converted values, not to the weight-axis ones it
+    # would be off by a Jacobian from, and it filters the data exactly as it
+    # does for a single spectrum plot.
     limit <- stats::median(compared_data[[2]])
     limited <- plotSpectra2(params, params, species = species,
                             resource = FALSE, power = 0, size_axis = "l",
                             ylim = c(limit, NA))
-    expect_true(any(limited$data[[2]] < limit))
+    expect_true(all(limited$data[[2]] >= limit))
     expect_equal(limited$scales$get_scales("y")$limits, c(log10(limit), NA))
+    single <- plotSpectra(params, species = species, resource = FALSE,
+                          power = 0, size_axis = "l", ylim = c(limit, NA),
+                          return_data = TRUE)
+    expect_equal(sort(limited$data[[2]][limited$data$Model == "First"]),
+                 sort(single[[2]]))
 })
 
 test_that("proportion_ylim widens to [0, 1] without ever hiding data", {
@@ -642,6 +650,125 @@ test_that("proportion_ylim widens to [0, 1] without ever hiding data", {
                      c(0, 1))
 })
 
+test_that("add_total_line sums series on a shared grid exactly", {
+    plot_dat <- data.frame(w = rep(c(1, 2, 4), 2), value = 1:6,
+                           Species = rep(c("a", "b"), each = 3))
+    out <- add_total_line(plot_dat)
+    total <- out[out$Species == "Total", ]
+    expect_identical(total$w, c(1, 2, 4))
+    expect_equal(total$value, c(1 + 4, 2 + 5, 3 + 6))
+    # The original rows are untouched
+    expect_equal(out[out$Species != "Total", names(plot_dat)], plot_dat,
+                 ignore_attr = TRUE)
+})
+
+test_that("add_total_line interpolates series that do not share a grid", {
+    # `b` lives on a shifted grid, so the total is evaluated on the union and
+    # each series is interpolated onto it, contributing nothing outside its own
+    # range.
+    plot_dat <- rbind(
+        data.frame(w = c(1, 4), value = c(1, 4), Species = "a"),
+        data.frame(w = c(2, 8), value = c(10, 40), Species = "b"))
+    total <- add_total_line(plot_dat)
+    total <- total[total$Species == "Total", ]
+    expect_identical(total$w, c(1, 2, 4, 8))
+    # At w = 1 only `a`; at w = 8 only `b`
+    expect_equal(total$value[[1]], 1)
+    expect_equal(total$value[[4]], 40)
+    # At w = 2, `a` interpolated linearly in log w between (1, 1) and (4, 4)
+    a_at_2 <- stats::approx(log(c(1, 4)), c(1, 4), xout = log(2))$y
+    expect_equal(total$value[[2]], a_at_2 + 10)
+    # At w = 4, `b` interpolated between (2, 10) and (8, 40)
+    b_at_4 <- stats::approx(log(c(2, 8)), c(10, 40), xout = log(4))$y
+    expect_equal(total$value[[3]], 4 + b_at_4)
+})
+
+test_that("add_total_line groups by the columns it is told to", {
+    plot_dat <- data.frame(w = rep(c(1, 2), 4), value = 1:8,
+                           Species = rep(c("a", "b"), each = 2),
+                           time = rep(c(1, 1, 2, 2), each = 2))
+    out <- add_total_line(plot_dat, by = "time")
+    totals <- out[out$Species == "Total", ]
+    expect_identical(nrow(totals), 4L)
+    # time 1 holds a = (1, 2) and b = (3, 4); time 2 holds a = (5, 6) and
+    # b = (7, 8)
+    expect_equal(totals$value[totals$time == 1], c(1 + 3, 2 + 4))
+    expect_equal(totals$value[totals$time == 2], c(5 + 7, 6 + 8))
+})
+
+test_that("the spectrum total is unchanged on a weight axis", {
+    # The total is now summed after the size axis is converted rather than
+    # before, which on a weight axis must make no difference at all.
+    for (power in 0:2) {
+        d <- plotSpectra(params, total = TRUE, power = power,
+                         return_data = TRUE)
+        total <- d[d$Legend == "Total", ]
+        fish_idx <- (length(params@w_full) - length(params@w) +
+                         1):length(params@w_full)
+        expected <- params@initial_n_pp
+        expected[fish_idx] <- expected[fish_idx] + colSums(params@initial_n)
+        expected <- expected * params@w_full^power
+        keep <- expected > 0 & params@w_full >= min(params@w) / 100
+        expect_equal(total[[2]], unname(expected[keep]))
+    }
+})
+
+test_that("the spectrum total appears on a length axis", {
+    on_l <- plotSpectra(params, total = TRUE, size_axis = "l",
+                        return_data = TRUE)
+    total <- on_l[on_l$Legend == "Total", ]
+    expect_gt(nrow(total), 0)
+
+    # All the species here share one weight-length relationship, so at each
+    # length the total is exactly the sum of the series drawn there
+    contributors <- on_l[on_l$Legend != "Total", ]
+    at <- total$l[[which.max(total$l)]]
+    expect_equal(total[[2]][total$l == at],
+                 sum(contributors[[2]][contributors$l == at]))
+
+    # The total is the total of everything the model holds, so it does not
+    # depend on what is drawn: neither hiding the resource nor selecting
+    # species changes it
+    small <- min(params@w)
+    total_at <- function(...) {
+        d <- plotSpectra(params, total = TRUE, return_data = TRUE, ...)
+        d[[2]][d$Legend == "Total" & d$w == small]
+    }
+    expect_equal(total_at(resource = FALSE), total_at(resource = TRUE))
+    expect_equal(total_at(species = species), total_at())
+})
+
+test_that("the comparison spectra keep their total on a length axis", {
+    params2 <- params
+    species_params(params2)$gamma <- species_params(params2)$gamma * 1.1
+
+    for (size_axis in c("w", "l")) {
+        compared <- plotSpectra2(params, params2, total = TRUE,
+                                 size_axis = size_axis)$data
+        expect_true("Total" %in% compared$Legend)
+        relative <- plotSpectraRelative(params, params2, total = TRUE,
+                                        size_axis = size_axis)$data
+        expect_true("Total" %in% relative$Legend)
+    }
+
+    # Each model in the comparison holds exactly what plotSpectra() draws for
+    # it on its own, total included, on either axis
+    for (size_axis in c("w", "l")) {
+        for (power in 0:2) {
+            compared <- plotSpectra2(params, params, total = TRUE,
+                                     power = power,
+                                     size_axis = size_axis)$data
+            single <- plotSpectra(params, total = TRUE, power = power,
+                                  size_axis = size_axis, return_data = TRUE)
+            expect_equal(sort(compared[[2]][compared$Model == "First"]),
+                         sort(single[[2]]))
+            expect_equal(sort(compared[[2]][compared$Legend == "Total" &
+                                                compared$Model == "First"]),
+                         sort(single[[2]][single$Legend == "Total"]))
+        }
+    }
+})
+
 test_that("validate_density_wrt accepts only the known measures", {
     expect_identical(validate_density_wrt(NULL), NA_character_)
     expect_identical(validate_density_wrt(NA), NA_character_)
diff --git a/vignettes/cheatsheet-analysis-and-plotting.Rmd b/vignettes/cheatsheet-analysis-and-plotting.Rmd
index f3caf46ba..2d3ac8744 100644
--- a/vignettes/cheatsheet-analysis-and-plotting.Rmd
+++ b/vignettes/cheatsheet-analysis-and-plotting.Rmd
@@ -237,7 +237,7 @@ dedicated `plot…()` functions below — share these optional arguments:
 | `wlim`/`llim` | numeric vector `c(min, max)` — restrict the size (x) axis |
 | `ylim` | numeric vector `c(min, max)` — restrict the value (y) axis |
 | `highlight` | character vector — draw named species with thicker lines |
-| `total` | logical — add a line for the community total |
+| `total` | logical — add a line for the community total. The total of *everything the object holds*, so it does not change when you select species or hide the resource; on a length axis it is summed at equal length |
 | `log_x`, `log_y` | logical — log-scale the x or y axis |
 | `size_axis` | `"w"` (default) or `"l"` — plot against weight or against length |
 
diff --git a/vignettes/upgrading.Rmd b/vignettes/upgrading.Rmd
index 85d24fa92..0017e5559 100644
--- a/vignettes/upgrading.Rmd
+++ b/vignettes/upgrading.Rmd
@@ -123,6 +123,50 @@ single `density_wrt` argument in place of `spectrum_power` and
 `spectrum_per_log_size`, and that the internal `array_spectrum_power()` is gone.
 The `power`-based interface of `plotSpectra()` and friends is unchanged.
 
+### The total is summed on the axis it is plotted against
+
+A total can only be formed once every line sits on the same coordinate. On a
+weight axis they do: every species shares the model's weight grid. On a length
+axis they do not, because each species — and now the resource — converts weight
+to length with its own allometric relationship, so at a given length the lines
+sit at different weights. That is why the `Total` line used to be dropped from
+length-based plots.
+
+It is now summed *after* the conversion, at equal length rather than at equal
+weight, interpolating each line onto the union of all the size coordinates
+(logarithmically in size, with a line contributing nothing outside its own
+range). Where the lines already share a grid — always on a weight axis, and on
+a length axis whenever the weight-length parameters agree — the union is that
+grid and the interpolation reproduces the values exactly. **The weight-axis
+total is unchanged**, for every power.
+
+`plotSpectra2()` and [`plotSpectraRelative()`](../reference/plotSpectraRelative.html) are fixed by the same change.
+They used to convert the size axis after assembling the two spectra, so the
+total they had been handed — already summed at equal weight — reached the
+conversion with no species to convert it by and was silently dropped. They now
+let `plotSpectra()` do the conversion, so the total they receive is the total
+on the axis being plotted.
+
+That also settles what `ylim` does there. `plotSpectra()` applies it both as
+the axis limits and as a filter on the data, with a hard floor at 1e-20.
+`plotSpectra2()` could not do the same on a length axis, because the values it
+was filtering were a Jacobian away from the ones the limits described, so it
+skipped the filter. Now that it converts first, the filter applies as it does
+everywhere else. The plot is unchanged — the axis limits hid those points
+anyway — but `return_data = TRUE` no longer hands back values outside the
+limits.
+
+One thing does change on the weight axis. `total = TRUE` now means the same
+thing everywhere: **the total of everything the object holds**, whatever is
+drawn. `plotSpectra()` always worked that way — the resource and every species
+count, whether or not the resource is shown and whichever species were
+selected — and it still does. The array plots did not: `plot(<array>,
+total = TRUE)` summed only the species selected for display, and only the sizes
+inside each species' own range. It now sums the whole array, so the total no
+longer moves when you change `species`, `all.sizes` or `background`, and a plot
+of two species can be read against the community total. If you were relying on
+the total of a selection, sum the selection yourself.
+
 ### Length and weight parameters follow the one you gave last
 
 A size can be given either as a weight (`w_mat`, `w_max`, …) or as the length it