A globally calibrated scheme for generating daily meteorology from monthly statistics : Global-WGEN ( GWGEN ) v 1 . 0

While a wide range of Earth system processes occur at daily and even subdaily timescales, many global vegetation and other terrestrial dynamics models historically used monthly meteorological forcing both to reduce computational demand and because global datasets were lacking. Recently, dynamic land surface modeling has moved towards resolving daily and subdaily processes, and global datasets containing daily and subdaily meteorology have become available. These meteorological datasets, however, cover only the instrumental era of the last approximately 120 years at best, are subject to considerable uncertainty, and represent extremely large data files with associated computational costs of data input/output and file transfer. For periods before the recent past or in the future, global meteorological forcing can be provided by climate model output, but the quality of these data at high temporal resolution is low, particularly for daily precipitation frequency and amount. Here, we present GWGEN, a globally applicable statistical weather generator for the temporal downscaling of monthly climatology to daily meteorology. Our weather generator is parameterized using a global meteorological database and simulates daily values of five common variables: minimum and maximum temperature, precipitation, cloud cover, and wind speed. GWGEN is lightweight, modular, and requires a minimal set of monthly mean variables as input. The weather generator may be used in a range of applications, for example, in global vegetation, crop, soil erosion, or hydrological models. While GWGEN does not currently perform spatially autocorrelated multi-point downscaling of daily weather, this additional functionality could be implemented in future versions.


Introduction
The development of the first global vegetation models in the 1970s (e.g., Lieth, 1975) brought about the demand for meteorological forcing datasets with global extent and relatively high spatial resolution, e.g., 1 • × 1 • .While a global weather-station-based monthly climate dataset was available at this time (Walter and Lieth, 1967), limitations in computers and storage allowed only the simplest treatment of these data.The first global simulations of the net primary productivity of the terrestrial biosphere (Lieth, 1975) thus used rasterized polygons of annual meteorological variables that had been crudely interpolated from the station-based climatology.The next decade saw the development of better computers and more sophisticated global vegetation models (Prentice et al., 1992;Prentice, 1989) that recognized the need for forcing at a subannual time step, and development of these models was done in parallel with the first global, gridded high-resolution (0.5 • ) monthly climatology (Leemans and Cramer, 1991).At the time, monthly meteorological data were the only feasible global data that could be produced in terms of the raw station data available to feed the interpolation process, the processing time required to produce gridded maps, and the data storage and transfer capabilities of contemporary computer systems and networks.Global gridded monthly climate data thus became the standard for not only large-extent vegetation modeling (Haxeltine and Prentice, 1996;Haxeltine et al., 1996;Kaplan et al., 2003;Kucharik et al., 2000;Woodward et al., 1995) but also for a wide range of studies on biodiversity and species distribution (e.g., Elith et al., 2006), vegetation trace gas emissions (e.g., Guenther et al., 1995), and even the geographic distribution of human diseases (e.g., Bhatt et al., 2013).
Published by Copernicus Publications on behalf of the European Geosciences Union.
Over subsequent years, the global gridded monthly climate datasets were improved (New et al., 1999(New et al., , 2002)), developed with very high spatial resolution (Hijmans et al., 2005), and expanded beyond climatological mean climate to cover continuous time series over decades (Harris et al., 2014;Mitchell and Jones, 2005;New et al., 2000).The latter was an essential requirement for forcing dynamic global vegetation models (DGVMs) (e.g., Sitch et al., 2003).However, despite increasing quality, spatial resolution, and temporal extent in these datasets, the basic time step remained monthly, partly for legacy reasons -models had been developed in an earlier era subject to computational limitations and therefore used a monthly time step for efficiency even if this was no longer strictly a constraint -and partly because of the challenge in developing a global, high-resolution climate dataset with a daily or shorter time step still presented a major data management challenge.
On the other hand, there was increasing awareness that accurate simulation of many Earth surface processes required representation of processes at a shorter-than-monthly time step.Global simulation of surface hydrology (Gerten et al., 2004), crop growth (Bondeau et al., 2007), or biogeophysical processes (Krinner et al., 2005) needed submonthly forcing to produce reliable results.To address this need for better forcing data, two main approaches were taken: either monthly climate data were downscaled online using a stochastic weather generator (e.g., Pfeiffer et al., 2013), or a subdaily, high-resolution, gridded climate time series was generated directly by merging high-temporal-resolution reanalysis data (e.g., NCEP, 6 h, 2.5 • ) with high-spatialresolution monthly climate data (e.g., CRU, 0.5 • ).The latter process resulted in the CRUNCEP dataset (Viovy and Ciais, 2016;Wei et al., 2014), which, while global, is large even by modern standards (approximately 350 GB), is not available at spatial resolution greater than 0.5 • , and covers only the period 1901-2014.Forcing data for global vegetation and other models with shorter-than-monthly resolution at higher spatial resolutions than 0.5 • , or for any other period than the last approximately 120 years, e.g., for the future or the more distant past, may therefore only be available through downscaling techniques.One approach to overcome the limitations of currently available datasets could be to use general circulation model (GCM) output directly; however, most GCM output currently available does not have greater than 0.5 • spatial resolution, with the current generation of GCMs typically approaching 1 • × 1 • .Furthermore, there is a general observation that daily meteorology produced by GCMs is not realistic, particularly for precipitation (Dai, 2006;Stephens et al., 2010;Sun et al., 2006).An alternative approach is, therefore, to perform temporal downscaling on monthly meteorological data using a statistical weather generator.
Statistical weather generators were first developed primarily for crop and hydrological modeling at the field to catchment scale (Richardson, 1981;Woolhiser and Pegram, 1979;Woolhiser and Roldan, 1982).The weather generator was parameterized using daily meteorological observations at one or more weather stations close to the area of interest, although some attempts were made to generalize the parameterization over larger, subcontinental regions (e.g., Wilks, 1999bWilks, , 1998;;Woolhiser and Roldán, 1986).Locally parameterized weather generators have been applied to a very wide range of studies (Wilks, 2010;Wilks and Wilby, 1999) and enhanced to include additional meteorological variables beyond the original precipitation, temperature, and solar radiation (e.g., Parlange and Katz, 2000).Applications of a weather generator at continental to global scales was still limited, however, because of the need to perform local parameterization.
The need to simulate daily meteorology in regions of the world with short, unreliable, or unavailable daily meteorological time series brought about the realization that certain features of weather generator parameterization might be generalized across a range of climates (Geng and Auburn, 1987;Geng et al., 1986).This ultimately led to the development of globally applicable weather generators (Friend, 1998) and their incorporation in DGVMs (Bondeau et al., 2007;Gerten et al., 2004;Pfeiffer et al., 2013).The original global parameterization (Geng et al., 1986) of these weather generators was, however, limited to seven weather stations, mostly in the temperate latitudes.Friend (1998) does not publish the parameters used in his global weather generator, but we assume these were the same as the original Geng and Auburn (1987) and Geng et al. (1986) models.Given the availability of (1) large datasets of daily meteorology and (2) computers powerful enough to process these data, we therefore decided that it would be valuable to revisit these parameterizations, perform a systematic and quantitative evaluation of the resulting downscaled meteorology, and potentially improve our ability to perform monthly to daily downscaling of common meteorological variables with a single, globally applicable parameterization.
In the following sections, we describe Global-WGEN (GWGEN), a weather generator parameterized using more than 50 million daily weather observations from all continents and latitudes.We demonstrate how updated schemes for simulating precipitation occurrence and amount, and for bias correcting wind speed, further improve the quality of the model simulations.We perform an extensive model evaluation and parameter uncertainty analysis in order to settle on a parameter set that provides the most accurate, globally applicable results.We comment on the limitations of the model and priorities for future research.
GWGEN is an open-source, stand-alone model that may be incorporated into any number of models designed to work at global scale, including, e.g., vegetation, hydrology, climatology, and animal distribution models.After smoothing the monthly input, the Markov chain is used to decide whether it is a dry or a wet day.If it is a wet day, we draw a random number from the gamma-GP (generalized Pareto) distribution.Furthermore, the other means of the variables (T min/max , c, w) are adjusted and their daily values are calculated using the estimated standard deviations and residuals.The wind speed furthermore undergoes a square root transformation before applying the cross correlation and in the end is corrected using the bias correction.A quality check in the end restricts our model to be within a 5 % range of the observed total precipitation and to replicate the number of wet days from the input.

Model description
GWGEN requires the following six monthly summary values as input: (1) total monthly precipitation, (2) the number of days in the month with measurable precipitation (i.e., wet days), (3, 4) monthly mean daily minimum and maximum temperature, (5) mean cloud fraction, and (6) wind speed.
The model outputs are the same variables at daily resolution.This section summarizes the basic workflow in the model which is also shown schematically in Fig. 1 and Algorithm 1.
The first approximation of the daily variables comes from smoothing the monthly time series using a mean-preserving algorithm (Rymes and Myers, 2001).
For precipitation, we then first use the Markov chain approach (Sect.3.2.1) to decide the wet/dry state of the day.If it is a wet day, we calculate the gamma parameters using Eqs.( 7) and ( 8).The resulting distribution allows us to draw a random number -the precipitation amount of the currently simulated day.If we are above the threshold µ, we draw a second random number from the generalized Pareto (GP) distribution parameterized via Eq.( 9) and the chosen GP shape.
The next step modifies the means of temperature, wind speed, and cloud fraction depending on the wet/dry state of the day (lines 11 and 15 in Algorithm 1).After that, we use the cross-correlation approach described in Richardson (1981) (lines 18-20 and Sect. 3.2.6)and calculate the daily values of these variables.Finally, we use the quantile-based bias correction described in Sect.3.4 to correct the simulated wind speed.Set T min, i = T min, wet , T max, i = T max, wet , c i = c wet , w i = w wet from Eqs. ( 10) and ( 12) and Tables 1-12: Set σ T min , i = σ T min , wet , σ Tmax, i = σ Tmax, wet , σ w, i = σ w, wet , σ c, i = σ c, wet from Eqs. ( 11), ( 13) and ( 14) and Tables 1-3 13: else 14: Set P i = 0 mm day −1 15: Set T min, i = T min, dry , T max, i = T max, dry , c i = c dry , w i = w dry from Eqs. ( 10) and ( 12) and Tables 1 and 3 16: Set σ T min , i = σ T min , dry , σ Tmax, i = σ Tmax, dry , σ w, i = σ w, dry , σ c, i = σ c, dry from Eqs. (11), ( 13) and ( 14) and Tables 1-3 Apply bias correction w (Eq.23) 22: j = j + 1 23: end for 24: end while 25: end for www.geosci-model-dev.net/10/1/2017/We restrict the weather generator to reproduce the exact number of wet days (±1) as the input and to be within a 5 % range of the total monthly precipitation (with a maximum allowed deviation of 0.5 mm).If the program cannot produce these results, the procedure described above is repeated (see line 4).

Model development
GWGEN is based on the WGEN weather generator (Richardson, 1981), using the method of defining the model parame-ters based on monthly summaries described by Geng et al. (1986) and Geng and Auburn (1987).GWGEN diverges from the original WGEN by using a hybrid-order Markov chain to simulate precipitation occurrence (Wilks, 1999a) and a hybrid gamma-GP distribution (Furrer and Katz, 2008;Neykov et al., 2014) to estimate precipitation amount.Temperature, cloud cover, and wind speed are calculated following Richardson (1981), using cross correlation and depending on the wet/dry state of the day.We further add a quantilebased bias correction for wind speed and minimum temperature, which improves the simulation results significantly.
In the following subsections, we first describe the global weather station database used to develop and evaluate the model, then describe the underlying relationships that we use to define GWGEN's parameters.

Development of a global weather station database
To parameterize GWGEN, we assembled a global dataset of daily meteorological observations.Precipitation and minimum and maximum daily temperature come from the daily Global Historical Climatology Network (GHCN-Daily) database (Menne et al., 2012a, b).The GHCN-Daily consists of observations collected at approximately 100 000 weather stations on all continents and many oceanic islands.As the GHCN-Daily stations are highly concentrated in some parts of the world, particularly in the conterminous United States, we selected stations for our study using a geographic antialiasing filter to avoid an especially strong geographic bias in the generation of the model parameters.Dividing the world up into a 0.5 • grid, we selected the single station with the longest record in each cell, if one was present.While the GHCN-Daily units for precipitation have a nominal precision of 0.1 mm, several of the stations in the US reported precipitation in fractions of an inch, which were later converted to mm.To ensure uniform precision across all of our calibration stations (this was particularly important when generating the probability density functions for precipitation amount), we selected only those GHCN-Daily stations where all precipitation amounts between 0.1 and 1.0 mm day −1 were reported in the record.This resulted in 9508 stations covering all continents, although the distribution was strongly heterogenous, with the majority of the stations in North America, despite our geographic filter (Fig. 2, top panel).For cloud cover, wind speed, and calculation of cross correlations between temperature, cloud cover, and wind speed, we used the Extended Edited Cloud Report Archive (EECRA) database (Hahn and Warren, 1999).The geographic distribution of the 6978 EECRA stations we selected is different than the GHCN-Daily, with more stations in Europe (Fig. 2, middle panel), but overall a relatively similar number of stations were used from both datasets.For the observations from both GHCN-Daily and EECRA, we made one additional filtering step, selecting only complete months, i.e., months with no days having missing observations, for further processing.Finally, we reserved some weather station records for model evaluation that were not used for model parameterization.These were individual stations or two stations separated by a maximum distance of 1 km, where all of the daily meteorological variables that GWGEN simulates (P , T min , T max , c, w) were recorded on the same dates in the EECRA database.This merged selection from EECRA and GHCN resulted in a set of 921 stations representing approximately 15 million daily records, with observations on all continents, although the geographic distribution is once again highly heterogenous, with a particularly high density of stations in Japan and Germany (Fig. 2, bottom panel).; the others were forced to (0, 0).The underlying data for the fits correspond to the means of the multi-year series for each month for each station.

Precipitation occurrence
Following Geng et al. (1986), we expect to find a good relationship between the fraction of days in a month with measurable precipitation and the probability that any given day will be wet.Following Wilks (1999a), we use a hybridorder model that retains first-order Markov dependence for wet spells but allows second-order dependence for dry sequences; this hybrid-order scheme has been shown to be a good compromise between performance and simplicity.To parameterize the precipitation occurrence part of the model, we thus calculated transition probabilities for a wet day being followed by a wet day (p 11 ), for a wet day being followed by a dry day being followed by a wet day (p 101 ), and for two dry days being followed by a wet day (p 001 ).We perform this analysis on a station-and month-wise basis: we first ex-tract each of the (complete) Januaries, Februaries, etc. for a given station and then merge all of the Januaries (Februaries, Marches, etc.) for this station into a single series representing each month.Merging months over several years is particularly important for stations that have relatively little precipitation in a given month; for example, it could take several years of observations to observe a single p 101 event.The final transition probabilities were then regressed against the fraction of days in the month with precipitation, which show the characteristic linear relationship described by Geng et al.
Because the transition probabilities (p 001 and p 101 ) must be zero by definition when the fraction of wet days (f wet ) is zero, i.e., a completely dry month, we force the linear regression between these quantities to pass through the origin.Likewise, we require the regression line for p 11 to equal 1 when f wet is 1.One has to note, however, that this methodol-ogy artificially increases the R 2 coefficient for the fit because we fix the intercept (see, for example Gordon, 1981).
The analysis results in the following relationships: In the weather generator (see line 6 in Algorithm 1), we determine if any given day will have precipitation by calculating the appropriate probability density function selected from Eqs. ( 1)-( 3) on the basis of the precipitation state of the previous day (or two).Comparing the calculated probability from the selected equation with a random number u ∈ [0, 1], a precipitation day is simulated if u is greater than its corresponding probability.

Precipitation amount
Following the original WGEN (Richardson, 1981), GWGEN disaggregates precipitation amount using a statistical distribution.A number of different probability density functions have been used to estimate precipitation amount in weather generators including, e.g., single exponential or mixed exponential, one-or two-parameter gamma, or Weibull distribution (Wilks and Wilby, 1999).The strong relationship between the gamma scale parameter and the mean precipitation on wet days noted by Geng et al. (1986) makes generation of precipitation amounts with only monthly input data feasible.It is based upon the fact that the expected value of a gamma random variable equals the product of its two parameters, i.e., E( ) = αθ .The gamma distribution, however, shows poor performance in simulating highprecipitation events consistent with observations.Furrer and Katz (2008) and Neykov et al. (2014) suggest that a hybrid probability density function, based on both gamma and the GP distribution, has superior accuracy in simulating extreme precipitation events when compared to gamma alone.Because of its superior accuracy and ease of implementation, we therefore adopt the hybrid gamma-GP distribution for simulating precipitation amount in GWGEN.
The probability density function (pdf) of the gamma distribution is defined as where α > 0 is the shape and θ > 0 the scale parameter.The pdf of the GP distribution is defined via with σ > 0 being the scale parameter and ξ ∈ R the shape parameter.µ is the location parameter.
Following Furrer and Katz (2008), we define the hybrid gamma-GP pdf as where F (µ) describes the cumulative gamma distribution function at the threshold µ.In our weather generator, however, we first draw a random number from the gamma distribution and, if we are above the threshold, we draw another random number from the GP distribution.Thus, the frequency of precipitation events larger than µ is determined by the gamma distribution, but the actual amount of precipitation simulated when above the threshold µ is determined by the GP distribution (Furrer and Katz, 2008).
To determine the parameters of the hybrid distribution for precipitation, we started with the simple strategy by Geng et al. (1986).As above, when calculating the Markov chain parameters, we created multi-year series for each of the parameterization stations for each month and extracted the days with precipitation.If a series contained more than 100 entries, we fit a gamma distribution using maximum likelihood to it in order to estimate the α and θ parameters.
Following Geng et al. (1986), we then fit a regression line of the gamma scale parameter against the mean precipitation on wet days p d (see Fig. 4) and found the relationship As proposed by Geng et al. (1986), we use this relationship in our model to estimate the scale parameter of the distribution.Using this approach, the gamma shape parameter α is a constant, given via The GP scale parameter σ on the other hand is calculated during the simulation following Neykov et al. (2014) via The other parameters of the GP distribution are obtained through a sensitivity analysis described in Sect.3.5.

Temperature
Following the standard WGEN methodology (Richardson, 1981;Geng et al., 1986), daily temperature is determined through two processes: first, the wet/dry state of the day and then the cross correlation (Sect.3.2.6).
In the weather generator, we know from the Markov chain (Sect.3.2.1)whether the current simulated day is a wet or dry day.Based upon the simple linear relationships Table 1.Fit results of temperature correlation for wet and dry days for Figs. 5, 6, 10, and 11.The coefficients c 0 to c 3 correspond to the coefficients used in Eqs. ( 10) and ( 14).
we adjust the monthly mean x of the variable x ∈ {T min , T max }.
To estimate the values of the parameters c 0 and c 1 in the above equations, we follow the same procedure as for the parameters of the Markov chain (Sect.3.2.1).We extracted the complete months for T min and T max from the GHCN-Daily dataset and created a multi-year series for each month and station.We then regressed the mean on wet and dry days separated against the overall mean of each month (Figs. 5  and 6).Through this procedure, we estimate the parameters necessary for Eq.(10) (see Table 1).
To estimate residual noise, we also need an estimate of the standard deviation of the variable (see Sect. 3.2.6). Figure 7 shows the correlation between standard deviation on wet and dry days and the corresponding mean.The means of the standard deviations (black bars in Fig. 7) indicate a strong but nonlinear relationship between the standard deviation and the corresponding mean.The correlation changes particularly at 0 • C. We therefore use two different polynomials of order 5 for the values below and above the freezing point.p 1 in Eq. ( 11) denotes a polynomial of order 1; p 5 a polynomial of order 5.The coefficients of the different polynomials are shown in Table 2.
These coefficients are based on the means of the standard deviation (black bars in Fig. 7).We chose this procedure to give the same weight to all temperatures.Otherwise, the fit would be dominated by the temperature values around the freezing points.

Cloud fraction
Monthly mean cloud fraction is disaggregated, as for temperature, using the standard WGEN procedure of adding statistical noise to a wet-or dry-day mean and accounting for cross correlation among the different weather variables.For the parameterization of the cloud fraction equations, we used Geosci.Model Dev., 10, 3771-3791, 2017 www.geosci-model-dev.net/10/3771/2017/   1.
the EECRA dataset.The original dataset contains eight measurements per day of the total cloud cover in units of octas, i.e., values ranging from 0 (clear sky) to 8 (overcast).Hence, to calculate the daily cloud fraction, those values were averaged and divided by 8 to produce a daily mean.
To adjust the monthly mean depending on the wet/dry state of the day, we could not use a simple linear relationship as we used for temperature because cloud fraction is bounded by a lower limit of 0 and an upper limit of 1. Furthermore, we observed that cloud cover on wet days is usually greater than or equal to the monthly mean cloud cover, whereas the cloud cover on dry days is usually less than or equal to the monthly mean cloud cover.This results in a concave curve for the wet case and a convex curve for dry days.We used a qualitative graphical analysis to develop "best guess" equa-tions that had the desired shape and propose the following formulae for the regression linking cloud cover on wet or dry days to the overall mean: with a c, wet < 0 and a c, dry > 0.
The standard deviation of cloud cover fraction becomes 0 when the mean monthly cloud fraction reaches both the minimum or maximum limits of 0 and 1.Hence, for c sd, dry and c sd, wet we have an concave parabola with the formula , not the multi-year monthly series as for, e.g., mean temperature; Figs. 5 and 6).Parameters of the fits are also shown in Table 1.
with a c, wet , a c, dry ≥ 0. Results of the fits can be seen in Figs. 8 and 9 and the parameters in Table 3.

Wind speed
The parameterization of the mean wind speed is based upon the same linear Eq. ( 10) as temperature.For the standard  3.   1.

Cross correlation
Following Richardson (1981) we use cross correlation to add additional residual noise to the simulated meteorological variables, which provides more realism in the daily weather result.This methodology, based on Matalas (1967), preserves the serial and the cross correlation between the simulated variables.It implies that the serial correlation of each variable may be described by a first-order linear autoregressive model.
Given the cross-correlation matrix M 0 ∈ R 4 × R 4 and the lag- The matrices A, B, M 0 , and M 1 are calculated using the stations from the EECRA database in Fig. 2. The results are   The columns and rows in the two matrices correspond to minimum and maximum temperature, cloud fraction, and square root of wind speed, respectively.
In the weather generator, the variables T min , T max , c, and w are then calculated using a combination of residual noise χ i (where i denotes the current simulated day) and the mean of the variables.χ i is determined by the other variables and the previous day using A and B from above (Richardson, 1981;Matalas, 1967).Hence, χ i is given via The daily values for the variables are then calculated via T max, i = χ T max • σ T max , wet/dry + T max, wet/dry with σ T min , wet/dry , σ T max , wet/dry from Eq. ( 11), σ c, wet/dry from Eq. ( 13), σ w, wet/dry from Eq. ( 14), T min, wet/dry , T max, wet/dry , w wet/dry from Eq. ( 10), and c wet/dry from Eq. ( 12).Since this procedure always requires the residuals from the previous day, χ i−1 , we initialize χ 0 with 0, simulate the month, and then simulate it again.
Note that, through the entire procedure, wind speed is subject to a square-root transformation (also when calculating M 0 and M 1 ) to account for the fact that it is not normally distributed.

Model evaluation
To evaluate GWGEN, we started with the daily meteorology at the evaluation stations described above and calculated monthly summaries.We used these monthly data to drive the model and simulate daily meteorology.The resulting daily series now has the same length as the observed meteorology from the GHCN and EECRA databases.Because we cannot expect the weather generator to reproduce the weather exactly as observed (for example, the number of rainy days in a month may be the same as observed but they may not occur in precisely the same order), our evaluation is restricted to comparing the statistical properties of the input observed versus the output simulated daily meteorology.
Figure 12 shows the comparison of simulated versus observed values for each of the five meteorological variables handled by GWGEN.For temperature, wind, and cloud fraction, the model does an excellent job of downscaling monthly input to daily resolution. 1 The comparison between precipitation amounts looks good when considering all of the data; however, a closer look into the results (Fig. 13) shows that while the higher-precipitation percentiles are well captured using the hybrid gamma-GP distribution, the lower percentiles show somewhat worse results.This observation of poor performance for very low values also holds true for wind speed (not shown here).The lower values of the two variables, however, are very close to the precision of the observation (0.1 mm for precipitation and 0.1 m s −1 for wind speed).Very small precipitation amounts and low wind speeds are also less biophysically and ecologically important compared to the higher percentiles.We therefore consider the results of the evaluation largely acceptable.
In Table 4, we also compare the simulated versus the observed frequencies for very light rain (≤ 1 mm), light rain (1-10 mm), heavy rain (10-20 mm), and very heavy rain (> 20 mm).As we can see, our model underestimates the occurrence of very light rain events (28.6 % instead of 36.4 %) and overestimates the light rain events (58.3 % instead of 48.6 %) but generally performs much better than GCMs (Dai, 2006;Sun et al., 2006), especially when it comes to the heavy rain events.

Bias correction
After evaluating the results of GWGEN for wind speed for the different quantiles (see Sect.  tematic bias between the simulated and the observed values.This observation led us to adopt a further measure to improve the quality of the model output by implementing a quantilebased bias correction. We use an empirical distribution correction approach (quantile mapping) (Lafon et al., 2012) to a posteriori correct the simulated data.In the quantile evaluation (Sect.3.3), we saw that the simulated wind speed is a linear function of the observed wind speed, i.e., w sim = intercept + slope • w obs (best fit line in Fig. 12).Therefore, we use two steps here: one is for the difference between simulation and observation (ideally 0); the other one is the fraction of observation and simulation (ideally 1).The first one corresponds to the intercept with the y axis in Fig. 12, the second one to the slope of the best fit line.The analysis is based on every second percentile between 1 and 100 (i.e., 1, 3, 5, and so on) and mapped to its corresponding random number u ∈ R from a normal distribution as it is used for the cross correlation in the weather generator (Sect.3.2.6;x axis in Fig. 14 and Richardson, 1981).
Regarding the intercept (Fig. 14, left panel), we see that it strongly follows an exponential function given through f exp (u) = e au+b , a, b, u ∈ R. (21) The slope (Fig. 14, right panel), on the other hand, can be described by a simple third-order polynomial given by Hence, given the best fit lines in Fig. 14, the simulated wind speed is corrected via with a = 1.1582, b = −1.3359,c 0 = 0.9954, c 1 = 0.8508, c 2 = 0.0278, and c 3 = −0.0671.

Sensitivity analysis
The generalized pareto part of the hybrid gamma-GP distribution, which we used to simulate precipitation amount, has two parameters: the GP shape and the threshold parameter.
Unlike the gamma parameters, we were unable to relate these GP parameters to any of the monthly summary data we use as input to GWGEN.Hence, we decided to set fixed values for these parameters, and determine them through a sensitivity analysis.
To select the "best" values of the GP parameters, we compared simulated with observed precipitation amounts, running GWGEN with a wide range of realistic parameter values.To quantitatively assess the model performance, we used .Basis for the wind bias correction.For the left plot, each data point corresponds to the difference of a simulated percentile to the observed percentile.For the right plot (wind speed), each data point corresponds to the fraction of simulated to the observed wind speed for a given percentile.The random number on the x axis represents the residual value from a normal distribution centered at 0 with standard deviation of unity, as it is used in the cross-correlation approach (Richardson, 1981).
two metrics: (1) direct comparison of the quantiles (see previous section) and (2) a Kolmogorov-Smirnov (KS) test that evaluates whether two data samples come from significantly different distributions.Our criteria were 1. the R 2 correlation coefficient between simulated and observed quantiles; 2. the fraction simulated precipitation observed precipitation from the slopes in Fig. 13 and its deviation from unity; 3. the fraction of simulated (station-specific) years that are significantly different (KS test) from the observation; and 4. the mean of the above values.
We tried two different approaches to select the gamma-GP crossover threshold: first, we tried a fixed crossover point; second, we used a quantile-based crossover point.For the latter, the model chooses to use the GP distribution if the quantile of the random number drawn from the gamma distribution is above a certain quantile threshold.This introduces a flexible crossover point in our hybrid distribution which, however, did not improve the results significantly.We therefore show here only the results using the fixed crossover point.
The values of the crossover point for our sensitivity analysis were 2, 2.5, 3, 4, and from 5 to 20 in steps of 2.5, and 20 to 100 in steps of 5. Furthermore, we varied the GP shape parameter from 0.1 to 3 in steps of 0.1 (810 experiments in total).The results of this sensitivity analysis are shown in the Supplement (Fig. A1).
In general, we found that the three criteria (1-3) could not be optimized all together at the same time.The R 2 is best for high thresholds and low GP shape parameters, the slope is best for low to intermediate thresholds and a low GP shape, and the KS statistic is best for low threshold and intermediate GP shape parameters.
However, R 2 did not vary that much (from 0.68 to 0.74), and from a visual evaluation of the corresponding quantile plots we saw that the higher quantiles (> 90) were much better represented for a better KS result.Hence, we chose to follow the KS test criteria, which is also the strictest of our evaluation methods but again compared the different quantile plots to get good results for the higher quantiles.Finally, we chose a threshold of 5 mm and a GP shape parameter of 1.5.For this setting, 81.7 % of the simulated years do not show a significant difference compared to the observation, the mean R 2 of the plots in Fig. 13 is 0.81, and the mean deviation of the slope from unity is 0.10 and for the upper quantiles (90 to 100) it is 0.017.
Nevertheless, in total, the results seem to be fairly independent of the two parameters since even the amount of years without significant differences varies from 73 % to only 83 %.It is, however, better than the gamma distribution alone which still has 78.6 % of station years not differing significantly but with a slope deviation from unity for the upper quantiles of 0.16.Thus, using the hybrid gamma-GP distribution improves the simulation of high-amount precipitation events by roughly a factor of 10 compared to a standard gamma approach.

Limitations
As demonstrated above, GWGEN successfully downscales monthly to daily meteorology with good correlation and www.geosci-model-dev.net/10/3771/2017/Geosci.Model Dev., 10, 3771-3791, 2017 low bias when compared to observations.However, there are a few limitations of the model as currently described that should be noted.Importantly, this version of GWGEN neither downscales all conceivable meteorological variables, nor does it provide a mechanism for generating daily meteorological time series across multiple points that are spatially autocorrelated.Concerning the former point, while GWGEN simulates daily precipitation, temperature, cloud cover, and wind speed, it does not currently handle other variables that might be important in land surface modeling, such as humidity or wind direction.On the latter point, the lack of explicit simulation of spatial autocorrelation may make GWGEN unsuitable for certain applications, e.g., regional high-resolution hydrological modeling in small catchments (<∼ 2500 km 2 ), where having the capability to simulate flood and other extremes is important.This is because the weather generator could, e.g., simulate rainfall on different days in different parts of the catchment, where in reality storm events would be highly autocorrelated in space and controlled by mesoscale meteorological conditions.
5 Discussion and outlook GWGEN successfully downscales monthly to daily meteorology for any point on the globe, in any climate, in any season, and in any time in recent Earth history and in the near future (e.g., next century).It extends the original Richardsontype weather generators to simulate wind speed along with precipitation, temperature, and cloud cover.The model requires only monthly values of the meteorological variables to be downscaled and does not rely on any other spatial information, e.g., whether or not the location is in the tropics.
In general, the results of our downscaled meteorology are excellent, with all simulated variables showing both very high correlation and limited bias when compared to observations.We improved the simulation of daily precipitation amount by replacing the gamma distribution used in the original Richardson-type weather generators with a hybrid gamma-GP distribution, which results in the improved simulation of heavy precipitation events.The GP distribution is based upon a globally fixed shape and location parameter, which may be an oversimplification, but is still 10 times more accurate than traditional methods that used gamma alone.Our extensive sensitivity analysis to determine the best coefficients for the shape and location parameters of the GP distribution suggest that further improvements might come through correlating the GP parameters to geographic region and/or seasonality (Maraun et al., 2009;Rust et al., 2009) or by introducing a dynamical location parameter (Frigessi et al., 2002).Finally, we introduced a step to correct for systematic bias in the downscaling of temperature and wind speed.
Despite the limitations noted above, GWGEN will be useful in a wide range of applications, from global vegeta-tion and crop modeling to large-scale hydrologic analyses, to understanding animal behavior, to forecasting of fire, insect outbreaks, and other ecosystem disturbances.GWGEN may even be envisaged as a potential replacement for very large and cumbersome gridded datasets of high-temporalresolution meteorology such as CRUNCEP (Viovy and Ciais, 2016), especially for models that use meteorological forcing at a daily time step.The weather generator is particularly suited for the incorporation into models that run on a spatial grid; for example, GWGEN can readily be incorporated into existing DGVMs such as LPJ-LMfire (Pfeiffer et al., 2013) or LPJ-ML (Bondeau et al., 2007) that already rely on a weather generator to provide daily meteorology for certain processes.
While GWGEN does not handle spatial autocorrelation, in most DGVMs there is no lateral connection between grid cells, and therefore an explicit representation of spatial autocorrelation in the driving daily meteorological data would have no effect on the model output.We further note that if the monthly data used to drive the model are spatially autocorrelated (this would be the case when using gridded climatology, for example) then the result of the weather generator will also preserve this autocorrelation, at least when the model results are analyzed on monthly or longer timescales.
The limitations present in this version of GWGEN could be addressed in future versions.Methods for simultaneous multi-site weather generation exist (Wilks, 1998(Wilks, , 1999b, c) , c) and could be adapted to GWGEN.However, even simpler methods to approximate spatial autocorrelation could be possible.Running GWGEN with gridded monthly meteorology (this is the primary application we foresee for the current version of the model) means that the input variables are already highly correlated in space, i.e., the monthly climate in one grid cell generally closely resembles neighboring cells outside of complex terrain containing sharp, monotonic climate gradients, e.g., rain shadows.Thus, one simple way of achieving a measure of spatial autocorrelation in GWGEN would be to impose a spatial autocorrelation field on the sequence of random numbers used to impose stochastic noise in the downscaling functions.If the random number sequence is similar between grid cells, then, e.g., rain is likely to fall on the same day, given that the transition probabilities will likely also be similar.Over moderate distances (< 50 km), it might even be sufficient to use the same random seed across all grid cells in a neighborhood.This would have the effect of producing strongly autocorrelated daily meteorology in space, with the only variations being imposed by the underlying input monthly climatology.
Furthermore, it would be straightforward to include additional meteorological variables in the model framework, handling, e.g., humidity in the same way that temperatures, cloud cover, and wind speed are disaggregated.Other variables, such as pressure and wind direction, might be more difficult using the basic GWGEN structure because of the importance of autocorrelation, particularly at high spatial resolution, and might benefit from a different approach towards weather generation.Finally, GWGEN only downscales meteorology from monthly to daily values; for models that require an even shorter time step, e.g., 6 hourly, some extension of the model functionality would be required.For certain variables, e.g., temperatures, subdaily downscaling could be easily implemented (Cesaraccio et al., 2001); for other variables, such as precipitation, a large literature on downscaling methods exists (e.g., Bennett et al., 2016), and global datasets of hourly meteorology for model calibration are available (e.g., the Integrated Surface Database; Smith et al., 2011).

Conclusions
Compiling a global database of daily precipitation, temperature, cloud cover, and wind speed measurements, we explored the relationship between daily meteorology and monthly summaries first described in the context of weather downscaling by Geng and Auburn (1987).Our analysis of more than 50 million individual records showed that daily to monthly relationships are relatively stable in space and time, and constant across a very wide range of stations from all latitudes and climate zones.With the resulting relationships, we parameterized a WGEN/SIMMETEO-type weather generator, with the intention of creating a generic scheme that could be applied anywhere over the Earth's land surface for the past, present, and (near) future.

Figure 1 .
Figure1.Schematic workflow of GWGEN.After smoothing the monthly input, the Markov chain is used to decide whether it is a dry or a wet day.If it is a wet day, we draw a random number from the gamma-GP (generalized Pareto) distribution.Furthermore, the other means of the variables (T min/max , c, w) are adjusted and their daily values are calculated using the estimated standard deviations and residuals.The wind speed furthermore undergoes a square root transformation before applying the cross correlation and in the end is corrected using the bias correction.A quality check in the end restricts our model to be within a 5 % range of the observed total precipitation and to replicate the number of wet days from the input.

Figure 2 .
Figure 2. Weather stations used for parameterization and evaluation of the weather generator.The uppermost panel shows the locations of the stations used for parameterizing precipitation and temperature; the middle panel shows the stations for cloud fraction and wind speed, as well as for calculating the cross correlations between temperature, cloud fraction, and wind speed.The lower plot shows the location of the stations used to evaluate the model, which were excluded from the parameterization stations.

Figure 3 .
Figure 3. Transition probabilities vs. wet fraction.The red density plot in the background shows the density of the observations, and the blue lines indicate the linear regression line of the probability against the wet fraction.The fit for the p 11 transition probability was forced to the point (1, 1); the others were forced to (0, 0).The underlying data for the fits correspond to the means of the multi-year series for each month for each station.

x
Figure 4. Mean precipitation-gamma scale relationship.The blue line represents the best fit line of the mean precipitation on wet days to the estimated gamma scale parameter of the corresponding distribution.Each data point corresponds to one multi-year series of 1 month for one station.
Furthermore, to account for the sparse data below −40 • C and above 25 • C for minimum temperature (or −30 and 35 • C for maximum temperature), we use an extrapolation for the extremes as indicated by the blue and violet lines in Fig. 7.The formulae for the standard deviations σ of minimum and maximum temperature are therefore a combination of four polynomials: σ Tmin, wet/dry = min, wet/dry , for T min, wet/dry ≤ −40 • C p 5 T min, wet/dry , for − 40 • C < T min, wet/dry ≤ 0 • C p 5 T min, wet/dry , for 0 • C < T min, wet/dry ≤ 25 • C p 1 T min, wet/dry , for 25 • C < T min, wet/dry σ Tmax, wet/dry = max, wet/dry , for T max, wet/dry ≤ −30 • C p 5 T max, wet/dry , for − 30 • C < T max, wet/dry ≤ 0 • C p 5 T max, wet/dry , for 0 • C < T max, wet/dry ≤ 35 • C p 1 T max, wet/dry , for 35 • C < T max, wet/dry .
Figure5.Correlation of minimum temperature on wet and dry days to the monthly mean.The y axes show the mean minimum temperature on wet or dry days, respectively; the blue line corresponds to the best fit line.Parameters of the fits are also shown in Table1.

Figure 6 .Figure 7 .
Figure6.Correlation of maximum temperature on wet and dry days to the monthly mean.The y axes show the mean maximum temperature on wet or dry days, respectively; the blue line corresponds to the best fit line.Parameters of the fits are also shown in Table1.

Figure 8 .Figure 9 .
Figure 8. Correlation of cloud fraction on wet and dry days to the monthly mean.The y axes show the mean cloud fraction on wet or dry days, respectively; the blue line corresponds to the best fit line.Parameters of the fits are also shown in Table3.
σ w, wet (w wet ) = c 1, w, wet w wet + c 2, w, wet w 2 wet + c 3, w, wet w 3 wet σ w, dry w dry = c 1, w, dry w dry + c 2, w, dry w 2 dry + c 3, w, dry w 3 dry .(14) This better resolves the complex behavior close to 0 m s −1 compared to a linear fit.The plots are shown in Figs. 10 and 11 and the parameters for the fits are shown in Table

Figure 10 .
Figure10.Correlation of wind speed on wet and dry days to the monthly mean.The y axes show the mean cloud fraction on wet or dry days, respectively; the blue line corresponds to the best fit line.Parameters of the fits are also shown in Table1.

Figure 11 .
Figure 11.Correlation of standard deviation of wind speed on wet and dry days to the corresponding monthly mean.The y axes show the standard deviation; the x axes the mean on wet or dry days, respectively.The blue line corresponds to the best fit line; a third-order polynomial corresponds to the underlying red density plot.The black bars have a width of 0.1 m s −1 , the accuracy of the input data, and indicate the mean standard deviations for the given interval range.Parameters of the fits are also shown in Table1.

Figure 12
Figure12.Q-Q plots for all variables with all quantiles(1, 5, 10, 25, 50, 75, 90, 95, and 99)  for µ = 5.0 mm and ξ = 1.5.The blue lines indicate linear regression from simulation to observation.The red line shows the ideal fit (the identity line).Blue shaded areas represent the 95 % confidence interval.The plots compare the simulated quantile from the list above one year of one station to the corresponding observed quantile of the same year and station.The plot for wind speed used the bias correction from Sect.3.4.

Figure 13
Figure 13.Q-Q plot for different quantiles for precipitation for µ = 5.0 mm and ξ = 1.5.The blue lines indicate linear regression from simulation to observation.The red line shows the ideal fit (the identity line).Blue shaded areas represent the 95 % confidence interval.The plots compare the simulated quantile of one year of one station to the corresponding observed quantile of the same year and station.
Figure14.Basis for the wind bias correction.For the left plot, each data point corresponds to the difference of a simulated percentile to the observed percentile.For the right plot (wind speed), each data point corresponds to the fraction of simulated to the observed wind speed for a given percentile.The random number on the x axis represents the residual value from a normal distribution centered at 0 with standard deviation of unity, as it is used in the cross-correlation approach(Richardson, 1981).

Table 2 .
Fit results of the correlation of temperature standard deviation with the corresponding mean on wet/dry days for Fig.7.The underlying equations are shown in Eq. (11).

Table 3 .
Fit results of cloud correlation for wet and dry days for Fig. 8. SE indicates standard error.

Table 4 .
Simulated and observed precipitation frequencies for certain ranges.The frequency is defined as the number of precipitation occurrences in the specified range divided by the total number of precipitation occurrences.