An option often overlooked by people is that you can pretty much draw anything in ggplot2 as long as you can express it in polygons. The downside is that it is a lot more work.
The following approach mimics my answer given here but adapted a little bit to work more consistently for multiple densities. It abandons the {ggridges} approach completely.
The function below can be used to slice up an arbitrary polygon along the y-position.
library(ggplot2)
library(ggridges)
library(polyclip)
fade_polygon <- function(x, y, yseq = seq(min(y), max(y), length.out = 100)) {
poly <- data.frame(x = x, y = y)
# Create bounding-box edges
xlim <- range(poly$x) + c(-1, 1)
# Pair y-edges
grad <- cbind(head(yseq, -1), tail(yseq, -1))
# Add vertical ID
grad <- cbind(grad, seq_len(nrow(grad)))
# Slice up the polygon
grad <- apply(grad, 1, function(range) {
# Create bounding box
bbox <- data.frame(x = c(xlim, rev(xlim)),
y = c(range[1], range[1:2], range[2]))
# Do actual slicing
slice <- polyclip::polyclip(poly, bbox)
# Format as data.frame
for (i in seq_along(slice)) {
slice[[i]] <- data.frame(
x = slice[[i]]$x,
y = slice[[i]]$y,
value = range[3],
id = c(1, rep(0, length(slice[[i]]$x) - 1))
)
}
slice <- do.call(rbind, slice)
})
# Combine slices
grad <- do.call(rbind, grad)
# Create IDs
grad$id <- cumsum(grad$id)
return(grad)
}
Next, we need to calculate densities manually for every month and apply the function above for each one of those densities.
# Split by month and calculate densities
densities <- split(lincoln_weather, lincoln_weather$Month)
densities <- lapply(densities, function(df) {
dens <- density(df$`Mean Temperature [F]`)
data.frame(x = dens$x, y = dens$y)
})
# Extract x/y positions
x <- lapply(densities, `[[`, "x")
y <- lapply(densities, `[[`, "y")
# Make sequence to max density
ymax <- max(unlist(y))
yseq <- seq(0, ymax, length.out = 100) # 100 can be any large enough number
# Apply function to all densities
polygons <- mapply(fade_polygon, x = x, y = y, yseq = list(yseq),
SIMPLIFY = FALSE)
Next, we need to add the information about Months back into the data.
# Count number of observations in each of the polygons
rows <- vapply(polygons, nrow, integer(1))
# Combine all of the polygons
polygons <- do.call(rbind, polygons)
# Assign month information
polygons$month_id <- rep(seq_along(rows), rows)
Lastly we plot these polygons with vanilla ggplot2. The (y / ymax) * scale does a similar scaling to what ggridges does and adding the month_id offsets each month from oneanother.
scale <- 3
ggplot(polygons, aes(x, (y / ymax) * scale + month_id,
fill = value, group = interaction(month_id, id))) +
geom_polygon(aes(colour = after_scale(fill)), size = 0.3) +
scale_y_continuous(
name = "Month",
breaks = seq_along(rows),
labels = names(rows)
) +
scale_fill_viridis_c()

Created on 2021-09-12 by the reprex package (v2.0.1)