Model Consistent Pseudo-Observations of Precipitation and Their Use for Bias Correcting Regional Climate Models

Lack of suitable observational data makes bias correction of high space and time resolution regional climate models (RCM) problematic. We present a method to construct pseudo-observational precipitation data by merging a large scale constrained RCM reanalysis downscaling simulation with coarse time and space resolution observations. The large scale constraint synchronizes the inner domain solution to the driving reanalysis model, such that the simulated weather is similar to observations on a monthly time scale. Monthly biases for each single month are corrected to the corresponding month of the observational data, and applied to the finer temporal resolution of the RCM. A low-pass filter is applied to the correction factors to retain the small spatial scale information of the RCM. The method is applied to a 12.5 km RCM simulation and proven successful in producing a reliable pseudo-observational data set. Furthermore, the constructed data set is applied as reference in a quantile mapping bias correction, and is proven skillful in retaining small scale information of the RCM, while still correcting the large scale spatial bias. The proposed method allows bias correction of high resolution model simulations without changing the fine scale spatial features, i.e., retaining the very information required by many impact models.


Introduction
Evaluation of climate models are best performed with gridded observational data sets as these better agree with the spatial average data produced by the models [1].However, the observational data sets are limited in their space and time resolution due to the underlying station density.Furthermore, as the precipitation field becomes increasingly heterogeneous with higher resolution, the assumptions and statistical methods of the gridding procedure have a larger impact on the results.With increasing space and time resolution of model outputs, this poses a problem for evaluation as well as bias correction applications for regional climate models (RCMs).
Bias correction of RCM data has over the last decade evolved from simple scaling to get the mean climatology corrected (e.g., [2]), to correcting several or all moments of the distribution (e.g., [3]).This increase in detail put higher demands on the quality of the observational data to correctly represent the distribution, and also regarding the length of the data set to have robust statistics for higher moments [4].If high enough time resolution data exists, it can be tempting to simply interpolate a spatially coarser observational data set to the high resolution of the model in order to perform the bias correction, however, commonly used bias correction methods will then remove fine scale spatial information from the RCM variables.
In large parts of the world, only coarse space and time resolution data are available, and only corrections of the mean are supported.Regionally, relatively high resolution precipitation data sets are available, e.g., the E-OBS data set over Europe (25 or 50 km and daily time step) [1], and locally even higher resolutions, e.g., PTHBV [5] over Sweden (about 4 km and daily resolution).Even in this case, the increasing spatial scales of RCMs now approaching a standard resolution of about 12.5 km for Europe are less and less supported by observational data, especially with an increased interest in sub-daily time steps.Furthermore, common for all these data sets is that the station density is not sufficient to alone describe the spatial details, and statistical methods and assumptions are applied to describe regional features of the precipitation field.Therefore, gridded observational data are increasingly becoming a product of the statistical model at finer resolutions.
Weedon et al. [6] presented a simple method to construct a merged data set of model reanalysis and observational precipitation data, known as the WATCH forcing data (WFD).They started from the assumption that the assimilation procedure of the reanalysis makes the simulated weather at single grid points synchronized with observations at a monthly time scale, such that similar events occur within that timeframe.Thus, the model simulated precipitation can be corrected to have the same monthly mean precipitation as the observed.Weedon et al. [6] calculate a simple scaling factor from the ratio of the model and observational data for each single month to correct the reanalysis.The scaling factor is then applied at the higher time output frequency for each month.The method allows the model to determine temporal scales below the monthly timescale, but is constrained to observations beyond that.This demands confidence in the model performance at the shorter timescales, thus the method is better applicable to climatological assessments, as opposed to studies of single events.
The main advantage of the WFD is that it provides a worldwide data set at high temporal resolution, which allows for coarser scale impact modeling also in regions of the world with sparse and coarse data availability.In this spirit, we have extended the Weedon et al. [6] method to a constrained RCM simulation that is corrected towards observations.To retain fine scale spatial details of the model, a spatial filter is applied, so that corrections are only performed at coarser spatial scales.This new high resolution pseudo-observational data set is evaluated and applied as reference data for bias correction of downscaling simulations with the same RCM.The paper outline is as follows: first a description of the study region (Section 2), a presentation of the model and data used (Section 3), a methods section where the creation of the pseudo-observations as well as the bias correction method are presented (Section 4), followed by the results (Section 5).The paper ends with a discussion (Section 6) and conclusions (Section 7).

Study Region
Initial setup and evaluation of the method is performed for the region of Sweden in northern Europe, see Figure 1.The main reason for restricting the analysis to this sub-domain of the model domain (see Section 3), is the availability of a high resolution gridded data set.Furthermore, this is a suitable region for initial setup and evaluation of the proposed method due to varying topographic conditions, with the Scandic mountain range along the western boundary to Norway and plains in the south (Figure 1), and conditions with dominantly large scale precipitation systems in winter, and a fair amount of convective small scale events in summer.The annual mean rainfall is highest along the Scandic mountain range, and in the inland regions from the south-western coastline (Figure 1).The rest of the country receives amounts around 600-900 mm•yr −1 , with local small scale variations.

Model and Data
Simulations over the Euro-CORDEX domain covering most of Europe at 0.11 degree (about 12.5 km) with the RCM RCA4 (Rossby Centre Atmosphere RCM version 4) are used for developing the proposed method.Details on an earlier version of RCA (version 3) are presented in Samuelsson et al. [7], and new developments consist mostly of re-coding and updating with respect to surface processes.In the standard setup, the model is run with prescribed sea surface temperatures, lateral boundary conditions from a driving model, and greenhouse gas concentrations.
The standard method for nesting RCA in a global climate model (GCM) is by relaxing the atmospheric variables of RCA at model levels to that of the GCM in a ten point relaxation zone at the lateral boundaries.An alternative strategy for imposing the lateral boundary conditions is to use spectral nudging, which has recently been introduced in the model as described in Berg et al. [8].When active, it imposes large scale features of the driving model on the interior domain of the RCM.This is carried out through a Fourier decomposition of 2-d fields of temperature and horizontal winds at model levels above the atmospheric boundary layer, and wavelengths above 800 km are constrained to the driving model.However, only a fraction of the driving model field is imposed on the RCM, and this fraction increases linearly with the model levels vertically from 0 at around 850 hPa to 10% at the top level at around 0.1 hPa, corresponding to e-folding times of about 10 h to 40 days [9].More details on the spectral nudging and the configuration for RCA4 can be found in Berg et al. [8].With a constrained circulation, the probability for the RCM to simulate observed precipitation time series at a given location increases, but it is also dependent on the model to actually produce a precipitation event at the correct time and place.Using a smaller domain imposes similar constraints on the circulation, which might be beneficial in some cases.
A set of experiments is carried out (see Table 1).First, a set of two downscaling simulations are carried out with ERA-Interim [10] reanalysis data at the boundaries; one (REI) with standard lateral boundary forcing, and the other with spectral nudging (REISN).For evaluation of the impact on bias correction, we make use of a control simulation with the EC-Earth GCM [11][12][13] from the CMIP5 ensemble applied at the boundaries (ECE).For construction of the pseudo-observational data, we make use of the E-OBS (v9.0) gridded precipitation data set [1], which includes around 70 to 100 gauges in Sweden, depending on the time period.It is available at both a 50 and 25 km resolution, but here we only use the former to better explore the possibilities of the proposed method.
For evaluation of the model at higher resolution (4 km), we use the PTHBV gridded data set [5].It is based on around 700 gauges.This means an average density of one station per 25×25 km 2 , however, the density is significantly lower in the mountainous regions in the north-west.Clearly, the high resolution of PTHBV is a product of the statistical interpolation model.
The observational data sets are remapped by bi-linear interpolation to the model 12.5 km grid for all purposes in this study.The difference in resolution between the original and target grids is fairly small, and the bi-linear approximation has therefore no significant impact on the results.For the coarser resolution E-OBS data set, there is a problem with missing data along the coastlines with the bi-linear interpolation method.The nearest neighbor with data was therefore replicated in order to fill out the data set.

Pseudo-Observational Data
The Weedon et al. [6] method requires the model to be synchronized with the observations such that the weather during each single month is similar, here in the sense of similar precipitation systems affecting the local regions, in the model and the observation.Reanalysis data are, through the assimilation of observational data, largely synchronized to the atmospheric circulation.However, here we are interested in much higher resolutions than current day reanalysis products provide, and therefore apply RCM downscaling.In practice, most RCM downscalings of reanalysis data produce internal variability and circulation biases for larger domains such that correlations between the RCM and the driving reanalysis data drops significantly.Assimilation is a time consuming task, and different simpler methods have been applied previously, such as re-initializations [14] and different versions of spectral nudging [9].Here we apply spectral nudging following Berg et al. [8] to constrain the large scales.
To correct the spectrally nudged downscaling simulation towards observations the number of dry days (days with less than 1 mm•day −1 of precipitation) are first corrected by setting excessive wet days with the least amount of precipitation in the model to zero.Then a simple scaling factor [6] is applied such that the pseudo-observational data, P P SEU DO (n, m), becomes . (1) MOD (OBS) denotes the model (observational) data, n and m are indices of the day in the month and the month, respectively.N(m) is the number of days in month m.This simple scaling will affect the higher intensities within a month more than the lower ones.A scaling taking also the variance into account could improve this, but when only monthly data are available (which is the premise investigated here), we are limited to scaling the mean value only.The scaling factor, i.e., the ratio in Equation (1), will impose the spatial variability of the observational data at the monthly timescale to the pseudo-observational data.In a gridded data set, this variability arises from statistical assumptions about precipitation amounts related to orography, and to the often low density gauge data, with station location often being biased to e.g., valleys.The dynamical model is more advanced in the simulation of spatial characteristics of the precipitation field, and is consistent, e.g., when it accounts for different flow directions which can affect orographic effects on precipitation amounts and distribution.Thus, the statistical and dynamical methods might differ, especially for smaller scale features.It is subjective whether one puts more trust in the statistical or the dynamical method as they are problematic to evaluate.However, best practice is to correct the model only when substantial information shows that the model is erroneous.Since generally no such substantial information exists for the finer spatial details, best practice suggests not to affect the fine scales of the model.Still, the coarser scales might need to be corrected.Thus, the pseudo-observational data should adapt large scale structures from the observations, but retain the small scale structures from the model.This can be performed by applying a spatial filter to remove small scale spatial features from the data sets prior to calculating the scaling factors in Equation ( 1).The smoothed scaling factors will then allow the small scale features of the model to be retained in the final product.
The spatial filter used here is a nine point spatial average which is applied to the model data for the monthly time step, such that the actual resolution of the model and the observations become more similar.A sensitivity study with a 16 point average was performed, however, with little actual differences in the results.
The pseudo-observational data (PSOBS) are constructed from REISN combined with E-OBS, as outlined above.

Distribution Based Scaling (DBS)
The DBS method [15] is a version of the commonly used quantile mapping bias correction methods in which meteorological variables are fitted to appropriate parametric distributions.For precipitation, DBS separates between wet and dry days and is applied in two steps: (1) remove drizzle generated by GCMs/RCMs to correct percentages of wet days and (2) transform the remaining precipitation to match the observed frequency distributions.In (1) a wet-day threshold for the model is identified by comparing the simulated daily precipitation to the observed daily precipitation sorted in an ascending order.Model values below the threshold are removed.In (2), for each data set two gamma distributions are fitted to the cumulative distributions of precipitation intensities below (above) the 95th percentile, and combined into one distribution.It allows a good fit to both low and high intensities.Correction factors are calculated separately for each season DJF (December-February), MAM (March-May), JJA (June-August) and SON (September-November), as the distributions can have significantly different shapes depending on the physical characteristics of the precipitation processes dominating the season.In this study, we are focusing on the effects of the reference data on the corrections.To emphasize these effects, we perform the analysis on the DBS calibration period, i.e., the full period, such that sporadic elements of the corrections on a validation period do not distract from the spatial features we are interested in.

Synchronizing Weather Events Using Spectral Nudging
A requirement of the proposed method is that the RCM can reproduce observed precipitation events from observations at the correct time and place, at least within a monthly time scale.Grid point time correlations between the model and the E-OBS data are calculated to investigate how well the model reproduces observations.
Starting with correlations at the monthly time scale, the standard model setup (REI), with only lateral boundary forcing, produces correlations of around 0.5 on average over Sweden for the different seasons, with higher values in winter than summer (Table 2), due to the type of precipitation systems present.With spectral nudging (REISN), the correlations increase significantly to around 0.8.These are high correlations, especially considering the highly complex atmospheric interactions to produce precipitation, which lends confidence in pursuing the construction of the pseudo-observations (PSOBS).Because PSOBS are corrected at the monthly time scale, its correlations are per definition close to perfect; with small deviations only due to the nine point spatial filter.Regionally, there are only small differences in the results, and the main features consist of generally lower correlations in southern Sweden in summer, which is probably due to stronger convective activity in this region.At the daily time step, the correlations drop to well below 0.5 for REI, but remain above that value for REISN (Table 2).The PSOBS only slightly increase the values of REISN, which is because it retains the time series of the model at sub-monthly time scales, with only corrected magnitudes.There are multiple reasons beyond model errors for the lower correlations at the daily timescale, e.g., we do not take into account smaller location errors, nor have we adjusted for the different timeframe for accumulating precipitation of the model (00:00-00:00) and the observations (06:00-06:00).

Evaluation
As presented in Section 5.1, DJF and JJA are the seasons with the largest contrasts and we therefore focus the presentations and evaluation to these seasons.MAM and SON were also investigated, but add negligible additional information to the analysis.
Figure 2 shows the mean precipitation for DJF and JJA for the two observational data sets as well as REISN and PSOBS.There are clear differences between PTHBV and E-OBS, mainly regarding the fine scale spatial features which are lacking in E-OBS due to the lower resolution, but also on larger scales with E-OBS having much lower amounts over the Scandic mountains to the north-west.REISN is, in general, in close agreement with PTHBV in winter, but overestimates amounts somewhat compared to E-OBS.In summer, REISN is too wet, especially compared to E-OBS.REISN and PTHBV seem to have similar fine scale features, however, a closer look in Section 5.3 will reveal some significant differences.PSOBS has corrected large scale bias to be close to E-OBS, but retains the smaller scale spatial features of REISN.
The probability distribution function (PDF) of daily precipitation intensities calculated from values throughout the domain is also slightly different between the two observational data sets (Figure 3).Clearly, E-OBS is underestimating intensities above about 20 mm•day −1 compared to PTHBV, especially in winter, but also has a higher dry day probability in both summer and winter.The former is probably to a large extent due to the coarser resolution of E-OBS which shifts the PDF towards lower intensities.The dry day (days with <1 mm of precipitation) differences could be related to the statistical method used.REISN has a wet bias in that it underestimates the number of dry days by a few percentage units, and overestimates intensities above around 20 mm•day −1 , for both winter and summer and compared to both PTHBV and E-OBS.Although PSOBS is only corrected for monthly means, the PDF is much closer to its E-OBS reference.This is due to a fortunate coincidence: An overestimation of the monthly values will give a scaling factor below 1, with the result that higher intensities are scaled down more than lower intensities, thus reducing the high intensity bias of the model.The larger wet bias in summer actually causes a reduction of the higher intensities below those of E-OBS.Although the dry days of REISN are corrected in the method to produce PSOBS, there are still smaller deviations from E-OBS.This is because the subsequent scaling of the data can shift precipitation amounts below the 1 mm•day −1 limit which defines the dry day threshold.
Differences in variance between model and observational data sets have been shown to affect climate change signals when bias correction methods include higher moments than the mean, such as quantile mapping does [4,16,17].The method to create PSOBS enforces the same standard deviation at monthly time scales, for which reason we instead investigate the standard deviation at intra-monthly timescales, i.e., the monthly mean was first subtracted from each daily value to produce monthly anomalies (Figure 4).PTHBV has higher standard deviation than E-OBS, which again is likely due to the difference in original resolution.REISN has on average a little higher values than E-OBS and lower than PTHBV in winter.However, in summer, REISN clearly overestimates the standard deviation by over 20% in a spatial average.PSOBS has similar standard deviation as E-OBS, although on average a small underestimation is found.Thus, the larger standard deviations of REISN is reduced by the method to produce PSOBS.This can be explained by the wet bias of the model, and the related stronger reduction in higher intensities (see Figure 3).This highlights the need for a good performance of the model to create the pseudo-observations from.

Bias Correction
An RCA4 downscaling simulation (ECE) driven by EC-Earth at the lateral boundaries, i.e., without spectral nudging, is bias corrected using the DBS method.Two different bias corrections are performed: one with E-OBS and the other with PSOBS as reference data.
ECE has a wet bias compared to E-OBS for both winter and summer (Figure 5).The smaller scale spatial details are clearly similar to REISN (Figure 2c,g), which is because the fine scale details are determined by the RCM rather than the driving model.Bias correction of ECE with E-OBS as reference (ECE∼E-OBS) corrects the large scale features, but also clearly shows the smoothing effect, which removes all small scales feature, from using a coarser reference data set (Figure 5).With PSOBS as reference (ECE∼PSOBS), the large scale features are corrected, however, the smaller scale information from the RCM are retained.
To better investigate that smaller scale information, a nine point spatial filter is applied to each of E-OBS, ECE, ECE∼E-OBS and ECE∼PSOBS, and subtracted from the respective original data.This removes larger scale differences between the data sets, and the small scale information emerges (Figure 6 as example for summer).Visual inspection of the different data sets indicates that PSOBS, ECE and ECE∼PSOBS have similar features, whereas PTHBV differs somewhat.ECE∼E-OBS lacks features at this scale, due to the coarser original resolution of E-OBS.To better quantify the differences, a Taylor diagram [18] is produced (Figure 6, bottom right), in which the spatial correlations and standard deviation are directly compared with PSOBS as reference data set.PTHBV clearly differs with correlations below 0.3 and lower spatial variance.In winter (not shown) the correlations are even slightly negative.ECE correlates with above 0.9 with PSOBS, but has larger standard deviation due to the wet bias.ECE∼PSOBS is close to perfectly correlated with PSOBS, i.e., the small scale information of the two fields is almost identical.ECE∼E-OBS has close to zero correlation and standard deviation, due to the complete lack of small scale features.

Discussion
The confidence in gridded observational data sets decrease with increasing spatial resolution, as the data becomes more and more a product of its statistical assumptions.One might therefore have more trust in a dynamical model, which is if not correct then at least consistent between different synoptic situations and different variables.
The pseudo-observations opens up possibilities for correction of RCM biases at higher time resolutions in a model consistent way.However, it is important to verify the RCM's performance at higher temporal resolutions as the pseudo-observations will only correct the magnitudes, and not the time line of events.This includes errors in the timing of the diurnal cycle of precipitation, and the method is generally not sufficiently detailed for studies of single events in the observed records.
The main limitation of the pseudo-observations is the reference RCM simulation which needs to be sufficiently close to observations regarding monthly means.If there are large biases, the scaling factor will significantly affect the intensity distribution as larger intensities are affected more by the scaling than lower ones.Also the temporal variance of the model will be affected.The reason for the bias needs to be investigated for each case.
When observations are too sparse and unreliable, or too coarse for the intended applications, the proposed method will create both a reference data set for impact modeling, and a reference for correcting RCM simulations in a consistent way.The method was here set up and investigated in a region with fairly good observational coverage at the daily time scale for reasons of evaluating the model, however, in such a region a more efficient way to retain the fine spatial scales would be to first separate the spatial scales and correct only the coarser scales of the model to the observational data, and then adding the un-corrected finer scales after the bias correction step.

Conclusions
A pseudo-observational data set was constructed from an RCM simulation with constrained large scale circulation to a driving reanalysis model.A scaling factor was then applied to the simulation to correct the monthly mean for each single month towards an observational data set, however, first introducing a spatial filter which makes sure that the fine spatial details of the model is retained after the correction.The pseudo-observations are shown to be consistent with the RCM at finer spatial and temporal scales, but shares large scale information as well as the intra-and inter-annual features of the observational data set.
The proposed method allows the user to decide whether the finer scale information should be determined by the model or the observations.For bias correction applications, this can reduce unnecessary corrections, and retain more of the model's own detailed features.This is especially interesting for higher resolution model applications, where no high resolution observations are available, or in regions where there is little trust in the observations.

Figure 1 .
Figure 1.Topographical map of Sweden and its near surroundings (left).Mean annual precipitation for the period 1980-2009 from the PTHBV data set at the RCA model resolution of about 12.5 km (right).

Figure 2 .
Figure 2. (a-d) Winter and (e-h) Summer Seasonal mean precipitation for (a,e) PTHBV, (b,f) E-OBS, (c,g) REISN and (d,h) PSOBS.The domain average value is given in the lower right box in each panel.

Figure 3 .
Figure 3. PDF of wet day precipitation intensities for PTHBV, E-OBS, REISN and PSOBS (symbols according to legend), for DJF (top) and JJA (bottom).The dry day probability (%) is given in the legend.

Figure 5 .
Figure 5. (a-d) Winter and (e-h) Summer Seasonal mean precipitation for (a,e) E-OBS, (b,f) ECE, (c,g) ECE∼E-OBS and (d,h) ECE∼PSOBS.The domain average value is given in the lower right box in each panel.

Figure 6 .
Figure 6.Small scale features of (a) PTHBV, (b) PSOBS, (c) ECE , (d) ECE∼E-OBS and (e) ECE∼PSOBS for summer, derived by subtracting a 9p spatial filtered field from each map of mean values.(f) A Taylor diagram shows the spatial correlation and standard deviation for each of the panels (a-e) with PSOBS as reference.

Table 1 .
List of model simulations, observational and merged data.The indicated time period is that used for the analysis, not that of the complete data set.

Table 2 .
Domain average results for correlations of precipitation towards E-OBS daily and monthly data.