Another workaround using patchwork. But still waiting for facet_wrap approach.

library(ggplot2)
library(dplyr)
library(tidyr)
library(patchwork)

windowsFonts(Palatino = windowsFont("Palatino Linotype"))

plot_histograms_patchwork <- function(data, columns = where(is.numeric), rename_xlab = NULL) {  
  # Select columns properly
  selected_cols <- if (is.numeric(columns)) {
    names(data)[columns]  # Select by index
  } else {
    names(select(data, {{ columns }}))  # Select by condition (e.g., numeric)
  }  
  plots <- list()  # Store individual plots  
  for (col in selected_cols) {
    x_label <- if (!is.null(rename_xlab) && col %in% names(rename_xlab)) rename_xlab[[col]] else "Value"    
    # Conditional y-axis title only for "Ca"
    y_label <- if (col == "Ca") "Frequency" else NULL    
    p <- ggplot(data, aes(x = .data[[col]])) +
      geom_histogram(aes(y = after_stat(count)), bins = 30, fill = "#69b3a2", color = "#e9ecef", alpha = 0.9, na.rm = TRUE) +
      theme_bw() +
      theme(
        text = element_text(family = "Palatino"),
        axis.title.x = element_text(size = 12),
        axis.title.y = if (col == "Ca") element_text(size = 12) else element_blank(),
        plot.title = element_blank()  # Remove titles
      ) +
      labs(x = x_label, y = y_label)    
    plots <- append(plots, list(p))
  }  
  # Arrange all plots using patchwork with 3 columns
  final_plot <- wrap_plots(plots) + plot_layout(ncol = 3)  
  return(final_plot)
}


plot_hist_all <- plot_histograms_patchwork(
  data = R_macro_rev, 
  columns = 7:15, 
  rename_xlab = c("N" = "Total N (mg/kg)", "P" = "Total P (mg/kg)", "K" = "Total K (mg/kg)", 
                  "Ca" = "Total Ca (cmol(+)/kg)", "Mg" = "Total Mg (cmol(+)/kg)", 
                  "NP" = "N:P ", "NK" = "N:K ", "PK" = "P:K ", "CaMg" = "Ca:Mg ")
)

plot_hist_all

the resulting plot here: