Un primo suggerimento che evita in qualche modo la soluzione che sto cercando è quello di combinare tutti i miei dati in una tabella di dati e utilizzare facet_grid sul mio variabile e simulazione
ggplot() + ... + facet_grid(variable~simulation, scales = 'free_y')
Questo produce una trama guardando bene che visualizza i dati in una figura, ma può diventare ingombrante se si considerano molte simulazioni. 
Per "hackerare" la trama per produrre ciò che voglio, ho prima determinato i limiti che desideravo per ogni variabile meteorologica. Questi limiti sono stati rilevati osservando le maggiori estensioni per tutte le simulazioni di interesse. Una volta determinato, ho creato una piccola tabella di dati con le stesse colonne dei miei dati di simulazione e l'ho aggiunta alla fine. I miei dati di simulazione ha avuto la struttura
'year' 'month' 'variable' 'run' 'mean'
1973 1 'rhmax' 1 65.44
1973 2 'rhmax' 1 67.44
... ... ... ... ...
2011 12 'windmin' 200 0.4
così ho creato una nuova tabella di dati con le stesse colonne
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))
che dà
'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
Mentre la scelta di anni e run è aribrary, la scelta di mese nee d essere qualsiasi cosa al di fuori 1:12. Poi allegate questo miei dati di simulazione
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')
Quando tracciare questi dati con i limiti y, limito i valori dell'asse x a quelle nei dati originali. La tabella di dati aggiunta con i limiti y è mese valori di 13. Poiché ggplot continua a ridimensionare gli assi sull'intero dataset, anche se gli assi sono limitati, ciò mi dà i limiti che desidero. È importante notare che se ci sono valori di dati superiori ai limiti specificati, ciò non funzionerà.
Prima: notare le differenze nei limiti y per ogni variabile meteo tra i pannelli.

Dopo: Ora i limiti y rimangono coerenti per ogni variabile meteo tra i pannelli. 

spero di modificare questo post nei prossimi giorni e aggiungere un esempio riproducibile per una migliore spiegazione. Per favore commenta se hai sentito qualcosa sull'aggiunta di questa funzionalità a ggplot.
Credo che se avessi intenzione di mostrare una linea di tendenza o una sorta di stat, dovresti usare scale_x_discrete (llmits = ...). xlim è l'abbreviazione di coord_cartesian che è essenzialmente solo uno zoom e quindi ggplot calcolerebbe anche una statistica con i dati aggiunti. forse qualcuno può confermare. – Dominik