Bias Correction of Climate Model’s Precipitation Using the Copula Method and Its Application in River Basin Simulation

During the last few decades, the utilization of the data from climate models in hydrological studies has increased as they can provide data in the regions that lack raw meteorological information. The data from climate models data often present biases compared to the observed data and consequently, several methods have been developed for correcting statistical biases. The present study uses the copula for modeling the dependence between the daily mean and total monthly precipitation using E-OBS data in the Mesta/Nestos river basin in order to use this relationship for the bias correction of the MPI climate model monthly precipitation. Additionally, both the non-corrected and bias corrected data are tested as they are used as the inputs to a spatial distributed hydrological model for simulating the basin runoff. The results showed that the MPI model significantly overestimates the E-OBS data while the differences are reduced sufficiently after the bias correction. The outputs from the hydrological models were proven to coincide with the precipitation analysis results and hence, the simulated discharges in the case of copula corrected data present an increased correlation with the observed flows.


Introduction
The integration of climate change projections, which are both obtained from General Circulation Models (GCMs) and Regional Climate Models (RCMs), in the hydrological models of rainfall-runoff has been conducted to investigate the response of the hydrologic cycle with an emphasis on the projected discharges under climatic variations [1][2][3].The most common implemented methodology in assessing the climate change impacts on water resources is the utilization of climate change series as input data to hydrological models.However, this specific approach presents various types of uncertainties both for GCMs [4,5] and RCMs [6,7] data.Kundzewicz et al. [8] present an overall scheme of the uncertainty that is not only driven by the non-correlated data (i.e., simulations versus observations) but also by the observed data itself, the modelling procedure and the future climatic projections.Nevertheless, the analysis of the uncertainty in projected hydrological changes under climate change in 12 large-scale river basins [9] demonstrated that hydrology models are the lowest contributors of uncertainty when dealing with small and mean flows.
The data from RCMs differ from the observed data due to its incapability in incorporating all the chemical and physical processes that take place in the atmosphere [10], especially at the regions where the topography is complex [11].Thus, the correction of biases between the climate models and observed values is mandatory in order to achieve greater accuracy in hydrological studies.In literature, various approaches have been presented, which includes the analysis of their implications in ameliorating the bias between GCMs or RCMs and the observed or reanalysis precipitation and temperature data [12,13].Watanabe et al. [14] evaluate five bias correction methods using monthly temperature and precipitation data from several GCM models while they propose a new method for bias correction with satisfactory results for the variation's coefficient, the mean and standard deviation values.Additionally, Räty et al. [15] evaluate nine different methods for the bias correction of ensemble RCM simulations and they conclude that the efficiency of each method is affected by different reasons, such as the location or the part of the study distribution.Consequently, they propose a combination of methods for an accurate bias correction.A recent study by Collados et al. [16] investigated some widely used methods for generating potential scenarios of the future climate using RCM simulations with consideration of the meteorological droughts in the region of Alto Genil Basin in Spain.Another research in the region of Eastern Spain (Serpis River Basin) by Pulido-Velazquez [17] considered hydrological scenarios for the future.Finally, an innovative propose analysis was conducted by Mao et al. [12] who used the copula method in order to correct the biases between the daily precipitation obtained by reanalysis and WRF model in Germany.Teutschbein and Seibert [18] present a thorough review of the applicable bias correction methods for RCMs.The derived results were used for the hydrological simulation of the selected river basin and the outputs were compared with observed datasets.
The copula method has also been used in some meteorological studies, mainly for the description of the dependence between two or more variables.Piani [19] used copulas for modelling the dynamical relationship between hourly values of temperature and precipitation for six stations in Germany.Using monthly values for the same parameters (temperature and precipitation), Cong and Brady [20] used the copula method in order to model their dependence for the region of Scania in Sweden.Additionally, an analysis for the dependence of temperature and precipitation between four stations located in the Eastern and Western Mediterranean area was conducted by Lazoglou and Angonostopoulou [21].However, apart from the studies that use copula for modeling the dependence between several climate parameters, there is a relatively limited number of studies that incorporate this method for bias correction [12,22,23].One of the most characteristic studies for bias correction was conducted by Mao et al. [12] who modelled the relationship between daily precipitation coming from reanalysis and WRF models before they used this function for the bias correction of the daily WRF precipitation.
The lack of hydrometereological data in various areas hinders the successful implementation of environmental modelling processes.O'Riordan [24] demonstrates that the lack of historic data or even the comprehensiveness of monitoring often contributes to a false impression of what is happening.A solution for this problem is provided by the use of climate models or reanalysis and daily gridded datasets, which namely produce the E-OBS data.Regarding the reanalysis and E-OBS data, the literature clearly demonstrates that the presented accuracy of this type of data is relatively high, particularly in the case of temperature [25], and they can be used as proxies for observations to force hydrological models, especially in the regions with limited meteorological stations [26].On the other hand, the back testing data (also known as hindcasts), i.e., historic data derived from predictive climate models, that will be used as forcing conditions for the climate change simulation demonstrate some substantial spatial biases [27].Consequently, bias correction techniques should be implemented.
Copula is a widely used method in the financial sector and in the fields of insurance [28][29][30][31].During the last decades, it mainly gained popularity in the hydrology.According to Zhang and Singh [31], the main advantage of the copula is the ability to analyze the dependence of several variables that have different marginal distributions.Additionally, copulas can model non-linear dependencies [32].Moreover, there are several studies in which the copula method is used for analyzing different characteristics of hydrological phenomena, such as drought.In 2010, Yang [33] uses the copula method for the investigation of drought characteristics with several drought indices.Several scientists [33][34][35] incorporate copulas for the investigation of the dependence between drought severity and duration.Apart from drought, the relationship with other climate variables, such as as precipitation [36] or modelling the interdependence between two rivers [37], has also been studied.
The present study aims to correct the biases between the E-OBS and MPI model's precipitation, which formulates a mathematical function (copula) that could also be used for studies that examine future data.The novelty of the present study is that the copula family is estimated using only the E-OBS data (monthly mean and totals) before it is applied to the MPI model for bias correction.Additionally, in this research, a large number of copula families are tested in relation to the ones proposed in the literature.The second scope and novelty is to investigate the performance of the proposed statistical method by (i) integrating the produced results to the simulation procedure of a transboundary river basin and (ii) validating the accuracy of the produced hydrological results.The study area is an important hydrological region both in Bulgaria and Greece and no similar studies (to the authors' knowledge) have been conducted in the past.The proposed methodology is believed to contribute to the bias correction of RCMs datasets and hence, the coupling of climate models with the hydrologic simulation models at the river basin scale could be assessed with increased accuracy.

Data and Case Study Area
In the present study, three different data sources are used and analyzed.Initially, E-OBS daily gridded (Version 17.0, 0.25 deg.regular grid) precipitation data produced by the European Climate Assessment & Dataset (ECA & D) [38] with a spatial resolution of 0.25 • × 0.25 • are exploited.Additionally, the daily precipitation records for the study area, which was namely the Mesta/Nestos river basin, were retrieved from the regional climate model of the Max Plank Institute (MPI model) with a spatial resolution of 0.5 • × 0.5 • .All the gridded datasets cover the studied river basin for a time period of 20 sequential years (1981)(1982)(1983)(1984)(1985)(1986)(1987)(1988)(1989)(1990)(1991)(1992)(1993)(1994)(1995)(1996)(1997)(1998)(1999)(2000).The specific transboundary river basin is located in southeastern Europe and is almost equally shared between Bulgaria (upstream country) and Greece (downstream area).It is a natural basin in terms of land uses since forest and semi-natural areas cover 77.34% of its 5619 km 2 while from its 1.62 × 10 9 m 3 of discharges, 58% and 42% accumulate in the Bulgarian and Greek parts of the basin, respectively [39].Finally, in order to evaluate the reanalysis dataset, the daily precipitation values from the three stations located close to the Greek part of the basin (Figure 1) were exploited, covering a time period of 20 years (1981-2000).
Water 2019, 11, 600 3 of 19 scientists [33][34][35] incorporate copulas for the investigation of the dependence between drought severity and duration.Apart from drought, the relationship with other climate variables, such as as precipitation [36] or modelling the interdependence between two rivers [37], has also been studied.The present study aims to correct the biases between the E-OBS and MPI model's precipitation, which formulates a mathematical function (copula) that could also be used for studies that examine future data.The novelty of the present study is that the copula family is estimated using only the E-OBS data (monthly mean and totals) before it is applied to the MPI model for bias correction.Additionally, in this research, a large number of copula families are tested in relation to the ones proposed in the literature.The second scope and novelty is to investigate the performance of the proposed statistical method by (i) integrating the produced results to the simulation procedure of a transboundary river basin and (ii) validating the accuracy of the produced hydrological results.The study area is an important hydrological region both in Bulgaria and Greece and no similar studies (to the authors' knowledge) have been conducted in the past.The proposed methodology is believed to contribute to the bias correction of RCMs datasets and hence, the coupling of climate models with the hydrologic simulation models at the river basin scale could be assessed with increased accuracy.

Data and Case Study Area
In the present study, three different data sources are used and analyzed.Initially, E-OBS daily gridded (Version 17.0, 0.25 deg.regular grid) precipitation data produced by the European Climate Assessment & Dataset (ECA & D) [38] with a spatial resolution of 0.25° × 0.25° are exploited.Additionally, the daily precipitation records for the study area, which was namely the Mesta/Nestos river basin, were retrieved from the regional climate model of the Max Plank Institute (MPI model) with a spatial resolution of 0.5° × 0.5°.All the gridded datasets cover the studied river basin for a time period of 20 sequential years (1981)(1982)(1983)(1984)(1985)(1986)(1987)(1988)(1989)(1990)(1991)(1992)(1993)(1994)(1995)(1996)(1997)(1998)(1999)(2000).The specific transboundary river basin is located in southeastern Europe and is almost equally shared between Bulgaria (upstream country) and Greece (downstream area).It is a natural basin in terms of land uses since forest and semi-natural areas cover 77.34% of its 5619 km 2 while from its 1.62 × 10 9 m 3 of discharges, 58% and 42% accumulate in the Bulgarian and Greek parts of the basin, respectively [39].Finally, in order to evaluate the reanalysis

Applied Methodological Framework
A formal definition of a copula was given by Nelsen in 1999 [40].Copula is a multivariate function with uniformly distributed marginal distributions between [0,1].If X and Υ are two studied random variables and F, G are their marginal distributions, respectively, there is function H, which is the joint distribution of the initial variables.The most important theorem for copulas was defined by Sklar in 1959 [41] and the general definition of this theorem for d dimensions says the following: given a d-dimensional distribution function H, there exists a d-copula C such that for all (x1, x2, . . ., xd) ∈ R n .
The copula C is uniquely defined on ∏ d j=1 ran Fj and is therefore unique if all the marginals are continuous.Conversely, if F1, F2, . . ., Fd are d (1-dimensional) distribution functions, the function H defined through Equation ( 1) is a d-dimensional distribution function.For the bivariate copula that is the simplest case, the dimension d is equal to two.This theorem is useful as it informs that every joint distribution can be composed by the marginal distributions of the studied variables and a copula function.
An important point that should be mentioned about the copula method is that the lack of independency among the studied variables is mandatory.Kendall's tau is a useful coefficient for calculating the dependency's level [42].Particularly, Kendall's tau estimates the concordance of data [43], meaning that it measures the cases when the two points (x1,y1), (x2,y2) have a positive slope ((x1 − x2)(y1 − y2) > 0).The existence of dependence between the studied parameters is important as the copula families describe this relationship.Consequently, copula families are mathematical functions that describe the joint distribution of the two or more studied parameters.The copula families can be categorized according to their behavior (symmetry, dependence on the upper or lower tail) in two main categories: the Archimedean and Elliptical copulas.Clayton, Gumbel and Frank are some widely used families belonging into the Archimedean category while Normal and Student-t are some of the most commonly used Elliptical ones.In general, Archimedean copulas are exchangeable and easily applied.Furthermore, they present asymmetric extensions.Elliptical copulas are mainly coming from elliptical distributions and are more appropriate for the description of parameters with symmetrical tail dependence.
The present study follows the next steps in order to achieve its aims: • Firstly, the accuracy of both E-OBS and MPI gridded precipitation data was tested.The daily precipitation values of three stations of the Hellenic National Meteorological Service that are located in and around of the study area were used as defaults for the comparisons (Figure 1).Particularly, the daily precipitation of both the closer E-OBS and MPI grid points was compared with the three Greek stations values (Thessaloniki, Chrysoupoli and Alexandroupoli) in order to test if the E-OBS grids can be used as the reference and if the MPI values need correction.

•
To reduce the autocorrelation of the monthly datasets, a sample covering the 75% of the initial values was used for the estimations, eliminating the periodicity of the data.The summary of the followed methodology is presented in Figure 2 while the detailed steps are: Firstly, for every modeled grid point (MPI), the closest continental reference (E-OBS) grid was selected [44][45][46][47] Using only the E-OBS data, the mean daily precipitation per month and the monthly total precipitation were estimated for the E-OBS data.
The dependence of the E-OBS mean daily and total precipitation for each month of the study period was modelled using the copula families.Particularly, more than 15 copula families (Gaussian, Student t, Clayton Gumbel, Frank, Joe, BB1, BB6, BB7, BB8 and their rotated versions) were tested in order to select the one, which can describe the dependence satisfactorily.Table 1 [48] presents the denotation and properties of the tested copula families.
Water 2019, 11, 600 The selection of the best fitted copula was based on the Akaike information criterion (AIC) [49] Apart from the copula family selection, the marginal distributions of the E-OBS monthly totals were estimated.Particularly, the normal, log normal, gamma, GEV and pareto distributions were tested in order to select the one that best fits the E-OBS monthly totals.The selection was made according to the AIC criterion.Before the use of copulas for the MPI bias correction, an evaluation of the used method was made for a different period as it also occurs in Mao et al. [12] who split the tested period into two datasets (fitting and testing).Accordingly, the period of 1981-2000 was used as the fitting period while the period 2001-2010 was used as the tested one.
The MPI mean daily precipitation values were estimated for each month of the study period.
Combining the MPI mean daily precipitation with the selected copula family, a new dataset was estimated.
The new dataset was fitted by the chosen marginal distribution and the bias corrected MPI total monthly precipitation was estimated (Figure 2).

•
The estimated bias corrected dataset was compared with the bias corrected datasets, which came from three other widely used statistical methods of bias correction (delta scaling and empirical quartile mapping).A synoptic analysis and comparison of the three tested bias correction methods is available by Yuanchao Xu [50].

•
Finally, the bias corrected values were used as inputs to a spatially distributed hydrological model.
Archimedean Families Name Function Parameter Range Tail Dependence (Lower, Upper) The rotated by 90 degrees version of the families: The calculations of the present research were conducted using the R programming language [51] and more specifically, the packages of "VineCopula" [52], the "copula" [53] and the "spcopula" [54].
from three other widely used statistical methods of bias correction (delta scaling and empirical quartile mapping).A synoptic analysis and comparison of the three tested bias correction methods is available by Yuanchao Xu [50].

•
Finally, the bias corrected values were used as inputs to a spatially distributed hydrological model.
The calculations of the present research were conducted using the R programming language [51] and more specifically, the packages of "VineCopula" [52], the "copula" [53] and the "spcopula" [54].

Copula Families Elliptical Families Distribution
Parameter Range Tail Dependence Archimedean Families Name Function Parameter Range Tail Dependence (Lower, Upper) The followed method.

Hydrologic Modelling
The three different datasets (i.e., E-OBS, MPI and bias corrected MPI precipitations) were used by the MODSUR (MODelisation du SURface) [55] hydrological model for the simulation of the discharges from the case study area.The model is a spatially distributed model and hence, this is based on a densely spaced grid and uses a progressive quadtree structure with varying cell sizes.The water budget in the MODSUR model is computed in each grid cell using a system of four reservoirs [56] that is responsible for the repartition of rainfall water into runoff, infiltration and evapotranspiration.All the parameters for the reservoirs are defined by the user on the basis of the physical and topographic characteristics of the basin, such as the geology, land cover, land use and topography.Knowing the daily precipitation, temperature and potential evapotranspiration over the basin, the model computes the surface discharges on preselective cells that are attributed the river network.The specific model has been coupled with atmospheric [57] and agricultural [58] models as well as having successfully represented the surface flow in many basins of varying scales [56,59], including the proposed case study basin [60].
As aforementioned, the MODSUR hydrological model has been applied to the Nestos basin with its calibration and validation periods of 1987-1993 (Coefficient of determination R 2 = 0.64) and 1994-1995 (coefficient of determination R 2 = 0.68), respectively [60].In the present research, the parameterization of the model remained exactly the same and the alterations were only focused on the input variables.Regarding the latter, the precipitation and temperature were derived by the E-OBS and MPI datasets while the evapotranspiration was calculated based on the Thornthwaite method utilizing the temperature time series of each dataset.Finally, the gridded variables were superposed over the hydrological model grid, which is composed of 9212 cells, and their nesting in the MODSUR grid was conducted with the use of spatial analyst tools for the geographic information system.

Comparison of Observed Data with E-OBS and MPI
In the first part of the present research, a comparison of both E-OBS and MPI precipitations with the observed precipitation in three Greek stations is attempted.Using daily data series from the three stations and from both E-OBS and MPI datasets, the monthly total precipitations have been estimated and compared.For the MPI and E-OBS precipitations, the grid points that are closest to the three Greek stations, which are not located in the sea, have been selected.
The analysis revealed that for the Thessaloniki station, the MPI model overestimates the observed precipitations after the 50th percentile.We determined that the MPI model overestimates the maximum precipitations by more than 200% while E-OBS datasets underestimate the extremes for only 20 mm (Table 2).Additionally, the correlation between observed and E-OBS datasets is high (0.94) while the correlation index of the observed and MPI precipitations is only 0.35.For the second station (Chrysoupoli), the results are similar.There is an overestimation of the Chrysoupoli precipitation from the MPI model, especially for the heights exceeding the median.However, E-OBS values are close to the observed values and the maximum recorded bias is almost 10 mm.The correlation between the observed and E-OBS data is lower than the respective results in Thessaloniki (0.8) but is still higher (almost double) from the correlation index between the observed and MPI precipitation (0.39).For the station of Alexandroupoli, Table 2 shows that the minimum and 25th percentile of both E-OBS and MPI values are close to the observed data while for the other statistical components (median, mean, 75% and max), the MPI overestimates importantly the real precipitation.On the other hand, E-OBS precipitations satisfactorily approach the real ones.Finally, the correlation between the observed and E-OBS precipitation in Alexandroupoli is equal to 0.97 while the respective value for the observed and MPI data is 0.2. Figure 3 presents the line plots of monthly precipitations for the observed and E-OBS datasets as well as for the observed and MPI datasets.precipitation from the MPI model, especially for the heights exceeding the median.However, E-OBS values are close to the observed values and the maximum recorded bias is almost 10 mm.The correlation between the observed and E-OBS data is lower than the respective results in Thessaloniki (0.8) but is still higher (almost double) from the correlation index between the observed and MPI precipitation (0.39).For the station of Alexandroupoli, Table 2 shows that the minimum and 25th percentile of both E-OBS and MPI values are close to the observed data while for the other statistical components (median, mean, 75% and max), the MPI overestimates importantly the real precipitation.
On the other hand, E-OBS precipitations satisfactorily approach the real ones.Finally, the correlation between the observed and E-OBS precipitation in Alexandroupoli is equal to 0.97 while the respective value for the observed and MPI data is 0.2. Figure 3 presents the line plots of monthly precipitations for the observed and E-OBS datasets as well as for the observed and MPI datasets.According to Figure 3 and Table 2, E-OBS data is significantly closer to the observed values while the MPI precipitation shows important biases.Based on the above results, the E-OBS dataset is used as default as it covers the whole study region and they are similar to the observed ones regarding  According to Figure 3 and Table 2, E-OBS data is significantly closer to the observed values while the MPI precipitation shows important biases.Based on the above results, the E-OBS dataset is used as default as it covers the whole study region and they are similar to the observed ones regarding both the precipitation amount and their variation.

Bias Correction Method Fitting and Evaluation
The E-OBS values for the period of 1981-2000 were used in order to model the dependence between the mean daily and total precipitation per month.For every grid point, the selection of the copula family describing the dependence between mean daily and total precipitation per month was made according to AIC criterion.Figure 4 shows the AIC values for a random grid (red point in Figure 1) in which the Student-t copula was selected as its value was the highest compared with the other tested families.Additionally, the percentiles of the selected copulas as well as for the selected marginal distributions in the grids of the study area are presented in Table 3.The Student-t family was selected for more than the half grids (60%) of the Nestos region.According to Demarta et al. [61] the Student-t copula (e.g., [62]) represents the dependence structure that is implicit in a multivariate t distribution.The second most frequently selected family was the Rotated Gumbel (20%) while the Frank and Gumbel copula were appropriate for the 11% and 9% of the grid points, respectively.As it concerns the marginals, the Weibull distribution was selected in the majority of grid points (87%) while in the grids that the Weibull was not appropriate, the gamma distribution was selected.

Bias Correction Method Fitting and Evaluation
The E-OBS values for the period of 1981-2000 were used in order to model the dependence between the mean daily and total precipitation per month.For every grid point, the selection of the copula family describing the dependence between mean daily and total precipitation per month was made according to AIC criterion.Figure 4 shows the AIC values for a random grid (red point in Figure 1) in which the Student-t copula was selected as its value was the highest compared with the other tested families.Additionally, the percentiles of the selected copulas as well as for the selected marginal distributions in the grids of the study area are presented in Table 3.The Student-t family was selected for more than the half grids (60%) of the Nestos region.According to Demarta et al. [61] the Student-t copula (e.g., [62]) represents the dependence structure that is implicit in a multivariate t distribution.The second most frequently selected family was the Rotated Gumbel (20%) while the Frank and Gumbel copula were appropriate for the 11% and 9% of the grid points, respectively.As it concerns the marginals, the Weibull distribution was selected in the majority of grid points (87%) while in the grids that the Weibull was not appropriate, the gamma distribution was selected.The proposed methodology was applied to a tested time period in order to evaluate its results.The period was split in two time periods: the fitted (1981-2000) and the tested (2001-2010).Figure 5 present the results of the evaluation for two randomly selected grid points (green points in Figure 1, grid 1: triangle, grid 2: circle) and for the mean grid of the study area.According to Figure 5, there are some grids where the bias corrected values approach almost perfectly the E-OBS (grid 1) while there are small deviations in others (grid 2).This fact is also presented by the RMSE values (9.01 and 17.2, respectively).The spatial mean grid combines all the results for the study grids.The results show that there is a small overestimation of the E-OBS data by the bias corrected values for the fitting period, which generally does not exceed the 10 mm.Finally, the RMSE for the mean grid is equal to  The proposed methodology was applied to a tested time period in order to evaluate its results.The period was split in two time periods: the fitted (1981-2000) and the tested (2001-2010).Figure 5 present the results of the evaluation for two randomly selected grid points (green points in Figure 1, grid 1: triangle, grid 2: circle) and for the mean grid of the study area.According to Figure 5, there are some grids where the bias corrected values approach almost perfectly the E-OBS (grid 1) while there are small deviations in others (grid 2).This fact is also presented by the RMSE values (9.01 and 17.2, respectively).The spatial mean grid combines all the results for the study grids.The results show that there is a small overestimation of the E-OBS data by the bias corrected values for the fitting period, which generally does not exceed the 10 mm.Finally, the RMSE for the mean grid is equal to 11.A comparison of the bias corrected precipitation (with copula method) with the bias corrected precipitation derived from three widely used methods (delta, scaling and empirical quantile mapping) was conducted, with the results presented for the mean grid of the study area as shown in Figure 6.According to Figure 6, the bias corrected precipitations with the copula method are closer to the E-OBS data compared with the other three methods.The results from the delta method overestimate the E-OBS values by almost 100 mm while the overestimation for the other two methods is higher.Consequently, the precipitation bias correction in the study area can be achieved satisfactorily following this study's methodology.Additionally, the RMSE for every tested method was estimated and the copula method presents the lower value (23.3).A comparison of the bias corrected precipitation (with copula method) with the bias corrected precipitation derived from three widely used methods (delta, scaling and empirical quantile mapping) was conducted, with the results presented for the mean grid of the study area as shown in Figure 6.According to Figure 6, the bias corrected precipitations with the copula method are closer to the E-OBS data compared with the other three methods.The results from the delta method overestimate the E-OBS values by almost 100 mm while the overestimation for the other two methods is higher.Consequently, the precipitation bias correction in the study area can be achieved satisfactorily following this study's methodology.Additionally, the RMSE for every tested method was estimated and the copula method presents the lower value (23.3).mapping) was conducted, with the results presented for the mean grid of the study area as shown in Figure 6.According to Figure 6, the bias corrected precipitations with the copula method are closer to the E-OBS data compared with the other three methods.The results from the delta method overestimate the E-OBS values by almost 100 mm while the overestimation for the other two methods is higher.Consequently, the precipitation bias correction in the study area can be achieved satisfactorily following this study's methodology.Additionally, the RMSE for every tested method was estimated and the copula method presents the lower value (23.3).

Evaluation of MPI and Corrected MPI Precipitation Data
The produced dataset contains the bias corrected MPI precipitations (C_MPI) and was evaluated temporally and spatially.It should be mentioned that the described procedure takes place for each grid point separately.

Spatial Analysis of Precipitations before and after the Bias Correction
After estimating the mean total precipitation for each grid point using the copula method, the mean total spatial and temporal corrected precipitation for the river basin are estimated.Figure 7 presents the seasonal and the annual boxplots of the MPI precipitation before the correction (MPI) and after the correction (C_MPI) as well as the data from the E-OBS grid points.As illustrated, the bias correction generally improves the precipitation values of the raw MPI data.According to the boxplots, there are grid points of MPI in which the seasonal total precipitation differs significantly to the E-OBS ones.On the other hand, for most of the C_MPI grid points, the seasonal total precipitation is closer to the E-OBS data.
In Figure 7, it is observed that in the winter and autumn, the MPI precipitation is overestimated with respect to E-OBS values by 40-45 mm.On the contrary, the C_MPI dataset is much closer to E-OBS as it presents a small overestimation of 5 mm (autumn) to 10 mm (winter).During the spring, the mean total precipitation amount of the grid points in the Nestos catchment is 40 mm (E-OBS).The MPI's estimation was almost 90 mm while the C_MPI precipitation amount is equal to 45 mm.Additionally, the range of the total precipitation in the MPI dataset is much greater compared to the other two datasets.Some of the MPI grid points have a lower precipitation total than that of other grid points.After the bias correction, the range of C_MPI total precipitation is between 30 and 45 mm, which is very close to the corresponding E-OBS (30 to 40 mm).For the summer, the results are similar but the data range is greater.Particularly, according to the MPI model, there are grids with a total summer precipitation that ranges from 10 to almost 50 mm.The corresponding E-OBS precipitation data vary from 20 to 30 mm while this range changes to 10-20 mm after the correction.The annual values summarizing all the results of the four seasons showed that after the copula bias correction, the total precipitation amount (C_MPI) is much closer to the E-OBS data for the majority of the grid points that cover the study period.
Figure 8 presents the Scatter plots of the E-OBS values with the MPI and C_MPI precipitation, for the four seasons and for the whole studied period (1981)(1982)(1983)(1984)(1985)(1986)(1987)(1988)(1989)(1990)(1991)(1992)(1993)(1994)(1995)(1996)(1997)(1998)(1999)(2000).The C_MPI values are closer to the E-OBS ones, while the MPI precipitations present important differences from the E-OBS respect values.The results in Figure 8 are similar with their respect boxplots (Figure 7).The C_MPI precipitation approach with high accuracy the E-OBS values for winter and autumn while the MPI precipitation differ importantly (more than 100 mm) from the E-OBS.For the seasons of spring and summer, the C_MPI precipitations are much closer to the E-OBS data compared with the MPI precipitation.However, the C_MPI for summer present higher differences from the E-OBS data compared with the other three seasons.Finally, the annual plots show that the E-OBS precipitations are very close to the C_MPI values while the MPI model overestimates them by more than 80 mm.
with respect to E-OBS values by 40-45 mm.On the contrary, the C_MPI dataset is much closer to E-OBS as it presents a small overestimation of 5 mm (autumn) to 10 mm (winter).During the spring, the mean total precipitation amount of the grid points in the Nestos catchment is 40 mm (E-OBS).The MPI's estimation was almost 90 mm while the C_MPI precipitation amount is equal to 45 mm.Additionally, the range of the total precipitation in the MPI dataset is much greater compared to the other two datasets.Some of the MPI grid points have a lower precipitation total than that of other grid points.After the bias correction, the range of C_MPI total precipitation is between 30 and 45 mm, which is very close to the corresponding E-OBS (30 to 40 mm).For the summer, the results are similar but the data range is greater.Particularly, according to the MPI model, there are grids with a total summer precipitation that ranges from 10 to almost 50 mm.The corresponding E-OBS precipitation data vary from 20 to 30 mm while this range changes to 10-20 mm after the correction.The annual values summarizing all the results of the four seasons showed that after the copula bias correction, the total precipitation amount (C_MPI) is much closer to the E-OBS data for the majority of the grid points that cover the study period.(1981)(1982)(1983)(1984)(1985)(1986)(1987)(1988)(1989)(1990)(1991)(1992)(1993)(1994)(1995)(1996)(1997)(1998)(1999)(2000).The boxplots created after the estimation of the mean seasonal values for each grid.
Figure 8 presents the Scatter plots of the E-OBS values with the MPI and C_MPI precipitation, for the four seasons and for the whole studied period (1981)(1982)(1983)(1984)(1985)(1986)(1987)(1988)(1989)(1990)(1991)(1992)(1993)(1994)(1995)(1996)(1997)(1998)(1999)(2000).The C_MPI values are closer to the E-OBS ones, while the MPI precipitations present important differences from the E-OBS respect values.The results in Figure 8 are similar with their respect boxplots (Figure 7).The C_MPI precipitation approach with high accuracy the E-OBS values for winter and autumn while the MPI precipitation differ importantly (more than 100 mm) from the E-OBS.For the seasons of spring and summer, the C_MPI precipitations are much closer to the E-OBS data compared with the MPI precipitation.However, the C_MPI for summer present higher differences from the E-OBS data compared with the other three seasons.Finally, the annual plots show that the E-OBS precipitations are very close to the C_MPI values while the MPI model overestimates them by more than 80 mm.Apart from the spatial analysis of the results (Section 3.3.1.),a temporal evaluation was also performed, the results of which are presented in Table 4 and Figure 9. Particularly, Figure 4 illustrates

Temporal Analysis of Precipitations before and after the Bias Correction
Apart from the spatial analysis of the results (Section 3.3.1),a temporal evaluation was also performed, the results of which are presented in Table 4 and Figure 9. Particularly, Figure 4 illustrates the line plots of the data series for the whole study region and for the three datasets (MPI, E-OBS, and C_MPI) while Table 2 statistically describes these datasets.The results derived after the calculation of the monthly total precipitation amount for the whole region (all grids).Table 4. Statistical analysis of the results temporally.For the estimation found in the table, the mean precipitation of all grid points monthly was firstly calculated before the statistical components of these data series were used for the table.

Assessment of MPI and Corrected MPI Precipitation on River Discharges
The hydrological simulation outputs are consistent with the results of the precipitation analysis.This means that the simulated discharges are overestimated and they do not adequately represent (coefficient of efficiency R 2 = 0.38) the river flows when the hydrological model is driven by the MPI non-corrected data as shown in Figure 10.On the other hand, the simulated discharges are well correlated (R 2 = 0.56) with the observed discharges when the MPI copula corrected datasets are used as the input.When the E-OBS data were used to trigger the hydrological simulation, the coefficient of efficiency was calculated to be 0.65.Moreover, while the observed mean flow for the time period of 1981-1995 is equal to 20.4 m 3 /s, the MPI non-corrected data discharges and the MPI copula corrected discharges are equal to 28.1 m 3 /s and 16.2 m 3 /s, respectively.The overestimation of the discharges in the case of the MPI non-corrected data is more obvious in the maximum flows where the corresponding value of 167.3 m 3 /s is more than double than the observed maximum flows of 74.0 m 3 /s and more than triple that of the MPI copula corrected maximum value of 53.7 m 3 /s.The increased frequency of maxima in the case of the MPI non-corrected data is also obvious in Figure 10.Moreover, the results show that RMSE is equal to 27.48, 13.34 and 11.98 for the MPI non-corrected, copula corrected and E-OBS datasets, respectively.This means that the implemented technique resulted in a reduced RMSE of the simulations triggered by the copula that is corrected in comparison to those of MPI non-corrected data.
A further analysis of the results demonstrated that: (i) the simulations driven by the MPI noncorrected data present larger bias when the river flows exceed 10.0 m 3 /s.while in the case of copula  According to Figure 9, the MPI data series overestimates the E-OBS equivalent series during the whole period while the C_MPI values satisfactorily simulate E-OBS precipitation.Moreover, a clear overestimation for all of the statistical amounts should be highlighted according to the results in Table 4. Particularly, the MPI minimum values are similar to the E-OBS ones as in the case of the C_MPI dataset.However, after the first quartile of the data series, the overestimation is much higher while the bias corrected values are close to the E-OBS.It can be concluded that MPI's median is almost double that of the E-OBS and the C_MPI values.This was observed for both the mean and 75% quartile.For the maximum values, MPI exceeds the E-OBS respective value by almost 150 mm while the MPI corrected values underestimate it by 40 mm.

Assessment of MPI and Corrected MPI Precipitation on River Discharges
The hydrological simulation outputs are consistent with the results of the precipitation analysis.This means that the simulated discharges are overestimated and they do not adequately represent (coefficient of efficiency R 2 = 0.38) the river flows when the hydrological model is driven by the MPI non-corrected data as shown in Figure 10.On the other hand, the simulated discharges are well correlated (R 2 = 0.56) with the observed discharges when the MPI copula corrected datasets are used as the input.When the E-OBS data were used to trigger the hydrological simulation, the coefficient of efficiency was calculated to be 0.65.
Water 2019, 11, 600 13 of 19 in Figure 11); (ii) MPI triggered simulations present the larger bias in winter (December-January-February (DJF)) and in spring (March-April-May (MAM)) than in summer (June-July-August (JJA)) and autumn (September-October-November (SON)) (the middle scatterplots in Figure 11); (iii) the copula corrected simulations are remarkably correlated with the observed ones in SON and DJF time periods; (iv) for the just aforementioned time periods, the E-OBS and copula corrected simulations present almost identical results; and (v) during summer periods, the copula corrected data do not successfully represent the actual low flows and the biases are more obvious.Moreover, while the observed mean flow for the time period of 1981-1995 is equal to 20.4 m 3 /s, the MPI non-corrected data discharges and the MPI copula corrected discharges are equal to 28.1 m 3 /s and 16.2 m 3 /s, respectively.The overestimation of the discharges in the case of the MPI non-corrected data is more obvious in the maximum flows where the corresponding value of 167.3 m 3 /s is more than double than the observed maximum flows of 74.0 m 3 /s and more than triple that of the MPI copula corrected maximum value of 53.7 m 3 /s.The increased frequency of maxima in the case of the MPI non-corrected data is also obvious in Figure 10.Moreover, the results show that RMSE is equal to 27.48, 13.34 and 11.98 for the MPI non-corrected, copula corrected and E-OBS datasets, respectively.This means that the implemented technique resulted in a reduced RMSE of the simulations triggered by the copula that is corrected in comparison to those of MPI non-corrected data.
A further analysis of the results demonstrated that: (i) the simulations driven by the MPI non-corrected data present larger bias when the river flows exceed 10.0 m 3 /s.while in the case of copula corrected data, the relatively small biases occur in flows larger than 25.0 m 3 /s (the upper scatterplot in Figure 11); (ii) MPI triggered simulations present the larger bias in winter (December-January-February (DJF)) and in spring (March-April-May (MAM)) than in summer (June-July-August (JJA)) and autumn (September-October-November (SON)) (the middle scatterplots in Figure 11); (iii) the copula corrected simulations are remarkably correlated with the observed ones in SON and DJF time periods; (iv) for the just aforementioned time periods, the E-OBS and copula corrected simulations present almost identical results; and (v) during summer periods, the copula corrected data do not successfully represent the actual low flows and the biases are more obvious.

Discussion and Conclusions
The proposed research demonstrates a methodological approach to precipitation biascorrection, which is derived from regional climate models.For that purpose, the precipitation data from regional climate models that correspond to a historic time period has been statistically improved with the use of E-OBS data.In the present study, E-OBS data have been used as defaults after a comparison with the observed precipitations of three Greek stations, which showed that the correlation differences were low and insignificant.In accordance with the present study, several researchers use E-OBS data as defaults [63,64] as they offer large spatial and temporal convergence, they are updated [65] and they present a lower degree of uncertainty in comparison to meteorological observations [66].The spatial location of the case study area could also play an important role in determining the accuracy of the E-OBS since the bias of the E-OBS dataset with data from a highdensity network of stations in the Czech Republic has demonstrated limitations with respect to its applicability for evaluating RCMs [67].
During the last few years, numerus studies proposed different methods for bias correction.However, the studies using the copula method for the bias correction of climate model's data using reanalysis or E-OBS data are limited [12].In 2015, Mao et al. used daily precipitation values from two datasets (reanalysis and WRF model) in order to model their relationship and to use it for the bias correction of the WRF model.Following Mao et al. [12], the present study tries to correct the biases of the MPI climate model using E-OBS data but on a monthly basis.Specifically, the present study models the dependence of the mean daily precipitation per month and total monthly precipitation

Discussion and Conclusions
The proposed research demonstrates a methodological approach to precipitation bias-correction, which is derived from regional climate models.For that purpose, the precipitation data from regional climate models that correspond to a historic time period has been statistically improved with the use of E-OBS data.In the present study, E-OBS data have been used as defaults after a comparison with the observed precipitations of three Greek stations, which showed that the correlation differences were low and insignificant.In accordance with the present study, several researchers use E-OBS data as defaults [63,64] as they offer large spatial and temporal convergence, they are updated [65] and they present a lower degree of uncertainty in comparison to meteorological observations [66].The spatial location of the case study area could also play an important role in determining the accuracy of the E-OBS since the bias of the E-OBS dataset with data from a high-density network of stations in the Czech Republic has demonstrated limitations with respect to its applicability for evaluating RCMs [67].
During the last few years, numerus studies proposed different methods for bias correction.However, the studies using the copula method for the bias correction of climate model's data using reanalysis or E-OBS data are limited [12].In 2015, Mao et al. used daily precipitation values from two datasets (reanalysis and WRF model) in order to model their relationship and to use it for the bias correction of the WRF model.Following Mao et al. [12], the present study tries to correct the biases of the MPI climate model using E-OBS data but on a monthly basis.Specifically, the present study models the dependence of the mean daily precipitation per month and total monthly precipitation for the E-OBS data before it applies the estimated function to the MPI data.Additionally, the present study evaluates more than 15 copula families while the majority of existing studies examine a smaller number.In accordance with Mao et al. [12], the bias correction results are closer to the E-OBS data for all seasons.
One of the main scopes of the study is also the use of bias corrected values that were calculated by the copula method as inputs in a hydrological model as the copula results have not been used in a hydrological model until now.The applied methodology is successfully verified by simulating a transboundary catchment and assessing the simulated discharges with gauged measurements.Particularly, the hydrologic simulation outputs showed that the copula method can be used for the bias correction of the MPI climate model's precipitations in the studied region.The satisfactory results of the copula method are due to its ability to describe mathematically the dependence between two or more variables [32] and therefore, it has been used by several scientists [21,33].For the studied area, the E-OBS precipitation were overestimated by the MPI values while after the bias correction, the new values are much closer to the E-OBS with a small underestimation during summer.
The bias correcting of climate change datasets is important since the overestimations or sub estimations of the under-simulation variables are presented in the cases where the raw climate change data grids are used as the input data for the hydrological models.Although recent methods used to bias correct precipitation from RCMs cannot overcome the limitations of the RCMs in simulating the precipitation sequence, which affects runoff generation [68], the uncertainty is being reduced when the statistical methods are applied.The correction of bias between climate model's historic data and reanalysis or E-OBS datasets has been thoroughly examined in the literature [69,70].The copula method has also been proposed for a dynamical bias correction of temperature and precipitation [19] as well as for bias correction of precipitations of the WRF model [12] based on reanalysis data.Nevertheless, the present study proposes the copula method for bias correction using the mean daily precipitation and the total precipitation of each month.
In terms of the E-OBS data and their interaction with hydrologic simulations, Photiadou et al. [71] used the aforementioned datasets to force a hydrological model, evaluating the annual mean and extreme discharges compared to observed discharges.After this, they aimed to assess the bias of the induced precipitation data.The output of the research showed that E-OBS performed well with respect to extreme discharge.E-OBS data are also routinely used as a reference for the bias correction of satellite-based precipitation products before the latter, using forcing conditions to hydrological models [72] while Ledesma and Futter [73] propose that the particular datasets could be more widely used as an alternative input to rainfall-runoff models, even in the cases where instrumental measurements are available.
To sum up, bias corrected climate model outputs are significant when investigating the impact of climate change on hydrological regimes.The areas where no monitoring facilities have been installed are usually isolated places and thus, the assessment of climate change in such places could provide outputs without requiring direct mankind interventions in the ecosystem.
The present study proposes a bias corrected method for the total monthly precipitation in a Bulgarian-Greek basin based on copula.Additionally, the evaluation of the method was achieved using the bias corrected data in a hydrological model for the aforementioned area.The dependence between the mean daily precipitation per month and the monthly total precipitation for the case study region was modeled in the present study and can be used for the model bias correction for future periods.Consequently, further work should evaluate the proposed methodology and the results for several different future periods, regions and climate models.

Figure 1 .
Figure 1.The studied area, which was the transboundary catchment of Mesta/Nestos.

Figure 3 .
Figure 3.Comparison of the observed, E-OBS and MPI precipitation for the three Greek stations.Figure 3. Comparison of the observed, E-OBS and MPI precipitation for the three Greek stations.

Figure 3 .
Figure 3.Comparison of the observed, E-OBS and MPI precipitation for the three Greek stations.Figure 3. Comparison of the observed, E-OBS and MPI precipitation for the three Greek stations.

Figure 4 .
Figure 4.The AIC values for a random grid.

Figure 4 .
Figure 4.The AIC values for a random grid.

Figure 5 .
Figure 5. Evaluation of the bias correction method for a period of ten years (2001-2010).

Figure 6 .
Figure 6.Comparison of the copula method with three other bias correction methods (delta, scaling, empirical quantile mapping).

Figure 7 . 19 Figure 7 .
Figure 7. Boxplots of the three studied datasets for the four seasons and for the whole studied period (1981-2000).The boxplots created after the estimation of the mean seasonal values for each grid.

Figure 8 .
Figure 8. Scatterplots of the three studied datasets for the four seasons and for the whole studied period (1981-2000).The scatter plots created after the estimation of the mean seasonal values for each grid.3.3.2.Temporal Analysis of Precipitations before and after the Bias Correction

Figure 8 .
Figure 8. Scatterplots of the three studied datasets for the four seasons and for the whole studied period (1981-2000).The scatter plots created after the estimation of the mean seasonal values for each grid.

Figure 9 .
Figure 9. Line plots of the data series for the whole studied region and for the three datasets (MPI, E-OBS, and C_MPI).

Figure 9 .
Figure 9. Line plots of the data series for the whole studied region and for the three datasets (MPI, E-OBS, and C_MPI).

Figure 10 .
Figure 10.Hydrological model simulated discharges at the borders of the two countries that share the basin.Observations (hashed area) and E-OBS simulation results (red dotted curve) versus MPI-copula corrected discharges (black curve) and MPI-no corrected discharges (grey dashed curve).

Figure 10 .
Figure 10.Hydrological model simulated discharges at the borders of the two countries that share the basin.Observations (hashed area) and E-OBS simulation results (red dotted curve) versus MPI-copula corrected discharges (black curve) and MPI-no corrected discharges (grey dashed curve).

Figure 10 . 19 Figure 11 .
Figure 10.Hydrological model simulated discharges at the borders of the two countries that share the basin.Observations (hashed area) and E-OBS simulation results (red dotted curve) versus MPI-copula corrected discharges (black curve) and MPI-no corrected discharges (grey dashed curve).

Figure 11 .
Figure 11.Simulated versus observed discharges at interannual and seasonal time scales.

Table 1 .
Denotation and properties of Elliptical and Archimedean copula families.

Table 1 .
Denotation and properties of Elliptical and Archimedean copula families.

Table 2 .
Statistical results for observed, E-OBS and MPI datasets in mm.The correlation of E-OBS and MPI was calculated between these datasets and observed datasets.

Table 2 .
Statistical results for observed, E-OBS and MPI datasets in mm.The correlation of E-OBS and MPI was calculated between these datasets and observed datasets.

Table 3 .
The percentiles of the selected copula families and marginal distributions.

Table 3 .
The percentiles of the selected copula families and marginal distributions.

Table 4 .
Statistical analysis of the results temporally.For the estimation found in the table, the mean precipitation of all grid points monthly was firstly calculated before the statistical components of these data series were used for the table.