Assessing atmospheric bias correction for dynamical consistency using potential vorticity

Correcting biases in atmospheric variables prior to impact studies or dynamical downscaling can lead to new biases as dynamical consistency between the ‘corrected’ fields is not maintained. Use of these bias corrected fields for subsequent impact studies and dynamical downscaling provides input conditions that do not appropriately represent intervariable relationships in atmospheric fields. Here we investigate the consequences of the lack of dynamical consistency in bias correction using a measure of model consistency—the potential vorticity (PV). This paper presents an assessment of the biases present in PV using two alternative correction techniques—an approach where bias correction is performed individually on each atmospheric variable, thereby ignoring the physical relationships that exists between the multiple variables that are corrected, and a second approach where bias correction is performed directly on the PV field, thereby keeping the system dynamically coherent throughout the correction process. In this paper we show that bias correcting variables independently results in increased errors above the tropopause in the mean and standard deviation of the PV field, which are improved when using the alternative proposed. Furthermore, patterns of spatial variability are improved over nearly all vertical levels when applying the alternative approach. Results point to a need for a dynamically consistent atmospheric bias correction technique which results in fields that can be used as dynamically consistent lateral boundaries in follow-up downscaling applications.


Introduction
The biases present in global climate model (GCM) output (for example see Rocheta et al 2014) require correction to improve the representation of variables before they are used in impact studies (Piani et al 2010, Teutschbein and Seibert 2012, Muerth et al 2013, Wilcke et al 2013. Recent bias correction of GCM output used to provide a regional climate model (RCM) with boundary conditions have been shown to improve the RCM simulation (Colette et al 2012, Bruyère et al 2014, White and Toumi 2013. Common practice involves bias correcting selected atmospheric variables separately without consideration of the final dynamical balance of those variables. These individual corrections lead to a modification in the spatiotemporal covariance structure in the climate variables, possibly creating a breakdown in the physical conservation principles (Liu et al 2013, Muerth et al 2013, Yu and Wang 2014 which might degrade the validity of downscaling and impact studies. The studies mentioned above have either ignored this effect or assumed the effect small, however White and Toumi (2013) who were interested in dynamical consistency, favoured the use of mean bias correction over quantile-quantile (QQ) mapping techniques as it at least maintains first-order spatial and intervariable dependencies while QQ does not.
A key strength of GCM and RCM outputs are that they ensure a close approximation of the dynamical consistency between variables; they produce auto-, spatial-, and intervariable correlations that are consistent with the model's atmospheric physics. This physical atmospheric state is important in driving lateral boundary conditions in dynamical downscaling or when using multiple variables simultaneously in an impact study. The focus of this paper is to quantify the effect on dynamical consistency of applying a standard bias correction technique, thus testing the assumption in previous studies that such an effect would be small. This approach may lead to bias corrected results which are more appropriate to use in dynamically sensitive applications.
The approach used here quantifies dynamical consistency though transforming the desired meteorological information into a single scalar which can be bias corrected and potentially inverted to produce corrected dynamically consistent variables. Ertel's isentropic potential vorticity (PV henceforth) is the ideal scalar variable to use in this case because of three inherent principles (Hoskins et al 1985): (1) PV contains both dynamical and thermodynamical components summarizing all relevant information about balanced flows (u, v, θ, and p fields) in a single scalar field; (2) PV is conserved in frictionless and adiabatic processes; and (3) under defined boundary conditions, PV can be inverted to recreate corrected variables in a dynamically coherent manner. The final aspect of PV which makes it suitable in this application is that atmospheric dynamical consistency is defined as the conservation of PV, angular momentum and total energy Bromley 1995, Byun 1999) and is regularly used as an indicator of atmospheric dynamic stability (Davis 1992, Whitehead et al 2014.
Here we use it to assess instabilities created during bias correction which are manifested as differences in the PV field compared with an ideal case and occur due to individual variable bias correction failing to consider intervariable relationships.
The rest of this paper is structured as follows. The next section present the data and methods used in the assessment. Section 3 presents results, while the following section discusses various factors that may affect the results. The last section formulates the main conclusions of this work and outlines the challenges ahead.

Data sources
Bias correction traditionally involves developing a model that maps a biased GCM simulation to an unbiased observed (or reanalysis, as is the case here) dataset. Once this model has been developed, it can be assumed to be applicable for other simulations of the same GCM, usually pertaining to a future climate. The GCM simulation used in this case was the Commonwealth Scientific and Industrial Research Organisation's Mk3.5 (CSIRO-Mk3.5) simulation of historical climate made available by the World Climate Research Program (WCRP) Coupled Model Intercomparison Project phase 3 (CMIP3) (Meehl et al 2007), and the reanalysis product was the ECMWF's ERA-Interim (ERA-I) reanalysis (Dee et al 2011) on 60 model levels.
The zonal (u) and meridional (v) winds, potential temperature θ ( ) and pressure (p) fields were extracted globally from each data source on model levels and interpolated to a fixed pressure level array containing 32 vertical pressure levels: 10,20,30,50,70,100,125,150,175,200,225,250,300,350,400,450,500,550,600,650,700,750,775,800,825,850,875,900,925,950,975, and 1000 hPa. CSIRO-Mk3.5's model top is specified at 4.5 hPa (Gordon et al 2010) limiting our results to a minimum pressure of 10 hPa in order to reduce the effect of differences in model top dampening schemes. Interpolation to a fixed cylindrical equidistant grid with a resolution of 1.5°latitude and longitude was performed on both models with areas below model topography masked. The assessment period presented in this paper was a randomly selected northern hemisphere winter month, January 2001, however similar results were obtained using other months.

Bias correction techniques
Commonly used correction techniques involve independent, individual variable corrections which was the baseline technique for this study. Before we describe the alternate correction approach, we identify commonalities between both methods. During both corrections the mean and standard deviation of the chosen GCM were corrected for each grid point separately in order to match the reanalysis. Equation (1) formulates the corrected variable ′ X based on the GCM variable X , the GCM and reanalysis variable means, X and Ȳ , and the GCM and reanalysis variable standard deviations σ X and σ Y .
( 1 ) X Y These corrections were applied using a single random split sample approach with 75% of the time period used in calibration (or estimation of the parameters X , Ȳ σ X , and σ Y ) and the remainder for result evaluation. The sensitivity of the split sample percentages (values chosen between 50% and 75%) and the selection of randomised samples were assessed with this framework and found to produce almost identical results.
The standard bias correction method, shown in figure 1, corrects each atmospheric variable individually to derive the corrected atmospheric fields: u′, v′, θ′, and p (in this case as all variables are on isobaric levels, p is not corrected). From these corrected variables PV′ is calculated for later comparison. Our proposed correction technique differs in that it involves transforming these variables into the scalar PV field first after which bias correction occurs on this single field producing PV″. The resultant bias corrected PV values (PV′ and PV″) were compared to PV E calculated from ERA-I reanalysis which was used here as an 'ideal' representation of PV. In all cases the results are displayed in PV units where 1 PVU = − − − 10 K m kg s 6 2 1 1 .

PV calculation
Ertel's isentropic PV is the product of absolute vorticity (dynamic component) and static stability (thermodynamic component), and is formulated as follows where g is the gravitational acceleration, ζ θ the relative isentropic vorticity, f the Coriolis parameter (planetary vorticity) and θ ∂ ∂p the partial differential of potential temperature as a function of pressure. The minus sign is convention to define PV as positive in the Northern hemisphere. The relative isentropic vorticity on pressure levels is defined as p where u, v are the zonal and meridional winds, and x, y the winds coordinate directions. R is the ideal gas constant, p represents pressure, T the temperature, θ the potential temperature, and σ the static stability.
Combining (2) and (3) gives Ertel's approximate isentropic PV θ ( ) PV on pressure coordinates for a fixed rectilinear grid, see Bluestein (1992) ( 4 ) p From (4) we see that any arbitrary change in the derivatives of u, v, θ, or p will lead to notable deviations in the resulting PV values.
As previously discussed, PV is useful in diagnosing anomalies produced by non-conservative processes. In this case we are relying on the inherent dynamical consistency of the ERA-I simulations to provide a baseline against which the corrected models are assessed for production of anomalies. This allows differences in the PV field to identify non-conserving processes which occur during each bias correction approach allowing for an objective comparison that determines which correction approach produces more realistic, dynamically coherent fields.

Results
Prior to bias correction, the correction parameters for individual variables were assessed to evaluate the spatial location of biases between ERA-I and CSIRO-Mk3.5 models. The zonal wind component (u) was found to have a greater average temporal and spatial bias when contrasted with the meridional wind (v). The top of the atmosphere was also found to contain larger biases over u, v, and, θ variables resulting from differences in the vertical dampening applied at the top of the atmosphere. In all cases PV was also several orders of magnitude greater in the stratosphere than the troposphere due to the reduced air density and increase in static stability above the tropopause (defined here as PV = 1.6 PVU (WMO 1985)). Additionally, ground levels were found to have larger biases due to the difference in model approaches to ground-surface interaction and differing resolution in surface topography.
The question we set out to answer was whether PV′ is quantifiably worse than PV″ when compared to PV E ? This was quantified in terms of the overall bias in the mean and standard deviation and properties of the spatial dependence (given by the variogram). Root mean square errors over time and other individual time-step correlations were not suitable evaluation metrics due to the lack of time correspondence between reanalysis and GCM simulations. By definition, the basic evaluation metrics of mean and standard deviation differences over the calibration period were negligible and are not presented here. However, there were specific but moderate differences in the ability of the two techniques to match the mean and standard deviation of PV E over the validation period. Figure 2 shows the zonally and temporally averaged bias, as well as the zonally averaged and temporal standard deviation bias, in the corrected models compared to ERA-I. Generally, over much of the vertical profile there is reasonable representation of PV resulting from the two correction techniques particularly throughout the troposphere. However, there is a moderate improvement when comparing PV″ to PV′, and in particular, specific biases in PV′ are observed above the tropopause in both the Northern and Southern hemispheres in both the mean and standard deviation differences which are reduced in PV″.

Variogram analysis
As well as the mean and standard deviation bias presented above, the spatial variability of the mean PV fields were assessed. The empirical variogram was used which calculates at each location the average squared difference of values separated by each lag distance, a measure of the separation of spatial regions. This technique is commonly used in geostatistics to describe how spatial data are correlated with distance. This metric assesses differences in the degree of spatial dependence within a field at varying spatial scales and was used here to indicate which correction technique produced spatial characteristics most similar to PV E . Semivariances were calculated for each vertical level in the zonal direction before being time averaged. Figures 3(a) through (d) shows the empirical variogram for four vertical levels while 3(e) and (f) show a heatmap of bias between the semivariance for each correction model and PV E , standardized by the maximum difference found on each level.
In general, the semivariance of PV″ is seen to be closer to the semivariance produced by PV E field with the following distinct trends: The magnitude of the semivariance reduces as we move down the atmosphere (the PV reduces substantially in magnitude from the stratosphere into the troposphere) until 750 hPa where it begins to slightly increase through to the 1000 hPa level, hence the need for standardization by the maximum difference in 3(e) and (f). The top levels down to 300 hPa show that PV″ performs better on average than PV′ as the latter consistently over represents the semivariance compared to PV E . Below 300 hPa PV′ semivariance sharply reduces in magnitude to be consistently below PV″ and PV E which are both aligned more closely. Three vertical levels (500, 600 and 900 hPa) show PV′ matching with PV E slightly better than PV″. In general we summarize from figures 3(e) and (f) that PV′ in the stratosphere over estimates spatial variance at a range of spatial scales compared to PV″ and PV E while in the troposphere for the majority of levels there is an under-estimation of spatial variance in PV′ with only a handful of levels showing an improvement over PV″. Generally PV″ provides a closer match to the spatial variance structure of PV E over the range of all spatial scales.

Effect of each variable
In addition to individually correcting all variables to produce PV′, correction was performed on subsets containing three of the four variables to evaluate the relative impact of bias correcting each variable on the mean and standard deviation of the PV′ field. (This individual variable correction is not possible with PV″ as bias correction occurs on the PV fields which are a composite of all four variables). When we compare the fields where u, v, and θ were not corrected, in figures 4(b) and (e), we see that the correction of v (meridional winds) provides the smallest improvement to the mean bias correction (seen in its similarity to figures 2(a) and (c)), with the lack of correction having little effect on the representation of PV′. This is due to the formation of meridional winds created from the meridional temperature gradient being similarly represented in both models and hence the correction factors for this wind being small compared to the zonal winds. Bias correction of both u and θ (figures 4(a), (d), and (c), (f)) provide more substantial improvements in the mean bias corrected PV shown by the magnitude of differences with PV E when these variables were uncorrected. We see in particular that θ (figures 4(c) and (f)) has a large effect on biases within the troposphere and changes the tropopause height whereas other biases generally occurred above the tropopause. Not correcting zonal wind, u, (figure 4(a)) creates biases mostly above the tropopause with some intrusions of PV bias descending into the troposphere. The standard deviation biases resulting from uncorrected u and v fields (figures 4(d) and (e)) are very similar to the correction of all variables meaning they are not providing a substantial contribution. In contrast to the mean bias where the correction of θ was very important to the proper representation of the PV field, the bias in standard deviation with uncorrected θ (figure 4(f)) is slightly reduced above the tropopause i.e., the correction of θ worsens the standard deviation bias in PV′.

Discussion
We have addressed two key aspects of PV that should be reproduced in dynamically consistent correction. In particular, biases in the mean and standard deviation of the PV field above the tropopause and the spatial variance in the PV field on all levels, both of which are improved with direct bias correction of the dynamical quantity over constituent variable bias correction. Now we discuss several factors that influence these results as well as a trajectory for future work in the field of dynamically consistent atmospheric bias correction.

Limitations
This paper presented results from a single month's randomised split sampled bias correction approach. Multiple split sampled analyses were undertaken and the results deemed to be indistinguishable and not dependant on the random selection of calibration and validation periods in the time period assessed and thus not shown here. The approach used in this paper could be improved through assessing longer term climatological differences between bias corrections techniques, including changes, if any, in the seasonality of the results. However in this case we were interested in identifying differences in the instantaneous or weather related representations of PV to validate the hypothesis rather than addressing the decadal mean states or climatological differences which would require longer time frames and is left as future work.
It is important to reiterate that the representation of PV is just one of the three measures of atmospheric dynamical consistency, with the other two being momentum and energy balances. In this study we address the first measure and demonstrate that standard bias correction leads to a measurable decrease in the dynamical consistency. To fully evaluate the improvement in dynamically coherent bias correction the PV fields will need to be inverted allowing the calculation of the other dynamical consistency measures.

Implications for dynamical downscaling
Bias correction of inputs intended for use in dynamical downscaling needs to consider the dynamical balance of the corrected variables. Although there were many regions of correspondence between the PV representations of the two correction techniques, individual variable corrections produce PV biases in the spatial variability on most vertical levels and in the mean and standard deviations above the tropopause. The individually corrected fields are in a physically anomalous state in these biased regions which would not otherwise be created in the model simulation. When using this anomalous state in the lateral boundary conditions of an RCM these inconsistencies must be damped through the relaxation zone. Within the RCM domain the model will find a new dynamically consistent result which may differ from that intended. Whether this difference is large enough to cause concern remains uncertain, however a method which allows bias correction while maintaining dynamical consistency would certainly be preferred.

Conclusions
ERA-I reanalysis and CSIRO-Mk3.5 GCM simulations are by definition dynamically consistent with the physical equations driving the models. This dynamical consistency is an important aspect of models and reanalysis products and should be maintained, where possible, during bias correction. Standard atmospheric bias correction techniques applied to boundary conditions for RCM simulations do not consider this dynamical consistency. The results presented here, which make use of the inherent dynamical consistency of the ERA-I reanalysis product, show that for much of the troposphere the inconsistencies introduced are small and may not be of concern, however, these inconsistencies can become relatively large near the tropopause and above. Applying bias correction directly on the PV field was better able to maintain dynamical consistency measured as smaller biases in the mean, standard deviation and spatial variance between the PV field from each corrected model and the PV from ERA-I. The next stage of this work is to undertake a PV inversion on the dynamically corrected fields in order to evaluate other components of atmospheric dynamical consistency: momentum and energy balances.