The first suggestion that wraps up the solution I'm looking for a bit is to combine all my data into one data table and use facet_grid for my variable and simulation
ggplot() + ... + facet_grid(variable~simulation, scales = 'free_y')
This creates a beautiful plot that displays the data in one drawing, but can become cumbersome when considering many simulations. 
To βcrackβ the idea of ββcreating what I want, I first determined what limits I wanted for each weather variable. These limits were found by looking at the greatest extent for all simulations of interest. After defining, I created a small data table with the same columns as my simulation data, and added it to the end. My simulation data was structured
'year' 'month' 'variable' 'run' 'mean' 1973 1 'rhmax' 1 65.44 1973 2 'rhmax' 1 67.44 ... ... ... ... ... 2011 12 'windmin' 200 0.4
So, I created a new data table with the same columns
ylims.sims <- data.table(year = 1, month = 13, variable = rep(c('rhmax','rhmin','sradmean','tmax','tmin','windmax','windmin'), each = 2), run = 201, mean = c(20, 100, 0, 80, 100, 350, 25, 40, 12, 32, 0, 8, 0, 2))
What gives
'year' 'month' 'variable' 'run' 'mean' 1 13 'rhmax' 201 20 1 13 'rhmax' 201 100 1 13 'rhmin' 201 0 1 13 'rhmin' 201 80 1 13 'sradmean' 201 100 1 13 'sradmean' 201 350 1 13 'tmax' 201 25 1 13 'tmax' 201 40 1 13 'tmin' 201 12 1 13 'tmin' 201 32 1 13 'windmax' 201 0 1 13 'windmax' 201 8 1 13 'windmin' 201 0 1 13 'windmin' 201 2
Although the choice of year and mileage is arbitrary, the choice of the month should be something outside 1:12. Then I added this to my simulation data.
sim1data.ylims <- rbind(sim1data, ylims) ggplot() + geom_boxplot(data = sim1data.ylims, aes(x = factor(month), y = mean)) + facet_wrap(~variable, scale = 'free_y') + xlab('month') + xlim('1','2','3','4','5','6','7','8','9','10','11','12')
When I draw this data with y constraints, I limit the x-axis values ββto those specified in the source data. The enclosed data table with y constraints has values ββof month 13. Since ggplot still scales the axes to the entire dataset, even if the axes are limited, this gives me the y constraints that I want. It is important to note that if the data values ββare greater than the limits you specified, this will not work.
Before: Note the difference in y values ββfor each weather variable between panels.


After: now the y limits remain the same for each weather variable between panels. 

I hope to edit this post in the coming days and add a reproducible example for a better explanation. Comment if you have heard anything about adding this feature to ggplot .