Benoît Gschwind

Lucien Wald

    Testing several methods for the fusion of the two HelioClim-1 and-3 data sets of daily downwelling solar radiation at surface-Application to the creation of a new consistent 1985-2016 data set

Introduction

The downwelling solar radiation incident on the surface of the Earth is known to be an essential variable in many domains. The solar irradiation is the density of energy received from the Sun on a horizontal surface at ground level and per unit surface, while the irradiance is the density of power. Other terms may be found in literature, such as solar exposure or solar flux, surface solar irradiation, downwelling shortwave irradiation, downwelling lar irradiation at the surface, or global surface irradiation. "Surface Solar Radiation Downward" (SSRD) is often used in reanalyses and is adopted in the following; it may refer either to irradiation or irradiance depending on the context.

Accurate assessments of the SSRD can be made from satellite images, see e.g. [START_REF] Amillo | A new database of global and direct solar radiation using the Eastern Meteosat Satellite, models and validation[END_REF], [START_REF] Chen | An efficient method of estimating downward solar radiation based on the MODIS observations for the use of land surface modeling[END_REF], [START_REF] Karlsson | CLARA-A2: the second edition of the CM SAF cloud and radiation data record from 34 years of global AVHRR data[END_REF], [START_REF] Müller | Digging the METEOSAT treasure-3 decades of solar surface radiation[END_REF], or [START_REF] Pfeifroth | Surface Radiation Data Set -Heliosat (SARAH) -Edition 2[END_REF]. The series of geostationary Meteosat satellites offer synoptic views of Europe, Africa and the Atlantic Ocean since 1983 with a spatial resolution of approximately 5 km at nadir and a temporal frequency of 30 min for the first generation and 3 km and 15 min for the second one. Several data sets of SSRD have been constructed from these images, such as the Surface Radiation Data Set -Heliosat (SARAH-2) [START_REF] Pfeifroth | Surface Radiation Data Set -Heliosat (SARAH) -Edition 2[END_REF] or the Copernicus Atmosphere Monitoring Service (CAMS) Radiation Service [START_REF] Qu | Fast radiative transfer parameterisation for assessing the surface solar irradiance: The Heliosat-4 method[END_REF].

Of particular interest here are the HelioClim-1 and -3 data sets of SSRD constructed from Meteosat imagery. The HelioClim Project is an ambitious initiative of MINES ParisTech launched in 1997 to increase the knowledge about SSRD and to offer data over a large geographical area and over a long period of time, to a wide audience [START_REF] Qu | Fast radiative transfer parameterisation for assessing the surface solar irradiance: The Heliosat-4 method[END_REF]. The project comprises the creation and management of several SSRD data sets, called HelioClim data sets and abbreviated hereafter as HC. Though they cover the same geographical area, the HC1 and HC3 data sets have different spatial and temporal resolutions, different accuracies, and different temporal coverages with an overlap of 23 months from February 2004 to December 2005. HC1 offers daily irradiation with a spatial resolution of approximately 30 km for the period 1985-2005. Given SEVIRI having 3 km of resolution at nadir, HC3 contains 15 min SSRD with a spatial resolution of approximately 5 km for the period from February 2004 onwards.

Benefiting from the fact that HC1 and HC3 share an overlap period, we have undertaken in this technical note to demonstrate the better accuracy of HC3 over HC1 and then improve the accuracy of HC1 using HC3 data set and a data fusion method. The main objective of this note is to study the merits and drawbacks of eight methods to perform this fusion of HC1 and HC3. Data fusion has been defined by the European Association of Remote Sensing Laboratories (EARSeL) as a formal framework in which are expressed means and tools for the alliance of data originating from different sources. It aims at obtaining information of greater quality; the exact definition of 'greater quality' will depend upon the application [START_REF] Wald | Some terms of reference in data fusion[END_REF].

Given a specific location, the fusion process adopted here locally adjusts the HC1 time series of daily SSRD onto coincident HC3 daily SSRD for the common period 2004-2005 to benefit from the better accuracy of HC3. The adjustment is performed at each pixel and is applied to the whole data set from 1985 to 2005. In other words, it creates an enhanced HC1 data set, abbreviated as HC1e, which, hopefully, is more accurate than the original HC1. The greater quality may be proven by a comparison with qualified ground measurements of the daily SSRD. The possible improvement in bias and other statistical indicators are assessed by comparing the performance of the original HC1 against these ground measurements on the one hand and those of the HC1e against these same measurements on the other.

The large use of the HC1 data by many researchers in various domains as reported by [START_REF] Lefèvre | The HelioClim-1 database of daily solar radiation at Earth surface: an example of the benefits of GEOSS Data-CORE[END_REF] supports this work and justifies providing some guidance to users of the HC1 data set so that they can implement themselves a fusion method.

As an application, a new data set, called SSRD-H3E, has been created by the concatenation of HC1e and HelioClim-3 using the quantile mapping of KT technique. This data set of consistent 1985-2016 time series of daily SSRD over Europe and Africa is publicly available and is called SSRD-H3E. It has been built with a spatial resolution of 0.25°. Its presentation and validation are another goal of this note.

Section 2 details the HC1 and HC3 data sets and presents the nine stations in Europe and Africa that are used for validation only. Section 3 presents the results of the comparison between the ground measurements and HC1 on the one hand, and HC3 on the other, for the 23 common months in 2004-2005. Such a comparison has not been published up to now to the best of the knowledge of the authors. It shows the greater accuracy achieved by HC3 thus highlighting the potential of a fusion.

The eight fusion methods are described in Section 4. In Section 5, they are applied to HC1 and validated at the 9 stations for the period 1985 to 2003 prior to the calibration period 2004-2005. Conclusions are drawn on the merits and drawbacks of each method.

Section 6 describes the SSRD-H3E data set and its validation for the same stations and the whole period 1985-2016 as well as for each year in this period, thus assessing the reliability of the performance with time. Section 7 concludes the note.

An annex provides a few lines of the Python code for implementing an adjustment based on quantile mapping.

Data sets

The HC1 and HC3 data sets

The HC1 data set has been constructed with images acquired between 1985 and 2005 by broadband imagers aboard the various Meteosat First Generation satellites. At the start of the HelioClim project in 1997, the large amount of data to be processed, the storage capacity, and the cost of the Meteosat images from the Eumetsat agency were obstacles to the production of a data set (Blanc et al., 2007). A data set of Meteosat images in reduced spatial and temporal resolutions has been purchased to overcome these obstacles. This data set was created for the International Satellite Cloud Climatology Project (ISCCP) and is called ISCCP-B2 [START_REF] Schiffer | ISCCP global radiance data set: A new resource for climate research[END_REF]. These images can be considered as having a spatial resolution of 30 km at nadir in a first approximation (Blanc et al., 2007). A great deal of effort was made to ensure the radiometric consistency of all these ISCCP-B2 images, so that the quality of the estimated SSRD is the same throughout the whole period at a given pixel (see details in Blanc et al., 2007). Because of the reduced temporal resolution, HC1 offers only daily SSRD.

The time sampling of the ISCCP-B2 images every 3 h induces limitations in the description of the daily mean of SSRD in HC1. A discrepancy between HC1 and HC3 is created by this effect, whose magnitude depends on the temporal variability of the SSRD as discussed by e.g. [START_REF] England | A study of the errors due to temporal sampling of the earth's radiation budget[END_REF]. Blanc et al. (2007) have studied the effect of the low number of ISCCP-B2 images usable per day on the difference between HC1 and HC3 for a given pixel; this number decreases as the latitude increases, and so does the quality of the HC1 retrieval as a whole. When investigating ocean pixels only in order to avoid effects of orographic origin, these authors found that this difference has a tendency to increase with the absolute value of the latitude, i.e. with the decreasing number of available ISCCP-B2 images, and that the difference decreases as the clearness index increases: i.e. for a given range of latitude, the error is smaller in case of clear skies. However, Blanc et al. (2007) underlined that outside the ocean pixels, the effect of the limited number of satellite observations within a day cannot be easily predicted. This shows in Section 3 where HC1 and HC3 are compared against ground measurements.

The HC3 data set is being constructed from images acquired by the multispectral imagers of SEVIRI aboard the operational series of Meteosat Second Generation satellites since their operational inception in February 2004. Contrary to the Meteosat First Generation images used for HC1, Eumetsat is embedding elements for the radiometric calibration in each SEVIRI image, ensuring the radiometric consistency of images. Once the images converted in radiances, and to improve the consistency between HC1 and HC3, the radiances of the two narrow visible bands of the SEVIRI at 600 nm and at 800 nm are combined to produce broadband radiances that are almost identical to those observed in the broadband channel of the imager aboard the Meteosat First Generation satellites [START_REF] Cros | Simulating Meteosat-7 broadband radiances at high temporal resolution using two visible channels of Meteosat-8[END_REF]. As those synthetic broadband radiances are similar to those used for HC1, the same method is applied to retrieve the SSRD, namely the Heliosat-2 method combined with the ESRA clear-sky model (see details in [START_REF] Lefèvre | Using reduced data sets ISCCP-B2 from the Meteosat satellites to assess surface solar irradiance[END_REF]. A clear-sky model is a model predicting the SSRD in cloud-free conditions. The volume of data to be processed, the storage capacity, and the cost were not any longer obstacles in 2004, thus HC3 is built with the full temporal and spatial resolutions of the images, respectively 4 images per hour and 3 km at nadir. The HC1 and HC3 data sets are live collection of data and are corrected a posteriori as flaws or drawbacks are discovered, yielding several versions. HC1v4 and HC3v5 are the most recent versions, the latter being obtained by the use of the Copernicus McClear clear-sky model [START_REF] Lefèvre | a new model estimating downwelling solar radiation at ground level in clear-sky condition[END_REF] following the method of [START_REF] Qu | Improving HelioClim-3 estimates of surface solar irradiance using the McClear clear-sky model and recent advances in atmosphere composition[END_REF]. The HC1 and HC3 data sets have been the subject of several comparisons against measurements of SSRD made at land-based sites or oceanic moorings, see references in [START_REF] Lefèvre | The HelioClim-1 database of daily solar radiation at Earth surface: an example of the benefits of GEOSS Data-CORE[END_REF] for HC1 and see e.g. [START_REF] Trolliet | Estimating downwelling solar irradiance at the surface of the tropical Atlantic Ocean: A comparison of PIRATA measurements against several re-analyses and satellite-derived data sets[END_REF] and its references for HC3.

For this work, the HC1, HC3 and the top-of-the-atmosphere time series of daily irradiation, respectively G hc1 , G hc3 and G0, have been downloaded from the SoDa Service [START_REF] Gschwind | Converting a successful research project into a sustainable service: the case of the SoDa Web service[END_REF] (http://www.soda-pro.com). The clearness indices KT hc1 and KT hc3 are respectively defined as the ratios G hc1 /G0, and G hc3 /G0. The clearness index is an indicator of the atmospheric transparency. It depends on several factors such as cloud presence and properties, content in water vapor or aerosol properties, and is a means to characterize the gross influence of clouds as the latter are the most influencing factor on the extinction of solar radiation by the atmosphere. It is typically close to 0.8 in cloud-free conditions, and close to 0.1 in overcast conditions with optically thick clouds.

Description of the stations and their measurements

Time series of daily irradiation G wrdc were downloaded from the web site (http://wrdc.mgo.rssi.ru) of the World Radiation Data Center (WRDC), a laboratory of the Voeikov Main Geophysical Observatory in Saint Petersburg, Russia, under the control of the World Meteorological Organization (WMO). G wrdc is given in J cm -2 as is usual in meteorology. If necessary, the daily mean of irradiance E wrdc is obtained by converting G wrdc in J m -2 and dividing the result by the number of seconds in 24 h:

E wrdc = G wrdc * 10000 / (24 x 3600) = G wrdc * 0.11574, in W m -2 .
(

The values of G wrdc were kept only when the quality flags set by WRDC were equal to 0 ("good value").

Otherwise, the corresponding day was declared as having invalid measurements. An additional visual screening was performed by the authors, looking for possible anomalies and rejecting very low values. This screening took into account the WMO guide on measurements of solar radiation (WMO, 2012) that says that the relative uncertainty (95 % probability) for "good quality" of measurements is 5 % if G wrdc > 800 J cm -2 ; otherwise the uncertainty is set to 40 J cm -2 .

Nine stations (Table 1) were used in this work spanning the period 1985-2016. Each station offers enough data for most of the years, i.e. i) more than 85 % of the days within a year bear valid daily SSRD, ii) no more than 5 days are missing within a month and iii) possible gaps cover less than 5 consecutive days. The constraints were somehow relaxed at the station Toravere, and to a lesser extent at Valentia, where the sun is very low above horizon in December and January yielding possible very low G wrdc that may fall below the uncertainty. The constraint on the number of missing days was relaxed at the station Toravere (up to 15 days), and to a lesser extent at Valentia (up to 10 days), where the sun is very low above horizon in December and January yielding possible very low G wrdc that may fall below the uncertainty. Using several periods is a means to assess the reliability of the performances with time. The filtered data set is available at (ftp://ftp.oielab.net/pub/hc1_hc3_article/ground_measurements_gschwind_wald.zip). 20) are those offering the least (Fig. 1). As a whole, El Arish exhibits the greatest means of G wrdc (Fig. 2) because it experiences many cloud-free days (KT wrdc >0.60 in Fig. 3). Then comes Maputo which is close to the Equator and thus it experiences large extraterrestrial daily irradiations G0. Carpentras and Nice are located in the Provence in France and experience a large number of cloud-free days with means of KT wrdc close to 0.55 (Fig. 3). In the comparison process against ground-based measurements, the two time series have different spatial supports of information: measurements are for a single point in space, whereas the estimated irradiances are for the area of a pixel (3 km for HC3v5 and 30 km for HC1 at nadir) or a grid cell (0.25° for SSRD-H3E). This approach is that used in many comparisons of satellite-derived quantities or meterological reanalyses. It means implicitly that the SSRD at coarse scale and at point-like scale can be compared and further that the SSRD is homogeneous at the coarse scale, i.e. does not offer any noticeable variation in space. This may reveal wrong depending on the site as cloud properties, including presence and optical depth, may vary within the pixel or grid cell. Such issues have been discussed long ago with respect to solar radiation, see e.g. [START_REF] Hakuba | Spatial representativeness of ground-based solar radiation measurements[END_REF], [START_REF] Pinker | Effects of spatial sampling of satellite data on derived surface solar irradiance[END_REF], or [START_REF] Zelenka | Effective accuracy of satellite-derived hourly irradiances[END_REF]. The magnitude of the effect due to the sole change in spatial support would depend on the actual underlying spatial variability of the SSRD. The lower this variability, the lower the magnitude of the effect. Most of the selected sites are in areas fairly homogeneous in terms of climate, orography and land cover and we may expect the effect of the differences in spatial support of information to be limited, because we are dealing with daily means. Using monthly averages is a means to further reduce the errors (Zhao et al., 2015). The Nice station is on the coast and may be affected though it does not clearly show on the performances when comparing with SSRD-H3E in Section 6. [START_REF] Zelenka | Effective accuracy of satellite-derived hourly irradiances[END_REF] studied the differences between the readings of ground-based pyranometers at different timescales and different distances in homogeneous areas. They found a standard deviation of 10-15% within 10 km relative to the hourly mean of SSI. Although it is difficult to predict because spatial variability is not a random variable and depends on each location and time because of the spatial and temporal properties of the cloudiness, we can expect an error of several per cent relative to the daily mean of SSRD. It should be noted that this error due to the mismatch in spatial support may not necessarily add to the error due to the sole modelling and that both errors can counterbalance, yielding more accurate results than expected.

Comparison of the accuracies of HC1 and HC3 in the period 2004-2005

Prior to the presentation of the fusion methods, this section presents the results of the comparison between the ground measurements and HC1 on the one hand, and the ground measurements and HC3 on the other, for the 23 common months in 2004-2005. To the best of the knowledge of the authors, such a comparison has not been published up to now and it is believed it was necessary to show the potential of a fusion of HC1 with the more accurate HC3.

Table 2 reports the bias, standard deviation of errors, correlation coefficient and slope of the fitting line for both HC1 and HC3 at each station. The bias is negative for HC1 and ranges from -314 (-24 % in relative value) to -22 J cm -2 (2 %) with the exception of El Arish where the bias is positive (101 J cm -2 , 5 %). The variation of the bias for HC3 is smaller as it ranges from -28 (-2 % in relative value) to 50 J cm -2 (3 %). The standard deviation of the errors ranges from 236 (13 %) to 335 J cm -2 (28 %) for HC1. It is much less for HC3 than for HC1 at all stations as it ranges from 109 (7 %) to 166 J cm -2 (13 %). The correlation coefficient is very large for HC1 but it is even greater for HC3. The slope of the fitting line for HC3 is similar to, or closer to 1 than that for HC1.

These figures demonstrate that HC3 is more accurate than HC1. Hence, this establishes the potential of an improvement of the accuracy of HC1 by a data fusion method.

Description of the eight fusion methods

The information entering a fusion process should present several properties [START_REF] Pau | Sensor data fusion[END_REF]. The data to be fused must be aligned, i.e. they must be represented within a common geographical time system. This is solved by using functions converting the image coordinates into latitude and longitude and using this canonical representation for geographical reference, and by using daily irradiation for both data sets.

Several methods were tested for the fusion. A brief review of other methods can be found in [START_REF] Polo | Preliminary survey on site-adaptation techniques for satellite-derived and reanalysis solar radiation datasets[END_REF].

The principle of the fusion is as follows. Given a method f, its parameters are established by comparing G hc1 and G hc3 , or KT hc1 and KT hc3 , for the common period February 2004 to December 2005 (23 months). Then, f is applied to G hc1 to yield the enhanced data set G hce . If f applies to KT hc1 , then it yields KT hce which after multiplication by G0 yields G hce .

Adjustment of means or medians (P50I and P50K)

One of the simplest methods consists of computing the means, or the medians, of G hc1 and G hc3 and then adding the difference of means or medians to G hc1 to yield G hce , here denoted as G hcP50I and G hcP50K . For the sake of simplicity, only the version with the median is presented as both versions give similar performances:

G hcP50I = G hc1 + (median(G hc3 ) -median(G hc1 )), (2)
G hcP50K = G0 [KT hc1 + (median(KT hc3 ) -median(KT hc1 ))]. (3)
It can be mathematically demonstrated that the standard deviation of errors, the correlation coefficient and the slope of the fitting line are the same for G hc1 and G hcP50I . No improvement can be expected for these indicators.

Ratioing methods (RatioI and RatioK)

The third and fourth tested methods compute a ratio of the means of G hc1 and G hc3 , or KT hc1 and KT hc3 , and then multiply G hc1 , or KT hc1 , by this ratio:

G hcRatioI = G hc1 (mean(G hc3 ) / mean(G hc1 )), (4)
G hcRatioK = G0 [KT hc1 (mean(KT hc3 ) / mean(KT hc1 ))]. (5)
It can be mathematically demonstrated that G hc1 and G hcRatioI have the same correlation coefficient. No improvement can be expected for this indicator.

Affine functions (AffI and AffK)

An affine function f = a G hc1 + b is adjusted onto G hc3 . It is then applied to G hc1 to yield G hcAffI . Same approach for KT to yield G hcAffK . The parameters a and b may be computed by a least square fitting or axis of inertia. For the least-square fitting, the parameters a and b are given by:

a² = cov(G hc3 , G hc1 ) / var(G hc1 ), (6)
b = mean(G hc3 ) -a mean(G hc1 ), (7)
where cov denotes the covariance, and var the variance.

The parameters a and b for the first axis of inertia are given by:

a = [[var(G hc3 )-var(G hc1 )] + [(var(G hc3 )-var(G hc1 ))² + 4 cov(G hc3 , G hc1 )]] / cov(G hc3 , G hc1 ), (8) b
= mean(G hc3 ) -a mean(G hc1 ). (9)
The results are fairly similar between both adjustment techniques. Hereafter the results for the first axis of inertia are presented. It can be mathematically demonstrated that G hc1 and G hcAffI have the same correlation coefficient.

No improvement can be expected for this indicator.

Quantile mapping (QMI and QMK)

The method of 'quantile mapping' consists of adjusting the cumulative distribution function of G hc1 onto that of G hc3 for the calibration period, thus yielding an abacus that is used to convert the G hc1 data set into G hcQMI , or G hcQMK if the quantile mapping is applied to KT hc1 and KT hc3 . A pseudo-code in Python is given in Annex A.

Quality of each fusion method

If the fusion is efficient, the following improvements in the indicators characterizing the quality of the fusion are expected:

 the absolute value of the bias should decrease for HC1e compared to HC1,  the standard deviation should decrease for HC1e compared to HC1,  the correlation coefficient should increase for HC1e compared to HC1,  the slope s of the fitting line should be closer to 1 for HC1e than for HC1. In other words, the absolute value of the difference between the slope and 1: ‖s-1‖, should be closer to 0 for HC1e than for HC1. This is the approach adopted here for a better visual analysis of the possible benefit of the fusion.

In the following, each fusion method is calibrated using the 23 common months of HC1 and HC3, and then applied to HC1 and validated at the 9 stations for the period 1985 to 2003 prior to the calibration period 2004-2005. Results are analyzed with respect to the expectations listed above.

Results of the tested fusion methods

Figure 4 exhibits the absolute value of the bias at each station for the original HC1 and the results each fusion method. A benefit is obtained by a fusion method if the absolute value of the corresponding bias is less than that for HC1, in red. The bias for HC1 is less than 100 J cm -2 at six stations in absolute value. It amounts to 240 J cm -2 at Wien (20 % of the yearly mean of the daily irradiation), 186 J cm -2 at Auxerre (16 %) and 240 J cm -2 at Payerne (14 %). Any of the fusion methods brings an improvement at Wien, Auxerre, Payerne, and Maputo. On the contrary, the bias is degraded at Valentia, Carpentras, Nice and El Arish whatever the method. Surprisingly, the P50I and P50K methods do not improve the bias at each station. The good adjustment obtained for the calibration period is not always appropriate to this former period. On the contrary, the bias is significantly increased by P50I at Carpentras and Nice. RatioK performs as one of the worst methods while RatioI, AffI, AffK, QMI, and QMK provide the best results among the eight methods as a whole. Figure 5 exhibits the standard deviation of errors at each station for the original HC1 and the results from each fusion method. A benefit is obtained by a fusion method if the corresponding standard deviation is less than that for HC1. The standard deviation mostly ranges between 220 and 270 J cm -2 for G hc1 . As written before, P50I has the same standard deviation than HC1: there is no benefit. P50K offers a status quo as a whole. RatioI and RatioK degrade the performance as the standard deviations are greater than those for HC1. Exceptions are found at El Arish and Maputo where a slight improvement and a status quo, respectively, are observed. AffI and QMI behave similarly. Except at El Arish, they exhibit a status quo or a degradation of performance. AffK and QMK behave similarly and tend to exhibit a standard deviation less than that of HC1. Figure 6 exhibits the correlation coefficient at each station for the original HC1 and the results from each fusion method. A benefit is obtained by a fusion method if the corresponding correlation coefficient is greater than that for HC1. The correlation coefficient for G hc1 is large and greater than 0.92 in all cases. As written before, P50I, RatioI, and AffI have the same correlation coefficient than HC1: there is no benefit. Actually, RatioK exhibits also the same correlation coefficients as HC1. QMI tends to slightly increase the correlation. P50K, AffK and QMK behave similarly and offer correlation coefficients greater than those of HC1. Exceptions are at Valentia where there is a status quo, and at Toravere where there is a status quo for P50K and a slight degradation for AffK and QMK. Figure 7 exhibits the absolute difference between the slope of the fitting line and 1 at each station for the original HC1 and the results from each fusion method. A benefit is obtained by a fusion method if this indicator is closer to 0 than that for HC1. The slope for G hc1 is close to 1.00±0.11 at all stations. As written before, P50I has the same slope than HC1: there is no benefit. The results for P50K are highly variable and offer either noticeable improvements (Wien, Payerne, Nice, El Arish) or noticeable degradation (Toravere, Valentia, Auxerre, Carpentras). RatioI and RatioK degrade the performance at all stations, except at El Arish and Maputo where an improvement is observed. AffI, AffK, QMI, and QMK behave fairly: the slope is improved at all sites, except at Auxerre and Maputo where degradation is observed.

Figure 7. Absolute difference between the slope of the fitting line and 1, at each station for the original HC1 and the results from each fusion method in the period 1985-2003. The closer the difference to 0, the better.

As a summary, the benefit of adjustment of the irradiance P50I is limited to the bias as the other statistical indicators are left unchanged. The adjustment of the clearness index P50K brings more improvement compared to P50I regarding also the other indicators. One may conclude that the fusion by adjustment of medians may bring a benefit in some cases but an improvement is not obtained at each station and each period. RatioI and RatioK offer both improvement and degradation for the bias depending on the station. They tend to degrade the performances of HC1 for the three other indicators or to offer a status quo. AffI and QMI on the one hand, and

AffK and QMK on the other tend to behave similarly. While AffI and QMI show a tendency to degrade the standard deviation or to offer a status quo, AffK and QMK tend to improve it. The four methods tend to improve the two other indicators or offer a status quo.

Discussion on the performance of the tested fusion methods

This assessment of performance shows that the results are not satisfactory for the simplest methods: P50I, P50K, RatioI and RatioK. There is no benefit brought by them, and even more one may observe degradation in performance. The other methods: AffI, AffK, QMI and QMK usually offer improvement to all indicators, though degradation may be observed for the bias. We have seen no clear picture of the relationship between the performances of a method and the station (climate, latitude, mean value, performance of HC1).

As a conclusion, one notes that an efficient fusion is hard to obtain for all periods and at all stations simultaneously. Reducing the bias is the greatest difficulty. Analysis of graphs (not shown) shows that the bias is more and more reduced from 1985 to 2004 as a whole as the year approaches the calibration period but it may still be degraded for the early years. A solution should be devised for that but it would not be simple given the short overlap between HC1 and HC3.

It was found that the methods applying to the clearness index behave better than those applying to the irradiation. QMK and AffK perform similarly for the bias and the standard deviation, QMK exhibits slightly greater correlation coefficients and AffK offers better results for slope. Hence, one may conclude that both methods are equivalent with respect to the studied indicators.

The affine model deals with mean, covariance and variances, while quantile mapping applies to the cumulative distribution functions. The latter could be a disadvantage since the original cumulative distribution function of HC1 at a pixel for the whole period 1985-2005 is forced to be close to that of HC3 for the period 2004-2005. However, the tests shown here demonstrate that there is no major drawback with this approach. Further tests applied to the cumulative distribution functions and quantiles (not shown) show that QMK tends to slightly better represent the cumulative distribution functions than the other methods. This is the reason why the QMK method was used to create the new SSRD-H3E which is described in the following section.

The new SSRD-H3E

The new SSRD-H3E data set is an attempt to offer consistent time series of daily mean of irradiance over the period 1985-2016 and after as updates of HC3 are available. This section describes this data set and its validation for the same stations and periods as above.

Construction of the SSRD-H3E data set and its access

To begin with, the HC3 data sets were aligned on the same temporal step, same geographical area, same spatial resolution and same grid cells of HC1 dataset. HC3 were computed on a latitude-longitude grid with a cell of 0.05° in size and summarized to daily irradiation. Then, a truncated and apodized sinus cardinal filters were applied to remove details less than 0.25° thus avoiding aliasing effects. Finally, we resampled HC3 in both latitude and longitude at a step of 0.25° for the geographical coverage [-60°, 60°] in latitude and [-60°, 60°] in longitude. The resampled HC1 data set was adjusted by the QMK method onto the resampled HC3 data set for the common period of the 23 months. The adjustment was then applied to the period January 1985 to January 2004 included, to yield the first part of the new data set. The second part was made of the resampled HC3 data set, from February 2004 up to December 2016. The SSRD-H3E, therefore, contains daily mean of irradiance in W m -2 contrary to the daily irradiation of HC1 and HC3.

The data are stored under the format HDF in a single file. Daily mean of irradiance at the top-of-the-atmosphere were calculated with the SG2 algorithm [START_REF] Blanc | The SG2 algorithm for a fast and accurate computation of the position of the Sun[END_REF] and are included in the data set. A DOI has been allotted to the new data set: www.doi.org/10.23646/385c3cd1-968e-4924-8c5c-17a098736aa6. Access to the data obeys the principles of data sharing in the Global Earth Observation System of Systems (GEOSS), and like HC1 [START_REF] Lefèvre | The HelioClim-1 database of daily solar radiation at Earth surface: an example of the benefits of GEOSS Data-CORE[END_REF], the SSRD-H3E data set is part of the GEOSS Data Collection of Open Resources for Everyone (GEOSS Data-CORE), which is a distributed pool of documented data sets with full, open and unrestricted access at no more than the cost of reproduction and distribution.

Validation of the SSRD-H3E data set

Two validations were performed by comparing the data G h3e against the same WRDC measurements G wrdc as before (Figs 123). In the first validation, all valid years were merged to yield the period 1985-2016, and in the second one, the comparison was performed for each single year. Table 3 reports the bias, standard deviation of errors, correlation coefficient and slope of the fitting line for the whole period 1985-2016 at each station. Values are given as irradiation or irradiance for the sake of comparison with the previous sections and other published works. Relative values in percent were computed by dividing the indicator by the mean daily irradiation. The bias ranges from -42 to 97 J cm -2 (-5 to 11 W m -2 ). It is low at Toravere, Wien, Auxerre, El Arish, and Maputo. It is noticeable at Valentia, Payerne, Carpentras, and Nice. The standard deviation of errors is fairly constant, around 209-229 J cm -2 (24-27 W m -2 ) with lows at 158, 179 and 191 J cm -2 (18, 21, and 22 W m -2 ). The correlation coefficient is large and greater than 0.94. The slopes are very close to 1, with a low at Maputo. Actually these indicators vary in time. The bias (Fig. 8), and its relative value (Fig. 9), offer variations in time that depend on the station. The relative values were computed by dividing the bias for a given year by the corresponding yearly mean of the daily irradiation. In the first years before 1989, the bias is positive (overestimation) at all stations except El Arish and Maputo where it is negative. There is a clear trend at all stations for the bias to converge towards 0 as time progresses. After 2004, this is clearly explained by the fact that SSRD-H3E is made from HC3 only, which is more accurate than HC1-QMK. Variations in time may also be observed for the standard deviation of errors (Fig. 10), and its relative value (Fig. 11), and depend on the stations. El Arish exhibits small variation from 1985 to 2007 (no data available after), around 140 J m -2 (approx. 16 W m -2 ). On the contrary, the other stations show two regimes: one before 2004 and one after. The period before 2004 is under the influence of HC1 even corrected by the fusion and there is more variability in the standard deviation than for the following period, which originates from HC3. Before 2004, the standard deviation is less than 250 J cm -2 (approx. 29 W m -2 ) in most cases at all stations. After 2004, the standard deviation is less than 140 J cm -2 (approx. 16 W m -2 ) in most cases at all stations and the variability is less for a given station. These two regimes appear clearly in relative values in Fig. 11. The correlation coefficient (Fig. 12) is large in all cases. The two regimes, before 2004 and after, are clearly visible. The correlation coefficients are highly variable before 2004. After 2004 they are greater than before and tend to approach unity. They are greater than 0.98 as a whole. Maputo is an exception with coefficients ranging from 0.94 to 0.97 for this period. A similar behavior can be observed for the slopes of the fitting lines (Fig. 13). Prior to 2004, a large variability may be observed in slopes from year to year and from station to station, though the range of variation of the slopes is limited between 0.9 and 1.1, except for Maputo where the slopes are less than 0.9 for this period. After 2004, the slopes are closer to 1 at all stations and the variability is less: the range of variation is limited between 0.96 and 1.04.

(J cm -2 ) (W m -2 ) % (J cm -2 ) (W m -2 )

Conclusions

This technical note has evaluated several methods aiming at improving the accuracy of the HC1 data set by fusion with the more accurate HC3 data set using an overlapping period of 23 months. The eight fusion methods tested here were fairly simple to implement with limited programming skills in Python, R or Matlab, or any language handling calls to web services, possibly in an automated way if building on the existing web services offering access to time series of the HC data sets.

The results of the several cases studied here shed light on the merits and drawbacks of each method and provided some guidance to users of the HC data sets so that they can implement themselves a fusion method. It was found that the simplest methods: P50I, P50K, RatioI and RatioK, bring limited benefit. Even more, instead of improving HC1, they may degrade its performance. Fusing the clearness index usually brings better results than fusing the irradiation. The methods AffK (affine function) and QMK (quantile mapping) perform the best, though degradation may be observed for the bias at several occasions. There are status quo or clear improvements of the standard deviation, correlation coefficient and slope of the fitting lines. As a whole, the authors note that the benefit brought by the two best methods is not always spectacular at all stations.

More sophisticated fusion methods may result in better performance. A simple example is the possible combination of the QMK and P50K methods in a cascade-type approach to avoid bias degradation. The task for an efficient fusion is complicated by the limited overlap of 23 months only between HC1 and HC3. If the fusion method calls upon statistical quantities, such as means, medians, cumulative distribution functions, which is the case of the eight methods here but AffI and AffK -and not on coincident time series of values like for AffI and AffK,-this limitation can be removed by using the whole data sets HC1 and HC3.

Other more complex mathematical tools may be used. One may also exploit the recent findings by [START_REF] Bengulescu | On the intrinsic time-scales of temporal variability in measurements of the surface solar radiation[END_REF] who reported that any time series of daily irradiation is composed of several characteristic timescales of variability, some having a deterministic origin, such as the yearly cycle, and others being realizations of various stochastic processes. One may devise a fusion method that comprises several sub-methods that are dependent on the timescale and the stochastic nature of the timescale.

The new SSRD-H3E data set has been created by the concatenation of the HC1e data set resulting from the QMK method between 1985 and January 2004 (included) and the HC3 data set for the following days. It offers 1985-2016 time-series of daily irradiation over Europe and Africa with a spatial resolution of 0.25°. The data set has been conceived as a simple and tractable means for the verification of models and for the adjustment of reanalyses and climate-change projections as discussed in e.g. [START_REF] Babst | Verification of NCEP reanalysis shortwave radiation with mesoscale remote sensing data[END_REF] with a spatial resolution close to those of the reanalyses and climate-change projections to ease the comparison and adjustment. Comparison with measurements made at nine stations for 15 periods of 2 years shows that the bias is variable with time for the early years. The standard deviation of errors is the greatest for the early years but is less than 250 J cm -2 (approx. 29 W m -2 for daily mean of irradiance). After 2004 (included), the standard deviation is less than 140 J cm -2 (approx. 16 W m -2 ). The correlation coefficient is always large and always greater than 0.92.

No direct comparison was made against reanalyses or other similar products. However, it can be speculated that the SSRD-H3E data set compares favorably to the reanalyses, given published results of comparisons between such reanalyses and measurements performed at terrestrial or marine stations. For similar stations, [START_REF] Boilley | Comparison between meteorological re-analyses from ERA-Interim and MERRA and measurements of daily solar irradiation at surface[END_REF] have reported standard deviations of errors, and correlation coefficients for the MERRA and ERA-Interim re-analyses and HC1. They found superior performances for HC1. For reanalyses, though there is not enough information to perform its detailed calculations, the standard deviation of errors is typically 350 J cm -2 (approx. 41 W m -2 ) as a whole; the reported correlation coefficients may be much less than 0.8. [START_REF] Trolliet | Estimating downwelling solar irradiance at the surface of the tropical Atlantic Ocean: A comparison of PIRATA measurements against several re-analyses and satellite-derived data sets[END_REF] published a comparison of measurements performed within the PIRATA network in the tropical Atlantic Ocean against several satellite-derived data sets, including HC3, and the ERA5 and MERRA-2 reanalyses. They found superior performances for the satellite-derived data sets. The standard deviation of errors is typically 320 J cm -2 (approx. 47 W m -2 ) for ERA5 and 420 J cm -2 (approx. 49 W m -2 ) for MERRA-2 whereas it is190 J cm -2 , approx. 22 W m -2 for HC3. The reported correlation coefficients may be much less than 0.6 for ERA5 and 0.5 for MERRA-2 and is always greater than 0.87 for HC3. Of course, reanalyses have their merits, including other meteorological variables, worldwide geographical coverage and several decades back in time.

The irradiation of a reanalysis may be corrected by the means of the SSRD-H3E, possibly using one the presented methods, though care must be taken when adjusting only one variable in a set of consistent and not independent meteorological variables [START_REF] Dekens | Multivariate distribution correction of climate model outputs: A generalization of quantile mapping approaches: Multivariate distribution correction of climate model outputs[END_REF].

Other sources of SSRD over the Meteosat coverage or broader are available. Among them are the data sets established by the Satellite Application Facility on Climate Monitoring (CM-SAF) of Eumetsat. The CLARA-A2 data record covers the period 1982 until 2015 and is derived from the AVHRR (Advanced Very High Resolution Radiometer) sensor carried by polar-orbiting, operational meteorological satellites [START_REF] Karlsson | CLARA-A2: the second edition of the CM SAF cloud and radiation data record from 34 years of global AVHRR data[END_REF]. It offers daily means of irradiance with a spatial resolution of 0.25°. As written above, the SARAH-2 data set is derived from Meteosat images [START_REF] Pfeifroth | Surface Radiation Data Set -Heliosat (SARAH) -Edition 2[END_REF]. It offers 30 min values of SSRD with a spatial resolution of 0.05° for the period starting in 1983. A slightly different version is also proposed, called SARAH-JRC [START_REF] Urraca | Extensive validation of CM SAF surface radiation products over Europe[END_REF]. Also available are HC3v5 of course and the Copernicus Atmosphere Monitoring Service (CAMS) Radiation Service [START_REF] Qu | Fast radiative transfer parameterisation for assessing the surface solar irradiance: The Heliosat-4 method[END_REF], abbreviated as CRS. These two latter data sets begin in 2004 and offer 15 min SSRD with the same resolution than the SEVIRI pixels, i.e. 3 km at nadir and degrading as the satellite viewing angle is increasing. Several works have compared the quality of HC3v5 and CRS and have reported a slightly better quality for HC3v5, see [START_REF] Trolliet | Estimating downwelling solar irradiance at the surface of the tropical Atlantic Ocean: A comparison of PIRATA measurements against several re-analyses and satellite-derived data sets[END_REF] and its references. [START_REF] Trolliet | Estimating downwelling solar irradiance at the surface of the tropical Atlantic Ocean: A comparison of PIRATA measurements against several re-analyses and satellite-derived data sets[END_REF] also performed a comparison between SARAH-2, CRS and HC3v5 for daily means of irradiance in the tropical Atlantic Ocean for the period 2012-2013 and reported fairly similar performances between the three data sets. It can thus be speculated that at the spatial resolution of 0.25°, the SSRD-H3E data set offers the same quality than the SARAH-2 data set for the period starting in 2004, and surpasses the CRS data set. One may not conclude for the previous periods but one may note that [START_REF] Pfeifroth | Trends and variability of surface solar radiation in Europe based on surface-and satellite-based data records[END_REF] reported better performances of SARAH-2 in the period 1992-2015 than in the period 1983-1991 in agreement with our findings for SSRD-H3E (graphs 8 to 13). [START_REF] Urraca | Extensive validation of CM SAF surface radiation products over Europe[END_REF] have assessed the quality level of the daily means of the SARAH-JRC and CLARA-A2 data sets against ground measurements in Europe for 2005 to 2015. Though comparison between their results and ours in Table 3 is difficult because of the differences in period, in stations and in spatial resolution for SARAH-JRC, one may note that performances are fairly similar between these products and SSRD-H3E for the bias taking into account that all products exhibit large changes from station to station. Table 3 in [START_REF] Urraca | Extensive validation of CM SAF surface radiation products over Europe[END_REF] reports the average of the biases and the standard deviation. If one assumes a normal law for the statistical distribution of the bias, 95 % of the bias are in the range [-8, 9] W m -2 for SARAH-JRC and [-7, 9] W m -2 for CLARA-A2, which are close to the range [0, 11] W m -2 reported in our own Table 3 for the European stations. A more thorough analysis is necessary to draw solid conclusions.

Figures 1 ,

 1 Figures 1, 2 and 3 provide the number of coincident valid daily irradiations per valid year between the WRDC measurements and HC1 and SSRD-H3E estimates, as well as the corresponding means of G wrdc and KT wrdc for each year at each station. Payerne, and Nice are the stations offering the greatest number of valid years (31); El Arish (13) and Valentia (20) are those offering the least (Fig.1). As a whole, El Arish exhibits the greatest means of G wrdc (Fig.2) because it experiences many cloud-free days (KT wrdc >0.60 in Fig.3). Then comes Maputo which is close to the Equator and thus it experiences large extraterrestrial daily irradiations G0. Carpentras and Nice are located in the Provence in France and experience a large number of cloud-free days with means of KT wrdc close to 0.55 (Fig.3).

Figure 1 .

 1 Figure 1. Number of coincident valid daily irradiations per year between the WRDC measurements and HC1 estimates from 1985 to 2016, or the SSRD-H3E estimates from 2004 to 2016 at each station

Figure 4 .

 4 Figure 4. Absolute value of the bias at each station for the original HC1 and the results of the each fusion method in the period 1985-2003. In irradiance on right axis.

Figure 5 .

 5 Figure 5. Standard deviation of errors at each station for the original HC1 and the results of the each fusion method in the period 1985-2003. In irradiance on right axis.

Figure 6 .

 6 Figure 6. Correlation coefficient at each station for the original HC1 and the results of the each fusion method in the period 1985-2003

Figure 8 .

 8 Figure 8. Bias in daily irradiation of the new data set SSRD-H3E for each year at each station in the period 1985-2016. In irradiance on right axis.

Figure 9 .

 9 Figure 9. Bias in daily irradiation of the new data set SSRD-H3E for each year at each station relative to the yearly mean of the daily irradiation in the period 1985-2016

Figure 10 .

 10 Figure 10. Standard deviation of errors in daily irradiation of the new data set SSRD-H3E for each year at each station in the period 1985-2016. In irradiance on right axis.

Figure 11 .

 11 Figure 11. Standard deviation of errors in daily irradiation of the new data set SSRD-H3E for each year at each station relative to the yearly mean of the daily irradiation in the period 1985-2016

Figure 12 .

 12 Figure 12. Correlation coefficient between WRDC measurements and SSRD-H3E estimates for each year at each station in the period 1985-2016

Figure 13 .

 13 Figure 13. Slope of the fitting line between WRDC measurements and SSRD-H3E estimates for each period at each station in the period 1985-2016

Figure 14 .Figure 15 .

 1415 Figure 14. Monthly means and standard deviations from the WRDC data and SSRD-H3E estimates at Wien Hohe Warte in the period 1985-2016

#

  Adjustment of clearness index KT funcK = matching_histogram_t(KThc3_2004-2005, KThc1_2004-2005, 1.0) ktQmK = funcK(KThc1_1985-2005) class matching_histogram_t: def __init__(self, x, y, maxValues): # adjust the values of a series y so that its histogram match that of the series x # x and y have the same length # return the adjusted y x = np.array(x, copy=True) y = np.array(y, copy=True) mask = np.isfinite(x)&np.isfinite(y) if np.all(~mask): yr1 = np.array([0.0, 1.0], dtype='f8') ya1 = np.array([np.NAN, np.NAN], dtype='f8') self.yr = np.linspace(0.0, maxValues, 100) self.ya = np.interp(self.yr, yr1, ya1) return x = x[mask] y = y[mask] # get the set of unique values and their corresponding indices and counts for x and y x_values, x_val_index = np.unique(x, return_inverse=True) x_freq = np.bincount(x_val_index) y_values, y_val_index = np.unique(y, return_inverse=True) y_freq = np.bincount(y_val_index) # take the cumsum of the frequencies, normalize by the number of values and get the empirical cumulative distribution functions x_quantiles = np.cumsum(x_freq).astype(np.float64) x_quantiles /= x_quantiles[-1] y_quantiles = np.cumsum(y_freq).astype(np.float64) y_quantiles /= y_quantiles[-1] # interpolate linearly to find the values in the x series that correspond most closely to the quantiles in the y series interp_x_values = np.interp(y_quantiles, x_quantiles, x_values) y_adjusted = interp_x_values[y_val_index] y_adjusted = np.clip(y_adjusted, 0.0, maxValues) # create f_adjust function that map y to y_adjusted. y_argsort = np.argsort(y) yr0 = y[y_argsort] # < store unadjusted y ya0 = y_adjusted[y_argsort] # < store adjusted y corresponding to yr0 # remove invalid values mask = np.isfinite(yr0)&np.isfinite(ya0) yr0 = yr0[mask] ya0 = ya0[mask]# add bound values, (0.0, 0.0) and (1.0, 1.0) yr1 = np.empty((yr0.shape[0]+2,)) # < yr0 with bound ya1 = np.empty((ya0.shape[0]+2,)) yr1, ya1) to the final (yr, ya) that have 100 points and make the f_adjust smoother. self.yr = np.linspace(0.0, maxValues, 100) self.ya = np.interp(self.yr, yr1, ya1) del yr1 # save few memory del ya1 # save few memory def __call__(self, y_to_adjust): return np.interp(y_to_adjust, self.yr, self.ya)

Table 1 .

 1 List of stations, their geographical coordinates and climate

	Station	Country	Latitude 1	Longitude 1 Altitude asl (m)	Climate (from Peel et al., 2007)
	Toravere	Estonia	58.250	26.467	70	Dfb: Cold climate without dry season and warm summer
	Valentia Obs	Ireland	51.933	-10.250	9	Cfb: Temperate climate without dry season and warm summer
	Wien Hohe Warte Austria	48.250	16.367	203	Dfb
	Auxerre	France	47.800	3.550	207	Cfb
	Payerne	Switzerland	46.815	6.944	491	Dfb
	Carpentras	France	44.083	5.059	100	Csb: Temperate climate with dry and warm summer
	Nice	France	43.650	7.200	4	Csb
	El Arish	Egypt	31.083	33.750	31	BWh
	Maputo	Mozambique -25.967	32.600	70	Aw: Tropical savannah

Table 2 .

 2 Comparison of the accuracies of HC1 and HC3 estimations at each station in the period[2004][2005]

	Station	Mean daily irradiation (J cm -2 )	Bias (J cm -2 )	Standard deviation errors (J cm -2 )	of	Correlation coefficient		Slope (unitless)
			HC1	HC3	HC1	HC3 HC1	HC3	HC1	HC3
	Toravere	1267	-205	-24	265	157	0.941	0.979	1.051	1.003
	Valentia	1109	-22	38	296	162	0.927	0.979	1.081	1.051
	Wien	1212	-285	-10	335	143	0.902	0.983	0.937	1.020
	Auxerre 1
	Payerne	1289	-314	-28	330	166	0.922	0.982	0.888	1.041
	Carpentras	1630	-162	28	264	109	0.954	0.993	0.978	0.967
	Nice	1613	-169	50	236	114	0.960	0.991	0.920	1.000
	El Arish	1838	101	46	237	153	0.951	0.974	1.180	1.067
	Maputo	1816	-59	-15	236	134	0.935	0.978	1.035	0.984
	1 Not enough data from the WRDC for this period

Table 3 .

 3 Accuracy of SSRD-H3E at each station, for the period 1985-2016

	Station	Mean daily irradiation (J cm -2 )	Bias	Standard deviation of errors	Correlation coefficient	Slope (unitless)

Coordinates follow the ISO standard, where positive latitude, respectively longitude, means Northern hemisphere and East of longitude 0°.

Acknowledgments:

The authors thank all ground station operators of the WMO network for their valuable measurements. They additionally thank the World Radiation Data Centre for hosting a website for downloading data. The authors thank the French company Transvalor, which takes care of the SoDa Service for the common good, thus providing an efficient access to the HelioClim databases. They are grateful to Philippe Blanc who designed the spatial low-pass filter used in the construction of the SSRD-H3E data set.


Shintaro Yashima

Satoshi Hirayama

Takayuki Kurokawa Cd

Thomas Salez De

Haruna Takefuji

Wei Hong

Jian Ping Gong
email: gong@sci.hokudai.ac.jp

Shearing-induced contact pattern formation in hydrogels sliding in polymer solution

The contact of a hydrogel during the rotational shearing on glass surface in concentrated polymer solution was observed in situ. Dynamic contact patterns that rotate in-phase with the rotational shearing of the gel were observed for the first time. The contact patterns with a periodicity in the circumferential direction appeared and became fine with the shearing time. The patterns appeared more quickly at elevated sliding velocity, polymer concentration, and normal pressure. Furthermore, the softness of the gel also substantially influenced the character of the patterns. The pattern formation was discussed in terms of the non-linear rheology of the polymer solution at the rotational soft interface.

Introduction

Biological systems display fascinating low friction properties, as find in human joints and eyes. The friction events in the biological systems take place between water-containing soft tissues intermediated by viscoelastic synovial fluid or mucus. [1][2][3][4][5] Hydrogels are a class of material similar to biological soft tissues, in terms of their water containing structure and softness. Hydrogels are therefore used as model substances to elucidate the low friction mechanism of soft tissues as well as biomaterials such as artificial cartilages and meniscus [6][7][8][9][10][11] . Previous researches have shown that hydrogels exhibit very rich friction behaviors, in which solid friction and fluid dissipation coexist. The surface chemistry, surface topology, bulk mechanical property of the hydrogels, and the polymer lubricants are all play roles in intertwined manner. [12][13][14][15][16][17][18][19][20] To elucidate the effect of synovial fluid, hydrogel friction in polymer solutions have been studied previously. 14,17 It has been found that the friction at low sliding velocity, at which the interfacial interaction is the dominant mechanism, decreases with the polymer concentration in the dilute and semi-dilute regimes, reaching a minimum value for a concentration around 10C*, and then the friction increases with the polymer concentration in the concentrated region. Here, C* is the overlap concentration of the polymer solution. Such an effect of polymer solution has been observed not only for sodium hyaluronate (HA) that is the main component of synovial fluid but also for poly(ethylene oxide), of various molecular weights, suggesting the generic effect of polymer solution as lubricants. 14,17 The polymer concentration dependence could be associated to two opposite effects to the hydrogel friction. The polymer at the interface screens the direct contact of the hydrogel to the substrate, which reduces the friction. On the other hand, the enhanced viscosity of the lubricant layer increases dissipation. As the viscosity of the polymer solution rapidly increases with the concentration above C*, the viscosity effect overwhelms the screening effect at high concentration. Direct observation of interfacial contact during the friction process is indispensable for elucidating the effect of polymer solution on the friction of hydrogels. In this work, we in-situ observed the macroscopic contact at the sliding interface of polyvinyl alcohol (PVA) hydrogels on glass substrates in concentrated HA solutions. The friction experiment was performed using a strain-controlled parallel-plate rheometer. The sliding interface was observed using a recently developed optical system based on the principle of critical refraction. Unique contact patterns with circular periodicity appeared with the sliding time. The effect of polymer concentration, sliding velocity, normal load, and modulus of the gel on the characters of the dynamic patterns were investigated. The probable mechanism for the contact pattern formation is discussed using the concept of gel/glass surface instability. [21][22][23][24] This journal is © The Royal Society of Chemistry 2019

Soft Matter

Please do not adjust margins Please do not adjust margins Methods Sample preparation Hydrogels: The physically crosslinked polyvinyl alcohol (PVA) gels were used in this study. The gels were prepared by the quenching method from PVA solution of mixture solvent (dimethyl sulfoxide:H 2 O=3:1, w/w). Unless specified, the PVA solution of 10 wt% concentration was prepared by heating PVA (Mw=88,000, Polymerization degree=2,000, Nacalai Tesque Co., Ltd) in the mixture solvent for 1 hour at 90˚C. After heating, the degassed PVA solution was poured into the reaction cell that consists of two glass plates with a silicone rubber spacer (3.0 mm) to form a hydrogel sheet, and this hydrogel sheet was quenched at -40˚C for 16 h. The hydrogel sheet was then allowed to warm up to room temperature and subsequently immersed in a large amount of distilled water for one week to extract the dimethyl sulfoxide. The water content, thickness, and Young's modulus E of the equilibrium-swollen gels were 19.5 wt%, 2.4 mm, and 40 kPa, respectively. The Young's modulus was measured by using a tensile-compressive tester (Tensilon, Orientec., Co.). The disc-shaped sample of 10 mm in diameter and 2.4 mm in thickness was compressed with a velocity of 10% thickness min -1 . Using the relation of strain and nominal stress, the Young's modulus was calculated at a strain of less than 10%. The modulus data was averaged over 3 samples. Polymer solution as a lubricant: HA with molecular weight of 1.9x10 6 g/mol was used. The HA aqueous solutions of various concentrations (10C* ~60C*) containing 0.1 M NaCl to inhibit aggregation of HA and 0.2 g/L sodium azide for suppression of bacterial growth were prepared. In 0.1 M NaCl, the overlap concentration, C*, of HA is 0.260 g/L. 25 In the studied concentration range of 10C* ~ 60C*, the HA concentrations were 0.26 wt% ~ 1.56 wt%. At the C* concentration, the osmotic pressure from the polymer is 1.27 Pa and from the Donnan effect is 2.58 Pa. 25 Accordingly, the overall osmotic pressures of the HA solutions were 31 Pa ~ 230 Pa, which were much lower than the modulus of the PVA gels. This ensures that the PVA gels were not dehydrated in the HA solution, and furthermore, it also ensures a large enough difference in the refractive indices between the gel and the HA solution, which is necessary for the optical observation of the contact interface.

Glass surfaces: Hydrophobically treated cover glasses (Micro cover glass, C050701; Matsunami Glass Ind. Ltd.) were used as rigid counter surfaces for the observation during sliding friction. For the surface treatment, these glass pieces were soaked in alkali solution for 1 hour and washed with water. Furthermore, an ozone-cleaning chamber (UV/O 3 ) was used for extra surface cleaning. To prepare the hydrophobic surface, the cleaned glasses were further exposed to binding silane vapor, 1H, 1H, 2H, 2H -Perfluorodecyltrichlorosilane (FDTS), at 20 kPa for 8 hours in a vacuum desiccator. The contact angle of water on the cover glass with FDTS coating was 110 ± 2˚.

Measurements

Friction measurement and in situ observation setup of the gel-glass interface: The friction measurement and in-situ observation of the gel-glass interface were performed at the same time. The friction of hydrogels against glass surfaces was measured using a rheometer (ARES, Rheometric Scientific F. E. Ltd.) that works in a constant compressive strain mode at 25 o C. The rheometer, with parallel-plate geometry, was equipped with a home-made optical system, which allows the in situ observation of the gel-glass interface (Fig. 1). 26 The equilibrium-swollen gels were cut into a disc shape of R=7.5 mm in radius and glued on the lower plate of the rheometer using cyanoacrylate instant adhesive agent (Toa Gosei Co., Ltd.). Small amount of polymer aqueous solution of prescribed concentration as a lubricant was poured on the gels to cover their whole surface. As the counter surface, the cover glass treated as described above was attached to the bottom of the prism surface, and the prism was attached to the upper plate of the rheometer. The counter surface was brought towards the gel surface at a constant speed of 2 m/s until the initial normal force reached the prescribed value. The rotation of the lower plate attached to the gel was started immediately after the normal force reached the prescribed value. The direction of rotation was clockwise in all measurements. The angular velocity  of the lower plate was varied in the range of 1.3 rad/s to 13.5 rad/s. The rotation lasted for a prescribed time at a constant angular velocity. The images of the frictional interface were recorded in situ by using a digital video camera (HDR CX550, Sony) set at an observation angle that is larger than the critical angle of the polymer solution but smaller than the critical angle of the gel. The video camera was located 1.5 m away so that the whole sample could be observed at the observation angle. In the contact images, the bright regions correspond to the gel in direct contact with the glass substrate, and the dark regions correspond to a gap with polymer lubricants. The friction force, F, was calculated as F = 4T/3R, where T was the frictional torque recorded during measurement and R (=7.5 mm) was the radius of the apparent contact area. 27 Here, we assumed that friction is linearly proportional to the velocity, and the torque T was the sum over the velocity range of 0 (center of axis) to R (periphery of gel). It should be noted that the "friction" described here refers to all the dissipative force at the gel/glass interface, including i) the dry friction at the gel/glass interface and ii) the viscous resistance of lubricant. Therefore, we represented friction stress the "shear stress at the gel/glass interface" for clarity. The average shear stress, , generated at the gel/glass interface was defined as the friction force per unit area,  = F/πR 2 .

Viscosity measurement:

The viscosity of the HA solutions was determined using a commercially available rheometer (ARES-G2) with cone-plate geometry of 50 mm in diameter at 25 o C, covering a shear rate range from 10 -3 to 10 3 s -1 .

Results

To reveal the rheological behaviors of the polymer solutions used as the lubricant, the shear rate dependency of the shear stress and viscosity of HA at various concentrations (10C* ~ 60C*) was investigated. As shown in Fig. 2a, the shear stress τ increases with the shear rate γ and saturates at high shear rate, indicating shear thinning. For the same shear rate, the shear stress increases with the HA concentration. As shown in Fig. 2b, all the samples exhibit shear thinning, where the viscosity begins to decrease as the shear rate is increased beyond a typical threshold value. With the increase of HA concentration, the threshold decreases, indicating that shear thinning occurs earlier for more concentrated HA solutions. We can extrapolate the rheograms to higher shear rate regions beyond the experimental limitation by using the Casson's equation [START_REF] Casson | Rheology of Disperse System[END_REF][START_REF] Graham | [END_REF] , .

(

Here k 0 and k 1 are constants, depending on the polymers and their concentrations. k 0 and k 1 are determined from the experimentally measured relationship between τ and γ at high shear rate as shown in Fig. 2c. Table 1 shows the k 0 and k 1 values for HA solutions of various concentrations. Fig. 3a shows the time evolution of the contact pattern at the gel-glass interface for the gel rotationally sliding on the hydrophobic glass in a 30C* HA solution. The images in Fig. 3a were captured from the video in the supplemental information (Video S1). At the very beginning of the sliding, a dark image is observed except at the periphery of the sample, indicating that the HA liquid is entrapped at the interface and the gel has no direct contact with the glass surface except at the edge of the sample. With time, a circular contact pattern with a periodicity in the circumferential direction appears gradually, starting from the periphery and developing towards the center of the sample. The contact pattern rotates in-phase with the rotation of the gel. Together with the contact pattern formation, the total contact area increases. This suggests that the squeezing of lubricants from the interface occurs during sliding. The contact pattern remains even 10 min after stopping the rotation. As shown in Fig. 3b, the corresponding frictional stress increases gradually with time. This increase of friction can be attributed to the increase of contact area during the pattern formation. As shown in Fig. 3a, initially there is no direct contact of the PVA gel on the substrate, and the dissipation is thus due to the viscous lubrication of the HA polymer. So, the frictional stress  originates from the shear stress τ of the fluid lubricant, that is,  at the very beginning of shearing (t ~ 0). This fact permits us to estimate the shear rate γ applied on the HA lubricant and then the average thickness h of the lubricant film from the measured frictional stress in Fig. 3b, using the rheograms in Fig. 2a. The frictional stress  at the very beginning of shearing (t ~ 0) at which the contact is hardly formed (Fig. 3b) is several times higher than the highest shear stress of the HA solution in Fig. 2a, and the average shear rate applied on the intermediated polymer lubricant is much higher than that of the upper limit of the experiments in Fig. 2a, even at the beginning of shearing. Accordingly, we can correlate the frictional stress with the shear rate γ through the Casson's equation ( 1), and the lubricant film thickness h is related to γ as, ,

k 0 [Pa 1/2 ] k 1 [Pa
which varies with the distance r from the center of the sample due to the rotational shearing. h should be a function of time and position. Assuming that h is independent of the distance r from the center, and that the frictional stress all comes from the viscous drag of the lubricant, we estimated the apparent lubricant film thickness (h) at The contact pattern of circular periodicity is only observed when the viscous HA solution is used as a lubricant. When water is used as a lubricant, the entire surface of the gel contacts with the substrate immediately, without showing any regular contact patterns even after certain sliding times (Supporting Information, Fig. S1). The frictional stress in water is much higher than that in the HA solution due to direct interaction of the PVA gel with the glass substrate. 17 The character and dynamics of the contact pattern strongly depends on the concentration of the HA solution and the sliding velocity. Fig. 4a shows a table of images of contact patterns for various HA concentration and angular velocities observed 90 seconds after the starting of the sliding motion. For a 10C* concentration, only two thick concentric circles are observed, and this pattern becomes clearer at high angular velocity. For a 20C* concentration, two concentric circles superimposed with periodic patterns in the circumferential direction appear. The circumferential pattern becomes finer with smaller domain size, but the concentric circles become less clear at high velocity. With the further increase of HA concentration to 30C*, the two concentric circles merge together and the circumferential pattern develops into even finer domain sizes. For 60C*, only two thin concentric circles with a circumferential pattern are observed. Apparently, the squeezing of the liquid from the interface is too slow to form the contact at the given observation time (90 seconds), due to the very high viscosity of the 60C* solution. The table of images of contact patterns with varied HA concentrations, and angular velocities observed at 150 seconds is shown in Fig. 4b. Essentially, Fig. 4b shows similar patterns as those observed in Fig. 4a for 90 seconds but the patterns display finer characters, especially for 60C*. Each contact pattern becomes finer with the increase of the angular velocity when compared at the same sliding time. As the concentration increases, that is, as the viscosity of the lubricant increases, the period in the circumferential direction tends to become smaller. Comparing Fig. 4a and Fig. 4b, it is apparent that the total contact area increases with the elapsed time, which also indicates the squeezing of lubricants from the interface as seen from Fig. 3a. The squeezing effect is more prominent for the low HA concentration and becomes less remarkable for high HA concentrations, which can be attributed to the increase of viscosity at high HA concentration.

Fig. 4c shows the friction behavior corresponding to the gel/glass interface images shown in Fig. 4a. The frictional stress increases weakly with the angular velocity for a fixed polymer concentration, which is a result of the coexistence of the solid friction generated at the contact region and the fluid resistance generated at the non-contact region. The slightly higher friction in 10C* concentration than that in 20C* at high velocity suggests that the solid friction dominates in the case of 10C*. The increase in the frictional stress with the polymer concentration from 20C* to 60C* indicates that the viscous drag of the polymer solution dominates the frictional stress at the observed conditions.

These results are consistent with Fig. 4a that shows large and homogeneous contact areas in the circular direction in the 10C* lubricant while small and patterned contact areas are observed in the 20C*, 30C*, and 60 C* lubricants. As shown in Fig. 4a, in the case of 60 C *, the gel is barely in contact with the glass substrate as compared with other concentrations. The frictional stress must come from the viscous resistance of the fluid. However, the frictional stress has weak velocity dependence as shown in Fig. 4c, one can deduce that strong shear thinning occurs in the 60C* system. Please do not adjust margins

Please do not adjust margins

To investigate the squeezing effect during sliding, the experiment was further performed at various normal pressures. Fig. 5 is the table of images showing the effect of normal pressure at an angular velocity of 5.3 rad/s and observed at 90 seconds. Patterns appear only at relatively large normal pressures (> 22 kPa) at the given observation time, and the patterns become finer with a smaller domain size when increasing the normal pressure. The result indicates that high normal pressure accelerates the squeezing of lubricant from the interface, and therefore, accelerates the contact pattern formation.

To study how the elasticity of the gels plays a role, we further performed the experiment using two gels of quite different moduli. Fig. 6 shows the time evolution of the contact images for soft and hard gels with Young's moduli of 27 kPa and 300 kPa, respectively. Under the same normal pressure (33 kPa), large part of the soft gel already formed contact, while the hard gel only had some contact at the periphery. For the soft gel, applying a shear rotation immediately induces contact patterns in all the pre-contacted area. While it took a longer time for the hard gel to develop the contact pattern towards the inner region. In addition, the contact patterns of the soft gel exhibit larger sizes with characters strongly twisted by the shear stress imposed on the soft gel. These results clearly show that with a soft gel, the contact pattern forms more quickly, displaying clear periodical characters.

Discussion

The dynamic contact pattern formation indicates the occurrence of an instability at the interface. Shearing is found to be essential for this instability, as we do not observe any pattern at a constant normal compression without any shear. This instability is not related to the Schallamach wave, in which the pattern travels in the opposite direction with the shearing, together with a stick-slip motion. 30 In the present case, the shearing and pattern rotations are in phase and there is no stick-slip like motion. Also, the shear velocity is too low for inertia to play any role. To see if the instability is governed by the local pressure variations, we compare the ratio of the frictional stress  to the normal pressure P, /P, to the ratio of the lubricant film thickness h to the wavelength of the pattern , h/, according to the Stokes equation in the lubrication approximation. [START_REF] Blossey | Thin liquid films: dewetting and polymer flow[END_REF] Reading from Fig. 3 and Fig. 4, /P is on the order of 0.1-0.01. The thickness h is roughly estimated as ~1 m from the Casson's plot, following the literature (see supporting information), while  reaches several millimeters, which gives h/ less than 0.001. This means that the patterns we observed are at least one order of magnitude larger than expected from the pressure variations.

Based on the experimental observations, we believe that there are two types of (interrelated) instabilities in the system, both involving the non-Newtonian rheology of the HA fluids. The primary mode of instability appears as annular contact patterns, as seen in the cases of low HA concentrations (Fig. 4a), and is believed to be a direct consequence of the shear thinning rheology of the fluid. For a non-Newtonian fluid with a rheology similar to the Casson's model, there exists a critical shear rate (or a yield stress) beyond which shear thinning occurs. The rotation in the rheometer provides a large enough circumferential shear such that the fluid in most regions (except the center) is in the shear thinning regime. In other words, the rotation yields the fluid, and enables significant gap flow in the radial direction, driven by the relatively small pressure gradient. Moreover, under a constant gap thickness, the circumferential shear rate is proportional to the distance to the rotation center, and thus the viscosity of the shear thinning fluid is smaller at the outer perimeter.

The resulting higher radial flow rate near the outer perimeter tends to deplete the HA solution in the intermediate-radius region while leaving a pocket of fluid in the center, leading to annular contact patterns. With the driving force for the radial flow being the radial pressure gradient and the resistance being related to the viscosity of the HA solution, the process has two dimensionless parameters, and , where is the initial thickness of the fluid, and  the angular rotation velocity. The pressure gradient buildup is related to those two dimensionless parameters, e.g. it is expected that solutions of higher (higher HA concentration) and/or rotating at a higher angular velocity tend to exhibit larger pressure gradients.

The secondary instability generates the radial-like channels through the annular contact region. The process is similar to the gap flow between a soft material initially in contact with a rigid wall: instead of uniformly opening the gap, the fluid tends to tunnel through the interface by locally deforming the soft solid. The characteristic length (circumferential wavelength)  of the instability results from the competition between the shear elastic modulus of the soft solid and the radial pressure gradient: Comparing to that caused by rotation, the shear rate induced by the radial flow is relatively small and thus the effect of shear thinning is minor in this secondary mode of instability. Nevertheless, the characteristic length  is indirectly related to the rheology of the fluid through the pressure gradient

. For example, solutions of higher HA concentrations, or under higher angular velocities, will lead to smaller feature sizes in the radial channels of the contact patterns. These trends agree qualitatively with the experiments, while the detailed model of the process will be addressed in a separate article due to the mathematical complexity of this highly nonlinear problem.

Conclusions

Dynamic contact patterns displaying periodicity in the circumferential direction are observed during the rotational sliding friction of hydrogels on solid substrate in concentrated polymer solutions. These patterns are formed in the region where the hydrogel and the glass are in contact in the polymer solution under shear stress. Although it belongs to the broad class of pattern formations under shear, it differs from phenomena such as shear banding in that the polymer solution itself forms a pattern. The appearing time and This journal is © The Royal Society of Chemistry 2019 Soft Matter Please do not adjust margins Please do not adjust margins character of the contact patterns are dependent on the polymer concentration, the sliding velocity, the normal pressure, and the modulus of the gel. More concentrated polymer solutions and higher sliding velocities produce contact patterns with finer characters. Higher pressures and lower moduli lead to a quicker pattern formation. The results suggest that the instability is related to the shear thinning of the HA solutions as lubricants of the soft hydrogels under rotational shearing.

Fig. 1

 1 Fig.1Schematic illustration of the setup for the in-situ observation of the frictional interface between a gel and a glass during sliding friction measurement. (a) A parallel-plate rheometer was used. The trapezoidal prism was fixed to the upper plate of the rheometer, and a cover glass was fixed to the prism as the counter surface of friction. Hydrogel was bonded to the rotating lower plate of the rheometer. (b) The hydrogel was compressed in the polymer solution by the upper plate with a prism. The images of the frictional interface were recorded in situ by using a digital video camera set at an observation angle between the critical angles of polymer solution and gel.

Fig. 2 Fig. 3 (

 23 Fig. 2 Dependences of (a) the shear stress and (b) the viscosity, with the shear rate in HA solutions, for various concentrations measured by rheometer at 25 o C, as indicated. (c) Casson's plots to determine constants k 0 and k 1 for HA solutions. C* is the overlap concentration (i.e. the onset of the semi-dilute regime). HA molecular weight: 2x10 6 g/mol.

Fig. 4 (

 4 Fig.4 (a, b) HA concentration and angular velocity effects on the contact pattern at the hydrogel/glass interface, and (c) the corresponding shear stresses. The results are for (a) 90 seconds and (b) 150 seconds after the beginning of the sliding motion. Gel modulus: 40 kPa; normal pressure: 22 kPa.

Fig. 5 Fig. 6

 56 Fig. 5 HA concentration and normal pressure effects on the contact pattern at the hydrogel/glass interface, observed 90 seconds after the beginning of the sliding motion. Gel modulus: 40 kPa; angular velocity: 5.3 rad/s.

.

Table 1 .

 1 k 0 and k 1 of Casson's equation for HA aqueous solutions (molecular weight 1.9x10

6

g/mol) at 25°C.

This journal is © The Royal Society of Chemistry 2019 Soft Matter

Acknowledgements

This work was supported by JSPS KAKENHI Grant Numbers JP24120003, JP17H06144, and JP17H06376. Institute for Chemical Reaction Design and Discovery (ICReDD) was established by World Premier International Research Initiative (WPI), MEXT, Japan.

Conflicts of interest

There are no conflicts to declare.

Supporting Information

Shearing-induced contact pattern formation in hydrogels sliding in polymer solution Shintaro Yashima, Satoshi Hirayama, Takayuki Kurokawa*, Thomas Salez, Haruna Takefuji, Wei Hong, Jian Ping Gong*

Estimation of lubricant film thickness h

We assume that lubricant film thickness h is independent of the distance r from the center of the gel disc, h(r)=h. The frictional force F is given by:

where  is the shear stress and R is the radius of the gel disc.

The shear rate  is written as:

Assuming the shear rate  to be high, by inserting Casson's equation (1) in equation (S1), the frictional force can be expressed as:

. (S3)

The friction stress can be written as:

. (S4)

Rewriting this quadratic equation (S4) by setting x=h -1/2 , one gets:

. (S5)

The lubricant film thickness h is then given as the solution of the latter:

. (S6)

Please do not adjust margins Video S1 Video image of contact pattern at hydrogel/glass interface under sliding motion at angular velocity 5.3 rad/s. Gel modulus: 40 kPa; lubricant: 30C* HA; normal pressure: 22 kPa. The images shown in Fig. 3 were extracted from the Video S1.

Video S2 Video image of contact pattern at hydrogel/glass interface under sliding motion at angular velocity 13 rad/s. Gel modulus: 40 kPa; lubricant: 30C* HA; normal pressure: 22 kPa. The images shown in Fig. 4a and Fig. 4b were extracted from the Video S2.