A Detailed Cloud Fraction Climatology of the Upper Indus Basin and Its Implications for Near-Surface Air Temperature*

Clouds play a key role in hydroclimatological variability by modulating the surface energy balance and air temperature. This study utilizes MODIS cloud cover data, with corroboration from global meteorological reanalysis (ERA-Interim) cloud estimates, to describe a cloud climatology for the upper Indus River basin. It has speciﬁc focus on tributary catchments in the northwest of the region, which contribute a large fraction of basin annual runoff, including 65% of ﬂow originating above Besham, Pakistan or 50km 3 yr 2 1 in absolute terms. In this region there is substantial cloud cover throughout the year, with spatial means of 50%–80% depending on the season. The annual cycles of catchment spatial mean daytime and nighttime cloud cover fraction are very similar. Thisregionaldiurnalhomogeneitybeliessubstantialspatialvariability,particularlyalongseasonallyvaryingverticalproﬁles(basedonsurfaceelevation). Correlations between local near-surface air temperature observations and MODIS cloud cover fraction conﬁrm the strong linkages between local atmospheric conditions and near-surface climate variability. These correlations are interpreted in terms of seasonal and diurnal variations in apparent cloud radiative effect and its inﬂuenceonnear-surfaceairtemperatureintheregion.Thepotentialroleofcloudradiativeeffectinrecognized seasonally and diurnally asymmetrical temperature trends over recent decades is also assessed by relating these locally observed trends to ERA-Interim-derived trends in cloud cover fraction. Speciﬁcally, reduction in nighttime cloudcoverfractionrelativetodaytime conditionsoverrecentdecades appears toprovide a plausible physical mechanism for the observed nighttime cooling of surface air temperature in summer months.


a. The upper Indus basin hydroclimatological and water resources context
The upper Indus River basin (UIB) covers a vast expanse of high-mountain Asia (164 867 km 2 ; Khan et al. 2014) and water from the Indus and its tributaries is the dominant supply source for Pakistan's irrigation, domestic consumption, and hydropower demands (Archer et al. 2010). Because of its importance for water supply, power generation, and, above all, food security, prospects for water resources in the UIB in the coming decades are a matter of great concern (Barnett et al. 2005), particularly in downstream areas (i.e., below the Tarbela reservoir). Figure 1 shows the main tributary catchments and locations of the meteorological stations used in this study.
Owing to its extremely mountainous nature, the dominant mechanisms generating runoff to feed the Indus are melting of glaciers and seasonal snowpack. These processes are constrained in their timing and progression by available energy inputs, which can be indexed as the fraction of catchment area where minimum temperatures remain above freezing (Forsythe et al. 2012). In the case of glacial regime tributary catchments, annual total runoff is governed by available energy, observed as mean summer temperature (Archer 2004). For this reason, climate processes including cloud cover, which influence local temperature, are of pragmatic interest for understanding present and future variability in water resources of the UIB.

b. Conceptual framework for radiative influence of UIB cloud fraction climatology
The precise and differentiated effects of cloud cover on the diurnal temperature cycle and surface have been and remain a key area of uncertainty for both climate scientists (Stephens 2005;Soden et al. 2008;Andrews et al. 2012;Gettelman et al. 2012;Betts et al. 2013Betts et al. , 2014 and hydrologists (Grundstein and Leathers 1998;Qian et al. 2007;Sherwood 2010). A detailed cloud climatology for the UIB is important for understanding both the present cloud radiative effects (CRE) on the local energy balance and for evaluation of the realism and coherence of the outputs of various numerical meteorological simulations.
The UIB, located between 348 and 378N, experiences considerable variation in incoming shortwave (SW) solar radiation over the course of the annual cycle (Fig. 2b) as attested by radiative flux climatologies from both ERA-Interim and local observations. These radiative flux data are described in the supplementary material to FIG. 1. Study area: the northwest upper Indus River basin (NW UIB). The focus study area is indicated by purple hashing. The inset map in the lower right corner shows the NW UIB (purple hashing) with respect to national boundaries in southern and central Asia.
this paper. The first physical mechanism of CRE is the reflection of this incoming solar radiation in the atmosphere, thus reducing the energy input reaching Earth's surface. The energy absorbed by the surface is then redistributed by a number of fluxes. First, depending on the relative temperature of the surface and the air immediately above it (Sade et al. 2011), energy can be redistributed in what is termed the sensible heat flux. This flux can operate in either direction (surface emitting or absorbing) and, in general, will follow regular patterns in the diurnal and annual cycles. Second, absorbed energy is also reemitted as (thermal) ''blackbody radiation'' with a rate dependent on the surface temperature according to the Stefan-Boltzmann law (Delany and Semmer 1998). This latter flux gives rise to the second physical mechanism for CRE as the outgoing thermal longwave (LW) energy may be limited by cloud cover analogously to a blanket or to the insulation of a building. Surface radiative fluxes and CRE mechanisms are illustrated conceptually in Fig. 2a.
Key elements radiative climatology of the UIB, shown in Fig. 2b (and presented in the supplementary material of this paper), may be summarized as follows: Net solar (SW) radiation, in the mountainous northwest UIB (NW UIB), reaches an annual maximum in early summer when 1) the ''top of atmosphere'' incoming SW radiation is also at a maximum; 2) premonsoonal prevailing ''mostly clear sky'' conditions allow much of the top-of-atmosphere radiation to reach the surface; and 3) in terms of surface area, the majority of ephemeral-seasonal snow cover has been ablated. The sign of net thermal (LW) radiation is opposite to net solar radiation as it quantifies the release of LW thermal energy from the surface to the atmosphere. The outgoing LW flux is dependent upon both surface temperature and overlying cloud and hence has a more complex annual cycle than net solar flux.
The interactions of the radiative and cloud climatologies of the UIB are expected to yield the following CRE influences on near-surface air temperature: d In the UIB, throughout the year, the magnitude of the daily net SW radiative flux is equal to or exceeds that of the oppositely signed net LW radiative flux with a minimum difference in the winter months and maximum difference in the summer.
d While both SW and LW radiative fluxes are modulated by instantaneous cloud cover, SW flux is limited to the daytime whereas LW flux is distributed through the diurnal cycle. The concentration of SW flux in the daytime exacerbates the magnitude difference with LW flux. Thus, the ''reflecting effect'' (reduced input energy) is expected to be the dominant daytime CRE mechanism, with direct influence on daily maximum temperature (Tmax).
d The absence of SW flux at night, however, means that the LW thermal insulating CRE mechanism is the controlling feature at night, with direct influence on daily minimum temperature (Tmin).
d The seasonal asymmetry of the radiative flux components implies that the daytime shortwave cooling mechanism is expected to dominate in summer months while the LW insulating effect is expected to be more important in the winter. ice cover, high surface elevations) with the UIB. A detailed study of cloud climatology and CRE at a low elevation mountain site in Germany by Ebell et al. (2011) found consistent negative (cooling) SW CRE, positive (warming) LW CRE, and net negative (cooling) CRE in a field observation campaign covering a period from April to December. Studying CRE in the Arctic, Shupe and Intrieri (2004) found that LW CRE had a roughly linear relationship with cloud cover fraction while shortwave CRE was influenced by total insolation (incoming SW) and surface albedo in addition to the cloud cover fraction during the daytime. More recently, Betts et al. (2013Betts et al. ( , 2014 have taken advantage of dataset of local cloud cover observations over the Canadian prairies with hourly temporal resolution to perform a comprehensive regional assessment of cloud cover fraction influence on the diurnal cycles of air temperature and relative humidity as well as SW and LW radiative fluxes. The Canadian prairies with substantial winter snow cover and large-amplitude annual temperature cycles provide an interesting analog to the UIB, albeit with the caveat of dramatically different topographic contexts. Key findings of Betts et al. (2013) with respect to CRE include the following: 1) in the warm season SW CRE dominates with negative correlations to Tmax and diurnal temperature range (DTR) but no substantial influence on Tmin; 2) in the cold season LW CRE dominates with both positive correlations to Tmax and Tmin; and 3) in the warm season, cloud cover-represented as effective cloud albedo-''is the primary control on the net surface radiative fluxes, which in turn drive the diurnal cycle of temperature and humidity'' (p. 8950). Betts et al. (2014) extended this work to investigate interactions between cloud cover and surface snow cover in governing the surface radiation balance finding that in transitional seasons ''snow cover acts as a fast climate switch'' (p. 1119), changing CRE mechanisms from warm season SW dominance to cold season LW dominance.

c. Rationale for investigation of UIB cloud fraction climatology
The larger strand of research of which this paper presents initial results is driven by two parallel working hypotheses: 1) the detailed spatiotemporal characterization of UIB cloud fraction climatology can usefully inform parameterization of energy inputs in distributed hydrological and land surface process modeling of basin tributary catchments; and 2) differentiated rates of change in diurnal components of cloud cover fraction could provide a partial causal mechanism to explain identified seasonally and diurnally asymmetrical patterns of change in UIB near-surface air temperatures in recent decades. To these ends, this paper examines currently available cloud cover datasets (MODIS and ERA-Interim) in order to describe in detail the present cloud climatology of the UIB and its potential implications for local air temperature variability. In ongoing and planned studies we intend to carry out intercomparisons of cloud fraction and related radiative flux climatologies for the UIB of an ensemble of recent reanalyses including NASA MERRA (Rienecker et al. 2011), NCEP Climate Forecast System Reanalysis (CFSR; Saha et al. 2010), and Japanese 55-year Reanalysis Project (JRA-55; Ebita et al. 2011) in addition to ERA-Interim. The purpose of this intercomparison will be to improve understanding of cloud influence on individual components of the radiative balance in the UIB.
The paper is structured as follows: Section 2 provides details of data sources and analysis methods used in this study. Section 3 presents the general cloud climatology of the UIB including a description of mean annual cycles and spatial variability and the vertical profile of annual cloud cover cycles. Section 4 addresses the issues of cloud radiative effect and apparent cloud-air temperature relationships in the NW UIB. Section 5 presents a discussion of some findings on the potential role of CRE in shaping temperature variability in the UIB and on difficulties in separating direct physical influences from apparent effects linked to day-night coupling. Section 6 concludes with a summary of the study findings and avenues for further research.

Data and analysis methods
Throughout this paper, regardless of data source, cloud cover fraction (CCF) is defined as the fractional area of a unit horizontal surface covered by all cloud within this footprint. This is for coherence with the variable returned by the MODIS data product MOD06L2 (Platnick et al. 2003). For data sources (e.g., meteorological reanalyses and some remote sensing data products) that differentiate or discretize cloud cover in a range of vertical levels the relevant total cloud fraction value has been extracted.
a. ERA-Interim estimates of cloud cover and nearsurface air temperature We use the meteorological reanalysis product ERA-Interim (Dee et al. 2011(Dee et al. ) for 1979 to the present at 0.758 horizontal resolution. In the ERA-Interim cloud cover is an ''analysis parameter'' available at 6-hourly synoptic time steps. The ERA-Interim was produced in a sequential manner for each time step increment. First, in the analysis phase, data from a wide range of observations were ingested to determine the values of a set of variables, designated ''analysis parameters,'' which were used to drive the second forecast phase determining the values of ''forecast parameters.'' Forecast parameters were subsequently used to initialize the analysis phase of the following time increment. In ERA-Interim cloud formation is determined to occur if relative humidity exceeds a pressure level-dependent threshold (Naud et al. 2014). Thus as cloud presence is determined by relative humidity, itself an analysis parameter, cloud fraction is also an analysis parameter in ERA-Interim (Dee et al. 2011). As a complement and counterpoint to local nearsurface air temperature observations, ERA-Interim 2-m air temperature, an analysis parameter in this reanalysis like cloud cover, was also analyzed. Spatial aggregates (means) for ERA-Interim in this study were calculated over the northwest UIB region as shown in purple (NW UIB) in Fig. 1. The reanalysis cloud cover data have the valuable attribute of a multidecadal (greater than 30 yr) record length, but the coarse spatial resolution limits its use for understanding spatial variability in the complex terrain of the NW UIB.

b. MODIS estimates of cloud cover (MOD06L2)
While reanalysis products estimate cloud formation/ presence through numerical simulation, including thresholds of relative humidity, data products from remote sensing use algorithms drawing on passive radiometric observations in a range of frequency bands to test for cloud presence. Thus, although the remote sensing algorithms may at times produce false positives or negatives, these data products can be considered as observations. This paper presents detailed analysis of a moderate resolution (5 km) cloud data product (MOD06L2, version 5.1) produced from MODIS imagery (Platnick et al. 2003;Li et al. 2005;Weisz et al. 2007) for a period covering 2000 to 2012. Data were downloaded from the NASA Level 1 and Atmosphere Archive and Distribution System (LAADS). Day and night cloud fraction subdatasets (SDS) were extracted from individual MODIS ''swath'' scenes (in native hdf file format) and reprojected to geographic coordinate system using the MODIS conversion toolkit plug-in (White 2013) for ENVI/IDL software. Using the RasterIO module (Holderness 2011) for Python, data were further regridded into a regular 0.058 resolution, clipped to a geographical window (27.58-40.08N, 70.08-85.08E) covering the western section of the Himalayan adjacent subregions, and temporally aggregated to produce daily and monthly means. Spatial aggregates (means) for MOD06L2 were calculated over the northwest UIB region as shown in Fig. 1. Additionally, MOD06L2 spatial statistics were calculated for each 100-m elevation section (derived from SRTM90 DEM) along the vertical (hypsometric) profile of the study area. In the present study, MOD06L2 data were only utilized from the MODIS instrument on the NASA Earth Observing System (EOS) satellite Terra. Terra is in a ''sun-synchronous polar orbit'' with local equatorial crossing times at approximately 1030 (day scenes) and 2230 (night scenes). The decision to use only MOD06L2 from Terra was made primarily due to data storage and processing constraints. Terra was selected because it entered service roughly two years prior to Aqua, thus providing an incrementally longer overlapping record with local observations and other historical datasets. Preliminary work comparing level 3 MODIS cloud data products (MOD08M3) from the instruments on both Terra as well as that on the EOS Aqua satellite (local equatorial crossing times at 1330 and 0130) showed that nighttime cloud fraction estimates over the NW UIB were essentially identical while daytime cloud fraction estimates differed by an offset ranging from 2% to 10% with greater cloud detected in the early afternoon (Aqua) than in late morning (Terra). This is illustrated in the supplementary material. MODIS daytime overpasses were compared to ERA-Interim means of the 0000 and 0600 UTC 6-h time slots (covering roughly 0500-1700 local time).
Nighttime MODIS overpasses were compared to ERA-Interim means of the 1200 and 1800 UTC time slots (covering roughly 1700-0500 local time). Preliminary work with NASA MERRA total cloud cover variable (not presented here), which is available at hourly time steps, showed very limited differences between ''single hour'' values extracted to correspond with MODIS overpass times and 12-h aggregates (corresponding to the ERA-Interim day and night means). The preliminary finding was that the differences between datasets dwarf any difference stemming from selection of time averaging.
While the good spatial resolution provided by MOD06L2 cloud cover fraction is a key strength for understanding spatial variability and topographically driven climate features in the complex terrain of the UIB, the short record length (,15 yr) of the MODIS instrument severely limits its use in assessing longerterm trends given the substantial underlying interannual variability of UIB climate. An ideal cloud cover dataset would provide equivalent (or better) spatial resolution to the MOD06L2 data product with a record length of at least that provided by ERA-Interim. At present, however, the two types of datasets, remote sensing and reanalyses, must be used in tandem to achieve the dual aims of detailed spatial characterization and assessment of longterm trends. An analysis of the agreement (correlation) for spatial mean estimates of NW UIB cloud cover from the two datasets is presented in the supplemental material.

c. Local observations of air temperature
Sparse, but spatially representative, point observations of daily maximum and minimum temperature and precipitation are available from valley meteorological stations operated by the Pakistan Meteorological Department (PMD) (see Fig. 1). Some of these stations have been operating for several decades and provide invaluable insight into the long-term climate variability of the region (e.g., Archer and Fowler 2004;Fowler and Archer 2006). Station records from Archer and Fowler (2004) and Fowler and Archer (2006) were extended to near present-day using updates provided by the Global Change Impact Studies Centre (Islamabad, Pakistan). Archer and Fowler (2004) validated the initial dataset and checked it for inconsistencies. The same methodology was applied to the updated records. Table 1 gives the elevation [in meters above mean sea level (MSL)] as well as details on the availability (record length, etc.) for each station. Although these valley stations are well below the catchment mean elevation, they provide a self-consistent record and an important record of interannual variability in climate parameters.
In the mid-1990s a network of automatic weather stations (AWS) was installed at consequential elevations (from 2100 to 4700 m MSL) by the Pakistan Water and Power Development Authority (WAPDA); they have since operated continuously. Although most AWS are below catchment mean elevations and key zones of accumulation and melt both for seasonal snowfall (nival) and perennial ice (glacial), the existence of these records provides useful data for extrapolation to higher elevations. This includes improving the vertical profiling of key parameters, including the development of realistic lapse rates, and for validating other data sources (e.g., remote sensing from satellite imagery or numerical simulation results from meteorological reanalysis or regional climate models). Unfortunately, the relatively short record length limits their applicability for assessment of climatological trends in the elevation bands covered. Furthermore, metadata for the instruments installed on these AWS units have not been made available to the authors, so data quality assessments are limited to corroboration of equivalent variables from independent data sources.

Spatial means and variability in annual cycle of diurnal cloud properties
a. Further context on the motivation for a detailed seasonal and diurnal cloud climatology for the upper Indus The general climatological interest of a detailed characterization of cloud cover in the upper Indus is further heightened by the search for causal mechanisms to explain seasonally and diurnally asymmetrical trends in local temperature observations over recent decades (Fowler and Archer 2006;Khattak et al. 2011), with cold months (September-May) experiencing rapid warming and warm months (summer; i.e., June-August) experiencing substantial cooling. Figure 3 shows that this seasonal asymmetry results from diurnal asymmetry. There have been increases in maximum temperature (Tmax) in cold months while minimum temperature (Tmin) has been nearly stationary. In contrast, summer temperature decreases have resulted from cooling Tmin while Tmax has remained stable. For simplicity, the trend rates presented in Fig. 3 are derived from linear regression. Confidence intervals (median and 5th and 95th percentiles) of temperature trend rates derived using the Theil-Sen slope method are presented in the supplementary material along with correlation coefficients and statistical significance estimates based on the Mann-Kendall and Spearman methods.
The observed summer cooling in the UIB is thought to play a role in the ''Karakoram anomaly'' of glacier stagnation and limited growth at the western end of the Himalayan arc in contrast to rapid glacial retreat in the central and eastern Himalayas (Hewitt 2005;Gardelle et al. 2012). The asymmetrical nature of cloud radiative effects over the diurnal cycle (i.e., predominantly cooling in daytime and insulating in nighttime) is a potential causal mechanism of asymmetrical temperature change for investigation. Validation of gridded cloud fraction datasets of any horizontal spatial resolution for the UIB is a challenge. In an ideal context, a substantial network of local cloud observations made either with active sensors (radar and lidar) or visually by trained experts, as in the Canadian dataset used by Betts et al. (2013), with an hourly frequency would enable collocated comparison of near-simultaneous observations. In reality, while visual cloud cover estimates have been made in past at UIB meteorological stations operated by PMD, no data are available to the authors for the period covering the MODIS observational record. Furthermore, the NW UIB ground-based radiometric observations available to this study are limited to daily cumulative shortwave radiation (converted from watt hours to mean watts per square meter).
While the absence of both manual and instrumentbased cloud observations in the study area inhibits the validation of gridded cloud fraction over the NW UIB, the physical context of the study area itself poses serious challenges for the performance of cloud detection algorithms based on passive radiometric observations. Frey et al. (2008) noted in their paper describing advances in the MODIS version 5 cloud algorithm that ''discriminating clear-sky from cloudy conditions is nowhere more difficult than in conditions of polar night'' (p. 1058). This statement arguably applies similarly to the cryosphere-dominated Himalayan mountain arc, often referred to as the ''third pole,'' as it does to the Arctic and Antarctic regions. The limitations of MODIS cloud detection at night over snow-covered areas was confirmed by Ackerman et al. (2008) who compared algorithm cloud mask results to cloud detection by the Geoscience Laser Altimeter System (GLAS) instrument which had a limited operational lifespan on the ICESat platform in the early years of the last decade.
The use of independent spaceborne cloud-detecting instruments for validation of MOD06L2 over the NW UIB does offer a potentially promising solution to the lack of local observations. In terms of ''active'' sensors to contrast the passive MODIS instrument, a merged data product combining cloud detection by the CloudSat millimeter-wavelength spaceborne weather radar and the CALIPSO lidar instrument (Mace et al. 2009) appears an initially logical choice for alternate observations. The CloudSat and CALIPSO instruments provide highly detailed vertical profiles of cloud presence and structure over the entire latitudinal range that are of great interest to scientists studying atmospheric processes. Unfortunately the detailed vertical resolution (,500 m) and ''along track'' horizontal resolution (,2 km) is counterbalanced by extremely limited area sample as the ''cross-track'' dimension sampled less than 2 km. This is in comparison with MODIS swath-widths measuring hundreds of kilometers. Thus comparing MOD06L2 with the merged CloudSat-CALIPSO product is analogous to comparing ''carpets to curtains.' ' Naud et al. (2014) did use both MODIS and CloudSat-CALIPSO data as observations for validation of ERA-Interim and NASA MERRA cloud climatologies over the Southern Ocean, but the two remote sensing datasets were used in different ways, with MODIS providing the horizontal dimension and CloudSat and CALIPSO offering vertical planes through identified warm season cyclonic events. Naud et al. (2014) also used data from an additional instrument, the Multiangle Imaging Spectroradiometer (MISR). MISR is one of the other instruments, in addition to MODIS, carried onboard Terra. As the instrument name suggests, MISR uses images captured a multiple cameras, albeit with a limited number of spectral bands, positioned at different viewing angles for cloud detection (Zhao and Di Girolamo 2004). This is in contrast to the MODIS cloud detection algorithm, which uses a single viewing geometry and relies upon the differences in ''brightness temperature'' between a substantial number of (thermal) spectral bands as well as ''reflectivity'' threshold tests in visible wavelengths . A drawback of the MISR instrument is that it functions solely in visible wavelengths and thus only produces imagery for daytime scenes. Of specific interest to this study is a level-3 MISR data product three-dimensionally aggregating ''cloud fraction by altitude'' (CFbA) to 0.58 grid boxes and 500-m vertical bins at a monthly time step (Di Girolamo et al. 2010). MISR CFbA includes a virtual bin representing total (horizontal-footprint) cloud fraction. This total cloud fraction is equivalent to the maximum value of the individual bins in the complete atmospheric column.
While the MISR CFbA data product is effectively independent of MOD06L2, there is no immediately evident basis to suppose that either of the two data products will be physically more reliable over the NW UIB. Thus comparing them is an exercise in confronting alternative estimates rather than assessing ''predictions against observations.'' In the spirit of exploring differences between alternative estimates, Fig. 4 compares daytime cloud cover fraction (CCFday) estimates from MOD06L2 to MISR CFbA as well as comparing MOD06L2 CCFday to daytime cloud climatology from three recent global reanalyses that are not otherwise examined in this paper (NCEP CFSR, NASA MERRA, and JRA-55). Each set of alternative CCFday climatology comparisons shows both commonalities and divergences. In all cases, winter [i.e., December-February (DJF)] CCFday values include annual maxima although the amplitudes of annual cycles differ markedly between data products. MOD06L2 and MISR CFbA (Fig. 4a) agree relatively closely in terms of period mean for the summer season (June-September) but disagree substantially on both the shape of the annual cycle and on absolute CCFday in winter and spring (November-May). MISR CFbA also shows much greater interannual variability-as measured by the range between the first (c10) and last (c90) decilesthrough the annual cycle than MOD06L2. MOD06L2 and NCEP CFSR (Fig. 4b) agree quite closely with the exception of midsummer to midautumn (July-October) when CFSR estimates marginally clearer conditions. This is to an extent coherent with Zib et al. (2012), who found CFSR and ERA-Interim to be the best performers among an ensemble of recent reanalyses in representing cloud climatologies of two observation sites in the Arctic. Although there are similarities in the shape of annual cycle of CCFday between MOD06L2 and MERRA (Fig. 4c), there are also consistent substantial divergences, with MERRA estimating consistently much clearer conditions than MOD06L2, particularly in midsummer to midautumn (July-October). This is coherent with Naud et al. (2014), who found MERRA to substantially underestimate cloud cover fraction over the Southern Ocean, particularly when shallow convective processes were likely to be present. Careful comparison of MOD06L2 to JRA-55 (Fig. 4d) shows that while the latter generally estimates substantially clearer conditions than the former, they agree relatively well on the shape of the annual cycle of CCFday albeit with JRA-55 presenting a much dampened (flattened) version. The conclusion we choose to draw from these four comparisons is that no presently available gridded cloud fraction climatology distinguishes itself as eminently more physically plausible than others for the NW UIB. There is general consensus on annual CCFday maxima occurring in winter with annual minima in early summer (June) or midautumn (October). There is, however, substantial disagreement on absolute values for both period means and on the range of interannual variability. In light of this uncertainty, as well as the advantages presented by MOD06L2 in terms of horizontal spatial resolution (5 km) and observational frequency (twice daily), we find that further exploration of MOD06L2 cloud cover fraction as an indicator of modulation of radiative inputs to the NW UIB climate and land surface processes is warranted.
c. Comparison of NW UIB catchment spatial mean annual cycle The MODIS and ERA-Interim spatially aggregated annual cycle of total CCFday and nighttime cloud cover fraction (CCFnight) for the northwest UIB (tributary catchments of Shigar, Hunza, Gilgit, and Astore Rivers plus the main channel from Yogo and Kharmong downstream to Besham) are shown in Fig. 4. CCF is the fraction of a given unit area (e.g., a MODIS cloud product 5-km pixel) that is covered by clouds. Even in periods with relatively clear prevailing conditions, on a mean-monthly catchment-wide basis MOD06L2 indicates substantial daytime and nighttime cloud cover (.35%). Figure 5 also shows that for much of the year the spatial mean total cloud cover from daytime and nighttime scenes are of quite similar magnitudes. The annual cycles described by ERA-Interim and the MODIS imagery match relatively well both in terms of mean values and interannual variability (spread between first and last decile, c10-c90). There are fractional discrepancies between the two datasets in some seasons. For period-mean daytime cloud (Fig. 5a), MOD06L2 shows a positive difference to ERA-Interim in from September to December. Nighttime cloud period-mean estimates (Fig. 5b) from the two datasets agree closely in winter months (November-March), but there is a consistent positive difference by MOD06L2 over ERA-Interim from April to October. In their assessment of reanalysis cloud cover estimates over the Southern Ocean, Naud et al. (2014) suggested that ERA-Interim underestimates of cloud cover fraction compared with MODIS might stem from parameterization of the cloud scheme in the underlying reanalysis model. Further investigation of this possible explanation is beyond the scope of the present work.

d. Spatial variability in annual cycles of diurnal cloud properties across the UIB and adjacent subregions
There is considerable spatial variability in diurnal cloud properties at the catchment scale (Fig. 6), which appears to be orographically defined. Daytime cloud appears to be concentrated over catchment boundaries (ridgelines). This suggests greater cloud over high-elevation surfaces, with implications for the amount of incoming radiation available to drive melt processes on glaciers and seasonal snowpack. Furthermore, there is a pattern of northward and upward migration of nighttime cloud concentrations [and night versus day excess (CCFday 2 CCFnight , 0); pink in Fig. 6] concentrated in the valleys, developing through spring, reaching its highest/northernmost extent in summer and retreating again during autumn. This highly complex spatiotemporal pattern demonstrates that some of the factors influencing local energy balance, and hence local temperature variation in the UIB, operate at a relatively restricted spatial scale, far smaller than the gauged tributary catchments. This in turn suggests that trends estimated from longrecord meteorological stations located in valley bottoms may be under the influence of physical mechanisms, such as cloud radiative effect, that are not spatially representative of wider catchment conditions. Further investigation, however, would be required to answer this question authoritatively. The absence of local active instrument cloud detection or of high-resolution (,10 km grid) regional dynamical downscaling of global meteorological reanalyses is a present obstacle to such analyses. With meticulous scene selection and geographic filtering, it should in theory be possible to investigate this using CloudSat and/or CALIPSO cloud-atmospheric profiles, but such an undertaking is beyond the scope of the present study.
Additional detail on the variability of NW UIB cloud cover at finer spatial and temporal resolutions (i.e., approaching the point scale at the daily time step) is presented in the supplementary material for those readers interested in processes operating at those scales.

e. Monthly vertical profiles of diurnal cloud properties, period means, and variability
The spatially aggregated catchment CCF for both day and night has a moderate interannual range (Fig. 5 . Spatial variability of CCF within each elevation band is described by the period medians (med) of the monthly spatial mean, maximum and minimum. Interannual variability is characterized by the period first (c10) and last (c90) deciles of the monthly spatial means.
whereas the spatial variability within the catchment, even for period means, is considerably greater (Fig. 6). In identifying the physical mechanisms that drive this spatial variability, the mean MOD06L2 CCF for each month from 2001 to 2012 was analyzed as a function of the elevation band of the underlying ground surface. This is shown in Fig. 7, which presents period medians of the monthly spatial mean for each calendar month and the 10th and 90th percentiles (i.e., first and last deciles) of the monthly spatial means within the period to indicate interannual variability. To characterize spatial variability the period medians of the monthly maximum and minimum are also shown.
This reveals a consistent vertical pattern of CCF distribution as a function of the underlying surface elevation. That is to say, while profile shapes vary considerably from month to month and between diurnal phases, for a given month and diurnal phase the vertical profiles of the statistics describing the distribution (mean, first and last deciles, etc.) have similar shapes; that is, they are separated by a roughly constant offset throughout the vertical range. CCFday (Fig. 7, columns 1 and 3) almost monotonically increases with elevation in every month. More specifically, in cold months (December-April) the vertical progression has a roughly stepwise shape whereas in warm months (May-November) the vertical gradient pattern is nearly linear. In terms of the annual cycle of CCFday spatial means presented in Fig. 5a, comparison of the individual months shown in Fig. 7 shows that decreases in CCFday in summer months compared to winter months are consistent throughout the range of catchment surface elevation.
For CCFnight (Fig. 7, columns 2 and 4), monthly vertical patterns are more complex, with relative profile maxima at, or near, the highest and lowest surface elevations in the catchment, suggesting that the most extensive cloud cover is found in the valleys and at the highest elevations. The vertical position of the CCFnight profile minima is also variable, specifically at lower elevations in winter months, moving upward to its highest elevation in summer months. The values shown in Fig. 7 are presented with respect to the fraction of surface area in that elevation band rather than the catchment total. Thus in the summer months the high CCFnight values corresponding to surface elevations below 3000 m, which make only a small percentage total catchment area, do not spatially compensate for incremental decreases in CCFnight in the middle elevation ranges (3000-5000 m) that comprise the bulk of the catchment. This is why the summer catchment spatial mean CCFnight is lower than in winter months despite the localized profile maxima.
While we have shown through Fig. 7 that vertical patterns of CCFday and CCFnight vary distinctly throughout the annual cycle, the question of the degree to which CCFday and CCFnight are linked also arises. Vertical and seasonal patterns of correlation between monthly diurnal cloud properties are explored in detail in the supplementary material.

Apparent cloud cover-air temperature relationships
In an ideal case, assessment of cloud cover influence on near-surface air temperature would be preceded by analysis of observed CRE on radiative fluxes in the study area. Unfortunately, in the case of the NW UIB only the AWS units are equipped to record radiometric fluxes, and then only cumulative SW radiation. Furthermore, AWS radiometer metadata including instrument characteristics, configuration, and calibration are not available to the authors of this study. Analysis of available SW flux observations is presented in the supplementary material in the context of MOD06L2 cloud cover and ERA-Interim radiation climatology for the NW UIB.
The lack of local NW UIB radiometric observations is analogous to the lack of local cloud fraction observations addressed in section 3a, and as such similar avenues to address this challenge-that is, comparison of an ensemble of radiative flux climatologies from meteorological reanalyses and remote sensing data products-present themselves. Logically then, the recurrent problem of the absence an effective ''ground truth'' in order to validate ensemble members remains. We intend to address these issues in ongoing and planned future studies, but the required analyses are beyond the scope of the present work.

a. Cloud cover-air temperature relationships
Because of the diurnal asymmetry of the SW and LW CRE components it is expected that maximum temperature (Tmax) from station observations would correlate negatively with CCFday. Similarly, minimum temperature (Tmin) would be expected to correlate positively with CCFnight. Furthermore, based on the variation of the relative magnitudes of the SW and LW fluxes throughout the year, CCFday-Tmax correlations would be expected to be most strongly negative in summer months, whereas positive CCFnight-Tmin correlations are anticipated to be strongest in winter months. However, results indicate that the relationship between CCF diurnal components and concurrent air temperature in the northwest UIB is more complex (Fig. 8). The spatial sampling of CCF used for these correlations is the sitespecific method described for Fig. S12a in the supplementary material (see also section S2.5 in the supplementary material). The correlations are carried out between daily time series of MOD06L2 CCF and local (PMD manned and AWS) temperature observations.
As expected, correlations between CCFday and Tmax (Fig. 8a) are nearly uniformly negative, with the strongest and most consistent negative correlations are in midspring (April) to early summer (June) and in late autumn (November). Positive CCFday-Tmax correlations are found only for December and January when incoming shortwave radiation is at its annual minimum. This is to be expected as these calendar months are the only phase of the year when longwave CRE is expected to exceed shortwave CRE during the day. The progressively stronger negative correlations in March and April compared to February may be linked to ablation of snow cover in the vicinity of the stations. While NW UIB catchment-wide ablation continues late into the spring, snow cover at the (lower) elevations where stations are located will melt significantly earlier. Symmetrically the abrupt transition from negative in November to positive in December could relate to expansion of seasonal snow cover to lower elevations. These rapid shifts in correlation sign and strength are coherent with the findings of Betts et al. (2014), who noted loss and gain of snow cover acting as a ''climatic switch'' shifting CRE effects between LW and SW dominance owing to changes in surface albedo. Figure 8 also indicates an apparent dependency on station type (manned or AWS), with the valley PMD stations showing consistently weaker correlations than those at the AWS units. This could indicate an elevation dependency, but comparison of the correlation strengths among AWS units with respect to their elevations does not show simple monotonic behavior. Another physical mechanism that could play a role here is ''hill-slopevalley circulations'' (e.g., katabatic winds and cold air drainage). Detailed investigation of these phenomena would likely require a nested, nonhydrostatic regional downscaling of a global meteorological reanalysis along the lines of that performed by Maussion et al. (2014) although perhaps to a higher spatial resolution (,5 km) to fully capture the influence of the extreme topographical relief of the study area. This, however, is beyond the scope of the present study.
Correlations between CCFnight and Tmin (Fig. 8b), however, show a greater degree of variation between stations in several months, in keeping with highly localized behavior of Tmin at valley stations (Bolstad et al. 1998;Lookingbill and Urban 1993). The expected, consistent positive CCFnight-Tmin correlations from autumn (October) to spring (April) are evident, but correlations in the remaining months yield both positive and negative values amongst the stations. Strong, statistically significant negative correlations are found in June through September. This may be partly explained in purely mathematical terms. There are strong positive correlations between CCFday and CCFnight (see the supplementary material). This means that in most cases daytime and nighttime cloud conditions are expected to be similar (in terms of sign and magnitude of anomalies), but given the physical nature of daytime (SW) and nighttime (LW) CRE the resultant temperature forcings will act in opposite directions. Which diurnal component dominates will depend upon the relative importance of SW and LW fluxes in the given season. The potential balance of diurnal CRE components can be gauged from the radiative flux climatology shown in Fig. 2. Therefore, in purely physical terms, in warm months the net cooling influence of daytime (shortwave) CRE is relatively large and carries over into the nighttime when the (LW) insulating influence of nighttime CRE is not sufficient to compensate for the daily net reduction in incoming energy. In summary it can be difficult to dissociate daytime and nighttime cloud effects in a correlation analysis. Indications of dependencies on station type and elevation in correlations between Tmin and nighttime cloud are limited. The strong negative (summer) correlations are predominantly, but not exclusively, found in AWS units. This is consistent, in the apparent carryover of SW cooling CRE, with the finding of tighter linking between day and night cloud cover at higher elevations (see the supplementary material). In a study of CRE influence on the diurnal cycle of temperature in the Canadian prairies, Betts et al. (2013) also found dominance of SW effects in the warm season, although in their results correlations between cloud cover and Tmin were weak or negligible rather than negative.
b. CRE as a potential explanation for seasonal and diurnal asymmetry in UIB air temperature trends Figure 9 presents a range of trends rates in CCF diurnal components over the northwest UIB from ERA-Interim (1979 and compares individual CCF components to corresponding components of near surface air temperature (Tmax, Tmin, and DTR) trends from PMD stations. Ideally trends would be assessed at a spatial resolution comparable to MOD06L2 rather than the coarse scale of global meteorological reanalyses. The short record length of the MODIS data product, however, prohibits its use in this way. Despite this limitation, correlations between monthly time series from ERA-Interim and MOD06L2 northwest UIB of spatial mean cloud cover for the common record period (2000-12) are respectively 0.857 and 0.797 (p , 0.001) for CCFday and CCFnight. Equivalent values from Mann-Kendall and Spearman rank correlation tests are 0.662 and 0.843 respectively for CCFday and 0.573 and 0.768 for CCFnight (in all cases p , 0.001). The Pearson correlations between the two datasets for CCFday and CCFnight by calendar month are shown in Fig. S4 of the supplementary material while Mann-Kendall and Spearman rank correlation values and estimates of statistical significance are given in Table S3 of the supplementary material.
Given the consistent patterns of correlation between MOD06L2 and locally observed near-surface air temperature presented in the preceding section, it is logical to look for causal mechanisms of diurnally asymmetrical temperature change in diurnally differentiated rates of cloud cover change. Figure 9a compares the estimated rate of change in CCFday and CCFnight along with the rate of change in CCFrange (i.e., the daily difference of CCFday minus CCFnight). The CCFrange trend is mathematically equivalent to the differential trend rate (i.e., CCFday trend minus CCFnight trend). In the ERA-Interim dataset, CCFday and CCFnight monthwise trends track very closely for much of the year. Both show substantial changes from autumn to spring (October-May). CCFrange trend in this period is near to zero. In the remaining (summer) months CCFrange is substantially positive (i.e., there is a relative decrease in CCFnight with respect to CCFday). This partial decoupling of CCF diurnal component trends, estimated from ERA-Interim, in summer months is consistent with the season variations in correlation between CCFday and CCFnight using the much higher spatial resolution of MOD06L2 (see the supplementary material). Figure 9b shows trends in diurnal temperature range (DTR) from PMD stations compared to trends in daily mean cloud cover (CCFavg) from ERA-Interim with the latter displayed on an inverted axis. The reason for the reversal of the CCFavg axis is the expected to negative correlation to DTR. By suppressing Tmax through SW CRE and enhancing Tmin through LW CRE, CCFavg dampens the diurnal temperature cycle, quantified as DTR. The Pearson's r values, and associated statistical significance, between the calendar month sequences of trends in CCFavg and DTR shown in Table 2 confirm that strong negative correlations exist between patterns of change in these variables at all three of the long record stations. Similar results, expected due to SW CRE, are found for the linkages between CCFday and Tmax as can be seen in Fig. 9c and Table 2. The expected linkage, due to LW CRE, between patterns of change in CCFnight and Tmin is not seen in Fig. 9d. Table 2 confirms weak correlations of variable sign lacking statistical significance. This absence of causal influence on Tmin is supported, albeit only for the warm season, by the findings of Betts et al. (2013) for the Canadian prairies. This ''non-result'' for CCFnight-Tmin correlation may be due the coarse resolution of ERA-Interim inadequately resolving spatial heterogeneity within the study area. Alternatively, it may be that SW CRE dominates temperature change throughout the diurnal cycle, with patterns of change in Tmin predominantly resulting as the sum of changes of Tmax and DTR. The dominance of SW CRE is suggested by large Tmax increases coincident with decreasing cloud cover (October-June). The strong warming concurrent with decreasing cloud cover is similar to findings for the central and eastern Tibetan Plateau (Duan and Wu 2006).
The Tmin decreases over a period of net stable cloud cover (July-September), however, require further explanation. This case could be explained by looking at the detail provided by Fig. 9a. Stable CCFday would imply stable daytime CRE. Decreasing CCFnight would imply decreasing LW CRE and greater radiative nighttime cooling (hence decreasing Tmin).
Spatially detailed detection of diurnally asymmetrical changes in CCF yielding persistent shifts in the surface energy balance would require a much longer dataset record than that available from MODIS. It is conceivable that customized multidecadal spatial climate data products including cloud cover estimates analogous to MOD06L2 could be developed for the Himalayan region using Advanced Very High Resolution Radiometer (AVHRR) imagery. This would satisfy the need recognized here for a long-record spatial cloud dataset with resolution comparable to MOD06L2 to substitute for the current use of coarseresolution reanalysis (ERA-Interim) in trend estimates. The authors are currently working to produce such a dataset (Hardy et al. 2015, unpublished manuscript) using snow-cloud discrimination methodology similar to that of Zhou et al. (2013).
An alternate hypothesis is that the locally observed air temperature trends are dominated, or at least substantially influenced, by advection of warm or cool air masses through changes in large-scale circulation. Mölg et al. (2014) found spatial patterns of statistically significant correlations between atmospheric flow (300-hPa wind speed) over the Tibetan Plateau and local summer air temperature. Changes in CCF could be driven by changes in circulation regime. Kotarba (2010) found strong links between atmospheric circulation types and CCF for a study of the Tatra Mountains in central Europe. Testing of this ''circulation hypothesis'' is beyond the scope of the present work, but further investigations are planned to examine the links between cloud cover variability and large-scale atmospheric circulation.

a. Factors beyond CRE mechanisms influencing air temperature variability in the UIB
It was expected that there would be a degree of ''diurnal carryover'' in apparent CCF air temperature correlations in any given month. Hence, in summer strong negative SW daytime CRE could be buffered by positive nighttime LW CRE while in winter strong positive LW CRE from seasonally high cloud cover could be dampened by increased negative SW daytime CRE. This appears to effectively be the case during the annual extrema (winter and summer).
Furthermore, while SW and LW CRE act locally via the surface energy balance to modulate air temperature, anomalies in cloud cover relative to reference climatology may in fact be a manifestation of a shift in the predominant regional circulation patterns. A tendency to prevailing clear or overcast conditions may be linked to specific source areas of incoming air masses. These air masses may bring with them, by advection, air temperature anomalies that either reinforce or counteract CRE. In the winter months, climate variability over the western Himalayas, including the UIB, is dominated by synoptic weather systems referred to as ''westerly disturbances'' (Yadav et al. 2009). Modeling studies have shown these systems to be dependent upon the subtropical westerly jet (Dimri et al. 2013). In the summer months regional climate variability is driven by the degree to which development of the South Asian summer monsoon displaces the westerly jet northward (Saeed et al. 2012). Improved understanding of the relationship between latitudinal positioning of the westerly jet and characteristic temperature and cloud cover anomalies could help to account for substantial air temperature variance not explained solely by CRE.

b. Applications to catchment hydrology in the northwest UIB
The fundamental cloud climatology described in this paper has important potential applications for understanding the processes which influence hydrological variability in the northwest UIB. Section 3e (Fig. 7) showed consistent, coherent monotonic gradients of CCFday with respect to elevation. This may have important implications for hydrological modeling in incorporating cryosphere dynamics to calculate runoff generation from melting snow and ice. Methods for snowmelt and glacier melt estimation range from relatively simple ''degree day'' temperature index methods (Singh et al. 2000) to full resolution of the surface energy balance (Marks and Dozier 1992) with intermediate levels of complexity (Pellicciotti et al. 2012). For methods using radiative fluxes it would be logical to use vertical gradients in CCFday to inform spatial modulation of incoming shortwave radiation (i.e., simulating spatially heterogeneous downward shortwave radiation due to topographically driven CRE).

Summary and conclusions
Clouds play a key role in hydroclimatological variability by modulating the surface energy balance and air temperature. This study has used MODIS MOD06L2 cloud cover data, with corroborating analysis of ERA-Interim cloud estimates, to produce a cloud climatology for the upper Indus River basin (UIB), with specific focus on tributary catchments in the northwest (NW UIB) contributing a large fraction of basin annual runoff. There is general agreement between MODIS and ERA-Interim on the annual cycle of CCF in the northwest UIB with spatial mean monthly values ranging from 50% to 80%. These means, however, conceal substantial spatial variability, particularly along seasonally varying vertical profiles. This study found that, despite a general spatial tendency for CCFday to exceed CCFnight for most locations within the catchment, there is a coherent progression upward and northward during the summer months of a regime where CCFnight exceeds CCFday. The relative balance between CCFday and CCFnight may have important implications for local net cloud radiative effect (CRE). This study also used the spatial resolution provided by MOD06L2 cloud data to investigate how the correlation between (concurrent) CCFday and CCFnight varies in the northwest UIB with surface elevation throughout the year.
Correlations between local observations of near-surface air temperature and MODIS cloud cover fraction confirm the strong linkages between local atmospheric conditions and near surface climate variability. The apparent CRE influence on near-surface air temperature in the northwest UIB was shown to exhibit substantial complexity. Correlations between daytime cloud cover and Tmax appeared to confirm the dominance of SW (cooling) CRE with strong negative values in all seasons except winter (annual SW minimum). Correlations between nighttime cloud cover and Tmin were less straightforward. Expected strong positive correlations, coherent with dominant LW (insulating) CRE, were found for winter months, but in summer correlations were mixed with some strongly negative values. This suggests that the summer month nighttime cloud cover-Tmin correlation reflects the net CRE on the average temperature (Tavg) with daytime cooling carried over into nighttime hours rather than indicating a physical cooling effect of nighttime cloud in the warm season. Alternatively, it could indicate that positive nighttime cloud anomalies in summer are linked to weather systems-circulation patterns that decrease Tmin through advection of cooler air masses.
Pursuing the consideration of possible causal influences of cloud cover fraction on northwest UIB air temperature, the potential role of CRE in recognized seasonally and diurnally asymmetrical temperature trends over recent decades was also assessed. Locally observed air temperature trends were first compared to a spatial mean analog from ERA-Interim, showing strong corroboration with asymmetrical patterns. Local temperature trends were also compared to ERA-Interim-derived trends in cloud cover diurnal components. Specifically, strong increases in maximum temperatures and diurnal temperature range occur in months when daytime cloud cover is decreasing. This further supports the role of shortwave CRE as the dominant element in cloud-air temperature linkage. Furthermore, there appears to be a plausible physical mechanism for decreasing minimum temperature in peak summer months. ERA-Interim trends in daytime and nighttime cloud cover suggest a differential change with relatively less nighttime cloud cover leading to reduced longwave CRE and hence greater overnight cooling. This summer nighttime cooling may be in part responsible for the recognized ''Karakoram anomaly'' of glacier stagnation or growth (Hewitt 2005;Gardelle et al. 2012) in sharp contrast to rapid glacial retreat-downwasting elsewhere along the Himalayan mountain arc.
While the cloud cover-temperature trend linkages suggested by the ERA-Interim cloud cover fraction variables are promising, the evidence would be more robust if the data were from a multidecadal record of remote sensing cloud observations of similar spatial resolution to the MOD06L2 data product. This underlines the value of CCF data products under development derived from the AVHRR instrument flown on successive generations of NOAA Polar-Orbiting Operational Environmental Satellites (POES).