1 Introduction

Local climate risk information is an essential component of an urban’s agglomeration comprehensive framework for responding to the risks of climate change and to its implementation of adaptation strategies (Bader et al. 2018). However, the available state-of-the-art climate change modeling experiments have a horizontal resolution of about 12 km (i.e. EURO-CORDEX) which can provide the climatological information needed to enable climate change impact studies and develop adaptation strategies on a national or regional level (Jacob et al. 2020a, b), but is still coarse to provide the adequate information at the local scale.

To overcome this drawback, while the latest “convection-permitting” kilometer-scale regional climate model (RCM) scenario simulations are in progress (Coppola et al. 2020; Pichelli et al. 2021), two different downscaling approaches are reported in the literature in order to generate high‐resolution regional climate models based on the large‐scale information (i) dynamical downscaling and (ii) statistical downscaling.

Dynamical downscaling usually refers to the use of RCMs, driven by General Circulation Models output or reanalysis data to produce regionalized climate information. Examples of dynamical downscaling in European cities can be found in the studies of Lauwaet et al. (2015; 2016). In particular, Lauwaet et al. (2015) used an urban boundary layer climate model (UrbClim) coupled to 11 global climate models contained in the Coupled Model Intercomparison Project 5 archive. The authors conducted 20-year simulations for the present (1986–2005) and future (2081–2100) climate conditions, considering the Representative Concentration Pathway (RCP) 8.5 climate scenario to examine the evolution of the urban heat island (UHI) in six European cities with a horizontal resolution lower than 1 km. UrbClim was also used to assess the climate change impact on UHI in Brussels between 2000–2009 and 2060–2069 under RCP4.5 and RCP8.5 climate scenarios at a resolution of 250 m (Lauwaet et al. 2016).

On the contrary, statistical downscaling is based on the use of empirical relationships that link local observations of a target variable (predictands) to a set of suitable large-scale variables representing the state of the atmosphere (predictors). The term statistical downscaling includes different methodologies such as analogue methods, weather types techniques, multiple linear regression, perfect prognosis and model output statistics (Smid and Costa 2018; Gutierrez et al. 2019). Recently, Manzanas et al. (2020) identified bias adjustment techniques as a good alternative for statistical downscaling techniques for climate change impact studies. Bias adjustment is the process of scaling raw climate model outputs to account for their systematic errors, in order to improve their fitting to observations. The aforementioned process can be achieved with various methods ranging from simple scaling of the mean and/or variance to more complex methods that scale the whole distribution (e.g. Casanuneva et al. 2020; Lazoglou et al. 2020, 2021). Furthermore, bias adjustment techniques can be categorized as univariate and multivariate. In univariate techniques, each variable is adjusted independently thus, intervariable dependence structure is ignored. In multivariate techniques on the other hand, two or more variables are jointly adjusted and the aforementioned structure is preserved (Cannon 2016; 2018).

Each approach has its strengths and limitations. In particular dynamical downscaling is based on climate models information with physical principles implemented in their code to reproduce local climates, but is computationally intensive. Statistical downscaling and bias adjustment are less computationally intensive and can be easily applied but they rely on the availability of high-quality long term observational records for a number of variables, both from station data and/or gridded observational datasets to establish robust predictor-predictant relationships and to adjust the raw model output, respectively (Casanuneva et al. 2020; Manzanas et al. 2020; Varotsos et al. 2021).

Regarding gridded observational datasets, E-OBS is considered the state-of-the art dataset for the European domain with a resolution of 0.1° (Cornes et al. 2018). In addition, regional/national gridded datasets with similar or higher resolutions have been compiled recently in Europe such as SAFRAN covering France, Spain and the Valearic islands (hourly dataset at 8 km grid spacing, Quintana-Seguí et al. 2017); Iberia01 for the Iberian Penisnula (daily datatest at 0.1°, Herrera et al. 2019a, b); PTHBV for Sweden (4 km daily dataset, Johansson, 2000); HYRAS (Frick et al. 2014; Razafimaharo et al. 2020) and Krähenmann et al. (2018) for Germany (daily dataset at 5 km grid spacing and hourly dataset at 1 km, respectively); seNorge2 for Norway (daily dataset with a 1 km resolution, Lussana et al. 2018); TabsD (MeteoSwiss, 2013a) and RhiresD (MeteoSwiss, 2013b) at 2 km for Switzerland; CY-OBS, a 1 km grid for Cyprus (Camera et al. 2014); CARPATCLIM, a 0.1° grid covering parts of nine countries along the Carpathian Mountains (Lakatos et al. 2013); a 0.11° grid for Poland (Herrera et al. 2019b); a ~ 1 km grid for Portugal (Fonseca and Santos, 2018) and a 1 km grid for Serbia (Sekulić et al. 2021). The vast majority of the aforementioned datasets are produced using dense networks of stations over the areas of interest whereas other studies have used a combination of station data with satellite and/or reanalysis data in areas with low density of stations (Doblas-Reyes et al. 2021 and references therein).

In this study, we present a framework in order to provide high resolution climate change projections, in the order of 1 km, in the Greater Athens Area (GAA) by constructing a gridded daily observational dataset for the period 1981–2000 which is then used as a reference dataset to statistically downscale RCM simulation data from the EURO-CORDEX initiative. To our knowledge this is the first study that presents both a gridded daily dataset for temperature and precipitation, as well climate change projections in the GAA under such a high horizontal resolution. Regarding the former, Mamara et al. (2017) and Gofa et al. (2019) presented a high resolution gridded dataset (about 1 km) for temperatures and precipitation, respectively for the period 1971–2000 based on mean monthly observational data from a number of stations within Greece while for the latter Georgoulias et al. (2022) used RCMs with a horizontal resolution of about 12 km with the future period simulations driven by RCP2.6, RCP4.5 and RCP8.5.

2 Data

2.1 Observational data

In the present work, we use available daily maximum (TX), daily mean (TG) and daily minimum (TN) temperatures as well as daily accumulated precipitation (RR) from 11 stations located within the GAA (Fig. 1) for the period 1981–2000. The observational data for the 10 meteorological stations (Fig. 1), namely: Tanagra (TAN, north of Attika in Viotia prefecture), Tatoi (TAT, northern Attika), Eleusis (ELE, western Attika), Nea Filadelfeia (NFIL, central Attika), Marathonas (MAR, northeastern Attika), Paiania (PAIA, eastern Attika), El. Venizelos airport (AIRP, eastern Attika), Helliniko (HEL, southern Attika), Piraeus (PIR, southwestern Attika), Aegina island (AEG, Saronic Gulf, southwestern Attika) are provided from the study of Soulis et al. (2018), from a dataset of daily observations of 140 stations in Greece, originally obtained from the Hellenic National Meteorological Service (HNMS) for the period 1971–2000. Regarding the last station, NOA (city of Athens) the daily observational data are provided, for the same period, by the Historical meteorological station of the National Observatory of Athens in Thissio (NOA,https://data.climpact.gr/en/dataset/1ce5d2ce-23df-412e-849e-ef2493319da9).

Fig. 1
figure 1

Locations of the stations with their altitude in parenthesis. In the background the elevation of the 3rd nested domain, as provided by the WRF simulation

2.2 WRF simulation

Moreover we apply the non-hydrostatic Advanced Research WRF (ARW) meso-scale meteorological model (Skamarock et al. 2008) version v3.9.1 for the high resolution (1 km 1 km) numerical simulation of the year 1995, using a two-way nesting approach. The Noah-Multiparameterization Land Surface Model (Noah-MP LSM) (Niu et al. 2011), coupled with the single-layer urban canopy model UCM (Tewari et al. 2007 and references therein) is used. The modified two-stream radiation transfer scheme (Yang and Friedl, 2003; Niu and Yang, 2004) is opted in the Noah-MP LSM. The Yonsei University (YSU) scheme (Hong et al. 2006) is applied for the planetary boundary layer parameterization and the Rapid Radiative Transfer Model for General Circulation Models (RRTM-GCM) (Iacono et al. 2008) is used for both long and short wave radiation parameterization. The revised MM5 similarity theory (Jiménez et al. 2012) is applied for the surface layer parameterization and the single-moment microphysics scheme (WSM) 3-class (Hong et al. 2004) for microphysics parameterization.

Regarding land use and soil types, we use the predefined datasets of Moderate Resolution Imaging Spectroradiometer (MODIS) with 21 land use classes and 16 soil categories. The numerical simulation is supplemented by high-resolution data on vegetation and urban land use derived from satellite image analysis. In particular, the land use of the inner and finer domain, covering the GAA, is remapped by extracting vegetation and urban-land use data from a 30-m resolution satellite image (Landsat Thematic Mapper) following the methodology described by Papangelis et al. (2012). Furthermore, a categorization of seven different urban land use categories is applied primarily based on building height and built-up density. Both codes of the Noah-MP and the single-layer UCM are adapted to account for the aggregated land-use information of the study area.

The initial, lateral and boundary conditions for the simulation of the year 1995 are obtained from the ERA-INTERIM Project (ds627.0), with spatial resolution ranging from 1.125° × 1.125° up to 0.703° 0.703° and are updated every 6 h. The WRF numerical simulation is performed by applying a triple nesting, leading to a high-resolution grid of 1 km x 1 km, which covers the GAA (Fig. S1 in Supplementary material). In the vertical axis, 38 full sigma levels resolve the atmosphere (model top at 50 hPa), with a finer resolution near the surface and above the urban-canopy layer.

As far as the year selected for the WRF simulation is concerned, the following criteria were examined: (a) the mean monthly annual cycle of the selected year should exhibit lowest deviations compared to the long term mean monthly annual cycle over the period 1981–2000 (mean absolute deviation from the climatological values calculated in a way similar to the Mean Absolute Error) for all variables examined and (b) the year should preferably be after the year 1990 since we wanted to implant in the WRF simulation the majority of the land use changes occurred in the area under study during the second half of the period 1981–2000. Regarding the first criteria, in Table S1 (Supplementary Material) the deviations, averaged for all available station data, between the mean (total for precipitation) monthly annual cycle for each individual year and the average mean (total for precipitation) monthly annual cycle over the period 1981–2000 are shown. The results indicate that the year 1995 exhibits for all temperature variables deviations lower than 1 °C and the second lowest deviation for precipitation (14 mm) among the years after the year 1990.

2.3 Regional climate models simulations

We employ daily maximum, minimum and mean temperatures as well as daily precipitation data from a transient realisation of the Regional Climate Model (RCM) RCA4 of the Swedish Meteorological and Hydrological Institute (SMHI; Strandberg et al. 2014 and references therein) driven by the Max Planck Institute for Meteorology model MPI-ESM-LR (Popke et al. 2013) with the simulations carried out in the framework of the EURO-CORDEX modeling experiment (http://www.euro-cordex.net). The horizontal resolution of the model is 0.11° (~ 12 km), while the simulated data in this study cover two periods: the 1981–2000 which is used as the reference period and the period 2081–2100 under RCP8.5 (Riahi et al. 2011). It should be mentioned that the specific climate change modelling system has been evaluated and used in previous studies examining climate change impacts in the area under study (Founda et al. 2019; van der Schriek et al. 2020).

3 Methods

3.1 Construction of the GAA high resolution gridded data set (GAA.HRES)

A crucial part in this study is the construction of a high resolution 1 km 1 km daily gridded data set for temperature and precipitation, hereafter named GAA.HRES, for the period 1981–2000 which is used as a reference dataset to bias-adjust the regional climate model simulations. To this aim, regarding the temperature variables, we blend the available observations with the WRF output through gridding techniques that have been used in previous state-of-the art daily gridded datasets in Europe. As it is shown in the forthcoming sections, this blending enhances the spatial variabilty of the gridded product in the complex orography of Attica and it was selected instead of using only the observational data, due to the low density of stations located in relatively low altitudes.

3.1.1 Spatiotemporal model for temperature variables

The steps to obtain the daily gridded values for temperatures are as follows:

Step 1: The mean monthly annual cycles over period 1981–2000 are calculated from the observations for each one of the temperature variables (M.o in Fig. 2). Consequently, the monthly means for the closest WRF grid point to the station location are calculated for the year 1995 (M.w in Fig. 2).

Fig. 2
figure 2

Flowchart of the methodology used for the construction of the high-resolution daily mean, maximum, and minimum temperature datasets over Attica (GAA.HRES). M.o indicates the observed mean monthly annual cycles over period 1981–2000, M.w the monthly means for the closest WRF grid point to the station locations, M.w* the perturbed WRF output, OK indicates Ordinary Kriging, KED indicates Kriging with External Drift while BA indicates Bias Adjustment

Step 2: The biases between the observed monthly means and the WRF ones (M.o – M.w) are interpolated on the 1 km WRF grid using Ordinary Kriging (OK) and then added to the WRF monthly means output. This local perturbation is performed in order to obtain a WRF product with long term climate temporal characteristics (M.w*) which will still maintain the spatial variability of WRF (Cornes et al., 2018; Fonseca and Santos 2018) since previous studies have shown that WRF exhibits a very good skill in capturing the spatial variability of temperature variables over the complex terrain of GAA (Papangelis et al. 2012; Giannaros et al. 2013; Politi et al. 2020; Dandou et al. 2017; 2021).

Step 3: The third step involves the interpolation of the observed temperatures to the WRF grid following a two-step process: (a) Kriging with External Drift (KED) is applied to the monthly values (12 values/year × 20 years) considering the elevation as a covariate to account for elevation dependencies (station’s elevation for modeling and WRF elevation for interpolation) and (b) the daily anomalies (from the monthly mean values, 365 values/year × 20 years) of only the observational data are interpolated by applying OK. The two interpolated products are superimposed by addition to obtain the interpolated daily values for all temperature variables, that is for each daily anomaly the corresponding mean monthly value is added. The aforementioned interpolation technique is similar to the ones used in other gridded daily datasets in Europe such as the early versions of E-OBS (version 16 and prior versions, Haylock et al. 2008) and Iberia01 (Herrera et al. 2019a, b).

Step 4: The final step to obtain the GAA.HRES daily temperatures is to transfer WRF spatial variability to the gridded daily datasets produced in the previous step based solely on the observations. This is achieved by using the unbiasing bias adjustment method (Déqué, 2007, BA Fig. 2). In particular, the mean monthly differences between the perturbed WRF simulation (M.w*) and the mean monthly observed values for the 20-year period (1981–2000) for each grid point are calculated. Consequently, to obtain the final daily gridded temperature data, the mean monthly differences are added to the daily gridded values. This method maintains the absolute trend as well as the temporal variability of the gridded data produced in Step 3 at all timescales (Hempel et al. 2013).

3.1.2 Spatiotemporal model for precipitation

For the daily precipitation dataset the method is based solely on the station data since the WRF simulation exhibited considerable biases when compared to the station data (not shown).The steps to obtain the daily gridded values for precipitation are as follows:

Step 1: The first step involves the interpolation of the observe monthly sums to the WRF grid using (KED), which is applied to the monthly values (12 values/year × 20 years) considering the elevation as a covariate (station’s elevation for modeling and WRF elevation for interpolation).

Step 2: For the daily precipitation the gridding technique is based on Generalized Additive Models (GAM; Wood, 2017). GAMs are a semi-parametric extension of Generalized Linear Models assuming that the underlying functions are additive and that the components are smooth, with the strength of GAMs being their ability to deal with highly non-linear and non-monotonic relationships between the response and the set of explanatory variables (Guisan et al. 2002). It should be mentioned that GAMs have been used in recently produced daily gridded datasets e.g. in Oregon USA (Parmentier et al. 2014) and in Europe (E-OBSv17 and later versions Cornes et al. 2018).

In particular, for the wet days daily accumulated precipitation amounts (RR), model fitting using GAM has the following form:

$$y \, = \, f_{1} \left( {{\text{lon}},{\text{ lat}}} \right) \, + \, f_{2} \left( {bg} \right) \, + e$$

where daily precipitation is modelled as smooth functions of (i) longitude (lon) and latitude (lat) and (ii) monthly background precipitation totals (bg) which are produced by Kriging with External Drift (KED) in Step 1. It should be noted that other variables such as elevation and distance from the coast were also examined however, no significant dependence between these variables and precipitation was found in our experiments on the daily timescale.

Step 3: The daily occurrence of rainfall (0 or 1 depending on whether RR > 0.1 mm) is also interpolated applying Indicator Kriging (IK) considering a threshold of 0.5 for assigning a wet day to a grid point (Herrera et al., 2019a, b). To obtain the final daily gridded dataset for precipitation, the two interpolated daily precipitation products obtained from Step 2 and Step 3 are superimposed by multiplication.

3.2 Statistical downscaling of regional climate model simulations

To downscale the RCM daily model variables to the GAA.HRES, grid a two-step post processing procedure is followed. Initially the RCM daily data are remapped on the GAA.HRES grid using bilinear interpolation and consequently the models’ output is bias adjusted using GAA.HRES as the reference dataset. This two-step approach is the reversed order of the bias correction and spatial disaggregation framework, which has been previously used to statistically downscale both regional and global models for climate change and seasonal forecast studies,respectively (Nilsen et al. 2022; Lorenz et al. 2021). In order to examine if the selection of bias adjustment method is important, we use two different bias adjustment techniques namely empirical quantile mapping (EQM) and quantile delta mapping (QDM). EQM works by constructing a transfer function calibrated over the reference period to map quantiles from the empirical cumulative distribution function of the model output onto the corresponding observed distribution (Iturbide et al. 2019; Casanueva et al. 2020). QDM works initially by detrending the quantiles of the RCM projections, afterwards quantile mapping is applied to the detrended series with the transfer function constructed in the reference period and finally the projected trends, absolute (for temperature) or relative (for precipitation), are added or multiplied to the bias adjusted quantiles (Cannon et al. 2015).

In the current study, average annual and seasonal means are examined for all temperature variables (sums for precipitation) as well as selected indices from the Expert Team on Climate Change Detection and Indices (ETCCDI) (Zhang et al. 2011). These are the number of days with daily TX > 25 °C (SU) and TX > 35 °C (SU35) for TX, the number of days with daily TN > 20 °C (TR) for TN and the number of days with RR > 1 mm (RR1) as well as the maximum monthly one day precipitation amount (RX1day) for RR. These indices are selected taking into account the main climatic characteristics of the examined area which exhibits Mediterranean-type conditions, with mild winters and warm to hot summers. In addition, the confidence ranges of the seasonal means (sums for precipitation) and the aforementioned indices calculated by the 95th percentiles confidence intervals as derived by bootstrap are calculated. More specifically, the confidence range (CR) is used as a measure of uncertainty where for the average value obtained for each grid box a \(\pm a\) value can be added or subtracted in order to get the upper and lower limits, respectively. To obtain α at the grid point scale, the parameter samples (20 values except for winter where the number of values is 19) are bootstapped 1,000 times with replacement. In each resample, the mean of each sample is calculated and the 95th percentile confidence intervals are then computed from the resulting series with α being the difference between the upper and the lower confidence intervals divided by two (Giannakopoulos et al. 2009, 2011, 2016).

4 Results and discussion

4.1 GAA.HRES average annual patterns and evaluation against observations over the period 1981–2000

In Figs. 3d, 4d and 5d the results of the average annual GAA.HRES patterns for daily TX, daily TN and daily TG are shown over the period 1981–2000, respectively as well as the results for the intermediate Steps 2 (Figs. 3b, 4b, 5b) and 3 (Figs. 3c, 4c, 5c) described in Sect. 3.1.1 In Figs. 3a, 4a and 5a the average annual WRF for daily TX, daily TN and daily TG are shown for the year 1995. In Fig. 6b the GAA.HRES average annual total precipitation is shown while in Fig. 6a the results for the intermediate Step 1.

Fig. 3
figure 3

Average annual TX a of the WRF simulation (1995), b of the WRF simulation after step 2 of the methodology for temperatures, c of the daily gridded dataset based only on observations and d of the final GAA.HRES for TX. In each panel, M denotes the spatial average over all grid points while in panels b, c, and Mo denotes the station mean values while Mc the mean values for the closest grid points to the stations locations

Fig. 4
figure 4

Average annual TN a of the WRF simulation (1995), b of the WRF simulation after step 2 of the methodology for temperatures, c of the daily gridded dataset based only on observations and d of the final GAA.HRES for TN. In each panel, M denotes the spatial average over all grid points while in panels b, c, and Mo denotes the station mean values while Mc the mean values for the closest grid points to the stations locations

Fig. 5
figure 5

Average annual TG a of the WRF simulation (1995), b of the WRF simulation after step 2 of the methodology for temperatures, c of the daily gridded dataset based only on observations and d of the final GAA.HRES for TG. In each panel, M denotes the spatial average over all grid points while in panels b, c and Mo denotes the station mean values while Mc the mean values for the closest grid points to the stations locations

Fig. 6
figure 6

Total annual RR average over the period 1981–2000 of a step 1 of the spatiotemporal model for precipitation and b the final GAA.HRTES for precipitation. In each panel, M denotes the spatial average over all grid points while in panels b, c and Mo denotes the station mean values while Mc the mean values for the closest grid points to the stations locations. KED indicates Kriging with External Drift

From these Figures the following can be extracted: (a) The use of the WRF simulation in the construction of the daily gridded temperature datasets enhance the spatial representation of these fields in the area under study, while the use of only, observational data that were available for this study, leads to a smoother temperature representation. This is more evident for TX (Fig. 3c) where the spatial correlation coefficients between the mean monthly maximum temperatures and the stations elevation is not found statistically significant or exhibit low positive values for the April-September period (the warming is uniform) leading to relatively low range of temperature values within the domain under study. (b) If the purpose of the study is not to produce a daily dataset but only to statistically downscale climate change projections, the output of step 2 for temperatures where the perturbed WRF simulation is obtained for the period 1981–2000 and the output of step 1 for precipitation could be sufficient, since methods that only use mean monthly values (total for precipitation) to bias adjust the climate change projections exist in the literature (e.g., the unbiasing method Deque 2007).

Regarding the results of the GAA.HRES temperatures, the average annual TX (Fig. 3d) over the period 1981–2000 and over the whole domain under study reaches about 20 °C ranging from about 12 °C at the highest mountain peaks to about 24 °C in the low elevation areas such as the Athens urban area as well as the industrial area in Elefsina (areas around NOA and ELE stations, respectively—Fig. 1). For the average annual TN (Fig. 4d) the results indicate higher spatial variability than TX (stable conditions occur under lower temperatures with less mixing and less uniform temperatures) with the highest (lowest) values shown for the areas with the highest (lowest) average annual TX, while intermediate values are shown for the rest of the domain. The average annual TN is about 10 °C (range 4–16 °C). Intermediate values are shown for the average annual TG (Fig. 5d), about 15.5 °C over the whole domain under study, with temperatures ranging from about 8 °C at the highest mountain peaks to about 20 °C in the low elevation areas. As far as the total annual RR is concerned (Fig. 6b), averaged over the twenty-year period, the spatial average over the whole domain is about 450 mm ranging from about 300 mm/year to about 1000 mm/year in the lower and higher altitudes, respectively.

In Fig. 7 the comparison between the observational data and the data from the closest grid point to the stations location from the GAA.HRES is shown for the average annual TX, average annual TN, average annual TG and the total annual RR as well as for selected indices. From the figure it is evident that GAA.HRES exhibits a very good performance when compared to the observations. In particular for the average annual TX the highest Mean Absolute Error (MAE) does not exceed 0.2 °C/year, with GAA.HREs exhibiting similar to the observed trends at all locations. Regarding the number of days with daily TX > 25 °C (SU) the maximum absolute difference 6 days/year is found only in one station (PIR) while for the rest of the stations the differences are about 3 days/year and less. Finally for the number of days with daily TX > 35 °C (SU35), similar results are shown for all station locations. Similar results are shown for average annual TN (and the number of days with daily TN > 20 °C (TR)) and TG. Regarding the annual total RR the highest MAE is about 35 mm/year in only one station location (NFIL) while for the rest of the locations is less than about 15 mm/year. In addition GAA.HRES captures the trends in each location as well as the average number of wet days (RR1) and the monthly maximum 1-day precipitation (RX1day). On the daily timescale the comparison between the observations and the data from the closest grid point to the stations location from the GAA.HRES indicates a perfectly linear positive correlation for all temperature variables, while for precipitation the lowest correlation coefficient is not lower than 0.98 (Figs. S3–S6 Supplementary Material). In addition the observed seasonal trends are maintained by GAA.HRES for all variables and at all locations.

Fig. 7
figure 7figure 7

Temperature average annual values and total annual precipitation for both the observational data (yellow color) and the data from the closest grid point to the station locations from the GAA.HRES (blue color)-from left to right TX, TG, TN and RR- and for all station location over the period 1981–2000. In each panel M indicates the average annual value of TX, TN and TG (in °C/year) and total annual precipitation (in mm/year), S is the annual trend as derived by the Sen's method (°C/year and mm/year for the temperature variables and precipitation respectively), MAE the mean absolute error between the the two datasets (in units similar to the y-axis of each panel). In addition the results for selected indices for each variable are shown. SU (number of days with daily TX > 25 °C) and SU35 (number of days with daily TX > 35 °C) for TX in days/year, TR (number of days with daily TN > 20 °C) for TN in days/year, RR1 (number of wet days) and RX1day (monthly maximum 1 day precipitation) for RR in days/year and mm/day, respectively

Moreover, a leave-one-out cross validation was performed for the daily TX, daily TN and daily TG as well as the daily precipitation sums for the period 1981–2000 (Table S2 Supplementary Material). In leave-one-out cross-validation each observed data point is removed and the observed value is recalculated separately from the remaining data by means of interpolation (Krähenmann et al. 2018). The results indicate that the average MAE over all months is about 0.7 °C, 1.2 °C and 0.8 °C for TX, TN and TG respectively while for RR is about 3 mm. Highest (lowest) MAEs for temperature variables are shown for the warmer (colder) months of the year while the opposite behavior is shown for RR.

4.2 GAA.HRES average seasonal patterns over the period 1981–2000

In Figs. 8, 9 and 10 (left columns) the average seasonal patterns for daily TX, daily TN as well as for total RR are shown, respectively. In addition, the confidence ranges of the above-mentioned variables calculated by the 95th percentiles confidence intervals as derived by bootstrap are provided (right columns). Moreover, in each figure apart from the average mean over all grid points in the area under study, the observed values as well as the values for the closest GAA.HRES grid point to the stations locations averaged over all locations is shown for each variable.

Fig. 8
figure 8

GAA.HRES average seasonal TX (winter (DJF), spring (MAM), summer (JJA) and autumn (SON)) for the period 1981–2000. In the right column the corresponding confidence ranges (CR) are shown. M denotes the spatial average over the grid points covering the area whereas Mo denotes the station mean values while Mc the mean values for the closest grid points to the stations locations. Units are the same as in the colourbar

Fig. 9
figure 9

GAA.HRES average seasonal TΝ (winter (DJF), spring (MAM), summer (JJA) and autumn (SON)) for the period 1981–2000. In the right column the corresponding confidence ranges (CR) are shown. M denotes the spatial average over the grid points covering the area whereas Mo denotes the station mean values while Mc the mean values for the closest grid points to the stations locations. Units are the same as in the colourbar

Fig. 10
figure 10

GAA.HRES total seasonal RR (winter (DJF), spring (MAM), summer (JJA) and autumn (SON)) for the period 1981–2000. In the right column the corresponding confidence ranges (CR) are shown. M denotes the spatial average over the grid points covering the area whereas Mo denotes the station mean values while Mc the mean values for the closest grid points to the stations locations. Units are the same as in the colourbar

In winter (DJF, Fig. 8), average TX is about 11.4 °C with the values ranging from about 3 °C at the highest mountain peaks to 15 °C at the low elevation areas such as the Athens urban area as well as the industrial area in Elefsina (areas around NOA and ELE stations, respectively—Fig. 1). CR indicates an average value over the whole domain of about 0.6 °C with values increasing with increased altitude. On the other hand, summer temperatures (JJA) depict mean values, over the whole domain, of about 30 °C with temperatures varying from approximately 22 °C up to 34 °C. Temperatures lower than about 27 °C are confined in areas with altitude higher than 500 m a.s.l. whereas the maximum overall temperatures are shown for the urban and industrialized areas (areas close to NOA and ELE stations), respectively where the UHI effect is developed. Another feature is the cooler temperatures (lower than about 30 °C) in the eastern coastal areas of the domain due to the impact of etesian winds and the stable sea surface temperatures prevailing at coastal sites (Koletsis et al. 2009; Founda et al., 2019). Average CR is about 0.4 °C with the highest values shown in the high altitude areas similar to the winter season. For the intermediate seasons, average TX vary from 10 to 21 °C in spring (MAM) and from 13 to 25 °C in autumn (SON), the latter thus being warmer on average.

For TN (Fig. 9) the DJF and JJA averages are about 4 °C (range from − 5 to − 9 °C) and 18 °C (range 9–25 °C), respectively. For the intermediate seasons, averaged temperatures are about 8 °C (range 0–14 °C) in MAM and about 12 °C (range 2–18 °C) in SON. CRs averaged over the whole domain do not exceed 0.6 °C for all seasons with the highest values mostly found high altitude areas.

Regarding precipitation sums (Fig. 10), the highest amounts of precipitation occur during DJF (about 190 mm/year) followed by MAM and SON (about 120 mm/year for both seasons) whereas for the dry season (JJA) precipitation does not exceed 27 mm/year on average. In all cases the highest uncertainty is found in the high altitude areas.

Concerning the spatial variability shown by the CRs, the highest CR values, for both TX and TN, are found in high altitude areas while the lowest values are found in metropolitan areas because of the stronger mixing that results in more stagnant conditions there. For similar reasons, the lower temporal variability, within the seasons, is found for TX in JJA where the warming becomes more uniform in the area. In addition in JJA the impact of the sea breeze and the etesian winds contribute to slightly higher variability in the southern and eastern coastal areas of the domain, compared to the inland areas. In contrast, higher values occur in MAM due to the more unstable meteorological conditions than in SON. For TN the highest CR values occur in JJA with the highest values shown in the higher altitude areas followed by the areas in the eastern coast of the domain due to the more effective cooling, under permitting conditions (windy days) in these areas during the night compared to the urban areas of the domain. Similar behavior is found in SON. Regarding RR, the orographic effect is clearly shown in the spatial distribution for the CRs in all seasons. The lowest averaged temporal variation is found in JJA when drying conditions prevail in the domain, while the highest ones are shown in DJF and MAM when most of the rain falls in the area.

It should be mentioned that the results presented in this study regarding the annual and seasonal temperature means as well as the corresponding precipitation sums are comparable to the results of previous gridded data sets produced for Greece (Mamara et al. 2017; Gofa et al. 2019) with the annual and mean monthly values for temperature and precipitation produced in the online tool developed by the Hellenic National Meteorological Service (http://climatlas.hnms.gr/).

4.3 Temperature and precipitation indicators

Besides the analysis of mean temperatures and precipitation sums on the annual and seasonal timescales, the examination of temperature and precipitation indicators is of high relevance, owing to the resulting impacts on many environmental and socioeconomic sectors (i.e. heat stress, floods, agriculture and tourism). The number of summer days (SU, Fig. 11a) is about 130 days/year (averaged over the whole domain) with the higher values (> than about 140–150 days/year) shown for the central, northern and eastern parts of Attica while lower values (< than about 80 days/year) are shown for the higher altitude areas. For the number of very hot days n (SU35, Fig. 11b) the values range from about 0–30 days/year (average about 8 days/year) with regional maximums (> than about 20 days/year) found in the industrialized area of Elefsina, the urban area of Athens as well as the northern coastal areas of the domain. The latter is attributed to the dry air descending from the mountain of Parnitha under southern direction winds. Similar regional maximums are also shown for TR with the highest values reaching about 120 days/year in the industrialized area of Elefsina (Fig. 11c).

Fig. 11
figure 11

GAA.HRES average annual a number of days TX > 25 °C (SU), b number of days TX > 35 °C (SU35), c number of days TN > 20 °C (TR), d maximum 1 day percipiation (RX1day) and e number of days RR > 1 mm (RR1) for the period 1981–2000. In each panel M denotes the spatial average over the grid points covering the area whereas Mo denotes the station mean values while Mc the mean values for the closest grid points to the stations locations. Units are the same as in the colourbar

Regarding precipitation indicators, monthly maximum 1 day precipitation amount (RX1day) as well as the number of wet days (RR > 1 mm) are shown in Fig. 11d, e, respectively. Regarding RX1day, the values averaged over the whole domain are about 18 mm/day with the highest values shown for the mountainous areas in the domain (maximum values reach 40 mm/day). In the lower altitude areas the values are in the range 10–16 mm/day In contrast to RX1day, the highest number of wet day/year are shown for the northern coastal part of the domain reaching 40–50 day/year while for the rest of the areas the values do not exceed about 35–40 days/year. This feature is expected due to the region’s coastal location and the influence of the weather fronts originating from the north east, passing over the Aegean Sea.

4.4 Raw and statistically downscaled climate projections

4.4.1 Temperature results

In Fig. 12 the results of climate projections are shown for both the raw model output (historical, future) as well the statically downscaled RCM model data under the EQM and QDM for the average annual TX (TN in the Supplementary Material, Fig. S7). Comparing the statistically downscaled climate projection for the historical period (top row, middle and right columns in both figures) with Fig. 3d it is evident that both methods adjust the model values towards the corresponding GAA.HRES ones for both variables. The statistically downscaled climate projections indicate higher spatial variability than the regridded raw model output following the spatial variability of the GAA observed gridded dataset. QDM maintains the raw climate change signal between the future and the reference period (increase of about 4.6 °C) while EQM indicates a higher climate change signal (about 4.9 °C). This alteration of the trend by EQM, with respect to the raw RCM signal has been reported in previous studies and it is, a by construction feature, related to the non-linear transformation of trends obtained by quantile mapping (Canon et al. 2015; Tong et al. 2021). Moreover, Maurer and Pierce (2014) demonstrated that underestimation of the observed variance by the raw model in the historical period leads to an amplified climate change signal by EQM while overestimation of the variance yields a dampened climate change signal. From the spatially averaged (over the whole domain) seasonal TX results presented in Fig. 13 (maps are shown in the Supplementary Material, Figs. S8 and S9), it is shown that for observed to raw model output CR ratios higher than 1 (averaged over the whole domain), EQM tends to amplify the climate change signal. This is clearly shown for DJF, MAM and SON (amplification of 0.3–0.6 °C depending on the season), while in JJA (CR ratio < 1) the climate change signal with respect to the raw model output is contracted by 0.2 °C. On the contrary QDM maintains the raw climate change signal in all seasons. In addition, for the absolute threshold based TX indices, SU and SU35 (Fig. 13 bottom row) EQM indicates a lower climate change signal compared to the raw model output, 53 days/year for EQM versus 56 days/year for raw model output and 45 days/year for EQM versus 47 days/year for raw model output, for SU and SU35, respectively. It should be mentioned that for SU the relationship between the CR ratio and the EQM climate change signal compared to the raw one is not consistent to the rest of results for TX, since it is not a robust estimator regarding the prediction of the effect of bias correction on trends as discussed in Maurer and Pierce (2014). Nevertheless, EQM performs better than QDM with respect to the raw climate change signal for both indicators. The latter can also be supported from the spatial distribution of SU35 (Fig. S10 Supplementary Material) where a very high number of days exceeding the 35 °C threshold are found in the industrialized area of Elefsina, the urban area of Athens as well as the northern coastal areas of the domain for the future period and under RCP8.5.

Fig. 12
figure 12

From left to right the results for the raw model output (left column) and the statistically downscaled simulation from EQM (middle column) and QDM (right column). Top row the average annual TX over the period 1981–2000 is shown, middle row the average annual TX over the period 2081–2100 and under RCP8.5 is shown while in the bottom row the absolute differences between the future and control simulations are shown. In each panel M denotes the spatial average over the grid points covering the area. Units are the same as in the colourbar

Fig. 13
figure 13

Top two rows: results for the spatially averaged seasonal (winter, DJF, spring MAM, summer JJA and autumn, SON) TX for all simulations. Bottom row: the average annual number of days TX > 25 °C (SU, left panel) and the average annual number of days with TX > 35°c (SU35, rights) are shown. In each panel RAW.c, EQM.c and QDM.c indicate the raw, EQM and QDM results for the 1981–2000 period while RAW.f, EQM.f and QDM.f ndicate the raw, EQM and QDM results for the 2081–2100 period under RCP8.5. M indicates the spatiotemporal average of each simulation while Δ indicates the absolute differences between the future and historical period of each simulation. CR ratio is the spatially averaged observed to raw model output confidence range (CR) ratios over the period 1981–2000

Regarding TN, the seasonal results are consistent with the ones found for the TX ones. From Fig. 14 it is evident that EQM amplifies the climate change signal compared to raw model output (CR ratios higher than 1) while QDM maintains the climate change signal in levels identical to raw model output for all seasons (maps are shown in the Supplementary Material, Figs. S11–S12). For TR (Fig. 14, bottom left panel and in the Supplementary Material, Figs. S13) both methods deliver changes close to the raw model output with the highest absolute deviation between the raw model output and the bias adjusted ones being less than 5 days/year.

Fig. 14
figure 14

Top two rows: results for the spatially averaged seasonal (winter, DJF, spring MAM, summer JJA and autumn, SON) TN for all simulations. Bottom row: the average annual number of days TN > 20 °C (TR, left panel) is shown. In each panel RAW.c, EQM.c and QDM.c indicate the raw, EQM and QDM results for the 1981–2000 period while RAW.f, EQM.f and QDM.f ndicate the raw, EQM and QDM results for the 2081–2100 period under RCP8.5. M indicates the spatiotemporal average of each simulation while Δ indicates the absolute differences between the future and historical period of each simulation. CR ratio is the spatially averaged observed to raw model output confidence range (CR) ratios over the period 1981–2000

4.4.2 Precipitation results

As far as the annual total precipitation is concerned (Fig. 15), the raw RCM output overestimates (on average over the domain under study) the GAA.HRES in the historical period (Fig. 6d). The projected relative decreases in RR are found slightly lower after bias adjustment: an average 31% reduction for the raw RCM data versus 26% and 22% after bias adjustment using EQM and QDM, respectively. For EQM the findings are consistent to what is found for TX, that is for observed to raw model output CR ratios higher than 1 (CR ratio 1.1) the climate change signal with respect to the raw model output is amplified. Nevertheless the net effect is relative small (about 5% on an annual basis) consistent to what is reported in Maurer and Pierce (2014). Moreover, QDM does not maintain the magnitude of the climate change signal as shown for the average annual TX and TN, which is a by construction feature since accounting for changes in all the quantiles does not guarantee the preservation in the mean value (Cannon et al. 2015).

Fig. 15
figure 15

From left to right the results for the raw model output (left column) and the statistically downscaled simulation from EQM (middle column) and QDM (right column). Top row the total annual RR over the period 1981–2000 is shown, middle row the total annual RR over the period 2081–2100 and under RCP8.5 is shown while in the bottom row the relative differences between the future and control simulations are shown. In each panel M denotes the spatial average over the grid points covering the area. Units are the same as in the colourbar

On the seasonal scale (Fig. 14 and Figs S14 and S15 in the supplementary material) both methods indicate precipitation relative changes close to the raw model output for DJF and MAM (observed to raw model output CR ratios > 1) where most of the rain falls in the area with EQM performing better than QDM. Nevertheless in both seasons the highest deviations for QDM do not exceed 7%. In JJA the highest relative decreases are shown: an average 71% decrease for the raw RCM data versus 65% and 52% after bias adjustment using EQM and QDM, respectively. In SON the average relative changes between the future and the reference period are not found statistically significant for the raw model output as well as for both the bias adjusted simulations. Regarding the precipitation indicators (Fig. 16 bottom row and Figs S16 in the Supplementary Material) for RR1 both methods underestimate the decreases in the number of days/year with QDM performing better than EQM. Finally, for the RX1day indicator both methods show absolute changes of similar magnitude to the raw model output (the relative decrease of − 13% is not relevant for QDM).

Fig. 16
figure 16

Top two rows: results for the spatially averaged seasonal (winter, DJF, spring MAM,summer JJA and autumn, SON) RR for all simulations. Bottom row: the average annual number of days wet days (RR1, left panel) and the average annual monthly maximum 1 day precipitation (RX1day, right panel) are shown. In each panel RAW.c, EQM.c and QDM.c indicate the raw, EQM and QDM results for the 1981–2000 period while RAW.f, EQM.f and QDM.f ndicate the raw, EQM and QDM results for the 2081–2100 period under RCP8.5. M indicates the spatiotemporal average of each simulation while Δ indicates the relative differences between the future and the historical period of each simulation except indicated otherwise. CR ratio is the spatially averaged observed to raw model output confidence range (CR) ratios over the period 1981–2000

5 Conclusions

In this study a methodological framework to obtain statistically downscaled high resolution climate projections over the Attica region was developed. The framework relies on the construction of a local daily gridded dataset for temperature variables (maximum, minimum and mean daily temperatures) and daily precipitation sums. To this aim, a mosaic of data that includes observations derived from station data and a high resolution WRF simulation for the year 1995 were blended using various gridding techniques to produce a 1 km 1 km high resolution daily gridded dataset for the period 1981–2000 and for the aforementioned variables.

The comparison of the gridded dataset against the observations revealed that the produced dataset maintains the long term statistical properties over the period 1981–2000 for both temperature and precipitation variables. Moreover, a leave-one-out cross-validation method was conducted to assess the methods performance, revealing an overall high skill of the datasets in reproducing the monthly average observed temperatures and precipitation sums. In addition, a qualitative comparison with an existing dataset of similar resolution showed relatively high coherency for the average annual temperatures and the annual total precipitation over the period 1981–2000.

The gridded dataset produced in this study, similarly to other gridded datasets produced in other areas, may be useful for any kind of research studies requiring high resolution climatic data such as agronomy, hydrology, forestry which have been found that are under pressure in the area of study (e.g. Giannakopoulos et al. 2011). In addition, the methodology proposed in this study may be useful in other areas with relatively low density of stations and/or complex orography.

To downscale the RCM daily model variables to the GAA.HRES grid the RCM daily data, for the 1981–2000 and the 2081–2100 under RCP8.5 periods, were remapped on the GAA.HRES grid using bilinear interpolation and consequently the models’ output was bias adjusted using GAA.HRES as the reference dataset. In addition, two different bias adjustment techniques, namely the empirical quantile mapping (EQM) and the quantile delta mapping (QDM), were used to examine the impact of bias adjustment in the downscaled high resolution climate projections. The analysis revealed that the selection of the bias adjustment method is important. In particular, the results showed that different bias adjustment methods can affect the simulated change signal in a different way and the final choice of method may depend on the specific goals of a given application. In particular, in our study QDM maintained the raw climate change signal on both the annual and seasonal timescales for both TX and TN. For the absolute threshold based temperature indices (SU, SU35, TR) EQM performed better than QDM, although a small modification of the raw climate change signal was found by the method. Regarding precipitation, the results indicated that EQM performed better than QDM when compared to the raw climate change signal for the annual and the seasonal timescales, with QDM performing better for RR1. Nevertheless, both methods performed quite well for the non fixed threshold precipitation indicator (RX1day). This is of great importance since statistically downscaled (through bias adjustment) datasets are produced widely for various sectors which in many occasions are then used by decision-makers in the public and private sectors to plan climate change related adaptation strategies.