A Lacour

H Chepfer

N B Miller

M D Shupe

V Noel

X Fettweis

H Gallee

J E Kay

R Guzman

And J Cole

How Well Are Clouds Simulated over Greenland in Climate Models? Consequences for the Surface Cloud Radiative Effect over the Ice Sheet

Using lidar and radiative flux observations from space and ground, and a lidar simulator, we evaluate clouds simulated by climate models over the Greenland ice sheet, including predicted cloud cover, cloud fraction profile, cloud opacity, and surface cloud radiative effects. The representation of clouds over Greenland is a central concern for the models because clouds impact ice sheet surface melt. We find that over Greenland, most of the models have insufficient cloud cover during summer. In addition, all models create too few nonopaque, liquid-containing clouds optically thin enough to let direct solar radiation reach the surface (21% to 23.5% at the ground level). Some models create too few opaque clouds. In most climate models, the cloud properties biases identified over all Greenland also apply at Summit, Greenland, proving the value of the ground observatory in model evaluation. At Summit, climate models underestimate cloud radiative effect (CRE) at the surface, especially in summer. The primary driver of the summer CRE biases compared to observations is the underestimation of the cloud cover in summer (246% to 221%), which leads to an underestimated longwave radiative warming effect (CRE LW 5 235.7 to 213.6 W m 22 compared to the ground observations) and an underestimated shortwave cooling effect (CRE SW 5 11.5 to 110.5 W m 22 compared to the ground observations). Overall, the simulated clouds do not radiatively warm the surface as much as observed.

Introduction

Currently, the Greenland ice sheet is losing mass because summer melt exceeds annual accumulation [START_REF] Rignot | Changes in the velocity structure of the Greenland ice sheet[END_REF][START_REF] Tedesco | The role of albedo and accumulation in the 2010 melting record in Greenland[END_REF]Tedesco , 2013;;[START_REF] Hanna | Increased runoff from melt from the Greenland ice sheet: A response to global warming[END_REF][START_REF] Krabill | Greenland ice sheet: Increased coastal thinning[END_REF][START_REF] Mernild | Greenland ice sheet surface melt extent and trends: 1960-2010[END_REF][START_REF] Vaughan | Observations: Cryosphere[END_REF]. Greenland mass loss has global sea level rise implications. Predicting future changes in Greenland ice sheet mass balance requires understanding regional atmospheric and oceanic processes, including their interactions with the ice sheet and how these interactions change with climate change.

Clouds affect the cryospheric mass budget of the Greenland ice sheet through their radiative effects and also through their impacts on precipitation (e.g., Shupe et al. 2013;[START_REF] Miller | Cloud radiative forcing at Summit, Greenland[END_REF][START_REF] Van Tricht | Clouds enhance Greenland ice sheet meltwater runoff[END_REF]). Yet, determining the magnitude and sign of the cloud radiative effect over Greenland is challenging as cloud radiative effects depend on cloud cover, cloud opacity and phase, cloud altitude, solar zenith angle, and surface albedo. At Summit, Greenland (3212 m above mean sea level), in situ measurements show that clouds radiatively warm the surface in all months of the year [START_REF] Miller | Cloud radiative forcing at Summit, Greenland[END_REF]. During winter, when shortwave cloud cooling is negligible, clouds warm the surface through their emission of longwave radiation. But even during summer, cloud warming through emission of longwave radiation exceeds shortwave cloud cooling [START_REF] Cawkwell | The impact of cloud cover on the net radiation budget of the Greenland ice sheet[END_REF][START_REF] Miller | Cloud radiative forcing at Summit, Greenland[END_REF]. During summer, the presence of opaque clouds prevents direct solar radiation from reaching the surface. Indeed, the high surface albedo of the ice sheet minimizes the shortwave cloud radiative cooling effect. These opaque clouds also absorb and reemit longwave radiation coming from the surface and the lower-tropospheric layers [START_REF] Guzman | Direct atmosphere opacity observations from CALIPSO provide new constraints on cloudradiation interactions[END_REF]).

In the Arctic, low-level liquid-containing clouds play a particularly important role in heating the surface and cooling the atmosphere (e.g., Shupe and Intrieri 2004). Liquid-containing clouds are common in the Arctic, occurring 30%-80% of the time, depending on location and season [START_REF] Cesana | How well do climate models simulate cloud vertical structure? A comparison between CALIPSO-GOCCP satellite observations and CMIP5 models[END_REF]Cesana and Chepfer 2013;[START_REF] Kay | Evaluating and improving cloud phase in the Community Atmosphere Model version 5 using spaceborne lidar observations[END_REF][START_REF] Morrison | Resilience of persistent Arctic mixedphase clouds[END_REF][START_REF] Shupe | Clouds at Arctic atmospheric observatories. Part II: Thermodynamic phase characteristics[END_REF]. Liquid-containing clouds occur frequently over Greenland with important impacts on radiation and surface energy budgets (e.g., [START_REF] Lacour | Greenland clouds observed in CALIPSO-GOCCP: Comparison with groundbased summit observations[END_REF]Shupe et al. 2013). Of particular importance, clouds can trigger surface melt over a large portion of the Greenland ice sheet [START_REF] Bennartz | July 2012 Greenland melt extent enhanced by low-level liquid clouds[END_REF][START_REF] Solomon | Cloud-atmospheric boundary layer-surface interactions on the Greenland ice sheet during the July 2012 extreme melt event[END_REF]. Greenland surface melting increases nonlinearly with increasing temperatures because of positive feedbacks between cloud microphysics, surface melting, and surface albedo [START_REF] Fettweis | Brief communication ''Important role of the mid-tropospheric atmospheric circulation in the recent surface melt increase over the Greenland ice sheet[END_REF]) and modulates the ice sheet mass balance [START_REF] Van Tricht | Clouds enhance Greenland ice sheet meltwater runoff[END_REF][START_REF] Hofer | Decreasing cloud cover drives the recent mass loss on the Greenland ice sheet[END_REF].

To assess the influence of clouds on Greenland surface melt, it is crucial to document cloud cover, opacity, and vertical distribution as well as surface conditions, such as the surface temperature and surface albedo. Often reanalysis datasets such as the European Centre for Medium-Range Weather Forecasts (ECMWF) interim reanalysis (ERA-Interim, hereafter ERAI; [START_REF] Dee | The ERA-Interim reanalysis: Configuration and performance of the data assimilation system[END_REF]) and the NCEP-NCAR reanalysis [START_REF] Kalnay | The NCEP/NCAR 40-Year Reanalysis Project[END_REF] are used to map surface and atmospheric properties over Greenland and assess clouds and precipitation. Reanalysis data are produced through a sequential data assimilation scheme using observations of atmospheric thermodynamic variables (temperature, pressure, humidity) combined with information from a forecast model used to estimate the atmospheric state at each time step. In this context, it is important to realize that clouds and precipitation in reanalyses are modeled quantities and are poorly represented over Greenland (e.g., [START_REF] Griggs | Assessment of cloud cover characteristics in satellite datasets and reanalysis products for Greenland[END_REF]. In other words, reanalysis datasets do not provide reliable cloud and precipitation observations. Direct observations with high spatial and temporal resolution are needed to understand cloud and precipitation processes. Thanks to enhanced observational networks put in place over the last decade, we have both in situ ground observations and satellite observations of Greenland clouds and precipitation. From the ground, the Integrated Characterization of Energy, Clouds, Atmospheric State, and Precipitation at Summit (ICE-CAPS; Shupe et al. 2013) project has collected cloud profiles, radiative fluxes, and precipitation at Summit since 2010. From space, the lidar CALIOP and radar CloudSat [START_REF] Winker | The CALIPSO mission: Spaceborne lidar for observation of aerosols and clouds[END_REF][START_REF] Stephens | The CloudSat mission and the A-Train: A new dimension to space-based observations of clouds and precipitation[END_REF] provide vertically resolved (30 and 500 m) and horizontally resolved (330 m and 1.7 km) cloud and precipitation properties over nearly the entire Greenland ice sheet since 2006. These active spaceborne sensors document the cloud properties and altitude with high confidence over cold and white surfaces (e.g., [START_REF] Van Tricht | Clouds enhance Greenland ice sheet meltwater runoff[END_REF][START_REF] Lacour | Greenland clouds observed in CALIPSO-GOCCP: Comparison with groundbased summit observations[END_REF]Morrison et al. 2018;[START_REF] Kay | Evaluating and improving cloud phase in the Community Atmosphere Model version 5 using spaceborne lidar observations[END_REF]. Of utility to constraining cloud opacity and thus cloud impacts on radiation, lidar observations also distinguish between optically thin clouds with a visible cloud optical depth less than '5 and opaque clouds with optical depth greater than '5 [START_REF] Guzman | Direct atmosphere opacity observations from CALIPSO provide new constraints on cloudradiation interactions[END_REF]Chepfer et al. 2013).

As climate models add interactive ice sheet components, the representation of cloud and precipitation processes over Greenland is becoming more critically important. Thus, here, we combine in situ and satellite observations available over the last decade to assess Greenland clouds and their radiative effects as simulated by current-generation climate models. Because the values for observed clouds depend on the instruments' properties and on the cloud detection algorithm [START_REF] Stubenrauch | Assessment of global cloud datasets from satellites: Project and database initiated by the GEWEX radiation Panel[END_REF], we use the well-known satellite simulator package the Cloud Feedback Model Intercomparison Project (CFMIP) Observation Simulator Package (COSP; [START_REF] Bodas-Salcedo | COSP: Satellite simulation software for model assessment[END_REF]. The simulators contained within COSP mimic the satellite observation scale and definition and thus enable a direct comparison of the modeled and observed clouds.

In this paper, we first assess cloud cover, cloud opacity, and cloud vertical distribution using the satellite records that provide data over almost the entire Greenland ice sheet (section 3a). Next, we evaluate the surface cloud radiative effects at Summit station (section 3c). Summit is unique because it is the only place where we have sufficient observations to make a robust assessment of cloud radiative effects over Greenland. As a result, it is important to asses both 1) if we can make robust links between cloud properties measured from space and cloud radiative effects at Summit and 2) if Summit is representative of the entire ice sheet for cloud biases (section 3b). Finally, we discuss (section 4) the joint constraint on Greenland cloud radiative effect provided by in situ measurements at Summit and spaceborne observations from CALIPSO. In the end, we provide guidance on the reliability of the current observational network to assess cloud radiative effects over the entire ice sheet.

Data and methods

a. Satellite observations

This study relies primarily on cloud observations from the lidar on the CALIPSO satellite. Specifically, we use the GCM-Oriented CALIPSO Cloud Product, version 3.0 (CALIPSO-GOCCP), an observational product that has been designed for climate studies and evaluation of climate models. More specifically, we use CALIPSO-GOCCP cloud cover maps, cloud fraction profiles [START_REF] Noel | A global view of horizontally oriented crystals in ice clouds from Cloud-Aerosol Lidar and Infrared Pathfinder Satellite Observation (CALIPSO)[END_REF], cloud phase mask (Cesana and Chepfer 2013), and cloud opacity mask [START_REF] Guzman | Direct atmosphere opacity observations from CALIPSO provide new constraints on cloudradiation interactions[END_REF].

In addition to the standard variables available in CALIPSO-GOCCP, we built two new gridded variables specifically for the current study: the opaque cloud fraction profile and the nonopaque cloud fraction profiles. We started from the orbit files containing the instantaneous CALIPSO scattering ratio [SR (CALIPSO-SR; [START_REF] Noel | A global view of horizontally oriented crystals in ice clouds from Cloud-Aerosol Lidar and Infrared Pathfinder Satellite Observation (CALIPSO)[END_REF]] profiles along the satellite orbits at a 480-m vertical resolution with a horizontal resolution of 330 m along track and 90 m cross track. Each SR profile contains the ratio between the attenuated backscatter (ATB) profile measured by the lidar and the ATB mol that would be measured by the lidar if there were only gas molecules and no particles (clouds, aerosols) in the atmosphere. For each SR profile, we use different SR thresholds to label each atmospheric layer as ''nonopaque cloud,'' ''opaque cloud,'' or ''clear'' (Table 1). Then, an ''opaque cloud fraction'' (''nonopaque cloud fraction,'' ''clear fraction,'' respectively) is computed for each vertical level in each 28 3 28 latitude-longitude grid box as the ratio of the number of ''opaque cloud'' layers encountered within this grid box over the total number of layers, Fraction cat (z) 5 å layers, cat (z) å all layers (z)

.

(1)

In doing this, we obtain the opaque cloud fraction profile, the nonopaque cloud fraction profile, and the clear fraction profile for each 28 3 28 grid box. For some profiles, the lidar laser beam is fully attenuated before reaching the ground. The portion of the atmosphere located between the surface and the altitude z_opaque, which corresponds to the lidar full attenuation, corresponds to the layers where no information could be retrieved by the lidar. For simplicity, following the approach by [START_REF] Guzman | Direct atmosphere opacity observations from CALIPSO provide new constraints on cloudradiation interactions[END_REF], we consider these layers to be in the ''opaque'' category because no direct solar radiation can reach these layers. Indeed, the direct solar radiation does not pass through the opaque clouds, so only diffuse shortwave radiation can reach the surface.

For the current study, we applied some specific processing to the CALIPSO-GOCCP dataset. In each 28 3 28 grid box, we adjusted the vertical profiles as a function of altitude above ground level (AGL) instead of above mean sea level, using the digital elevation model TerrainBase (National Geophysical Data Center 1995). TerrainBase does not exactly reproduce the complex Greenland topography [START_REF] Arabelos | Intercomparisons of the global DTMs ETOPO5, TerrainBase and JGP95[END_REF] but since the CALIPSO-GOCCP data resolution is coarse (28 3 28 3 480 m) the elevation error would need to be very large to affect our results (see online supplemental material Fig. S3).

Our analysis uses eight years of CALIPSO-GOCCP data from January 2008 to December 2015. Observations before 2008 were discarded because of the different pointing angle of the CALIPSO lidar at the time [START_REF] Noel | A global view of horizontally oriented crystals in ice clouds from Cloud-Aerosol Lidar and Infrared Pathfinder Satellite Observation (CALIPSO)[END_REF], which affects cloud retrievals. The data covers all of Greenland except the area north of 828N that is not sampled by the satellite. To characterize clouds near the Summit ground station in the middle of Greenland, we extract the 28 3 28 cells with surface elevation above 3000 m above mean sea level. This area covers roughly 100 000 km 2 . The CALIPSO-GOCCP observations above this area have been compared with the local ground observations at Summit in [START_REF] Lacour | Greenland clouds observed in CALIPSO-GOCCP: Comparison with groundbased summit observations[END_REF].

Lidar attenuation can prevent cloud detection. In CALIPSO-GOCCP, the cover of opaque clouds (no surface echo detected) is equivalent to the percentage of profiles that are fully attenuated at the surface. Based on the surface echo detection, we found that only 5% of the space lidar profiles are fully attenuated over Summit throughout the year except in July and August (see supplemental material Fig. S1). In July and August, up to 50% of the profiles are fully attenuated, probably because of optically thick liquid-containing clouds frequent at Summit during that period [START_REF] Miller | Cloud radiative forcing at Summit, Greenland[END_REF]. This attenuation occurs below 2 km AGL on average and is coincident with the presence of liquid-containing clouds (see supplemental material Fig. S2).

b. In situ observations of CREs at Summit

We use observations of radiative fluxes and clouds from a comprehensive site at Summit. The total surface cloud radiative effect (CRE tot ) is defined as the difference between the surface net flux in all sky conditions and the surface net flux in cloud-free (clear) conditions in the longwave and shortwave domains,

CRE surf tot 5 (F net surf all sky 2 F net surf clear sky ) LW 1 (F net surf all sky 2 F net surf clear sky ) SW , (2)
and with the net fluxes defined as the difference between the downward and upward fluxes,

F net surf LW=SW, all sky 5 F Ysurf LW=SW, all sky 2 F [surf LW=SW, all sky . (3)
At the surface, the longwave and shortwave radiative fluxes are retrieved from the ground-based observations at Summit at 30-min temporal resolution [START_REF] Miller | Cloud radiative forcing at Summit, Greenland[END_REF](Miller et al. , 2017) ) [START_REF] Miller | Cloud radiative forcing at Summit, Greenland[END_REF]. Clear sky fluxes are estimated using a radiative transfer model (Rapid Radiative Transfer Model;[START_REF] Clough | Atmospheric radiative transfer modeling: A summary of the AER codes[END_REF] constrained by observations of surface temperature, surface albedo, and profiles of temperature, humidity, and ozone. The temperature and humidity profiles come from twice daily radiosonde measurements and the ozone profiles from weekly ozonesondes measurements.

From space, we examined two surface radiative flux products: CERES EBAF-surface [START_REF] Kato | Surface irradiances consistent with CERES-derived top-of-atmosphere shortwave and longwave irradiances[END_REF] based on the CERES instrument on board MODIS and level 2B fluxes and heating rates with lidar (2B-FLXHRlidar; [START_REF] Henderson | A multisensor perspective on the radiative impacts of clouds and aerosols[END_REF] based on CloudSat and CALIPSO observations. However, we found significant disagreements between these datasets over Greenland. The CERES EBAF-surface product shows a negative bias of the surface shortwave cloud radiative effect (CRE SW ) relative to models all year long and the 2B-FLXHR-lidar product shows a positive bias relative to models of the surface CRE SW (not shown). We concluded these products cannot currently be used to evaluate biases in models at this location. Another study would be necessary to properly understand the differences between those two datasets, which is outside the scope of this paper.

c. Model evaluation method

The definition of a cloud differs between climate models and satellite observations and between satellite instruments (e.g., [START_REF] Stubenrauch | Assessment of global cloud datasets from satellites: Project and database initiated by the GEWEX radiation Panel[END_REF]) and between climate models. To use a consistent definition of cloud to allow a more consistent comparison between clouds simulated by climate models and observed by satellite, CFMIP has developed COSP [START_REF] Bodas-Salcedo | COSP: Satellite simulation software for model assessment[END_REF], which is composed of satellite simulators that mimic clouds as they would have been observed by different satellites, if they observed the atmosphere predicted by a climate model. Within COSP is a CALIPSO simulator [START_REF] Chepfer | Use of CALIPSO lidar observations to evaluate the cloudiness simulated by a climate model[END_REF][START_REF] Chiriaco | The ability of MM5 to simulate ice clouds: Systematic comparison between simulated and measured fluxes and lidar/radar profiles at the SIRTA atmospheric observatory[END_REF]) that mimics the unique information on the cloud vertical distribution observations at a resolution of 480 m that is provided by CALIPSO-GOCCP. The current paper uses the following standard gridded COSP, version 1.4 (v1.4), outputs: the cloud cover, the cloud fraction profiles, and the height-intensity histograms [START_REF] Chepfer | Use of CALIPSO lidar observations to evaluate the cloudiness simulated by a climate model[END_REF][START_REF] Noel | A global view of horizontally oriented crystals in ice clouds from Cloud-Aerosol Lidar and Infrared Pathfinder Satellite Observation (CALIPSO)[END_REF]. COSP outputs from the GCMs have been compared against observations in various regions of the globe to evaluate clouds simulated by climate models. The models present common documented biases. They tend to create too many highly reflective clouds in the tropics [START_REF] Nam | The 'too few, too bright' tropical low-cloud problem in CMIP5 models[END_REF][START_REF] Klein | Are climate model simulations of clouds improving? An evaluation using the ISCCP simulator[END_REF] but not enough in the Southern Ocean (Bodas-Salcedo et al. 2014). They also have difficulties correctly representing the cloud vertical distribution [START_REF] Cesana | How well do climate models simulate cloud vertical structure? A comparison between CALIPSO-GOCCP satellite observations and CMIP5 models[END_REF]) and the cloud thermodynamic phase (Cesana and Chepfer 2013). Between 2408 and 08C, the models struggle to simulate coexisting liquid and ice clouds, resulting in insufficient liquid cloud cover in the Arctic [START_REF] Cesana | How well do climate models simulate cloud vertical structure? A comparison between CALIPSO-GOCCP satellite observations and CMIP5 models[END_REF][START_REF] Kay | Evaluating and improving cloud phase in the Community Atmosphere Model version 5 using spaceborne lidar observations[END_REF]. While many regions have been studied with the help of satellite simulators and observations, clouds over Greenland have not been studied yet. To do this comparison, we use COSP outputs from the models participating in the CFMIP, phase 2 (CFMIP-2), which was part of phase 5 of the Coupled Model Intercomparison Project (CMIP5; [START_REF] Taylor | An overview of CMIP5 and the experiment design[END_REF]). Our study is limited to the eight models that provided COSP/lidar outputs (Table 2). The eight models provided the cloud cover and the cloud fraction profiles through the COSP/lidar simulator. Four of these eight models also provided the height-intensity histograms (SR profiles). From these four, we computed the simulated opaque and nonopaque cloud fraction profiles, following the same method as for the observations (section 2a). Cloud opaque and nonopaque categories are classified the same way as in the satellite observations from the COSP outputs. Nevertheless, unlike CALIPSO-GOCCP, which looks for the surface echo to detect cloud opacity, in climate models, lidar attenuation is defined when SR , 0.01. This difference means the opaque fraction could be underestimated in the climate models by 5% (supplemental material Fig. S4). To evaluate surface CRE and surface melting, we use climate model outputs of the all sky and clear sky longwave and shortwave surface fluxes and the near-surface temperature.

Results

a. Model-observation comparisons over the entire Greenland ice sheet

We begin by comparing total cloud cover over Greenland in the spaceborne lidar observations with the climate models. Results are presented both for winter [December-February (DJF)], and for summer [June-August (JJA)] using maps (Figs. 1,2) and averages over the entire Greenland ice sheet (Table 2). The mean observed total cloud cover in winter (DJF; Fig. 1a) is 56% over Greenland. In summer (JJA; Fig. 2a), the mean observed total cloud cover increases to 64%, locally reaching 80% in the central region, largely because of the frequent occurrence of liquid-containing clouds [START_REF] Lacour | Greenland clouds observed in CALIPSO-GOCCP: Comparison with groundbased summit observations[END_REF]. Many models underestimate cloud cover over Greenland during both winter (Fig. 1) and summer (Fig. 2); however, the spread in both the sign and the magnitude of the biases among the eight models that were evaluated is large. In winter (Figs. 1c-j We next evaluate the cloud fraction profiles in the models using the spaceborne lidar observations (Fig. 3). Four models (CanAM4, HadGEM2-A, MRI-CGCM3, and CAM5; full lines in Fig. 3) have small cloud fraction profile biases (,5%) in both winter and summer. Below 2 km AGL in winter and below 4 km AGL in summer, most models underestimate the cloud fraction compared to lidar observations. Despite the diversity of cloud fraction biases identified, most models have more clouds in summer than in winter as in the observations. In the midtroposphere, the cloud fraction difference between winter and summer is consistent with observations, except for the Meteorological Research Institute (MRI) model, simulating more high clouds in winter than in summer. At low altitudes, CanAM4 and the MPI model only simulate more low clouds in winter than in summer. Overall, the comparisons suggest that models produce too much high-level cold clouds and not enough lowlevel clouds.

Next, we compare observed and modeled cloud opacity diagnostics. We remind the reader that because the cloud opacity diagnostics we are using are relatively new [START_REF] Guzman | Direct atmosphere opacity observations from CALIPSO provide new constraints on cloudradiation interactions[END_REF]), we can only analyze cloud opacity diagnostics in a small subset of climate models in which they are currently available (CanAM4, HadGEM2-A, MRI-CGCM3, and CAM5). We compare mean clear sky, opaque cloud, and nonopaque cloud fractions near the surface, at 480 m AGL, over all Greenland (Fig. 4). During winter (Figs. 4a,c), on average the near-surface clear sky fraction is above 70% in both models and observations. CanAM4 has the most underestimated near-surface clear sky fraction, but also has the biggest standard deviation, because of its coarse 2.88 3 2.88 grid and smaller sample size. In the Greenland-wide observations, the near-surface opaque fraction is larger than nonopaque cloud fraction. On average, the GCMs tend to underestimate the opaque fraction, although given the uncertainties, the opaque fraction values in CanAM4 and HadGEM2-A are closer to the observed values. Generally CAM5 produces far too few opaque clouds. In JJA (Figs. 4c,d), the cloud fraction (nonopaque and opaque) is higher than in DJF. The observed increase in near-surface opaque fraction from DJF to JJA is reproduced in each model for all of Greenland for both. In the observations, we also observe an increase in the nonopaque fraction in summer, but this increase is generally not reproduced in the models.

After looking at the cloud fraction at the ground level, we evaluate the vertical distribution for each cloud category in CALIPSO-GOCCP and in climate models (Fig. 5). To have more details on how well models represent clouds, we divide the nonopaque clouds category into two new categories: those with SR , 30 related to optically thinner nonopaque clouds and those with SR . 30 related to optically thicker nonopaque clouds. This distinction gives more insight into the implications of the cloud opacity model biases. At low levels, both models and observations show a larger opaque fraction Only one model (CanAM4) creates optically thicker nonopaque clouds (SR . 30) in DJF. In JJA, only CanAM4 creates optically thick nonopaque clouds at all levels, while the other models create some near the ground. However, all models underestimate the thick nonopaque cloud fraction at all levels.

Finally, we evaluate cloud thermodynamic phase (liquid/ice) in the models (Fig. 6). We remind the reader that cloud phase detection for our diagnostics in both the observations and the models is based on the polarization of the lidar signal (Cesana and Chepfer 2013). As expected, liquid clouds are more often opaque than nonopaque (Figs. 6a,e,b,f). The liquid clouds are located mostly below 2 km AGL in winter and below 4 km AGL in summer. A distinction between the optically thinner (SR , 30) and thicker (SR . 30) nonopaque clouds gives more information on the distribution of the ice and liquid clouds in the nonopaque cloud category (Figs. 6c,g,d,h). For nonopaque clouds, those with SR . 30 are much more likely to contain liquid water in both seasons, and almost all cases with SR . 30 below 2 km have liquid nonopaque clouds. The difficulty of the models to create both opaque clouds and the thicker nonopaque clouds (SR . 30) over Greenland tends to confirm the known difficulty for the models to create liquid clouds in the Arctic [START_REF] Cesana | How well do climate models simulate cloud vertical structure? A comparison between CALIPSO-GOCCP satellite observations and CMIP5 models[END_REF][START_REF] Kay | Evaluating and improving cloud phase in the Community Atmosphere Model version 5 using spaceborne lidar observations[END_REF].

b. Limitations of model-observation comparisons over the entire Greenland ice sheet

Before going forward with the rest of the results, we want to introduce and discuss two limitations. First, the (2008) only. This limitation occurs because the climate model data that are available through CFMIP-2 are only available for one year. The second limitation is more fundamental. There are no direct observations of the surface cloud radiative effects over the entire Greenland ice sheet. As a result, it is complicated to connect the lidar-observed and lidar-simulated cloud properties directly with quantities of most interest to the climate: cloud influence on radiation and surface energy budgets. We next address each of these concerns in turn.

Starting with the first limitation-we acknowledge that using only a single year for comparisons could be a very serious limitation. In fact, the seriousness of this limitation depends on if the model biases are larger than the interannual variability. To evaluate the seriousness of this single-year limitation, we compare our single-year evaluation with a multiyear evaluation. Specifically, we compare winter and summer cloud fraction data for eight years (2008-15 for observations, 2001-08 for models) with those that were done for one year (Fig. 3). The results of this comparison Figs. 7 and 8 reveal that the interannual variability of cloud fraction profiles over Greenland is less than 5% at all levels for all models except IPSL-CM5A-MR in winter. In the observations, the cloud fractions for the year 2008 are representative of the mean cloud fraction between 2008 and 2015; they are within 62% or better at all heights. Thus, the model biases are larger than the interannual variability, at least for the models that we could assess.

Moving to the second, more fundamental limitationwe are not able to compare lidar cloud properties with direct surface radiation observations over all of Greenland. To address this observational limitation, we assess the relationship between cloud influence on surface radiation and cloud properties at the one location where both are available: Summit.

To assess if the relationships between cloud properties and radiation at Summit are relevant for the rest of the Greenland ice sheet, we compare cloud property biases over all of Greenland with those over just Summit (Tables 3 and4). In winter, in many models, the sign of the cloud cover biases match over Greenland and Summit, but in two models (CanAM4 and HadGEM2-A) the sign of the biases differ. Compared to other models, the biases in CanAM4 and HadGEM2-A are small and can be due to regional differences in the two models. In summer, every model underestimates the cloud cover as well as the opaque and nonopaque cloud fraction. Only the cloud elevation biases are different; most of the time, the models underestimate the mean cloud altitude above Greenland and overestimate it at Summit. The consistency of most of the cloud property biases between the Greenland average and Summit fosters the use of observations at Summit only. In the next section, we look in more detail at how the simulated cloud properties and radiative effects compare to the observations in the Summit environment only.

c. Evaluation of surface CRE in models using observations from Summit

We begin with winter because the lack of shortwave radiation during winter makes it simpler to relate cloud properties to cloud radiative effects. Indeed, the total CRE is equivalent to longwave cloud radiative effect (CRE LW ) during winter. The longwave flux at the surface depends mostly on cloud and atmosphere properties with cloud cover being the main property affecting CRE LW . During winter, the monthly spaceborne-lidar cloud cover at Summit ranges between 40% and 80% (Fig. 9a). The corresponding monthly mean CRE LW measured at Summit ranges from 130 to 150 W m 22 (Fig. 9b).

During winter at Summit, there are strong connections between model biases in the cloud properties (cover, opacity, altitude) and model biases in CRE LW . We find models with the smallest cloud cover have the smallest CRE LW . For example, CAM5 is the model with the least cloud cover, and it has a correspondingly small value of CRE LW . Cloud cover bias and CRE LW bias are not always so simply related to each other. In fact, HadGEM2-A and CanAM4 have cloud cover between 50% and 80% but underestimate the CRE LW by ;210 W m 22 compared to the ground observations, suggesting cloud properties influence the CRE LW bias in these two models.

Having investigated winter, we next move to summer. Evaluating model CRE and connecting model cloud property biases to model CRE biases is more complex in summer because we must consider both shortwave and longwave radiation. The shortwave CRE is influenced by the surface albedo and thus is measuring the cloud cooling but relative to the surface albedo. At Summit, the observed cloud cover and observed CRE LW are both larger in summer than in winter. During summer, the cloud cover ranges between 70% and 90% (Fig. 10a) and CRE LW is 60 W m 22 (Fig. 10c).

During summer, clouds in the models do not radiatively warm the surface enough compared to the observations. In summer, this bias is larger than in winter for CanAM4 and HadGEM2-A because these two models do not reproduce the seasonal cycle of the cloud cover at Summit. In contrast, only CAM5 produces more cloud cover in summer than in winter. While CAM5 better reproduces the seasonal cycle, it still simulates too little summer cloud cover and thus underestimates the surface warming compared to observations. Cloud cover biases explain most of the CRE SW biases: The models with the smallest cloud cover, such as CAM5 (Fig. 10a), also produce the smallest shortwave cloud radiative cooling (Fig. 10d). While MRI-CGCM3 simulates just a bit more clouds than CAM5 in JJA (Fig. 10a), it overestimates the CRE SW at the surface compared to CAM5 and to ground observations. The two models produce distinct shortwave forcing for the same cloud cover, with average values of CRE SW, CAM5 5 23.7 W m 22 and CRE SW, MRI-CGCM3 5 221.5 W m 22 at Summit (Fig. 10b).

In some cases as described above, the cloud cover alone does not explain all the CRE SW biases. Thus, next we look at the other factors influencing the CRE SW . 2) The surface albedos in summer are very different between the two models. In MRI-CGCM3, the surface albedo reaches values as low as 0.73, whereas it stays above 0.8 in CAM5 (Fig. 11d). In the Summit ground observations, the average measured surface albedo is 0.86. A higher albedo leads to less shortwave radiative cooling due to the clouds, with all else constant. Thus, the very low albedo in MRI-CGCM3 helps to explain why its CRE SW is so large even though the model underestimates clouds. Hereafter, we looked at the contribution of these last two factors on the CRESW. Using the surface shortwave flux outputs, we investigate how the average CRE in MRI-CGCM3 would change if the surface albedo and cloud opacity were changed (Fig. 11b). If the surface albedo in MRI-CGCM3 is increased to be similar to CAM5 (0.8 in July instead of 0.73), then CRE SW, MRI-CGCM3 5 212 W m 22 in JJA indicating that the difference in cloud opacity between the models explains (23.7) 2 (212) 5 8.3 W m 22 of their difference in CRE SW . On the other hand, by theoretically decreasing the cloud opacity such that MRI-CGCM3 matches CAM5, the downward shortwave fluxes become 400 W m 22 in July instead of 350 W m 22 , and the remaining difference in CRE SW is 7.8 W m 22 . Thus, cloud opacity and surface albedo both contribute about half of the CRE SW difference between CAM5 and MRI-CGCM3.

Combining the CRE LW and CRE SW components reveals that the deficit of modeled longwave radiative warming compared to the ground observations is generally not compensated by the deficit of shortwave radiative cooling (Fig. 10b). Every model included in this study underestimates the net cloud radiative surface warming in summer. For the models with the best distribution of net CRE, this underestimation arises from compensating biases in longwave and shortwave FIG. 9. Cloud cover and surface longwave radiative effect over Summit in winter. (a) PDF of the cloud cover from models and from satellite observations and (b) PDF of the surface longwave CRE from models and from ground-based observations. PDF are built from 28 3 28 gridbox monthly mean values for a bin size of 10% for the cloud cover and 10 W m 22 for CRE surf LW . The number of sampled grid boxes (time and grid points) used for the calculation of each PDF is given in brackets in the legend. This number depends on the model's spatial resolution; the higher the resolution is, the higher is the number of sampled grid boxes. Cloud cover for the models comes from the model 1 COSP/lidar. components. In MRI-CGCM3, the clouds typically radiatively cool the surface as a result of its low surface albedo.

Discussion

Since we only have surface radiative flux measurements over Summit, the model surface radiation can only be evaluated over the Summit area (section 3c) and not over all Greenland. On the other hand, satellite observations provide an evaluation of the cloud cover, cloud fraction profile, and cloud opacity over all of Greenland; all these cloud variables influence the surface radiation.

Greenland cannot be seen as a homogeneous icecovered area, and Summit might not be representative of other regions. For instance, in the ablation zones at the margin of Greenland, the ice albedo is lower compared to Summit [START_REF] Goelles | Albedo reduction of ice caused by dust and black carbon accumulation: A model applied to the K-transect, west Greenland[END_REF]. Over such low-albedo surfaces, decreasing cloud cover enhances the surface melting [START_REF] Hofer | Decreasing cloud cover drives the recent mass loss on the Greenland ice sheet[END_REF]), when at the opposite, clouds' presence can enhance the ice sheet melt in the middle of Greenland [START_REF] Bennartz | July 2012 Greenland melt extent enhanced by low-level liquid clouds[END_REF]. These examples prove the necessity of having validated satellite observations above ice-covered surfaces. In the meantime, the Summit observations can be useful to evaluate models with the largest biases observed all over Greenland, including Summit.

The evaluation of the models showed that despite some constant biases, results may vary a lot from one model to another. Along with the clouds, the parameterization of the Greenland ice sheet drives the variability of the surface albedo and thus the CRE. Yet only few general circulation models are able to represent the surface of the Greenland ice sheet [START_REF] Cullather | Evaluation of the surface representation of the Greenland ice sheet in a general circulation model[END_REF]. In future experiments, dynamic ice sheet models will be integrated in climate models [START_REF] Eyring | Overview of the Coupled Model Intercomparison Project Phase 6 (CMIP6) experimental design and organization[END_REF]. This reflects that before the fifth assessment of the IPCC, Greenland received little attention from modeling centers. Further, cold clouds are hard to This number depends on the model's spatial resolution; the higher the resolution is, the higher is the number of sampled grid boxes. Cloud cover for the four CMIP5 models comes from the COSP/lidar outputs.

model, and we lack observational constraint in the arctic regions.

To help constrain the models' observational simulator are powerful tools. To improve the model evaluation of the cloud properties, the last CALIPSO-GOCCP opaque clouds diagnostics used in this study should be added in more climate models [START_REF] Guzman | Direct atmosphere opacity observations from CALIPSO provide new constraints on cloudradiation interactions[END_REF]. Hence, the evaluation of the cloud properties will be improved and help to understand the interactions between the clouds and the surface fluxes and thus help to predict future changes in the Greenland ice sheet mass balance.

Conclusions

In this study, we have evaluated climate model outputs compared to CALIPSO and ground-based observations in Greenland.

The comparison with CALIPSO observations revealed that the group of climate models examined here has some spread in their representation of vertical cloud fraction, cloud cover, and cloud opacity, with some constant biases in each model. Among the four models containing the height-intensity histogram outputs, they all underestimate the nonopaque clouds with large SR (SR . 30), and most models tend to underestimate opaque clouds observed in CALIPSO-GOCCP. These clouds are frequent in the observations, especially below 2 km AGL and are mostly liquid or mixed-phase clouds.

Biases in these cloud categories confirm the difficulty of the models to create liquid-containing clouds at temperatures below 08C [START_REF] Cesana | How well do climate models simulate cloud vertical structure? A comparison between CALIPSO-GOCCP satellite observations and CMIP5 models[END_REF]). On the other hand, some models create too many opaque clouds, especially in summer. This may be a compensating error that reduces biases in Earth's radiative budget [START_REF] Nam | The 'too few, too bright' tropical low-cloud problem in CMIP5 models[END_REF].

The comparison to Summit observations revealed that the models underestimate the CRE at the surface especially in summer. The CRE biases may be a consequence of the underestimated cloud cover in the models in summer (ranging from 246% to 221%). Fewer clouds causes an underestimation of the longwave radiative warming effect (235.7 to 213.6 W m 22 ) and an underestimation of shortwave cooling effect (11.5 to 110.5 W m 22 ) even though the model albedos are lower than that measured at Summit ('0.8 in the models and '0.9 in the observations). Overall, the simulated clouds do not radiatively warm the surface as much as in the observations.

The effect this bias has on the CRE over Greenland during summer depends on the surface albedo. CanAM4 produces the most opaque clouds, but it creates them over high-albedo surfaces, while HadGEM2-A creates opaque clouds above surfaces with a low albedo (,0.7). Thus, for the same amount of simulated opaque clouds, the shortwave radiative cooling effect is higher in HadGEM2-A than in CanAM4. Since the overall cloud radiative warming is underestimated in the models, we may expect an underestimation of Greenland surface melting. However, misrepresentation of clouds is not the only contributor to biases in the modeled surface melting. Previous works have shown that near-surface temperature is colder in the models than in the reanalysis independently of cloud presence [START_REF] Chapman | Simulations of Arctic temperature and pressure by global coupled models[END_REF]. Indeed, the surface area over Greenland with a near-surface temperature around zero is two times smaller in the models than in the reanalysis (not shown). The extent to which cloud biases impact surface temperature biases requires further study.

In this paper, we focused on the role of clouds on the surface melting, which is one side of the ice sheet mass balance; on the other side, the precipitations act as a source of mass. Precipitations are obviously connected to cloud presence so that satellites such as CloudSat and CALIPSO can contribute to the documentation of precipitating clouds. Currently, studies based on model prediction suggest that more precipitation will fall in a warmer climate [START_REF] Gregory | Ice-sheet contributions to future sea-level change[END_REF][START_REF] Fettweis | Brief communication ''Important role of the mid-tropospheric atmospheric circulation in the recent surface melt increase over the Greenland ice sheet[END_REF]. In future work, spaceborne radar combined with spaceborne lidar observations will enable a full assessment of the influence of clouds on surface melt and precipitation over Greenland.

  between January 2011 and October 2013. Both downward and upward fluxes are collected by instruments maintained by the Swiss Federal Institute of Technology [Eidgenössische Technische Hochschule (ETH); Zurich]. The downward and upward longwave fluxes come from two Kipp and Zonen CG4 pyrgeometers, sensitive to the spectral range 4.5-40 mm. The downward and upward shortwave fluxes come from two Kipp and Zonen CM22 pyranometers, sensitive to the spectral range 200-3600 nm. The uncertainties are 65 W m 22 for the longwave and shortwave fluxes

  ), four models [MRI-CGCM3, CAM5, MIROC5, and BCC_CSM1.1(m)] underestimate cloud cover, three models (CanAM4, HadGEM2-A, and MPI-ESM-LR) overestimate and underestimate cloud cover depending on the area, and the last model (IPSL-CM5A-MR) overestimates cloud cover. During summer (Figs. 2c-j) most of the models have insufficient cloud cover, but some [HadGEM2-A, MPI-ESM-LR, and BCC_CSM1.1(m)] have excessive cloud cover in the west and southwest of Greenland compared to the observations.

FIG. 1 .

 1 FIG. 1. Maps of mean cloud cover in winter 2008 for (a) observed, (b) multimodel mean bias, and (c)-(j) in each of the CMIP5 models. The observations are from CALIPSO-GOCCP. The multimodel mean is calculated from eight CMIP5 models with COSP/lidar: CanAM4, HadGEM2-A, MRI-CGCM3, CAM5, MIROC5, IPSL-CM5A-MR, MPI-ESM-LR, and BCC_CSM1.1(m). To calculate the multimodel mean, each model is interpolated into the CAM5 grid. The green contour lines indicate the 1000-, 2000-, and 3000-m surface elevation.

FIG. 3 .

 3 FIG. 3. Mean cloud fraction profiles over Greenland from observations and models in (a) winter 2008 and (b) summer 2008. The observations are from CALIPSO-GOCCP. The models means' are calculated from eight CMIP5 models with COSP/lidar: CanAM4, HadGEM2-A, MRI-CGCM3, CAM5, MIROC5, IPSL-CM5A-MR, MPI-ESM-LR, and BCC_CSM1.1(m). The shaded area is the envelope of the models' mean profiles.

FIG. 5 .

 5 FIG. 5. Profiles of opaque cloud fraction (gray), nonopaque cloud fraction (green), and clear sky fraction (blue) over Greenland in (a)-(e) DJF 2008 and (f)-(j) JJA 2008 observed from (a), (f) CALIPSO-GOCCP-OPAQ and simulated by four CMIP5 models 1 COSP/ lidar: CanAM4 in (b),(g); HadGEM2-A in (c),(h); MRI-CGCM3 in (d),(i); and CAM5 in (e),(j). The nonopaque cloud fraction (green) is split into two categories: the optically very thin clouds (light green) and the optically less thin clouds (dark green).

FIG. 7 .

 7 FIG. 7. Mean cloud fraction profile over Greenland calculated in DJF for (a) CALIPSO-GOCCP for the years 2008-15 and for (b)-(i) eight CMIP5 models for the years 2001-08.

Figure 11

 11 Figure 11 suggests two other factors explain the difference of CRE SW between CAM5 and MRI-CGCM3: 1) More downward shortwave fluxes reach the surface in CAM5 compared to MRI-CGCM3 (150 W m 22 in summer; Fig. 11c). It is likely because more opaque clouds are produced by MRI-CGCM3 compared to CAM5 (Fig. 4). 2) The surface albedos in summer are very different between the two models. In MRI-CGCM3,

FIG. 8 .

 8 FIG.8. As in Fig.7, but for JJA.

FIG. 10 .

 10 FIG. 10. PDFs of cloud cover and surface radiation over Summit in summer for (a) cloud cover, (b) net cloud radiative effect, (c) LW cloud radiative effect, and (d) SW cloud radiative effect. PDFs are built from monthly mean values in 28 3 28 grid boxes located in the Summit region (z . 3000 m) in JJA 2008. The bin size of the PDF are values 10% for the cloud cover, 10 W m 22 for CRE surf LW , and 10 W m 22 for CRE surf SW . Values in the brackets of the legend indicate the number of sampled grid boxes (time and grid points) used for the calculation of each PDF.This number depends on the model's spatial resolution; the higher the resolution is, the higher is the number of sampled grid boxes. Cloud cover for the four CMIP5 models comes from the COSP/lidar outputs.

  FIG. 11. Simulated annual cycle of the (a) cloud cover, (b) surface SW CRE, (c) all sky downward shortwave flux, and (d) surface albedo as simulated by two CMIP5 climate models (CAM5 and MRI-CGCM3) over the Summit region (z . 3000 m) in 2008. The cloud cover is obtained from models 1 COSP/lidar. Solid lines are standard models and models 1 COSP/lidar outputs, while the dotted lines in (b) are perturbed ones.

TABLE 1 .

 1 Cloud classification based on the lidar SR and groundecho detection[START_REF] Chepfer | Use of CALIPSO lidar observations to evaluate the cloudiness simulated by a climate model[END_REF][START_REF] Guzman | Direct atmosphere opacity observations from CALIPSO provide new constraints on cloudradiation interactions[END_REF]).

	Category	Fulfilled condition
	Clear sky	0.06 , SR , 1.2
	Nonopaque cloud	5 , SR , 80
	Opaque cloud	No ground detection or SR , 0.01 (COSP)

TABLE 2 .

 2 CMIP5 models included in the study.

			Available model outputs
			Surface all sky and clear	Cloud	Histogram	Temporal	Spatial
	Model name	Reference	sky radiative fluxes	cover	of SR	resolution	resolution
	CanAM4	Von Salzen et al. (2013)	Yes	Yes	Yes	Daily, monthly	64 3 128
	HadGEM2-A	Martin et al. (2011)				Monthly	145 3 192
	MRI-CGCM3	Yukimoto et al. (2012)					160 3 320
	CAM5	Hurrell et al. (2013)				Daily, monthly	192 3 288
	MIROC5	Watanabe et al. (2010)			No	Monthly	128 3 256
	IPSL-CM5A-MR	Dufresne et al.

Acknowledgments. The authors would like to acknowledge the support of NASA, CNES, and ClimServ for access to the CALIOP data. This work was supported in part by the CNES/EECLAT project. JEK acknowledges funding support from NASA Grant 15-CCST15-0025. NBM and MDS acknowledge funding support from the U.S. National Science Foundation Grants PLR1303879 and PLR1314156. The authors thank the three anonymous reviewers who greatly helped to improve the manuscript.


Henri Casanova

Julien Herrmann

Yves Robert

Computing the expected makespan of task graphs in the presence of silent errors

Applications structured as Directed Acyclic Graphs (DAGs) of tasks occur in many domains, including popular scientific workflows. DAG scheduling has received an enormous amount of attention. Many of the popular DAG scheduling heuristics make scheduling decisions based on path lengths. At large scale compute platforms are subject to various types of failures with no longer negligible probabilities of occurrence. Failures that have recently received increasing attention are "silent errors," which cause data corruption. Tolerating silent errors is done by checking the validity of compute results and re-executing a possibly task from scratch. The execution time of a task then becomes a random variable, and so do path lengths in a DAG. Unfortunately, computing the expected makespan of a DAG (and equivalently computing expected path lengths in a DAG) is a computationally difficult problem. Consequently, designing effective scheduling heuristics in this context is challenging. In this work we propose an algorithm that computes a first order approximation of the expected makespan of a DAG when tasks are subject to silent errors. We find that our proposed approximation outperforms previously proposed approaches both in terms of approximation error and of speed.

Introduction

This paper introduces a new algorithm to approximate the expected makespan of a workflow application, i.e., an application structured as a Directed Acyclic Graph (DAG) of tasks, in which tasks can fail and must be re-executed. A key question when executing workflows on parallel platforms is the scheduling of tasks on the available compute resources, or processors. When not considering task failures, list scheduling algorithms are the de-facto standard [START_REF] Kwok | Benchmarking and comparison of the task graph scheduling algorithms[END_REF], and tools are available that use such algorithms for scheduling workflow applications onto large-scale platforms in practice [START_REF] Wieczorek | Scheduling of scientific workflows in the askalon grid environment[END_REF][START_REF] Deelman | Workflows and escience: An overview of workflow system features and capabilities[END_REF]. The prevalent list scheduling heuristics prioritize tasks with large bottom levels. The bottom-level of a task is defined as the longest path from that task to the end of the execution, assuming unlimited processors. This heuristic is known as CP-scheduling (Critical Path scheduling [START_REF] Brucker | Scheduling Algorithms[END_REF][START_REF]Introduction to Scheduling[END_REF][START_REF] Pinedo | Scheduling: Theory, Algorithms, and Systems[END_REF]), and it has been extended to handle heterogeneous environments (see the HEFT algorithm [START_REF] Topcuoglu | Performance-effective and lowcomplexity task scheduling for heterogeneous computing[END_REF]).

It is well-documented that large-scale platforms are increasingly subject to errors that can cause task failures. In particular, the occurrence of "silent errors" (or SDCs, for Silent Data Corruptions) has been recently identified as one of the major challenges for Exascale [START_REF] Cappello | Toward exascale resilience: 2014 update[END_REF]. Silent errors can be caused by external causes, including cosmic radiations or packaging pollution. In addition, silent errors can occur when using DVFS (Dynamic Voltage Frequency Scaling) to reduce energy consumption. For instance, at low voltage the probability that a task produces incorrect results increases [START_REF] Dixit | The impact of new technology on soft error rates[END_REF][START_REF] Zhu | The effects of energy management on reliability in real-time embedded systems[END_REF].

Regardless of their causes, an effective approach for avoiding propagating results corrupted by silent errors is to use a verification mechanism after executing each task [START_REF] Hérault | Fault-Tolerance Techniques for High-Performance Computing[END_REF]. When an error is detected, the task is then re-executed. The verification mechanism can be general-purpose (e.g., based on replication [START_REF] Fiala | Detection and correction of silent data corruption for largescale high-performance computing[END_REF], with re-execution only when the two outputs do not match) or applicationspecific. Many application-specific error detection methods are available for classical High performance Computing (HPC) applications, such as approximate re-execution for ODE and PDE solvers [START_REF] Benson | Silent error detection in numerical time-stepping schemes[END_REF], orthogonality checks for Krylovbased sparse solvers [START_REF] Chen | Online-ABFT: An Online Algorithm Based Fault Tolerance Scheme for Soft Error Detection in Iterative Methods[END_REF][START_REF] Sao | Self-stabilizing iterative solvers[END_REF]), or Algorithm-Based Fault Tolerance (ABFT) for numerical linear algebra [START_REF] Bosilca | Algorithm-based fault tolerance applied to high performance computing[END_REF].

Silent errors make it challenging to define efficient list scheduling algorithms to schedule task graphs. After the first execution of a task, the error detector is used to check the result. If the result is correct, the task's execution is marked as successful, and its successor tasks are marked as ready. If the result is incorrect, then the task must be executed (and, given that the first execution has failed, the second execution will typically succeed with very high probability). While this scheme is conceptually simple, it greatly complicates the computation of the bottom-level of a task. And yet, computing the expected length of the longest path in a DAG with unlimited processors (or equivalently the expected bottomlevel of a task in the DAG), is key to designing silent-error-aware versions of effective list scheduling heuristics (CP-scheduling, HEFT).

It is known that computing the expected length of the longest path in a DAG whose task weights are probabilistic is a difficult problem [START_REF] Pinedo | Scheduling: Theory, Algorithms, and Systems[END_REF]. Even in the case in which each task is re-executed at most once, i.e., when task weights are random variables taking only two discrete values, the problem remains #Pcomplete [START_REF] Hagstrom | Computational complexity of pert problems[END_REF] (see Section 2 for a detailed discussion). In this work we develop an algorithm to compute an accurate first-order approximation of the expected length of the longest path in a general DAG in which tasks are subject to silent errors, and whose execution lengths can take two different values, depending upon there is a re-execution or not. More specifically, our contributions are:

• We develop an exact first-order approximation of the expected makespan of a general DAG, which can be computed in polynomial time;

• We compare our approach to two previously proposed approximations for DAGs for three sets of task graphs;

• We quantify approximation errors via comparison to a brute-force Monte Carlo approach; and

• We show that our proposed approximation (i) leads to lower or similar error than previously proposed approximations (and, importantly, to much lower error when the probability of task failure is low) and (ii) can be computed in practice much faster than previously proposed approximations.

This paper is organized as follows. Section 2 reviews related work. Section 3 formalizes the problem that we address and states assumptions. Section 4 describes our proposed approximation. Section 5 presents evaluation results. Section 6 concludes with a brief summary of results and with perspectives on future work.

Related Work

In this section we first review previous works on estimating the makespan of probabilistic DAGs, and then review relevant literature on silent errors.

Expected makespan of probabilistic 2-state DAGs

Computing the expected makespan of a DAG whose task execution times obey arbitrary probability distributions is known to be a difficult problem, even with unlimited processors. When task weights are fixed (deterministic), the makespan is the length of the longest path in the graph (also called "critical path"). But assume instead a probabilistic 2-state DAG where task weights obey a simple 2-state probability distribution: task T i has weight a i,1 with probability p i and weight a i,2 with probability 1 -p i . Assume also that all the probability distributions of all tasks are independent. The makespan of the DAG is now a random variable. It is known that computing its probability distribution, or even just its expected value, is a #P-complete problem [START_REF] Hagstrom | Computational complexity of pert problems[END_REF]. Recall that the class of #P problems is the class of counting problems corresponding to NP decision problems [START_REF] Valiant | The complexity of enumeration and reliability problems[END_REF][START_REF] Provan | The complexity of counting cuts and of computing the probability that a graph is connected[END_REF][START_REF] Bodlaender | A note on the complexity of network reliability problems[END_REF], and that #P-complete problems are at least as hard as NP-complete problems.

An informal explanation of why computing the expected makespan of probabilistic 2-state DAGs is a difficult combinatorial problem, is as follows. The main intuitive reason is that the expected value of the maximum of two random variables is not the maximum of their expectations. As a result, when computing the length of a path in the DAG, one must keep track of all possible values for the starting time of each task, together with their probabilities, and there may be an exponential number of such values [START_REF] Möhring | Scheduling under uncertainty: Bounding the makespan distribution[END_REF].

In practice, there are three standard methods to compute the expected makespan of a probabilistic 2-state DAG, as described hereafter.

Monte Carlo simulations

The Monte Carlo approach works as follows [START_REF] Mitzenmacher | Probability and Computing: Randomized Algorithms and Probabilistic Analysis[END_REF][START_REF] Van Slyke | Monte carlo methods and the pert problem[END_REF]: For each task in the DAG, a value of its weight is sampled from its probability distribution. Once this is done, the DAG is deterministic and its longest path can be computed as explained in Section 3. One then repeats this operation for a large number of trials, generating a new value at each trial. The set of these values empirically approaches the actual distribution of the DAG makespan as the number of trials increases.

An interesting question is that of determining the number of trials to obtain a high confidence level in the result. We refer the reader to the relevant discussion in [START_REF] Canon | Correlation-aware heuristics for evaluating the distribution of the longest path length of a DAG with random weights[END_REF]. A key drawback of the Monte Carlo approach is that it is computeintensive since the necessary number of trials is typically high. In this work, we only use Monte Carlo to compute a ground truth so as to assess the accuracy of our and previously proposed algorithms that compute approximations of the expected makespan. Hence, instead of determining a minimal number of trials, we conservatively use a very large number of trials so as to guarantee that our ground truth is accurate.

Approximation by a series-parallel graph

Basic probability theory tells us how to compute the probability distribution of the sum of two random variables (by a convolution) and of the maximum of two random variables (by taking the product of their cumulative density functions). This simple consideration leads to an exact method to compute the expected makespan when the DAG is series-parallel (see [START_REF] Valdes | The recognition of series parallel digraphs[END_REF] for a definition of series-parallel graphs). The problem with probabilistic 2-state series-parallel graph remains NP-complete in the weak sense and admits a pseudo-polynomial solution [START_REF] Möhring | Scheduling under uncertainty: Bounding the makespan distribution[END_REF].

When the DAG is not series-parallel, one approach is to approximate it by a series-parallel graph, which is constructed iteratively, first by a sequence of reductions and then by duplicating some vertices. Dodin's method [START_REF] Dodin | Bounding the project completion time distribution in PERT networks[END_REF] constructs such an approximated series-parallel graph, whose expected makespan is used to estimate that of the original DAG. See [START_REF] Möhring | Scheduling under uncertainty: Bounding the makespan distribution[END_REF][START_REF] Canon | Correlation-aware heuristics for evaluating the distribution of the longest path length of a DAG with random weights[END_REF] for a detailed description of Dodin's method. We include this method in our quantitative experiments in Section 5.

Approximation with normality assumption

The central-limit theorem states that the sum of independent random variables tends to be normally distributed as the number of variables increases. The expected makespan of the DAG is a combination of sums and maximums of the original task weights, so a popular approach proposed by Sculli [START_REF] Sculli | The completion time of PERT networks[END_REF] is based on the normality assumption:

• Approximate the distribution of each task by a normal distribution of same mean and variance. This step has constant cost per task for probabilistic 2-state DAGs.

• Use Clarke's formulas in [START_REF] Clark | The greatest of a finite set of random variables[END_REF] to compute the mean and variance of the sum and maximum of two (correlated) normal distributions, and then assuming that they also follow normal distributions.

• Traverse the original DAG and compute the mean and variance of the makespan.

See [START_REF] Canon | Correlation-aware heuristics for evaluating the distribution of the longest path length of a DAG with random weights[END_REF] for a full description of Sculli's method, which we we also include in our quantitative experiments in Section 5.

Silent errors

Considerable efforts have been directed at verification techniques to handle silent errors. A guaranteed, general-purpose verification is only achievable with expensive techniques, such as process replication [START_REF] Fiala | Detection and correction of silent data corruption for largescale high-performance computing[END_REF][START_REF] Ni | ACR: Automatic Checkpoint/Restart for Soft and Hard Error Protection[END_REF] or redundancy [START_REF] Lyons | The use of triple-modular redundancy to improve computer reliability[END_REF][START_REF] Elliott | Combining partial redundancy and checkpointing for HPC[END_REF]. However, application-specific information can be exploited to decrease the verification cost. Algorithm-based fault tolerance (ABFT) [START_REF] Huang | Algorithm-based fault tolerance for matrix operations[END_REF][START_REF] Bosilca | Algorithm-based fault tolerance applied to high performance computing[END_REF][START_REF] Shantharam | Fault tolerant preconditioned conjugate gradient for sparse linear system solution[END_REF] is a wellknown technique to detect errors in linear algebra kernels using checksums. Various techniques have been proposed in other application domains. Benson et al. [START_REF] Benson | Silent error detection in numerical time-stepping schemes[END_REF] compare a higher-order scheme with a lower-order one to detect errors in the numerical analysis of ODEs. Sao and Vuduc [START_REF] Sao | Self-stabilizing iterative solvers[END_REF] investigate self-stabilizing corrections after error detection in the conjugate gradient method. Heroux and Hoemmen [START_REF] Heroux | Fault-tolerant iterative methods via selective reliability[END_REF] design a fault-tolerant GMRES capable of converging despite silent errors. Bronevetsky and de Supinski [START_REF] Bronevetsky | Soft error vulnerability of iterative linear algebra methods[END_REF] provide a comparative study of detection costs for iterative methods.

Recently, detectors based on data analytics have been proposed to serve as partial verifications [START_REF] Berrocal | Lightweight silent data corruption detection based on runtime data analysis for HPC applications[END_REF][START_REF] Gomez | Detecting silent data corruption through data dynamic monitoring for scientific applications[END_REF][START_REF]Detecting and Correcting Data Corruption in Stencil Applications through Multivariate Interpolation[END_REF]. These detectors use interpolation techniques, such as time series prediction and spatial multivariate interpolation, on scientific dataset to offer large error coverage at the expense of a negligible overhead. Although not perfect, the accuracy-to-cost ratios of these techniques tend to be very high, which makes them attractive alternatives at large scale.

The rate of silent errors also depends upon the execution framework. Temperature and power consumption are known to have a high impact [? ]. Furthermore, as mentioned in Section 1, lowering the voltage/frequency is believed to have an adverse effect on system reliability [START_REF] Dixit | The impact of new technology on soft error rates[END_REF][START_REF] Zhu | The effects of energy management on reliability in real-time embedded systems[END_REF]. In particular, many papers (e.g., [START_REF] Zhu | The effects of energy management on reliability in real-time embedded systems[END_REF][START_REF] Zhao | Reliability-aware dynamic voltage scaling for energy-constrained real-time embedded systems[END_REF][START_REF] Aupy | Energy-aware scheduling under reliability and makespan constraints[END_REF][START_REF] Das | Combined DVFS and mapping exploration for lifetime and soft-error susceptibility improvement in MPSoCs[END_REF]) have assumed the following exponential error rate model:

λ(s) = λ 0 • 10 d(smax -s) smax-s min , (1)
where λ 0 denotes the average error rate at the maximum speed s max , d > 0 is a constant indicating the sensitivity of error rate to voltage/frequency scaling, and s min is the minimum speed. This model suggests that the error rate increases exponentially with decreased processing speed, which is a result of decreasing the voltage/frequency and hence lowering the circuit's critical charge (i.e., the minimum charge required to cause an error in the circuit). This is a concern for resilience, because these studies suggest that minimizing energy via DVFS techniques will also lead to an increased number of silent errors.

A silent error must be detected for each task via some verification mechanism [START_REF] Hérault | Fault-Tolerance Techniques for High-Performance Computing[END_REF]. A naive but general-purpose verification technique would consist in re-running the task once so that an output mismatch indicates a silent error [START_REF] Fiala | Detection and correction of silent data corruption for largescale high-performance computing[END_REF]. In this case, the task can be re-executed one extra time so that the correct output can be picked via majority voting. Many application-specific, and thus efficient, error detection methods have been proposed [START_REF] Benson | Silent error detection in numerical time-stepping schemes[END_REF][START_REF] Chen | Online-ABFT: An Online Algorithm Based Fault Tolerance Scheme for Soft Error Detection in Iterative Methods[END_REF][START_REF] Sao | Self-stabilizing iterative solvers[END_REF][START_REF] Bosilca | Algorithm-based fault tolerance applied to high performance computing[END_REF]. Because silent errors are not detected when their occur, but only at the end of the execution of the task (by the verification mechanism), the task must be fully re-executed in case an error is detected. Assuming that the task re-execution allows a correct output to be produced, we then need to consider probabilistic 2-state DAGs, in which each task has a weight for the no-error case, and another (longer) weight for the error case. For simplicity, in this work we assume that the weight of a task T is a if there is no silent error, and 2a if a silent error is detected (which requires a re-execution). However, our work is immediately extensible to arbitrary such two weights to account, for instance, for the overhead of the silent error detection mechanism.

Finally, this work assume that task failure arrival times are Exponentially distributed with Mean Time Between Failure (MTBF) 1/λ. Extensions to other probability distributions such as those considered in [? ? ] apply whenever using the distribution MTBF as the (average) inverse of the failure rate and neglecting second order terms (i.e., terms proportional to the square of that failure rate) would lead to accurate approximations of the corresponding failure probabilities.

Problem statement

We consider a general model of computation in which an application is structured as a Directed Acyclic Graph, in which vertices represent tasks and edges represent task precedence. More formally, let G = (V, E) be a DAG, with V a set of tasks, and E ⊂ V × V a set of edges. For each task i, let a i be its weight, i.e., its failure-free execution time. Let pa(i, j) denote the length (as a sum of task weights) of the longest path from task i to task j, if such a path exists, otherwise let pa(i, j) = -∞. The longest path length in G is then defined as We consider that tasks fail independently and task failure arrival times are exponentially distributed with Mean Time Between Failure (MTBF) 1/λ. Therefore, the probability that task i fails during its first execution attempt is 1 -e -λai , in which case the task must be re-executed from scratch.

d(G) = max i,j∈V {pa(i, j)}.
Our objective is to compute an approximation of the expected makespan of G, i.e., the longest path length in G taking into accounts that tasks can fail and must be re-executed. The failure-free makespan defined above is a clear lower bound on the expected makespan.

Approximating the Expected Makespan

In this section we compute a first-order approximation of the expected makespan of a DAG G, which we denote as E(G). Our approximation relies on the fact that in practice λ is close to zero, which allows us to neglect O(λ 2 ) terms.

The probability that the first execution attempt of task i succeeds is:

p i = e -λai = 1 -λa i + O(λ 2 ) .
The probability that the first execution attempt of task i fails but that its second execution attempt succeeds is :

(1 -e -λai )e -λai = λa i + O(λ 2 ) = 1 -p i + O(λ 2 ) .
Neglecting the O(λ 2 ) terms leads to the approximation that a task either takes time a i , with probability 1 -λa i , or time 2a i , with probability λa i . In other terms, our first-order approximation consists in assuming that a task never fails more than once. Hereafter when we say that a task fails, we mean that its first execution attempt fails and that its second execution attempt succeeds.

For any S ⊂ V , let P (S) denote the probability that all tasks in S fail and that no task in V \ S fails. Let also L(S) denote the length of the longest path in G when all tasks in S fail and no task in V \ S fails. E(G) is thus defined as:

E(G) = S⊂V P (S) × L(S) .
We note that:

P (∅) = i∈V (1 -λa i + O(λ 2 )) = 1 - i∈V λa i + O(λ 2 ) , P ({i}) = (λa i + O(λ 2 )) × j∈V \{i} (1 -λa j + O(λ 2 )) = λa i + O(λ 2 ) , and
P (S) = O(λ 2 ) if |S| > 1 .
Therefore:

E(G) = P (∅) × L(∅) + i∈V P ({i}) × L({i}) + O(λ 2 ) .
By definition, L(∅) = d(G) (G's failure-free makespan). Similarly, L({i}) = d(G i ), where G i a DAG identical to G but such that task i has weight 2a i instead of weight a i . We thus obtain:

E(G) = (1 - i∈V λa i ) × d(G) + i∈V (λa i ) * d(G i ) + O(λ 2 ) = d(G) + λ i∈V a i (d(G i ) -d(G)) + O(λ 2 ) .

Evaluation

Makespan Approximation Techniques

In this section we evaluate three expected makespan approximations techniques:

1. First Order -The approximation described in Section 4; 2. Dodin -The bound in [START_REF] Dodin | Bounding the project completion time distribution in PERT networks[END_REF], which is computed by transforming any general DAG into an approximately equivalent series-parallel graph and then computing the exact expected makespan of this graph using the approach explained in Section 2.1.2. 3. Normal -The approach that consists in approximating the discrete execution time of each task (i.e., a i with probability p i and 2a i with probability 1 -p i ), by a continuous Normal distribution of same mean and standard deviation. The overall expected makespan in the approximated as explained in Section 2.1.3.

Application DAGs

We measure the error of the three makespan approximation techniques using 3 DAG datasets.

Linear Algebra

This dataset consists of DAGs used in numerical linear algebra computations. More specifically, we consider 3 classical factorizations of a k × k tiled matrix: Cholesky, LU, and QR factorization. Each tile has size b × b, where b is a platform-dependent parameter. Hence the actual size of a k × k tiled matrix is N × N , where N = kb. For each factorization, the number of vertices in the DAG depends on k as follows: the Cholesky DAG has 1 3 k 3 + O(k 2 ) tasks, while the LU and QR DAGs have 2 3 k 3 + O(k 2 ) tasks (but the tasks in QR entail, on average, twice as many floating-point operations as in LU).

Figures 1,2, and 3 show examples for k = 5. The tasks in these DAGs are labeled by the corresponding BLAS kernels [START_REF] Choi | The design and implementation of the ScaLAPACK LU, QR, and Cholesky factorization routines[END_REF], and their weights are based on actual kernel execution times as reported in [START_REF] Augonnet | StarPU: A unified platform for task scheduling on heterogeneous multicore architectures[END_REF] for an execution on Nvidia Tesla M2070 GPUs with tiles of size b = 960. For simplicity, in all that follows we call k the DAG size (i.e., the larger k the more tasks). For each DAG class we perform experiments with k = 4, 6, 8, 10, 12, for a total of 3 × 5 = 15 DAGs with up to 650 tasks.

Assembly Trees

This dataset consists of assembly trees for a set of sparse matrices obtained from the University of Florida Sparse Matrix Collection (http://www.cise. ufl.edu/research/sparse/matrices/). These matrices satisfy the following assertions: not binary, not corresponding to a graph, square, with a symmetric pattern, with between 20,000 and 2,000,000 rows, with an average number of non-zero elements per row at least equal to 2.5, and with a total number of non-zero elements at most equal to 5,000,000. Among these matrices, for each group of matrix in the dataset we pick the one that has the largest number of non-zero elements, resulting in 76 matrices (given the content of the dataset at the time this article is being written). We first order the matrices using MeTiS [START_REF] Karypis | MeTiS: A Software Package for Partitioning Unstructured Graphs, Partitioning Meshes, and Computing Fill-Reducing Orderings of Sparse Matrices Version 4[END_REF] (through the MeshPart toolbox [START_REF] Gilbert | Geometric Mesh Partitioning: Implementation and Experiments[END_REF]) and amd (available in Matlab), and then build the corresponding elimination trees using the symbfact routine of Matlab. We also perform a relaxed node amalgamation on these elimination trees to create assembly trees. We thus obtain a large set of instances by allowing 1, 2, 4, and 16 (if more than 1.6 × 105 nodes) relaxed amalgamations per node. Among these assembly trees, we selected four of them for our experiments: 'dump.1.1.amd.Goodwin.rim-447' (a.k.a. 'Goodwin'), 'dump.1.1.amd.Li.li-732' (a.k.a. 'amd.Li.li'), 'dump.1.1.amd.FEMLAB.sme3Dc-932' (a.k.a. 'FEMLAB') and 'dump.1.3.metis.Li.li-732' (a.k.a. 'metis.Li.li').

Random task graphs

This dataset consists random DAGs generated using the Directed Acyclic Graph GENerator (DAGGEN) 1 . DAGGEN uses four popular parameters to define the shape of a DAG: size, width, density and jumps.

• The size determines the number of tasks in the DAG (tasks are organized in levels). • The width determines the maximum parallelism in the DAG, that is the number of tasks in the largest level. A small value leads to "chain" graphs and a large value to "fork-join" graphs. • The density denotes the number of edges between two levels of the DAG, with a low value leading to few edges and a large value to many edges. • Random jump edges are added that go from level l to levels l + 1 . . . l + jumps. The second and third parameters take values between 0 and 1. Our DAG generation procedure is similar to the one used in [START_REF] Suter | Scheduling delta-critical tasks in mixed-parallel applications on a national grid[END_REF]. Specifically, we generate three sets of random DAGs:

• Daggen1: 100 randomly generated DAGs with size = 1000, width = 0.3, density = 0.5 and jumps = 5. A sample Daggen1 DAG is depicted in Figure 4.

• Daggen2: 100 randomly generated DAGs with size = 1000, width = 0.3, density = 0.9 and jumps = 7.

• Daggen3: 100 randomly generated DAGs with size = 1000, width = 0.9, density = 0.5 and jumps = 10.

Experimental Methodology

For each DAG, and for a given failure rate λ (see Section 3), we compute the First Order, Dodin, and Normal approximations of the expected makespan.

To compute approximation errors one would ideally use the exact expected makespan (the computation of which is a P#-complete problem). Instead, we resort to the brute-force Monte Carlo approach described in Section 2.1.1. We use 300,000 random trials. For each task in a trial, the task succeeds or fails as determined by sampling a random time-to-next-failure value from an Exponential distribution with parameter λ. We then approximate the expected makespan as the average makespan over the 300,000 samples. This method is prohibitively expensive in practice, but provides us with a reasonable ground truth in our experiments. In all results hereafter we report on the relative error between the approximations and this ground truth.

To allow for consistent comparisons of results across different DAGs (with different numbers of tasks and different task weights), in our experiments we simply fix the probability that a task of average weight fails, which we denote as p fail , and compute the failure rate. Formally, for a given DAG G = (V, E) and a given p fail value, we compute the average task weight as ā = i∈V a i /|V | and pick the failure rate λ such that

p fail = 1 -e -λā .
We evaluate the performance of our proposed approximation and compare it to its competitors for a range of p fail values. We present results for a wide range of values, noting that the higher values correspond to error rates much higher than those observed or expected on current and future large-scale computing platforms. The average execution time of a task in our experiments is ā = 0.15 seconds. So, for instance, p fail = 0.01 leads to an error rate λ = 0.067. The corresponding MTBF is µ = 1/λ = 14.9 seconds. For a platform with 100, 000 processors, this corresponds to an individual MTBF (per processor) of 17.3 days, quite an unrealistic value since individual processor MTBFs are typically estimated to be several years. Table 1 shows the p fail values used in our experiments and corresponding platform MTBF values (assuming 100, 000 processors). The results in the next section show that the lower p fail , the lower the error incurred by our proposed approximation. In particular, our approximation outperforms its competitors for all realistic p fail values (i.e., p fail ≤ 0.001), and sometimes even for unrealistically higher values.

Approximation Error Results

In this section we show relative error results, relative to expected makespans computed using the Monte Carlo method. Negative values denote an underestimation, while positive values denote an overestimation. All figures use a logarithmic scale on the vertical axis.

Linear algebra DAGs

Figures 5,6, and 7 show relative error vs. graph size for Cholesky graphs and for p fail = 0.01, 0.001, and 0.0001, respectively. The data set has 5 graphs representing the Cholesky factorization of a k × k tiled matrix for k = 4, 6, 8, 10, 12, as described in Section 5.2.1, which correspond to graphs with 20, 56, 120, 220 and 364 nodes. For p fail = 0.01, we see that First Order leads to the lowest relative error for graphs sizes below 10, and that Normal leads to the lowest relative error for larger graphs. Dodin leads to the largest relative error, but for the smallest graph size. Considering the largest graph size, First Order has 1.9% relative error while Normal has 0.3% relative error and Dodin has 8.7% relative error. For p fail = 0.001, First Order leads to dramatically lower error than its competitors (at least one order of magnitude lower). For instance, for the largest graph, First Order has 0.03% relative error, compared to 5.8% for Dodin and 0.9% for Normal. These results are even more striking for p fail = 0.0001. In this case, still for the largest graph, First Order has -0.0006% relative error, compared to 3.7% for Dodin and 0.4% for Normal. In this case, First Order is an underestimation of the expected makespan.

Figures 8,9, and 10 show similar results for LU DAGs. The data set has 5 graphs representing the QR factorization of a k × k tiled matrix for k = 4, 6, 8, 10, 12, as described in Section 5.2.1, which correspond to graphs with 30, 91, 204, 385 and 650 nodes. The overall message from these results is the same, meaning that First Order leads to drastically lower error than its competitors as p fail decreases. Here again Dodin leads to the largest errors overall. For p fail = 0.01, out approach leads to similar errors of the same order of magnitude as the errors for Normal (and in this case Normal is most often For p fail = 0.01, First Order leads to similar or lower error than Normal. Here again, Dodin leads to the highest errors across the board.

Assembly tree DAGs

Figures 14 to 18 show results for our Assembly tree DAGs described in Section 5.2.2. These results confirm the general observations from the previous set of results. Here again, Normal always outperforms Dodin. For the low values of p fail First Order leads to an overestimation of the makespan, while its competitors lead to underestimations, but the magnitude of the estimation errors are comparable. As soon as p fail ≤ 0.01, First Order leads to better results that its competitors, and for p fail ≤ 0.0001 this improvement is by orders of magnitude (but for the metis.Li.li trees for which its margin of improvement over its competitors is smaller). Overall, we conclude that First Order leads to similar or lower error than Normal (and Dodin), with drastic improvements for lower, more realistic, p fail values.

Randomly generated DAGs

Figure 19 to 23 show results for randomly generated DAGs described in Section 5.2.3. We observe similar trends as for the two previous sets of DAGs. Normal outperforms Dodin across the board. For low values of p fail , First Order leads to error comparable to that of Normal and Dodin. For p fail ≤ 0.0001 error, it leads to orders of magnitude improvements over its competitors. Our overall conclusion is that our proposed approximation approach improves over previously proposed approaches in terms of accuracy, and drastically so for realistically low p fail values. To assess the scalability of the three approximation methods, we run experiments with k = 20 for the LU DAG, i.e., 2, 870 tasks, and p fail = 0.0001. Results for other DAGs show similar relative ranking of the three approximation results, and we picked this particular DAG for our scalability analysis because it is from a real-world application, because it is large enough for a scalability study, and because it is small enough that the Monte Carlo method (which provides our ground truth) can run in a reasonable amount of time. For this large graph, we ran the Monte Carlo for ten hours, so as to ensure the accuracy of the ground truth. Error (normalized difference with Monte Carlo) and execution times are shown in Table 2. We see that for this large DAG, Dodin exhibits very large error. First Order is roughly two orders of magnitude more accurate than Normal. We also see that First Order can be computed in under a second, while Normal requires about 20 minutes, i.e., about three orders of magnitude longer. All these experiments are performed on one core of a 2.1GHz AMD Opteron(TM) Processor 6272.

Scalability Results

Summary of Results

Our results show that our proposed approximation, First Order, is accurate as long as the probability of task failure is sufficiently low. In particular, the realistically low probability of task failures it is more accurate that its two competitors by orders of magnitude. When the task failure probability is high, its accuracy is comparable to that of the Normal approximation. Across the board the Dodin approximation leads to the largest error. This is because the DAGs that we consider are far from being series-parallel. As a result, the seriesparallel graph constructed by Dodin is a poor approximation of the original DAG, hence the large errors. Finally, not only is First Order more accurate, but it is also faster. Its computation for a large DAG with 2, 870 tasks requires more than three orders of magnitude less time that Dodin and Normal.

Conclusion

We have proposed an algorithm to compute a first-order approximation of the expected makespan of a Directed Acyclic Graph (DAG) of tasks in which tasks are subject to silent errors that may require task re-execution from scratch due to corrupted results, which we term a failure. Our approximation can be computed in O(V 2 + V.E) time for a DAG with V vertices and E edges. It is exact in the first order and neglects second order O(λ 2 ) terms, where λ is the exponential failure rate. This amounts to assuming that each task may need to be re-executed at most once. The problem of computing the expected makespan of a DAG of tasks with this assumption is actually a computationally difficult problem (#P-complete). As a result, techniques to approximate the expected makespan have been proposed in previous works [START_REF] Dodin | Bounding the project completion time distribution in PERT networks[END_REF][START_REF] Canon | Correlation-aware heuristics for evaluating the distribution of the longest path length of a DAG with random weights[END_REF]. We have evaluated our proposed approximation and these previously proposed techniques for three sets of DAGs (linear algebra DAGs, assembly trees, and randomly generated DAGs). In our evaluations we quantify the approximation error via comparison to a ground truth computed using a brute-force Monte Carlo method.

Our results show that our proposed approximation is more accurate than previously proposed approximations by several orders of magnitude for realistically low failure rates. In addition, it can be computed much more quickly that these previously proposed approximations, which is crucially important for solving problems at scale. Overall, we have proposed a novel and improved approximation of the expected makespan of probabilistic DAGs in which tasks are subject to silent errors that mandate task re-executions.

A possible future research direction would to use our general approach to compute a more complicated, but still tractable, second order approximation (i.e., an approximation that does not neglect the O(λ 2 ) terms. While the improvement due to including the second order would be likely insignificant for low failure rates, it may be significant for relatively high failure rates, such as those observed for several recent systems [? ? ? ]. These higher failure rates, however, may not be relevant for the scale of current platforms and the failure probabilities of their processors. A broader and more promising future direction is to adapt existing list scheduling scheduling algorithms, or to develop novel such algorithms, that rely on our proposed approximation to make sound scheduling decisions.

  Let us call d(G) the failure-free makespan of G. Because G is acyclic, we can compute d(G) in O(|V | + |E|) time as follows: add two zero-weight vertices v 1 and v 2 to G, where v 1 represents a unique source task and v 2 a unique sink task. Also add an edge from v 1 to any entry task in G (a task without predecessor), and an edge from any exit task in G (a task without successor) to v 2 . Then d(G) = pa(s 1 , s 2 ) can be computed in time O(|V | + |E|) [42, Section 24.2].

  For a DAG G = (V, E), d(G) can be computed in O(|V | + |E|) time. Therefore, the above approximation can be computed in O(|V | 2 +|V |.|E|) time. Lower complexity can be achieved by exploiting the fact that G and the G i 's differ in only the weight of one task.

Figure 1 :

 1 Figure 1: DAG of a Cholesky factorization on a 5 × 5 tiled matrix

Figure 2 :Figure 3 :

 23 Figure 2: DAG of a LU factorization on a 5 × 5 tiled matrix

Figure 4 :

 4 Figure 4: Sample Daggen1 DAG

Figure 5 :

 5 Figure 5: Cholesky (k = 4, 6, 8, 10, 12), p fail = 0.01

Figure 6 :

 6 Figure 6: Cholesky (k = 4, 6, 8, 10, 12), p fail = 0.001

Figure 7 :

 7 Figure 7: Cholesky (k = 4, 6, 8, 10, 12), p fail = 0.0001

Figure 8 :

 8 Figure 8: LU (k = 4, 6, 8, 10, 12), p fail = 0.01

Figure 9 :

 9 Figure 9: LU (k = 4, 6, 8, 10, 12), p fail = 0.001

Figure 10 :

 10 Figure 10: LU (k = 4, 6, 8, 10, 12), p fail = 0.0001

Figure 11 :

 11 Figure 11: QR (k = 4, 6, 8, 10, 12), p fail = 0.01

Figure 12 :

 12 Figure 12: QR (k = 4, 6, 8, 10, 12), p fail = 0.001

Figure 13 :FEMLABFigure 14 :FEMLAB

 1314 Figure 13: QR (k = 4, 6, 8, 10, 12), p fail = 0.0001

Figure 15 : 1 - 10 - 2 -Figure 16 :FEMLAB

 15110216 Figure 15: Results on the 4 assembly trees, p fail = 0.001

Figure 17 : 10 - 1 - 10 - 2 -Figure 18 :

 1710110218 Figure 17: Results on the 4 assembly trees, p fail = 0.00001

Figure 19 :

 19 Figure 19: Results on the 3 sets of randomly generated graphs (of size 1000 and various shapes), p fail = 0.01

Figure 20 :Figure 21 :

 2021 Figure 20: Results on the 3 sets of randomly generated graphs (of size 1000 and various shapes), p fail = 0.001

Figure 22 :Figure 23 :

 2223 Figure 22: Results on the 3 sets of randomly generated graphs (of size 1000 and various shapes), p fail = 0.00001

Table 1 :

 1 p fail values and corresponding processor MTBF values assuming a 100, 000-processor platform.

	p fail	processor MTBF
	0.01	17.3 days
	0.001	173.5 days
	0.0001	4.7 years
	0.00001	47.5 years
	0.000001 475.6 years

Table 2 :

 2 Results for a LU DAG with k = 20 (i.e., 2,870 tasks) and p fail = 0.0001. Execution time ∼ 2 minutes ∼ 20 minutes < 1 second

		Dodin	Normal	First Order
	Error	-0.97	954 × 10 -6	7 × 10 -6

Code publicly available at https://github.com/frs69wq/daggen.

Acknowledgments. This work was supported in part by the European project SCoRPiO, by the LABEX MILYON (ANR-10-LABX-0070) of Université de Lyon, within the program "Investissements d'Avenir" (ANR-11-IDEX-0007) operated by the French National Research Agency (ANR), by the PIA ELCI project, and by NSF award SI2-SSE:1642369. The authors would like to thank the reviewers, whose comments and suggestions have helped improve the final version.