Multivariate bias correction of regional climate model boundary conditions

Improving modeling capacities requires a better understanding of both the physical relationship between the variables and climate models with a higher degree of skill than is currently achieved by Global Climate Models (GCMs). Although Regional Climate Models (RCMs) are commonly used to resolve finer scales, their application is restricted by the inherent systematic biases within the GCM datasets that can be propagated into the RCM simulation through the model input boundaries. Hence, it is advisable to remove the systematic biases in the GCM simulations prior to downscaling, forming improved input boundary conditions for the RCMs. Various mathematical approaches have been formulated to correct such biases. Most of the techniques, however, correct each variable independently leading to physical inconsistencies across the variables in dynamically linked fields. Here, we investigate bias corrections ranging from simple to more complex techniques to correct biases of RCM input boundary conditions. The results show that substantial improvements in model performance are achieved after applying bias correction to the boundaries of RCM. This work identifies that the effectiveness of increasingly sophisticated techniques is able to improve the simulated rainfall characteristics. An RCM with multivariate bias correction, which corrects temporal persistence and inter-variable relationships, better represents extreme events relative to univariate bias correction techniques, which do not account for the physical relationship between the variables.


Introduction
The direct use of Global Climate Models (GCMs) is limited at regional or hydrological catchment scales because of the coarse spatial and temporal scales used (Caldwell et al. 2009). Regional Climate Models (RCMs) forced with GCM simulations through initial, lateral, and lower boundary conditions are more commonly used to simulate detailed information at finer scales.
Although the GCM-driven RCM simulations show better performance than the raw GCM simulations (Diffenbaugh et al. 2005;Leung and Qian 2009;Di Luca et al. 2016;Li et al. 2018), systematic biases that can be propagated from GCM data into the input boundary conditions of the RCM still remain. These improper boundary conditions can affect model outputs, leading to biases when compared to observation (Caldwell et al. 2009;Xu and Yang 2012;Rocheta et al. 2017;Kim et al. 2020). Applying bias correction, sometimes called bias adjustment but herein as bias correction, to climate model outputs is therefore an important step in improving model performance. Different approaches for correcting and characterizing systematic bias have been applied to either surface variables of the RCM outputs before use for impact assessment, or the GCM simulations prior to downscaling to remove biases in the RCM input (lateral and lower) boundary conditions. The second approach above forms the focus of this study.
For the case where bias is removed as a post-processing step following the climate model simulation, past studies have evaluated various bias-correction approaches that vary in complexity in terms of the mathematical operations used, as well as the nature of the attributes that are corrected. Most common approaches use quantile mapping as the basis for adjusting the distribution of a simulated variable to match 1 3 that of the observations to correct both current and future climate simulations (Wood et al. 2004;Li et al. 2010;Piani et al. 2010). Alternatives that address the model biases in temporal persistence are also widely used given their relevance in hydrological applications (Johnson and Sharma 2012;Ojha et al. 2013). Recently, bias correction methods have evolved toward more advanced techniques that correct inter-variable relationships among the atmospheric variables of the model outputs, such as precipitation and temperature, which possess a strong co-dependence (Mehrotra and Sharma 2015;Sharma and Mehrotra 2016;Cannon 2016Cannon , 2017François et al. 2020). The abovementioned studies, however, were specifically designed for correcting the surface fields of the climate model output and have not been applied to the full atmospheric fields required for RCM lateral boundary inputs.
The method for correcting the RCM input boundary conditions has evolved from simple to more complex techniques that attempt to mimic observed multi-scale relationships in simulations. Xu and Yang (2012) investigated mean and variance bias correction of the RCM boundary conditions, and showed that the bias corrections improved the model performance for climatological means and extremes. Bruyere et al. (2014) corrected the RCM boundary conditions to investigate simulations of hurricanes, which can produce extreme precipitation. The results showed that simple mean bias correction produced the greatest improvement.
The bias correction approaches that have been used to correct the RCM input boundaries in previous studies, however, often overlooked the biases in the persistence-related attributes in the model outcomes, and corrected the biases at a single time scale (Bates et al. 2008;Nguyen et al. 2016). To overcome this limitation, Rocheta et al. (2017) investigated correcting low-frequency rainfall variability representation using nested bias correction, which includes correcting lag1 auto-correlation correction at multiple time scales. They showed that this method produced an added improvement compared with the simple scaling approach. Kim et al. (2020) also used nested bias correction to improve RCM simulations of rainfall extremes. They showed that even simple bias correction techniques could significantly reduce bias in terms of extreme rainfall events. They also showed that nested bias correction performed better than simple corrections for seasonal extremes, suggesting that correcting persistence characteristics is essential for the simulations.
While bias corrections of RCM input boundaries have shown consistent improvement in mean and extreme events in previous studies, these approaches have addressed only the correction of individual variables under the assumption that inter-variable bias is not of key importance. However, if the physical relationships between the variables in dynamically linked fields are not considered, the errors introduced can influence the outputs, including precipitation, temperature, and humidity (Chen et al. 2011;Rocheta et al. 2014). The misconception of interacting physical processes across multiple temporal and spatial scales can potentially lead to the underestimation of extreme events, such as heavy rainfall and droughts (Zscheischler et al. 2018).
In a recent study, Kim et al. (2021) investigated whether the RCM could reproduce spatio-temporal and inter-variable dependence with respect to observed variables. They showed that the univariate bias correction techniques improved temporal and spatial dependence but not multivariate dependence. They highlighted that the mismatch in the physical relationships among the atmospheric variables was not sufficiently adjusted through the relaxation zone, suggesting further investigation to address all systematic biases in the boundary conditions.
In brief, RCM outputs with univariate bias-corrected boundaries still contain large biases in the multivariate dependence even though they exhibit adequate performance for simple statistics, leading to substantial anomalies in simulated extremes. (Kim et al. 2020).
Considering this, the focus of this study is the impact that multivariate bias correction can have on lateral boundary conditions (LBCs) and the lower boundary condition, which represents the sea surface temperature (SSTs) field, particularly for rainfall characteristics. Along with addressing biases in the cross-dependence attributes of an inter-variable field, we addressed biases across multiple time scales using the nesting approach. We also used mean, mean and variance, and univariate nested bias correction methods to compare the model performance between the simple scaling and the more complex approaches.
The paper is organized as follows. In Sect. 2, the datasets and methods which include bias correction approaches and the Weather Research and Forecasting (WRF) model setup are described. Results are presented in Sect. 3. The discussion and limitations of this work are presented in Sect. 4. Conclusions are described in Sect. 5.

Models and data
This study used the Australian Community Climate and Earth System Simulator Earth System Model Version 1.5 (ACCESS-ESM1.5) GCM simulation made available by Commonwealth Scientific and Industrial Research Organization (CSIRO) for the purpose of contributing to the internationally coordinated Coupled Model Intercomparison Project Phase 6 (CMIP6) (Ziehn et al. 2020). For RCM simulations, the Weather Research and Forecasting model (WRF) with dynamical core (ARW), version 4.2.1 (Skamarock et al. 2019) was used. The resolution of the ACCESS-ESM1.5 is based on the N96 Gaussian grid (approximately 1.875°EW × 1.25°NS) with 38 vertical levels extending from the surface to 40 km.
The reanalysis data used for correcting GCM biases was ERA5 which is the fifth-generation model reanalysis of the global climate from the European Centre for Medium-range Weather Forecasts (ECMWF) (Hersbach et al. 2020). It has a resolution of 31 km with 37 pressure levels. The ERA5driven RCM simulation was considered a "perfect" regional climate simulation in this study.
The atmospheric variables [specific humidity q (g/kg), temperature T (K), zonal wind u (m/s), and meridional wind v (m/s)] and the sea surface temperature (SST) in the lateral and lower boundary conditions of RCM were corrected towards those of ERA5. The ERA5 boundary variables were first regridded using both the conservative (for the specific humidity) and bilinear (for the other four variables) remapping methods and then linearly interpolated to match the ACCESS-ESM1.5 horizontal and vertical resolution. The pressures from the surface to the top were recalculated using the bias-corrected fields based on the hypsometric equation. Therefore, the pressure field, including the surface pressure, is also adjusted during the bias correction procedure. All other variables remained identical, and the boundary conditions were built after bias correction. The model was simulated over 31 years from 1 January 1982 to 31 December 2012. The first year of the model simulation was ignored as a spin-up period in order to remove issues associated with the equilibrium state for various soil types (Cosgrove et al. 2003;Chen et al. 2007).

Bias correction approaches
In this study, four different types of bias correction methods were implemented to correct the atmospheric variable (specific humidity, temperature, zonal and meridional winds) and surface fields (sea surface temperature) at the boundaries. Simple scaling methods, which correct climatological mean and mean and variance with respect to those of reanalysis data, were used to correct the magnitude of events independently for each variable at a single time scale. In this study, these corrections were on a daily basis, and the daily statistics were calculated using the data falling within a centered moving window of 31 days (Sharma and Lall 1999). A more sophisticated alternative that also corrected lag1 autoand lag0 cross-correlation across variables with respect to those of reanalysis data, was applied at multiple time scales using a nesting approach which can correct the biases in the persistence related attributes in the model outcomes (Johnson and Sharma 2012). The bias corrections were implemented at 6-hourly time scale over the 31 years GCM dataset. Here we denote X g () is the GCM data, X o () is the reanalysis data, and X g () is the transformed GCM variable. The subscripts 6h , d , m , s , y indicate the time scale for 6-h, day, month, season, and year, respectively. The superscripts g , o indicate the model and observation, respectively.

Mean bias correction
This correction was implemented to calculate climatological means in which each of the GCM simulations was transformed toward that of the reanalysis data. The transformation function can be expressed as follows: where X g d and X o d are the climatological means of the original GCM and the reanalysis data at the daily time scales, respectively. The daily statistics were calculated using a centered moving window of 31 days.

Mean and variance bias correction
Another simple approach for correcting biases in boundary variables was the mean and variance bias correction. This approach is similar to the mean bias correction, but it is extended to include the variability of the atmospheric variables and surface fields. This can be expressed as below by applying a standard deviation correction to the mean correction: where s g d and s o d represent the standard deviation of the original GCM and the reanalysis variables at the daily time scale, respectively. After this bias correction, the mean and variance of bias-corrected GCM simulations ( X g 6h ) are similar to those of reanalysis data. These statistics were estimated on 6-hourly GCM data in a daily group over 31 years, and the centered 31 days moving window was applied.

Nested bias correction (NBC)
A more comprehensive alternative extended to include lag1 auto-correlation attributes on multiple nested timescales was used. This approach has shown an improvement in the rainfall mean and variability at multiple time scales (monthly, seasonal, and annual time scale) in the previous studies (Johnson and Sharma 2012;Kim et al. 2020). Based on a standard autoregressive lag1 model, lag1 auto-correlation of GCM was replaced by the observed lag1 auto-correlation as follows: where z g t is standardized GCM variable for time step t, r o t and r g t are lag1 auto-correlation of the reanalysis and GCM variables, respectively. The corrections were applied at daily, monthly, seasonal, and annual time scales, which are extended to the daily time scale with respect to the previous studies. These different sets of bias-corrected time series are then combined with the nesting approach. This means that the bias-corrected daily GCM data is aggregated to a longer timescale (month, season, and year) and follows the same procedures for other timescales as well. The daily GCM data thus incorporates the effect of bias correction at longer time scales, which include nested bias-corrected monthly, seasonal, and annual time scales. The general process of the nesting approach is described in (Johnson and Sharma 2012;Mehrotra and Sharma 2012). From Srikanthan and Pegram (2009), the corrected variables at different time scales can be applied to the raw 6-hourly GCM data as a form of weighting factor: where X g 6h is the 6-hourly corrected value and X g 6h is the original 6-hourly value, X d,m,s,y is the daily GCM data for day d, month m, season s, and year y. X m,s,y is the aggregated monthly value, X s,y is the aggregated seasonal value, and X y is the aggregated yearly value. Finally, X g h exhibits the same persistence attributes as X o () across the nesting time scale. These bias corrections were carried out at each grid cell over the 31 years GCM dataset. The full details of the nested bias correction can be found in Johnson and Sharma (2012).

Multivariate bias correction (MBC)
The univariate bias correction methods as described above are designed to correct variables independently. Although univariate distribution features can be adjusted toward reference data, they can generate inappropriate cross-dependence structures among multiple variables simulated by GCMs, leading to substantial anomalies in the simulation.
To correct the inter-variable relationships, we investigate the multivariate recursive nesting bias correction method which is extended to also include lag0 and lag1 auto-and cross-correlations among the variables. This approach was presented by Mehrotra and Sharma (2015) and has shown that it can correct for variability and persistence biases in the atmospheric variables on a range of timescales. A large number of parameters, however, may cause large estimation errors and possibly hide (overestimate) the true quality of the model (Ehret et al. 2012;Mehrotra and Sharma 2021). In this study thus some simplifications in the model structure were applied and lag1 cross-correlations are assumed to be unimportant and are ignored. We investigate the distributional attributesmean and variance and the dependence attributes-the lag1 auto-and lag0 cross-correlation coefficients were calculated at four selected time scales, daily, monthly, seasonal, and annual. Three iterations (for the recursive scheme) were adopted to reduce biases in the nesting procedure across all timescales (Mehrotra and Sharma 2012). We denote the vectors Z g () and Z o () are the standardized GCM and reanalysis variables with zero mean and unit variance, respectively.
A multivariate transformation for the day t is based on a standard Multivariate AutoRegressive order 1 (MAR1) model and can be expressed as follows (Salas 1980): and where , and , are coefficient matrices of the lag1 and lag0 cross-correlations of the reanalysis and GCM data, respectively. The random vector ε t is assumed to be mutually independent, indicating that where is the identity matrix, and T denotes the transpose of the matrix.
Since the elements of ε t are mutually independent having a zero mean and unit variance, the coefficient matrices and or and can be expressed as follows (Matalas 1967): which are the lag0 and lag1 cross-correlation matrices. The elements of these matrices corresponding to variables i and j are obtained from: where N is the total number of samples. The coefficient matrices and are obtained by singular value decomposition, which is a method of decomposing a matrix into three component matrices that are easy to manipulate and analyze.
These procedures are applied by rearranging the Eq. (5b) for ε t : where ε t represents a standardized vector that has removed the lag0 and lag1 auto-and cross-correlations obtained from the Z g t series. This vector is then used to modify Z g t to Ź g t that has the lag0 and lag1 attributes of reanalysis data.
Finally, by adding back the means and standard deviations of reanalysis data, the bias-corrected output provides appropriate attributes in means, standard deviations, lag1 auto-, and lag0 and lag1 cross-correlations.
Similarly, MAR1 with periodic parameters can be calculated as: where Z g y, represents the vectors of the standardized GCM periodic series in time interval , and year y , indicating that the series all contain the seasonal cycle. The matrix parameters of this model can be obtained following as: where 0, and 1, are the periodic lag0 and lag1 matrix correlation matrices. The elements of these matrices are obtained as similar to those with constant parameters. As abovementioned, some simplifications in the MAR1 model structure are applied in this study to avoid a large number of parameters, which may hide the true quality of the model simulations. The lag1 cross-correlations are thus assumed to be unimportant, and the parameters and are considered as diagonal matrices in which the elements outside the main diagonal are all zero. This simplification is also applied to the periodic model. The elements of these matrices corresponding to variables i and j can be obtained as: These different sets of bias-corrected time series are then combined with the nesting approach as abovementioned in the Eq. (4).

Performance assessment
Three statistics were used to quantify the impact of bias correction: mean absolute error (MAE), Pearson correlation coefficient (R), and bias. The MAE is defined as where N is the total number of grid cells, X mod n and X obs n represent the climatological WRF outputs and observation data at each grid cell, respectively.
The differences of the means between the model simulation and observation can be defined by mean bias as: The biases are calculated for each vertical level at each grid cell.
To assess the ability of the model simulations in terms of multivariate aspect, Fisher z-transformation was used here to convert the correlation coefficient to a transformed variable (Fisher 1915(Fisher , 1921. The test statistic ( z ) can be represented as a standard normal deviate where z M r and z O r are Fisher's z transformation of the model and observed correlation coefficient ( z r = 1 2 ln 1+r 1−r ), and s M and s O are sample sizes of the model outputs and observation, respectively. This study assumed a 5% significance is a criterion for rejecting the null hypothesis, which means that there is a significant difference in the correlation between the two simulations if the computed z-value is not in the z -value distribution table. Full details of this formulation can be found in Kim et al. (2021).

Results
This section will first evaluate the statistics of the biascorrected atmospheric variables in the boundary conditions with respect to the ERA5 datasets. RCM simulations are then compared to the ERA5-driven RCM outputs over the Australasia CORDEX domain to reveal the effectiveness of bias correction in terms of the climatological mean, standard deviation, lag1 auto-correlation, and multivariate dependence measures. The ERA5-driven RCM simulation is treated as a simulation with "perfect" boundary conditions. To evaluate the model performance, the outermost five grid cells from the boundary were trimmed off as the relaxation zone to avoid a bias due to boundary effects.

Evaluation of bias-corrected GCM datasets
This section evaluates the impact of bias correction on the GCM outputs that have been subsetted along the lateral boundary conditions. Four different bias correction approaches introduced in Sect. 2 were applied to the atmospheric variables (u, v, T, q) and sea surface temperature (SST) to the raw GCM datasets. From the results presented here, we can evaluate whether the bias correction approaches reduce a bias before downscaling. Figure 1 shows a scatterplot for specific humidity q (g/kg) comparing ERA5 to the four bias correction techniques: GCM(M), GCM(MSD), GCM(NBC), and GCM(MBC). These indicate the raw GCM datasets, GCM with mean bias-corrected boundary conditions, GCM with mean and standard deviation bias-corrected boundary conditions, GCM with nested bias correction, and GCM with multivariate bias-corrected boundary conditions, respectively. The results show the effect of each bias correction approach at daily, monthly, seasonal, and annual time scales for the three statistics over 38 vertical levels along the western boundary. Details regarding the other variables, T and w, along the eastern, northern, and southern boundaries can be found in the online supplemental material: Fig. A1 to A10. Each statistic presented in the figure shows a separate point for each day, month, season, and year; hence there is a single point for each grid cell at each time scale.
It is clear that the bias corrections produce an improvement as expected by the construction of each formulation that is specified to adjust in the mean, standard deviation, and lag1 auto-correlation fields. We see that there is a large difference between before and after bias correction for lag1 auto-correlation, particularly for the annual time scale, which can affect the model performance in simulating persistence characteristics and the temporal variability of the rainfall. Despite the GCM being different in all statistics, the biases for the simple statistics, mean and standard deviation, are noticeably reduced compared to the old versions of simulations (Rocheta et al. 2017;Kim et al. 2020). This indicates that the GCM performance has been improved in CMIP6 as it represents state-of-the-art climate modeling. Figure 2 shows a scatter plot comparing the bias-corrected and uncorrected atmospheric variables to ERA5 in terms of multivariate relationships among the three atmospheric variables over all vertical levels. Specific humidity above the 24th level has been ignored during multivariate bias correction as it is always close to zero which can lead to an unstable output. Unlike the model performance shown in Fig. 1, there are large biases even after the bias corrections, GCM(M), GCM(MSD), and GCM(NBC), particularly in the annual time scale. This indicates that the univariate bias corrections are not capable of correcting physical relationships, which may cause unrealistic behavior of the atmospheric variables. On the other hand, GCM(MBC) performs well, showing most points are spread around the 45-degree line, meaning that bias in the relationships among the variables is well corrected. Details along the eastern, northern, and southern boundaries can be found in the online supplemental material: Fig. B1 to B3. Figure 3 shows a bias map of seasonal SST over 31 years comparing the bias-corrected and uncorrected SST to ERA5 covering the research domain for three statistics. The results show that the raw GCM contains a large positive bias for mean and standard deviation over the domain and a negative bias in the tropics. The results show that each bias correction performs well with regard to the statistics as they are constructed.

Evaluation of RCM simulations over the Australasia CORDEX domain
This section presents the RCM outputs to assess whether the influence of the bias-corrected boundary conditions is preserved even after the RCM simulations in terms of four aspects: mean, standard deviation, extremes, and the multivariate relationship.

Climatological mean
We first evaluate the RCM simulations with regard to the climatological mean. Figure 4 shows a bias map of RCM simulations compared with ERA5-driven RCM outputs for three surface variables for monthly mean statistics over 30 years.
The results show that RCM(GCM) produces a large positive bias for the specific humidity and temperature over the domain. It can also be seen that the RCM output with the raw GCM boundary conditions are significantly different for precipitation showing a significant negative (dry) bias over the Western Pacific Warm Pool and positive (wet) bias along the South Pacific Convergence Zone. In contrast, the RCM with bias-corrected boundary conditions shows a noticeable improvement across the variables. While negative biases are shown in western and southern Australia for the specific humidity, a bias is substantially reduced with the 0.1 mean absolute error over the domain. We see that the specific humidity and temperature are quite insensitive to the bias correction approaches with respect to a monthly climatological mean. It is clear that correcting variability of the atmospheric variables along the boundaries can better represent the monthly mean precipitation in the tropics compared to the simple mean bias correction technique.

Inter-annual variability
We then evaluate the model performance in terms of the standard deviation aspect. Figure 5 shows a bias map of the coefficient of variation (CV) of annual precipitation. The higher the CV, the more variable the inter-annual precipitation of a location is. The mean absolute error over the domain is represented in the bottom right of each map. We see that RCM(GCM) produces a large positive bias compared to the ERA5-driven RCM outputs in the tropics, Indian Ocean, and South Australia. This rainfall variability can be influenced by several global and local climate drivers, such as the Indian Ocean Dipole, El Niño Southern Oscillation, and northwest cloudbands (Risbey et al. 2009). Whereas the RCM outputs with bias-corrected boundary conditions show improvement, biases still remain in the tropics and on Australia's northwest coast. The results show that RCM(M) can improve rainfall variability even though the boundary variables were corrected without regard to the standard deviation. From the results we see that RCM(MBC) performs well with the lowest mean absolute error, 4.5, over the domain. Figure 6 compares the models and RCM(ERA5) for standard deviation of the three surface variables at the presented here shows some degradations in MAM for precipitation, suggesting that correction on mean statistics only may increase the model uncertainty in capturing output variability. The figure shows that RCM(MSD), RCM(NBC), and RCM(MBC) improve most time scales for the mean absolute error compared to RCM(GCM). While the complex approaches often show a marginal improvement, and can be worse than the simple approaches for some variables in some seasons (DJF for specific humidity and SON for precipitation), we see that they generally perform better over the domain. This means a corresponding improvement in the statistics of RCM outputs over the domain, indicating that the rainfall simulation can be changed as a result of changes in the atmospheric variables in the boundary conditions.

Extreme events
In order to evaluate the effect of bias-corrected boundary conditions on extreme events, we used the 99th percentile determined for each grid cell and each variable at a daily time scale individually. Figure 7 shows a bar plot similar   Fig. 6, we broadly see that the better representation of relationships between the atmospheric variables in the boundaries translates effectively to the fine-scale extreme events. RCM(NBC) shows worse performance for precipitation in DJF than the RCM(MSD), even though it includes the lag1 auto-correlation attribute, indicating that using more complex techniques without regard to the physical relationship between the variables does not guarantee the better performance for rainfall extremes.

Multivariate relationship
This section investigates the effect of bias correction with regard to multivariate dependence in the RCM outputs to assess whether or not the relationship among the variables is preserved even after passing through the boundary relaxation zone.
The impact of lateral boundary conditions is linearly reduced through the relaxation zone (from the outermost specified zone to the inner zone). One can expect that the mismatch in the multivariate relationship will be increased when compared to that of ERA5. Here pairwise correlation of the atmospheric variables was computed for each RCM simulation and then compared against RCM(ERA5) with a 5% significance level. Grid cells in the specified zone and the sixth grid cells, which are just inside the relaxation zone, were assessed along the boundaries as suggested in Kim et al. (2021). Table 1 shows the results of the multivariate cross-correlation between the RCM simulations and RCM(ERA5) at the 5% significance level as a criterion for rejecting the null hypothesis that the correlation is equal. For example, one can test the hypothesis that the correlation between temperature 'a' and specific humidity 'b' of an RCM simulation is the same as that for temperature 'c' and specific humidity 'd' of observation (i.e., H 0 ∶ ab = cd ). The percentage indicates a number of grid cells showing a significant difference in multivariate relationship compared to RCM(ERA5). As mentioned above, specific humidity close to zero was ignored during calculation. The results show that 79.9% of RCM(GCM) grid cells are significantly different compared to RCM(ERA5) for the variable wind speed (w) and temperature (T) at a daily time scale along the specified zone at western boundary. Although the RCM(M), RCM(MSD), and RCM(NBC) marginally improve the multivariate aspects of the variable combination, large differences for all three relationships are still shown ranging from 62.7 to 78.3% along the western boundary. This indicates that those bias corrections lack the ability to capture the inter-variable relationship presented in RCM(ERA5). On the other hand, RCM(MBC) reduces the bias by showing 21.6-72.0%, and substantial improvements are shown in correlations with wind speed. While the improvement falls into 50.0-63.5% after passing through the relaxation zone as expected, RCM(MBC) still shows better correlations between the variables. Although there are some exceptions where the effect of multivariate bias correction is reduced in the interior model domain, w and q at 6th grid line in the north and T and q at 6th grid

Discussion
The work presented in this study demonstrates the influence of bias correction techniques on the input boundary conditions of an RCM. The focus here is whether multivariate bias correction can reduce bias with respect to the intervariable relationships and, if so, whether this improvement is preserved through the relaxation zone where the model is nudged or relaxed towards the driving model. This section is divided into three parts. Section 4.1 discusses the performance of the different versions of the raw GCM simulations. Section 4.2 evaluates the effect of the relaxation zone, Sect. 4.3 discusses the implications of the bias correction techniques used in this study, and Sect. 4.4 evaluates the advantages of the multivariate bias correction approach. Finally, Sect. 4.5 discusses the limitations and future work.

Comparison with previous research
This study used state-of-the-art GCM simulations and the WRF model to assess the performance of RCM simulations using various bias-correction alternatives. The ACCESS-ESM1.5 model used here represents Australia's contribution in CMIP6 and has shown significant improvements over the historical period against observations over the previous version, ACCESS-ESM1 (Ziehn et al. 2020). One can expect that the systematic bias will decrease as the model is updated; the newer the model, the lower the bias. The results of this study also agree with this expectation by showing that the RCM simulation without bias correction over the domain has improved when compared to previous studies that have used older generation GCMs (Rocheta et al. 2017;Kim et al. 2020Kim et al. , 2021.
However, it is also clear that the newer model still contains biases even though it has been updated. Although the GCM used here shows a relatively good performance in the climatological mean, there are significant differences in other statistics: standard deviation, lag1 auto-, and lag0 cross-correlation, leading to considerable anomalies in simulation of extremes.
It should be noted that correcting RCM boundary conditions before downscaling is important to reduce the biases in the outputs, and more complex techniques, which extend to lag1 auto-and lag0 cross-correlation, appear to be improving details that are important in capturing persistence and maintaining physical consistency (Fig. C1 to C4). This is expected to provide a complete picture of the efficacy of the ACCESS models in simulating the historical climate, as well as a better estimation of future climate change as noted in .

Can the impact of the bias correction along the outermost zone be preserved inside the RCM domain?
The results in Table 1 showed that improvements are certainly evident, and MBC reduces the bias related to physical inconsistency by up to 53.9% in the western boundary compared to those obtained from RCM(GCM). We see that MBC is capable of preserving the inter-variable relationships even after passing through the relaxation zone, showing up to 24.5% improvement. This indicates that the bias-corrected dependent structure in the boundaries of the RCM can be preserved inside the domain despite the influence of the internal dynamics of WRF that propagate into the atmospheric fields in complex ways. This result is in line with a previous study ) that has shown that a large portion of bias-corrected information infiltrating into the model interior is lost in the process of generating the lateral boundary conditions and through relaxation zones where GCM data is passed into the RCM. We see that although MBC still shows the greatest improvement at the specified zone and adjacent grid cells to the relaxation zone, the effect of bias correction, which was almost perfectly corrected toward ERA5, immediately deteriorates after generation of the boundary conditions. This indicates that the multiple interpolations for the purpose of bias correction and the model configuration may introduce additional bias in the creation of the lateral and lower boundary conditions. These are generally unavoidable limitations of a regional numerical model that has boundaries where the model could suffer from changes in resolution with differing parameterizations and physics schemes derived from the GCM, and further assessment is needed to address such issues (Warner et al. 1997;Wu et al. 2005).

Can the sophisticated bias correction of the GCM outputs effectively improve the RCM rainfall simulation?
Improvements in the model performance relative to simple bias correction techniques, which do not account for the dependence structure, are shown in the RCM(NBC) and RCM(MBC) outputs over the Australian CORDEX domain, particularly for precipitation. We see that the RCM performs better in capturing the rainfall coefficient of variation and standard deviation (Figs. 5 and 6) with the increasing complexity of the bias correction techniques. In contrast, RCM with uncorrected boundary conditions introduces a large bias across the variables used in this study. Although RCM with complex bias-corrected boundary conditions, RCM(NBC) and RCM(MBC), show degradation in comparison with RCM(MSD) for precipitation in SON, they broadly perform better at different time scales. We see that the RCM simulations tend to produce more bias in the ocean than inland areas. RCMs with simple bias correction, RCM(M) and RCM(MSD), present larger bias in rainfall climatological mean and variability in the ocean, indicating that they may be weak to capture the seasonal SST variation affected by a large-scale circulation. (Table A1 in the online supplemental material). Previous studies have shown that complex correction techniques affect the model's internal process and can correct details for rainfall extremes as well as its magnitude (Kim et al. 2020). This work represents that the effectiveness of increasingly sophisticated techniques, which aim to correct temporal persistence and physical consistency, adds improvement to rainfall characteristics, particularly in the simulation of rainfall variability and extremes.

What are the advantages of multivariate bias correction over the univariate bias correction?
This work focuses on the inter-variable relationships in RCM input variables and the expected improvement in simulations by correcting multivariate dependence at the boundaries. A combination of interacting physical processes across multiple climate variables, such as relative humidity, temperature, and wind speed, often leads to floods, wildfires, heatwaves, droughts, and even compound events (Zscheischler et al. 2018;Kim et al. 2023). Poor representation of these dependencies makes them challenging to foresee, leading to a significant impact on the risk assessment. These implications could be understood through the results presented in this study. Although simple techniques show comparable performance to the complex ones in the mean-field, they produce large biases in the tropics and the Indian Ocean.
Although RCM(NBC) shows good performance in the attributes used here (Fig. 1), it should be noted that more complex techniques that do not take into account inter-variable correlations do not guarantee better performance for extreme events. This identifies the importance of multivariate bias correction, especially when persistence biases are accounted for. It is observed that RCM(MBC) better preserves the physical relationships between the atmospheric variables (Table 1) in the specified zone and after passing through the relaxation zone, generally showing a better representation of rainfall variability and extreme events (as shown in Figs. 5, 6, and 7).

Limitations and future work
Several limitations may influence the results presented here. While this work has used newer models to assess model performance with several bias correction approaches, we only focused on one RCM (WRF) driven by one GCM (ACCESS-ESM1.5). Hence, coupling different GCM and RCM may change the results. Based on previous work that investigated the models' performance for CORDEX dynamical downscaling, even the updated models still contain biases (Di Virgilio et al. 2022). This indicates that the impact of the bias correction on the RCM boundaries is not significantly different if alternate GCMs are used.
The RCM used in this study was WRF4.2.1 which has complex internal dynamics, alternative options and an elaborate setup. Although the parameterizations for the WRF simulations were chosen following past evaluations, different setups would affect the model outputs.
In this study, the bias corrections were applied towards ERA5 fields, meaning that the ERA5-driven WRF outputs were used as ideal cases to assess the model performance against. Although reanalysis produces the most accurate outputs overall compared to the in-situ data (Kim et al. 2020), the model performance can be limited by the bias that comes from the reanalysis data, in particular for precipitation (Moalafhi et al. 2016(Moalafhi et al. , 2017. It should also be noted that specific humidity at higher levels where it is always close to zero was ignored during the multivariate bias correction as it makes the models unstable and creates physically unrealistic values. The parameters in the correlation structure, for example, including the lag1 auto-and lag0 cross-correlation in a contemporaneous model, are applied to a normalized variable. These changes, even if minor, influence the original value that was very close to zero to become unrealistic. The sea surface temperature was also corrected by nested bias correction only for RCM(MBC) as it shows high seasonal and annual variability but relatively much lower daily variability than the atmospheric variables. This large difference in variability between the daily and the longer time scales can cause unrealistic monthly variation during the multivariate correction.
Furthermore, the structure of the bias correction model makes implicit assumptions where biases are significant. For instance, correcting bias at daily and longer time scales implicitly assumes that diurnal patterns are simulated accurately. While these assumptions may not impact the overall results, there can be situations where their implications become important.
As mentioned above, for a better understanding of the underlying physical consistency and interaction between the variables under climate change, this study can be extended to future simulations. We note that the performance of MBC applied to future climate projections will be affected by any non-stationarity in the inter-variable correlation Guo et al. 2019Guo et al. , 2020Nguyen et al. 2020). Future work should investigate the impact of this effect. Additionally, it would be of interest to extend the bias correction approaches presented here to consider the spatial dependence across grid cells to reduce the uncertainty in the bias estimates (Kim et al. 2021;Switanek et al. 2022).

Conclusions
Despite the efforts to understand single drivers of extremes, most major weather and climate-related catastrophes are caused by a joint occurrence of different types of extreme events (Zscheischler et al. 2018). Correcting single variables independently, as many previous studies have done under the assumption that there is no (or negligible) bias in the dependence structure, limits the model performance in the simulation of extreme events. Although several studies have used sophisticated approaches that also corrected lag1 auto-and lag0 cross-correlation across variables to deal with inter-variable relationships, only the surface variables of the model outputs, such as precipitation and temperature, have been addressed.
This study investigated whether correcting the RCM input boundary conditions could reduce bias inside the domain. Several bias correction techniques, from simple scaling, which has been used in previous studies, to sophisticated techniques, which can correct persistence and physical relationships between the variables compared to the reanalysis data, have been applied. The corrections were applied to all the vertical levels of GCM based on four statistics: mean, standard deviation, lag1 auto-correlation, and lag0 cross-correlation. Substantial improvements in model performance are shown after applying bias correction to the boundaries of the RCM in the statistics used here. This work shows that the effectiveness of increasingly sophisticated techniques substantially improves rainfall characteristics. The RCM with multivariate bias-corrected input boundary conditions represents extreme events better than univariate bias correction techniques, which do not account for the physical relationship between the variables, and considerably better than the approach of using uncorrected GCM simulations to serve as the lateral or lower boundary inputs into an RCM.