MODIS observed increase in duration and spatial extent of sediment plumes in Greenland fjords

The freshwater flux from the Greenland Ice Sheet (GrIS) to the North Atlantic Ocean carries extensive but poorly documented volumes of sediment. We develop a suspended sediment concentration (SSC) retrieval algorithm using a large Greenland specific in situ data set. This algorithm is applied to all cloud-free NASA Moderate Resolution Imaging Spectrometer (MODIS) Terra images from 2000 to 2012 to monitor SSC dynamics at six river plumes in three fjords in southwest Greenland. Melt-season mean plume SSC increased at all but one site, although these trends were primarily not statistically significant. Zones of sediment concentration> 50 mg L−1 expanded in three river plumes, with potential consequences for biological productivity. The high SSC cores of sediment plumes ( > 250 mg L−1) expanded in one-third of study locations. At a regional scale, higher volumes of runoff were associated with higher meltseason mean plume SSC values, but this relationship did not hold for individual rivers. High spatial variability between proximal plumes highlights the complex processes operating in Greenland’s glacio–fluvial–fjord systems.


Introduction
The Greenland Ice Sheet (GrIS) is the largest ice mass in the Northern Hemisphere, containing the equivalent of 7 m of global sea level (Bamber et al., 2013).High rates of seasonal mass turnover, coupled with decadal-scale ice sheet mass imbalances (142 ± 49 Gt yr −1 between 1992 and 2011), result in hydrologic fluxes from the GrIS exceeding 300 km 3 yr −1 (Hanna et al., 2008;Shepherd et al., 2012).
Mass loss from the GrIS primarily occurs through iceberg calving and meltwater runoff.These losses occur at two distinct types of outlets: (1) approximately two-thirds of mass loss occurs at marine terminating outlets where both processes operate, and (2) one-third of mass loss occurs via meltwater runoff at land terminating outlets (Bamber et al., 2012).This paper focuses on rivers draining the landterminating ice sheet margin.Total ice sheet imbalance has increased four-fold over the past two decades, increasing from 51 ± 65 Gt yr −1 in the 1990s to 211 ± 37 Gt yr −1 between 2000 to 2011 (Shepherd et al., 2012).This imbalance is the result, in part, of an increase in the ablation area of the GrIS (McGrath et al., 2013) and increased rates of meltwater production (Ettema et al., 2009).
Runoff exits the ice sheet at more than 400 distinct outlets (Lewis and Smith, 2009) and typically contains large volumes of sediment produced beneath the ice sheet.Sediment is produced through two key mechanisms beneath temperate based glaciers: first, subglacial water pressure can create and/or exploit cracks in bedrock that allow ice to quarry and remove it and second, rocks embedded in basal ice (e.g., those quarried) abrade bedrock over which it flows (Anderson and Anderson, 2010).These processes are both partially controlled by the presence and volume of water, which varies on temporal scales of hours to months (McGrath et al., 2011).For instance, variations in surface melt entering the ice sheet can modulate the effective pressure at the ice sheet's bed, thus enhancing quarrying and can also modify the basal sliding velocity, which in turn can increase abrasion (Anderson and Anderson, 2010;Zwally et al., 2002).Sediment produced subglacially may be stored until removed by water or ice (1).Once exported by the glacier (2), sediment can be transported by rivers (3).While transported by the river, it can potentially be stored or passed through lakes ( 4), stored on the river's outwash plain (5) until it is discharged into the ocean (6).Sediment stored on the river's outwash plain may be remobilized at a later time by the river (5).Additional sediment can be sourced locally from the landscape, for example through bank erosion.Photograph courtesy of Dorthe Petersen, ASIAQ.
A complex and highly variable relationship exists between the glacier hydrologic network and the transport of sediment (Fenn et al., 1985;O'Farrell et al., 2009;Østrem, 1975).Spatio-temporal variability in sub-glacial hydrologic networks can drive variability of sediment emerging from glaciers (Collins, 1990) and transient flushes of sediment from glaciers in Greenland have been observed (Stott and Grove, 2001).Water accessing areas of the bed not previously reached may cause increases in SSC per unit discharge (Alley et al., 1997).On the other hand, areas of the bed that the sub-glacial hydrologic network often reaches may be depleted in sediment (Alley et al., 1997).
Once water and sediment exit the glacier's hydrologic system and enters a river, fluvial processes control its transport.As Fig. 1 illustrates, sediment produced subglacially may be stored until removed by water or ice (1).Once exported by the glacier (2), sediment can be transported by rivers (3).While transported by the river, it can potentially be stored or passed through lakes (4), stored on the river's outwash plain (5) until it is discharged into the ocean (6).Sediment stored on the river's outwash plain may be remobilized at a later time by the river (5).Additional sediment can be sourced locally from the landscape, for example through bank erosion, but this is not treated explicitly here (Mernild and Hasholt, 2009).
Upon entry into the ocean, river water decelerates rapidly and entrained sediment settles.In the ocean sediment travels as a function of the velocity of the water that moves it away from the river mouth and the rate it sinks from the surface layer (and satellite's view).Initially, the rivers velocity at its mouth (u 0 ) sets down fjord velocity (Syvitski, 1990).The extent of fluvial control typically ends ∼ 5.2 river widths seaward of its mouth (Albertson et al., 1950).Beyond this point, the fjord's estuarine circulation conveys the sediment seaward (Inall and Gillibrand, 2010).Sediment settling velocity (removal rate) is a function of its grain size, density, and biological process such as flocculation that alter the particle's characteristics (Syvitski and Murray, 1981).A sediment plume's centerline SSC can be written as a function of distance from river mouth (x) as (modified from Syvitski et al., 1988) where SSC 0 is the SSC of the river as it enters a fjord and λ is the size dependent particle removal rate constant (Peckham, 2008;Syvitski et al., 1985).
A number of additional processes operating in the fjord environment can influence sediment plume dynamics.The large-scale Coriolis effect can deflect the plume, although this is dependent on river mouth geometry, fjord width, and water column stratification in the fjord (Cushman-Rosin et al., 1994).Wind can cause vertical mixing of the water column in sea-ice free conditions (Cottier et al., 2010).Often the topography of the fjord helps steer winds up or down fjord.During summer, sea breezes caused by differential heating of the land and water surface can cause strong up fjord winds.In winter, strong katabatic winds can cause strong down fjord winds (Seidenkrantz et al., 2007).This bi-modal wind flow can help aid or retard estuarine circulation and cause localized upwelling (Skogseth et al., 2007).
The relationship between water discharge and river SSC is largely unstudied in Greenlandic rivers, although a power law relationship often exists between a river's SSC and its discharge, Q: river SSC = aQ b (Morehead et al., 2003;Syvitski et al., 2000), where a and b are rating parameters.This relationship is often valid for braided rivers with glacierized catchments.Church (1972) found this relationship for the braided Lewis River (90 % glacierized catchment including the Barnes ice cap) and the upper south river, Ekalugad fjord (50 % glacierized catchment on Baffin Island).Two examples from Greenland support this simple relationship.First, in situ observations on the Watson River show increased sediment load during high discharge years (Hasholt et al., 2013) and second, Chu et al. (2012) found that Polar MM5 modeled positive degree-days, a proxy for meltwater production, was the most significant driver of regional fjord SSC variability.
While Greenland fjord SSC is understudied it offers the potential to be an important indicator of environmental change and is relevant to subglacial sediment dynamics and coastal biological productivity.Sediment delivers essential nutrients to the coastal ocean (Statham et al., 2008), although high levels of suspended sediment also reduces the availability of light (Retamal et al., 2008).Sediment may impact biological productivity, such as pelagic zooplankton feeding efficiency and production rates (Arendt et al., 2011).
Widespread in situ observations of river and fjord SSC in Greenland are logistically infeasible; in fact Hasholt (1996) and Hasholt et al. (2006) identified < 15 locations where river sediment dynamics had been studied, even for a short period of time.Instead, recent efforts have focused on using orbital remote sensing to monitor fjord sediment dynamics (Chu et al., 2009(Chu et al., , 2012;;McGrath et al., 2010;Tedstone and Arnold, 2012).
Here, we build on previous efforts using orbital remote sensing to monitor fjord sediment dynamics along the coast of the GrIS (Chu et al., 2009(Chu et al., , 2012;;McGrath et al., 2010;Tedstone and Arnold, 2012).We develop an improved SSC retrieval algorithm based on a robust set of in situ SSC observations.This algorithm is applied to the entire 13 year cloud-free MODIS visible satellite imagery record at three sites along the western margin of the GrIS to better characterize SSC variability at these sites.Further, we explore the complex relationship between fjord SSC and GrIS runoff and test the hypothesis that increased freshwater discharge due to increased ice sheet melt has elevated melt-season SSC.

Study areas
Our study sites consist of three fjords in southwest Greenland (Pakitsuup, Kangerlussuaq, and Ameralik), fed by six proglacial river systems (Fig. 2).GrIS catchment areas (Fig. 3) were calculated from the local hydrostatic pressure field of the ice sheet, which combines surface elevation data with basal topography (Lewis and Smith, 2009;Cuffey and Patterson, 2010).Surface elevations are from the Greenland Mapping Project (GIMP) digital elevation model (30 m posting, Howat et al., 2014).Bedrock topography of the Pakitsuup region is from the 500 m resolution University of Kansas Center for the Remote Sensing of Ice Sheets (CRESIS) Jacobshavn basal topography data set (https://data.cresis.ku.edu/data/grids/).For the Kangerlussuaq and Ameralik regions a coarser 5 km resolution basal topography from Bamber et al. (2001) was used.Flow routing and catchment delineation was performed following a D-infinity approach using RiverTools 3.0 (Peckham, 2009).Catchment areas describe only the portion of the watershed covered by glacial ice, as these areas are the primary sources of water and sediment to the systems studied.

Pakitsuup fjord
Pakitsuup fjord (also spelled Pakitsoq) is located near the town of Ilulissat and receives melt from two small river systems.The Pakitsuup north river (69 • 31 45.5 N, 50 • 13 23.7 W) has two branches, respectively, 8 and 10 km long.As Fig. 3a shows, it drains a catchment of the GrIS ∼ 302 km 2 (Table 1), which reaches a maximum height of ∼ 1169 meters above sea level (m a.s.l.).The Pakitsuup south river (69 • 26 5.9 N, 50 • 27 53.3W) flows 7 km from the margin of the GrIS that feeds it along its southern tributary.A northern tributary flows 20 km through a series of lakes, before combining with the southern tributary approximately 4 km from the river mouth.It is the smallest catchment studied (∼ 143 km 2 ) with the lowest maximum elevation (∼ 1020 m a.s.l.).Mernild et al. (2010) estimate that the equilibrium line altitude (ELA) for the region is approximately 1125 m a.s.l.(Fig. 3a).

Kangerlussuaq fjord
Three major river systems discharge water and sediment into Kangerlussuaq fjord (Søndre Strømfjord).The Watson River (66 • 57 53.7 N, 50 • 51 49.8 W) enters at the northeast head of Kangerlussuaq fjord.A discharge gauging station has operated on the river since 2007 at the town of Kangerlussuaq (asterisk on Fig. 2b), 10 km from the river mouth (Hasholt et al., 2013).The Watson River receives water from two tributaries, each ∼ 34 km in length and has a combined GrIS catchment area of 3639 km 2 , (Fig. 3b).This is smaller than previously published values (e.g., Hasholt et al., 2013) where the catchment area is based on the surface topography.This catchment has a maximum elevation of 1860 m a.s.l.(Table 1) and a mean ELA of ∼ 1550 m a.s.l.(van de Wal et al., 2012).
The Umiiviit River (66 • 50 01.6N, 50 • 48 37.3 W) flows into the southeastern head of Kangerlussuaq fjord and drains a catchment area of 6320 km 2 with a maximum elevation of ∼ 2390 m a.s.l.In comparison, the Umiiviit River is > 85 km in length, more than double that of the Watson River branches.
The Sarfartoq River (66 • 29 29.6 N, 52 • 01 29.5 W) discharges into Kangerlussuaq fjord ∼ 79 km down fjord from the northeast fjordhead.From its mouth, the river stretches 53 km in a narrow valley to a lake, Tasersiaq (T on Fig. 2b).Tasersiaq is ∼ 71 km long and 14 m deep (Goldthwait et al., 1964), and bordered by the Sukkertoppen ice cap, lobes of the GrIS, and an 18 km long lake, Tasersiaq Qalia (TQ on Fig. 2b).The Sarfartoq River has 5385 km 2 of glaciated catchment with a maximum height of 2157 m a.s.l.Approximately 919 km 2 of this area is covered by minor ice caps and glaciers not part of the GrIS (e.g., the Sukkertoppen ice cap).It is the only river in the study that receives melt from both peripheral glaciers and the GrIS.

Ameralik fjord
Ameralik fjord (also known as Lysefjord, or Ameralgda fjord if considering only the innermost part of the fjord) is located close to Greenland's capitol city, Nuuk.The Naujat Kuat River (64 • 12 37.5 N, 50 • 12 31.0W) flows 26 km from the margin of the GrIS, receiving melt from a catchment 460 km 2 , with a maximum height of 1342 m a.s.l.Mote (2000) estimated the regional ELA near 1450 m a.s.l.Meltwater discharge originates from two lobes of the GrIS.The southern lobe discharges water into an 8 km long lake, Isortuarssuk (Fig. 2), a short distance from its terminus.The northern lobe, Kangaussarssup Sermia (KS on Fig. 2), discharges water to the river uninterrupted by lakes.On the north side of Kangaussarssup Sermia ice margin lakes exist.These lakes historically drained down a valley (the Austmannadalen) to a point ∼ 7 km from the river mouth.Since 2004, the lakes drain north into Kangiata Nunaata Sermia (KNS on Fig. 2c, Weidick and Citterio, 2011).As the MODIS record straddles this drainage realignment, Fig. 3c and Table 1 provide two delineations, one with the lakes (Naujat Kuat) and one without (Naujat Kuat-Small).At its smallest, the on-ice catchment area is 356 km 2 , at largest it is ∼ 460 km 2 .This uncertainty in drainage basin delineation does not affect the maximum elevation of the drainage or the general shape of the hypsometric curve.Unless otherwise noted, the Naujat Kuat-Small catchment was used in subsequent analysis.3 Methods

Suspended sediment samples
We collected surface water samples as part of oceanographic surveys during the 2008, 2010, 2011, and 2012 summer seasons (Table 2).In total, 140 samples, covering 12 days in all major melt months were used to develop the SSC retrieval algorithm.Three additional water samples were collected from Orpigsoq fjord (68 • 38 30.9 N, 50 • 54 45.5 W) on 2 July 2011 which are used in algorithm development, although this fjord was not included in the present study.Surface water samples were collected along variably spaced transects in each fjord (Fig. 2) by dipping a 1000 mL bottle into the near-surface layer while a handheld GPS noted the geographic location.Samples were then vacuum filtered through pre-weighed Whatman GF/F filters, rinsed with distilled water to remove salt, freeze dried, and reweighed to determine SSC.Fjord SSC ranged between 1.2 and 716 mg L −1 with a mean SSC of 73 mg L −1 .One sample collected from the shallow river mouth of the Watson river on 9 July 2011 had an SSC of 1581 mg L −1 but was excluded from analysis because it could not be determined if the sampling boat had disturbed bottom sediments.

SSC retrieval algorithm based on MODIS reflectance
MODIS Terra MOD02 250 m imagery and the MOD03 geolocation data product for each swath image were downloaded using NASA's Reverb Earth Science Data Discovery Tool (http://reverb.echo.nasa.gov/reverb/).Image processing was completed using Exelis ENVI 4.8 software, including the MODIS Conversion Toolkit (http://www.exelisvis.com) and included georeferencing, a bow tie correction, and an atmospheric correction using dark object subtraction.The latter is completed using a cloud-free open ocean pixel from each scene (Kirk, 2011;Miller and McKee, 2004).Georeferenced images were compared against the Danish National Survey and Cadastre (DNSC) Greenland coastline vector file and were manually adjusted to agree with the DNSC coastline if a discrepancy existed.SSC values were compared to co-located band one (620 to 670 nm) and two (841 to 876 nm) reflectance values.Barbieri et al. (1997) defined the MOD02 reflectance data product as the bidirectional reflectance factor of the Earth (ρ ev ) multiplied by cosine of the solar zenith angle at the earth scene (θ ev ) ρ ev × cos θev . ( With ρ ev defined as where L ev is the earth view spectral radiance (W m −2 sr −1 nm −1 ), and E i is solar irradiance (W m −2 nm −1 ).Combing these two equations, reflectance is calculated as Ocean color remote-sensing literature often define R rs as radiance reflectance (Kirk, 2011) or remote-sensing reflectance (Mobley, 1999), though technically E i should be downwelling irradiance directly above the surface of the water (E d ), not from the top of atmosphere (E i ).Despite this difference, we use the R rs notation.Following Miller and McKee (2004), variations in scene solar zenith angle were corrected on an individual pixel basis to arrive at We developed a fjord SSC retrieval algorithm using both MODIS band one and two and the co-located in situ sediment samples (locations shown on Fig. 2; relationship between variables on Fig. 4).Unlike previous SSC algorithms developed using MODIS imagery for Greenland, which became saturated above 80 mg L −1 (e.g., Chu et al., 2009)  infrared, allows the algorithm to be sensitive to higher SSC.The final retrieval algorithm was The in situ SSC data set included samples from fjords with variable salinity.Using a Seabird 19 conductivity, temperature, and depth probe, we measured salinities as high as 25 PSU (Ameralik fjord) and as low as 1 PSU (Kangerlussuaq fjord) in the top meter of the surface layer coincident with surface water sample collection locations.The variable salinity during our oceanographic surveys, and the independence of Eq. ( 6) on this parameter, suggests that this retrieval technique can be applied to disparate locations with variable oceanographic properties.Following algorithm development, we processed ∼ 300 sea ice free swath scenes per year (spanning day of the year 135 to 295) to map fjord SSC values over the entire melt season.Appendix A describes the full processing details.

Error from sampling irregularities caused by clouds
Since the year 2000, MODIS Terra satellite observations have provided an ongoing daily record in this region, although, as common in polar regions, the presence of clouds compromise the imagery's utility as an observational tool.In total, clear-sky conditions existed for between 5 and 21 % of all available fjord pixels during a given summer season.
Computational experiments following Gregg and Casey (2007) were used to assess this potential sampling bias.A simulated data set based on the entire 13 year data record was created and then masked with MODIS observed cloud cover to understand how cloud obscuration impacted SSC statistics.
To create the simulated data set, we calculated a 13 year mean SSC map on a pixel-by-pixel basis as 13 year SSC =  where N was the number of times the pixel passed all tests and was used in the analysis, i is pixel location in the horizontal direction and j is pixel location in the vertical direction.The cloud masking scheme set SSC values to equal nota-number if a cloud was detected, and thus excluded cloud influenced reflectance values from the statistics.
We then used the 13 year mean SSC map as a spatial template to mimic seasonal SSC fluctuations.For each scene in the season, the 13 year mean SSC map was multiplied by a sinusoidal factor (the scene intensity, I s ) that peaked near day 190 of each melt-season simulation (Fig. 5).SSC i,j = 13 year SSC i,j × I s , where x was a number between two and four randomly generated for each scene, and y began at one and increased by one with each scene used.This range was selected by tuning simulated SSC values to the range found in the observed MODIS data set.As an example, Fig. 5 plots how mean SSC was allowed to vary seasonally during computational experiments.

Optimizing the quantification of plume characteristics
Three metrics were selected to describe plume dynamics and trends over the 13 year record: melt-season mean SSC (SSC msm ), percent available fjord pixels in a melt season above 50 mg L −1 (P >50 ), and above 250 mg L −1 (P >250 ).
The melt season was defined as day of year 135 through 295 (∼ 15 May through 22 October).SSC msm was selected because cloud masked simulated values were closest to noncloud masked simulated values.P >50 was selected because it matches the threshold identified by Arendt et al. (2011) as having a significant impact on copepods and secondary production near Nuuk in Godthåbfjord.P >250 was selected because SSC above 250 mg L −1 is a proxy for the plume's zone of flow establishment (Albertson et al., 1950;Syvitski et al., 1988).P >50 and P >250 were also the most consistent metrics between years in that they were similarly influenced by clouds each year.
Computational experiments explored cloud bias in a number of annual metrics (mean, median, etc.) by imposing the MODIS derived 13 year record of cloud cover on the simulated data set (Fig. 6).In addition, these experiments evaluated how metrics calculating the percent of available fjord pixels in a year observed above a specified SSC threshold (50 or 250 mg L −1 ) performed.These P >t metrics describe the duration and the spatial extent of a river plume for a given  region of interest (ROI).They were calculated with two steps.First, for each scene a threshold mask (STM) was created, where all pixels in a ROI greater than the given SSC threshold t were assigned a value of one and all others a value of zero.Second, STM values were summed for an entire melt season and corrected for pixels classified as fjord (F p ) and the total number of scenes in a year (S) with Eq. ( 10): Variable plume area and fjord shape caused each ROI to vary in size and fjord area.The Pakitsuup north ROI contained a fjord area of 4.4 km 2 and Pakitsuup south 4.3 km 2 .Watson, Umiiviit, and Sarfartoq ROIs contained 41. 4, 42.8, and 65.1 km 2 of fjord area, respectively.The Naujat Kuat ROI contained a fjord area of 13.0 km 2 .Tables 3 through 5 report ROI metrics normalized by each ROI 13 year mean to account for this variability.
As Fig. 7 illustrates, SSC msm deviated the least from the simulated data set.At most, it underestimated simulated values by 16 % and overestimated by 38 %.On average it underestimated by 16 % and overestimated by 12 %.The maximum range of percent errors found for SSC msm was 54 % (Naujat Kuat), with a mean of 36 %.P >50 estimates were between 67 and 91 % below simulated values.The maximum range in the deviation from the simulated using P >50 for a plume was 22 % (Umiiviit), with an average range of 15 %.As Fig. 10d-f show, P >250 estimates were similar, between 63 and 92 % below simulated values.The maximum range in the deviation from the simulated was 20 % (Sarfartoq), with an average range of 16 %.Cloud obscuration prevents P >t metrics from correctly calculating the true number of times fjord pixels exceed a given value.Thus, absolute values reported may not be reliable.However, we argue this bias did not compromise the usefulness of the metric to detect trends.
We assigned error bars that symmetrically divided the mean range of variability found (SSC msm ± 18 %, P >50 ± 7.5 %, P >250 ± 8 %).The large bias in P >t values should also be remembered when interpreting results.We then added error associated with the SSC retrieval to bring overall error assigned to melt-season metrics to ±19 % for SSC msm and ±10 % for P >t metrics.

RACMO2 runoff data
To extend our analysis beyond the in situ record, we used yearly catchment scale RACMO2 modeled ice sheet runoff volume (Bamber et al., 2012, hereafter referred to as RACMO2 runoff) as a substitute for observed discharge.Catchments geographically match those in Fig. 3, though differences in area exist due to RACMO2's 5 km grid.For the four year period of overlap, Watson River observed discharge and RACMO2 runoff were positively correlated (r 2 = 0.5), although RACMO2 runoff overestimates observed discharge values by 23 to 113 % per year (0.79 to 2.95 km 3 ) perhaps due to underestimation of englacial storage (Overeem et al., 2014).We did not use RACMO2 runoff for either the Sarfartoq catchment or the Umiiviit River.For the former, no geographically appropriate outlet could be identified and for the latter, the misfit between RACMO2 runoff values and in situ field discharge measurements was too great to be reliable.

Trends in mean suspended sediment concentration
We hypothesize that increased freshwater discharge due to increased ice sheet melt has elevated the melt-season SSC of Greenlandic fjords.Overall, SSC msm increased for five out of six plumes, although this increase is only statistically significant for the Sarfartoq River plume (Fig. 8).The Pakitsuup south river decreased a statistically insignificant 0.6 ± 3.9 mg L −1 yr −1 .The Sarfartoq River plume SSC msm increased 2.1 ± 1.8 mg L −1 yr −1 (p = 0.013).We offer two plausible explanations for why this river's trend is statistically significant.First, the Sarfartoq ice sheet catchment may be most sensitive to atmospheric warming because it has the largest below ELA catchment area studied (3291 km 2 , 61 % of its total catchment), and hence large volumes of meltwater production are likely to occur, and are more likely to exit the catchment area, rather than refreeze in the firn (Harper et al., 2012).Second, the Sarfartoq watershed contains two long lakes, Tasersiaq and Tasersiaq Qalia, that may influence plume SSC.Tasersiaq's length (∼ 71 km) suggests it has a high sediment trapping efficiency for all but the finest particles (Brune, 1953) and that sediment particles that do pass through the lake have a longer transit time than particles that are only transported in rivers.These two characteristics may dampen the magnitude of SSC variability from year to year.
High interannual variability in this natural system, coupled with a relatively short study period, impede the detection of statistically significant trends in other plumes.Further, GrIS mass balance estimates (Shepherd et al., 2012;van den Broeke et al., 2009;Rignot et al., 2008) have shown acceleration in GrIS mass loss beginning in the 1990s.Consequently, the MODIS imagery record only captures the latter half of the response to recent mass balance changes.

Trends in plume area metrics
SSC P >t metrics showed that both areas potentially detrimental to fjord biological productivity (P >50 ) and the core of plumes (P >250 ) have grown larger and/or more persistent.P >50 and P >250 significantly increased for all river plumes except for the Pakitsuup south river (Figs. 9 and 10, respectively).Within a given plume ROI, P >50 increased between 0.08 and 0.25 % yr −1 when normalized by its 13 year mean and P >250 increased between 0.05 and 0.08 % yr −1 when normalized by its 13 year mean.Expressed as the sum of the area of pixels exceeding these thresholds for all images in a melt-season (where 1 pixel = 0.0625 km 2 )P >50 zones increased between 1.4 and 43.2 km 2 yr −1 .P >250 zones increased between 0.9 and 15.5 km 2 yr −1 .

Comparison of fjord SSC to river discharge and ice sheet runoff
Implicit to our initial hypothesis was that plume SSC and river discharge would be positively correlated, as river SSC and discharge often are.However, at the individual plume level, we find plume SSC msm does not directly correlate to the yearly volume of river discharge to the ocean.Total annual river discharge did not closely correspond to SSC msm using the Watson gauge data (Hasholt et al., 2013).For example, moderate yearly discharge (3.5 +1.68 −0.52 km 3 yr −1 ) in 2007 corresponded with the highest SSC msm recorded (216 ± 41 mg L −1 ), while in 2011 a comparable discharge (3.9 +1.76  −0.56 km 3 yr −1 ) produced one of the lowest SSC msm values (127 ± 24 mg L −1 ).The poor correlation between these variables also exists when RACMO2 annual runoff totals were used to extend the analysis in both time and to additional rivers.Of note, the Naujat Kuat plume SSC msm was negatively correlated to RACMO2 runoff (r 2 = −0.6,p = 0.05).However, when data from multiple river plumes were aggregated, we found larger discharge catchments predicted by the RACMO2 model tended to have higher SSC msm values (Fig. 11; r 2 = 0.64, p < 0.001, n = 44).This result agrees with Chu et al. (2012) who found melt (as determined from a positive degree day model) was correlated to fjord SSC for large, spatially aggregated 100 km × 100 km grid cells.

Trends between rivers
During a given year, SSC exhibits a complex spatial relationship between the individual fjord sites.At times, meltseason metrics were similar between individual river plumes (Tables 3 through 5).For example, all river plumes in 2009 displayed SSC msm at or below 2000-2012 mean conditions, while during the warm 2010 year all plumes displayed SSC msm at or above 2000-2012 mean conditions.However, just as commonly, plume conditions exhibited limited spatial coherence (Figs. 12 and 13).All fjords studied with multiple rivers displayed complex SSC responses in certain years, even for plumes with nearby catchments that presumably experienced similar surface melt conditions.During the extreme melt year of 2012 (Nghiem et al., 2012;McGrath et al., 2013), the Watson River plume SSC msm was at or below its 2000 to 2012 mean, while the Umiiviit River plume, with a catchment directly south of the Watson's, experienced a SSC msm 43 ± 17 % higher than its 13 year mean.This asynchronicity is especially notable as the rivers discharge into the same fjord, and thus background oceanographic conditions are very similar.Differences in SSC for the core of the plume appear to drive this asynchronicity.Both the Watson and Umiiviit have P >50 metrics elevated above there 13 year mean, but Watson's P >250 metric was close to mean conditions, while the Umiiviit metric was 83 ± 19 % above.(Hasholt et al., 2013) and SSC msm plotted for comparison, but these data were not included in yearly mean SSC-runoff relationship.Figure 13d shows that the distal plume was above mean values, while pixels near the river mouth were below mean values.
River plumes displayed asynchronicity in other years.During another reportedly high melt year of 2002 (Hanna et al., 2008), the Umiiviit River plume SSC msm was close to its 13 year mean SSC, while the Watson River plume SSC msm was 28 ± 14 % below its 13 year mean SSC.In 2005, the Umiiviit River plume SSC msm was at or above its 13 year mean SSC, while the Sarfartoq River plume SSC msm , with a catchment abutting Umiiviit's GrIS catchment, was 23 ± 15 % below its 13 year mean SSC.
Observed asynchronicity is indicative of the complex sediment dynamics expected in systems controlled by glaciers.Collins (1990) hypothesized that variability in sub-glacial hydrologic networks drove variability in SSC emerging from glaciers.Stott and Grove (2001) recorded transient flushes of increased SSC without increases in discharge for a glaciated (non-GrIS) catchment in Greenland.As mentioned in Sect. 1, river sediment dynamics often show a relationship between discharge and SSC.However, the variable response found may indicate in certain years that ice sheet sediment dynamics, not river sediment dynamics are the dominant control on suspended sediment delivery to the coastal ocean.Fjord-oceanographic processes, such as the strength of estuarine circulation, may influence plume dynamics.However, data limitations and lack of long-term fjord monitoring restrict quantitative characterization of these processes.The plume SSC observed is consistent with river SSC variability seen in Landsat imagery (Hudson et al., 2014a).Thus, we argue that glacial and river processes dominate the observed variability.

Conclusions
We developed a retrieval algorithm capable of using MODIS reflectance to measure fjord SSC around Greenland.The SSC retrieval algorithm utilized the largest, most geographically widespread fjord SSC data set for Greenland.Data were collected over four melt seasons between 2008 and 2012, three melt months, and four fjords with highly variable oceanographic conditions.Most of the geologic terrain eroded by the GrIS that drain into study fjords is similar, consisting of Archean and Paleoproterozoic age gneiss (Escher and Pulvertaft, 1995) and the geologic terrain in southwest Greenland sourcing river plumes is generally similar.As a result we argue that the retrieval algorithm is applicable to a larger area than just the study fjords.
The retrieval utilized both the red-visible (band one, 620 to 670 nm) and near-infrared (band two, 841 to 876 nm) 250 m resolution bands of MODIS.The combination of these bands provided accurate retrievals of SSC concentrations much higher than previous single band retrieval algorithms (e.g., Chu et al., 2009).In addition, a rigorous cloud detection scheme was developed for differentiating between turbid water and clouds, which minimizes cloud contamination of the retrieval (Appendix A).
An analysis using this retrieval on 13 years of MODIS imagery showed that SSC msm has increased for most plumes.Intra-annual and regional variability is large in Arctic regions though, and this makes it difficult to arrive at statistically significant trends over the 13 year record.
Metrics describing the duration and the spatial extent of river plumes responded more pronouncedly to the increased melt of the GrIS.High-concentration plumes that can impair biological productivity expanded in most fjords and one-third of the river plumes expand its high-concentration core zone in time and extent (95 % confidence interval).
While the relationship between the SSC msm of a river plume and its volume of water discharged is not correlated interannually, there appears to be a relationship between yearly modeled runoff totals and SSC msm when rivers are studied as a regional aggregate.We found that larger rivers tended to have higher plume SSC msm .The behavior of fjord sediment plumes over the MODIS record was complex, and variable from catchment to catchment, even in catchments that are proximal to one another.This high spatial and temporal variability suggests that approaches that treat the GrIS as a whole, or even composed of distinct regions, may not be applicable to the study of sediment dynamics.Numerous and frequently non-linear processes operating in the glacio-fluvial-fjord systems preclude a simple relationship between these parameters.Because of this complexity, orbital remote sensing of sediment dynamics offers a unique and important perspective for observing change in the GrIS system and should be continued.

B.Figure 1 .
Figure 1.Conceptual diagram of processes that influence delivery of sediment from ice sheet to fjord, plotted on top of aerial photograph of the Pakitsuup south river and Pakitsuup fjord.Sediment produced subglacially may be stored until removed by water or ice(1).Once exported by the glacier (2), sediment can be transported by rivers(3).While transported by the river, it can potentially be stored or passed through lakes (4), stored on the river's outwash plain(5) until it is discharged into the ocean(6).Sediment stored on the river's outwash plain may be remobilized at a later time by the river(5).Additional sediment can be sourced locally from the landscape, for example through bank erosion.Photograph courtesy of Dorthe Petersen, ASIAQ.

Figure 2 .
Figure 2. Overview map of study fjords, with surface suspended sediment concentration (SSC) sample locations used in retrieval overlaid on 13 year mean SSC values derived from MODIS algorithm.Background images are Landsat 7 ETM+ from (a) 7 July 2001, (b) 20 August 2002, and (c) 3 August 2001.Inset shows location of fjords relative to Greenland Ice Sheet.An asterisk shows the location of the Watson River gauge.

Figure 3 .
Figure 3. Land surface digital elevation model (DEM), ice sheet topographic contours, and glacier catchment area for rivers and fjords discussed in paper.Land DEM and ice sheet contours were derived from the GIMP DEM.Equilibrium line altitudes (ELA) discussed in Sect. 2 also plotted on ice sheet for reference, as are 13 year mean SSC concentration.For (c), Naujat Kuat-Small is delineated with sold pink color.The entire Naujat Kuat catchment combines both solid and diagonal striped pink regions.Latitude and longitude grid is in 0.25 decimal degree intervals.

Figure 4 .
Figure 4. Relationship between MODIS band one and two corrected remote-sensing reflectance (R rsc ) and surface water sample SSC values used in retrieval algorithm.

Figure 5 .
Figure 5. Example of how computational experiment allowed simulated scene suspended sediment concentration (SSC) to vary durring a simulated melt season.For each scene a scene intensity was randomly generated that controls the simulated scene mean SSC (plotted).Three melt seasons are plotted to show a wider range of SSCs allowed.

Figure 6 .
Figure 6.Example of (a) simulated image and (b) cloud masked simulated image for Watson River region of interest (ROI).

Figure 7 .
Figure 7. Percent deviation of annual plume metrics from simulated data set was used to constrain error budget.P >50 values (not plotted) are similar to P >250 values.

Figure 11 .
Figure 11.Relationship between yearly RACMO2 modeled runoff and melt-season mean suspended sediment concentration (SSC msm ) for four river plumes.Watson River observed discharge(Hasholt et al., 2013) and SSC msm plotted for comparison, but these data were not included in yearly mean SSC-runoff relationship.

Figure 12 .
Figure 12.Maps of 2010 melt-season mean suspended sediment concentration (SSC msm ) and its anomaly from its 13 year mean for plume regions of interest.Landsat 7 false color image from 7 July 2001 underlaid for reference.

Figure 13 .
Figure 13.Maps of 2012 melt-season mean suspended sediment concentration (SSC msm ) and its anomaly from its 2000 through 2012 mean for plume regions of interest.Landsat 7 false color image from 7 July 2001 underlaid for reference.

Table 1 .
Various measures of glacier catchments in paper.

Table 3 .
Melt-season mean suspended sediment concentration (SSC msm ) percent of 13 year mean.Underlined values are statistically below the mean, bold values statistically above.

Table 4 .
Percent available fjord pixels in a melt season above 50 mg L −1 (P >50 ) as a percent of statistics' 13 year mean.Underlined values are statistically below the mean, bold values statistically above.

Table 5 .
Percent available fjord pixels in a melt season above 250 mg L −1 (P >250 ) as a percent of statistics' 13 year mean.Underlined values are statistically below the mean, bold values statistically above.