Greenland ice sheet rainfall climatology, extremes and atmospheric river rapids

Greenland rainfall has come into focus as a climate change indicator and from a variety of emerging cryospheric impacts. This study first evaluates rainfall in five state‐of‐the‐art numerical prediction systems (NPSs) (CARRA, ERA5, NHM‐SMAP, RACMO, MAR) using in situ rainfall data from two regions spanning from land onto the ice sheet. The new EU Copernicus Climate Change Service (C3S) Arctic Regional ReAnalysis (CARRA), with a relatively fine (2.5 km) horizontal grid spacing and extensive within‐model‐domain observational initialization, has the lowest average bias and highest explained variance relative to the field data. ERA5 inland wet bias versus CARRA is consistent with the field data and other research and is presumably due to more ERA5 topographic smoothing. A CARRA climatology 1991–2021 has rainfall increasing by more than one‐third for the ice sheet and its peripheral ice masses. CARRA and in situ data illuminate extreme (above 300 mm per day) local rainfall episodes. A detailed examination CARRA data reveals the interplay of mass conservation that splits flow around southern Greenland and condensational buoyancy generation that maintains along‐flow updraft ‘rapids’ 2 km above sea level, which produce rain bands within an atmospheric river interacting with Greenland. CARRA resolves gravity wave oscillations that initiate as a result of buoyancy offshore, which then amplify from terrain‐forced uplift. In a detailed case study, CARRA resolves orographic intensification of rainfall by up to a factor of four, which is consistent with the field data.


| INTRODUCTION
While solid precipitation is the main source of Greenland ice sheet mass input (e.g., Box et al., 2004;van den Broeke et al., 2016), climate warming has caused an increasing rainfall fraction of total precipitation (Niwano et al., 2021). Rain water delivery can accelerate ice sheet flow (Doyle et al., 2015), destabilize tundra snowpacks (Abermann et al., 2019) and alongside surface heating, initiate the snowmelt-albedo feedback .
Rainfall extremes can occur as part of concentrated poleward transport of moisture and heat occurring in 'atmospheric river' (AR) episodes (Neff, 2018). Neff et al. (2014) identified AR episodes as promoting high Greenland ice melt by advection of air masses over the ocean including upstream development over the 2012 record-setting summer North American heatwave. Two July 2012 AR episodes were responsible for record-high ice ablation in Greenland Fausto, van As, et al., 2016;Mattingly et al., 2018). Observations during AR episodes include liquid water clouds near the ice sheet summit (Bennartz et al., 2013), 100 mm rainfall at 1500 m in the accumulation area of the northern ice sheet (Niwano et al., 2015) and rainfall above 2850 m elevation . The occurrence of moist ARs reaching Greenland has been increasing (Mattingly et al., 2016), driven by more frequently occurring highly amplified jet-stream patterns (Francis & Skific, 2015).
Here, we present a rainfall climatology for the Greenland ice sheet and its peripheral ice masses through application of the new and fine-resolution Copernicus Arctic Regional Reanalysis (CARRA) dataset. We compare rainfall from CARRA and four other state-of-the-art numerical weather prediction systems with a 4-year set of independent in situ rainfall measurements, mainly from the ice sheet. A CARRA rain and snowfall climatology for 1991 to 2021 is presented. CARRA rainfall and rain fraction of total precipitation maps are presented for the extreme high melt year 2012 and for the 1991-2021 period, including a difference mapping with ERA5 to illustrate added value from CARRA. We examine the atmospheric dynamics and thermodynamics and surface energy transfers during an extreme rainfall event in September, 2017 using the CARRA data and in situ ice sheet meteorological/glaciological observations. 2 | DATA AND METHODS 2.1 | Rainfall observations 2.1.1 | On-ice rainfall data Since 2016, the Geological Survey of Denmark and Greenland (GEUS) has maintained precipitation measurements on glaciated areas of Greenland. These include the Qagssimiut ice lobe of the southern ice sheet (QAS) and the K-transect east of Kangerlussuaq (KAN) (Figure 1, Table 1). The rainfall data obtained from the precipitation measurements are not assimilated by any of the numerical weather prediction systems evaluated in this study.
The rain gauge employed in this work is the Texas Electronics TR-525I, having a 15.2-cm (6-in.) diameter orifice. The instrument is sold by the Onset Corporation as the HOBO RG3-M rain gauge. This gauge has been mounted on the PROMICE automatic weather station (AWS) mast 0.3 m below the level of wind speed measurements ( Figure 1). The gauge is unheated and unshielded. The justification of its use is to focus on rainfall only, in which catch efficiency is much higher than for snow (Førland et al., 1996), and to focus on high-magnitude rainfall events in which the signal-to-error ratio is highest.

| Narsaq rainfall data
A meteorological station at the Narsaq heliport has been operated since 1996 by Asiaq, the Greenland Survey. We use data from July 2016 onwards, overlapping with the time period of the GEUS on-ice rainfall measurements. For this period, the sensor in operation is a Pluvio2 precipitation gauge with a Tretyakov shield.

| Rainfall gauge errors
Precipitation gauge undercatch errors and correction efforts are widely documented (e.g., Førland et al., 1996;Goodison et al., 1998;Sevruk et al., 2009). In Greenland rainfall studies, various undercatch corrections have been applied to data from land-based precipitation gauges operated by the Danish Meteorological Institute (Yang et al., 1999;Mernild et al., 2015;Koyama & Stroeve, 2019;Huai et al., 2021;Niwano et al., 2021). In addition to precipitation undercatch, which results from distortion of the windfield around the rain gauge by the measurement platform and the gauge itself (Sevruk et al., 1991), error sources include uncertainty about precipitation liquid versus solid phase. Here, we simplify the problem by considering only rainfall by excluding cases with hourly air temperatures under 0 C. The gauges may, nonetheless, accumulate snow, which can lead to errors from delayed snow melt. Therefore, we utilize the data only after mid-June and until mid-September under prolonged above-freezing temperatures, when any such snow would have ablated. Incidentally, we found no evidence that the TR-525I tipping mechanism is triggered under high wind speed conditions from vibration effects.

Cohesion and evaporation losses
Precipitation undercatch also results from water droplet adhesion above the measurement device, that is, 'wetting loss' (Yang et al., 1999), which includes evaporation of the droplets before the mass flux can be registered. In the case of the TR-525I data, rain water that is insufficient to drain through to reach the tipping mechanism counts as trace precipitation. The evaporation and wetting losses increase in importance as the total precipitation decreases (Yang et al. 1999), that is, at the relatively dry locations KAN_L and KAN_B.

Rain gauge tilt
The PROMICE AWS, and thus the TR-525I instrument, has variable tilt as the surface ablates from beneath the AWS by up to 6 m each year. When not level, the TR-525I tipping function can be reduced, leading to a dry bias. Based on this consideration, we limit our comparison to those measurements having a tilt not exceeding 5.6 , a maximum value obtained during the 14-15 September 2017 extreme rain event studied here. Under this condition, we found no obvious tilt dependence in the bias between the field data and the numerical weather prediction data. For the relatively small sample of QAS_U data (because the record started later, in 2018, and under lower temperatures), a large tilt ($18 ) prevents the use of most QAS_U rain data.

| Undercatch corrections
On-ice TR-525I undercatch correction Given that the TR-525I gauge has a cylindrical shape, we apply an undercatch correction for an unshielded Hellmann-type precipitation gauge under liquid-only precipitation conditions with a corresponding catch efficiency (k) correction after Yang et al. (1999): with wind speed (U) in m s À1 at gauge height. We use the PROMICE AWS hourly average wind speed observations interpolated linearly in time to each TR-525I datum, which are logged as tipping events with seconds precision. While the correction was developed using daily averages in the U ≤ 6.5 m s À1 range, we assume that the equation is also valid for hourly wind speeds. An extrapolated correction is used when the wind speed exceeds 6.5 m s À1 , occurring in 7%-12% of the cases. While Yang et al. (1999) concluded that wetting and evaporation losses correspond with a daily measurement error of less than 0.1 mm, which is below the 0.2 mm resolution of the TR-525I gauge measurement, we coarsely account for trace precipitation by setting the lowest possible rainfall correction to 1.01 (+1%). The applied undercatch correction averaged +10%.
As with unattended precipitation measurements anywhere, substantial uncertainty remains with the corrected values: for example, the undercatch is not entirely explained by temperature and wind speed (Køltzow et al., 2020). Indeed, the explained variance of the correction of the Helmann gauge relative to the Double Fence Intercomparison Reference (DFIR) is only 48% and for a moderate sample size (N = 223). Yet, for liquid precipitation, Yang et al. (1999) computed the undercatch to be under 20% for Greenland locations, which was also confirmed by Niwano et al. (2021).

Narsaq undercatch corrections
The Tretyakov shield liquid precipitation correction scheme [Equation 7 in Yang et al. (1995)] is applied to the Narsaq data. During the extreme rain event between 18:00 UTC on 13 September 2017 to 12:00 UTC on 15 September, Narsaq wind measurements were under 5 m s À1 , producing a moderate (under +22%) undercatch correction.

| Rainfall from numerical prediction systems
A wide variety of precipitation estimates from numerical prediction systems (NPSs) ( Table 2) are here compared with the independent field data.

| HARMONIE-AROME CARRA
The Copernicus Climate Change Service (C3S) Arctic Regional ReAnalysis (CARRA) reanalysis system (Yang et al., 2020) applies HARMONIE-AROME, a nonhydrostatic, convection-resolving weather forecast model (Bengtsson et al., 2017), to assimilate an extensive collection of observations within the Greenland domain at 2.5 km horizontal grid spacing and 65 vertical levels. The CARRA reanalysis is run with a 3-h assimilation interval.
T A B L E 1 Greenland locations with liquid precipitation gauge data used in this study.  Huai et al. (2022) database were updated with ArcticDEM data (Porter et al., 2018), which were further improved manually (pers. comm. Palmason, 2018). Glaciated area including a sub-grid fractional mask corresponds to 1,804,032 km 2 for the ice sheet and 42,944 km 2 for ice masses peripheral to the ice sheet. Peripheral glaciated areas were obtained using watershed segmentation ('scikit-image'). Model physics improvements include the assimilation of the Box et al. (2017) and Wehrlé et al. (2021) satellite-derived glacier albedo data in combination with a snow albedo increase to 0.85 in cases of modelled snowfall (Nielsen, 2019).

| Liquid phase of precipitation
Here, external to the HARMONIE-AROME simulations, the CARRA precipitation phase is estimated from an approach (Rohrer, 1989;Hock & Holmgren, 2005) assumed to be accurate given that it lies in the middle of seven phase identification approaches (Huai et al., 2021) with total precipitation with the transition air temperature from snow to rain set to vary linearly 1 C above and below 1.5 C.

| ERA5
ERA5, the fifth-generation ECMWF global reanalysis [Copernicus Climate Change Service Climate Data Store (CDS); Hersbach et al., 2020] is based on the Integrated Forecasting System (IFS) Cy41r2, operational at ECMWF in 2016. ERA5 therefore benefits from a decade of developments in model physics, core dynamics and data assimilation compared to its predecessors, for example, ERA-Interim and ERA-40. ERA5 includes finer (31 km) horizontal grid and temporal (hourly) output compared to ERA-Interim. By assimilating historical observations into a numerical weather prediction model, ERA5 provides more spatially and temporally consistent estimates of atmospheric and surface variables. ERA5 improves on ERA-Interim for precipitation estimates (Hersbach et al., 2020). For 2-m air temperature and surface energy balance components on the west Greenland ice sheet, Huai et al. (2020) found that ERA5 better represents the observations than ERA-Interim, although the improvement was only statistically significant for albedo, consistent with Delhasse et al. (2020).

| NHM-SMAP
The polar Non-Hydrostatic atmospheric regional climate Model with the Snow Metamorphism and Albedo Process (NHM-SMAP) (Niwano et al., 2018) was applied to Greenland at a grid spacing of 5.0 km and produced 1-h output. NHM-SMAP was forced at its boundaries every 6 h by the Japanese 55-year reanalysis JRA-55 dataset (Kobayashi et al., 2015). Daily JRA-55 atmospheric profile initialization prevents large deviations between JRA-55 and NHM-SMAP atmospheric fields. A daily simulation is performed starting at 18:00 UTC the previous day, thus including a 6-h spinup. NHM-SMAP initialization uses JRA-55 sea ice extent and sea surface temperature data.

| MAR
The regional climate polar model MAR (version 3.11.5) run at a grid spacing of 6 km is used here to downscale the reanalysis ERA5, which forces MAR at its lateral boundaries as well as at the surface of its ocean (sea surface temperature and sea ice cover) every 6 h from 1950 to 2021. With respect to the version 3.9 of MAR, intensively evaluated over the Greenland ice sheet , the main improvements in MARv3.11 include updates in the cloud scheme and in the bare ice albedo with the aim of getting better fit with observations. The cloud scheme of MAR was originally based on Meyers et al. (1992), with some improvements on the basis of the work done in the Mixed-Phase Arctic Cloud Experiment (Fridlind et al., 2007). We refer to Hofer et al. (2019) for an evaluation of this scheme over Greenland where rainfall is automatically converted to snowfall for near-surface temperature below À1 C. The MAR output is downscaled to a 1-km grid .

| RACMO
The Regional Atmospheric Climate Model (RACMO) was developed by the Royal Netherlands Meteorological Institute (KNMI) (van Meijgaard et al., 2008). The version used here (RACMO2.3p2) has 40 vertical atmospheric levels and 5.5 km horizontal grid spacing, and is forced by a combination of climate reanalyses including ERA-40 (1958ERA-40 ( -1978 and ERA-Interim (1979-1989) on a 6-hourly basis and by ERA5 (1990-2021) on a 3-hourly basis (Noël et al., 2019). A polar version was developed to simulate the surface mass balance over the ice sheets of Greenland and Antarctica . The cloud waterto-snowfall conversion coefficient remains constant for liquid (above 0 C) and mixed-phase clouds (À23 C to 0 C), whereas it decreases with temperature for ice clouds (less than À23 C), resulting in changes in snowfall production (Van Wessem et al., 2014). Precipitation computation described in Noël et al. (2015) has liquid/solid phase determined by a cloud phase scheme in the 5.5-km product. In the 1-km product, precipitation is assumed solid when the air temperature is <À1 C. The precipitation phase is further corrected as a function of elevation through statistical downscaling. The scheme allows for ice supersaturation, resulting in improved relative contributions of rainfall and snowfall to total precipitation .
Here, we use two RACMO2.3p2-based daily rainfall products: one at the native 5.5 km resolution, and a rain product that is statistically downscaled from the outputs of RACMO2.3p2 to a 1-km grid, which much better resolves, for example, the low-lying marginal glaciers . To improve the partitioning between solid and liquid precipitation, the statistical downscaling technique uses daily elevation gradients of the snowfall fraction (SF frac ), which strongly correlates with altitude on both icecovered and surrounding tundra regions. As a result, a corrected daily snowfall (SF) and rainfall (RA) product at 1 km is estimated for the ice sheet and tundra regions as where PR is the daily total precipitation (including snowfall and rainfall) from RACMO2.3p2 at 5.5 km and bi-linearly interpolated to the 1-km grid. SF frac is the daily snowfall fraction, that is, SF/PR statistically downscaled to 1 km resolution (Huai et al., 2022).

| Uncertainty of precipitation estimates in numerical prediction systems
The accuracy of precipitation in NPSs (Table 2) depends on the description of physical processes and the NPS initialization. Different NPSs therefore often agree on large-scale precipitation patterns but differ substantially on precipitation amount and phase (liquid versus solid) on local scales (Barrett et al., 2020;Edel et al., 2020;Køltzow et al., 2019). However, even if field instrumentation and the NPSs were perfectly accurate, the models and observations will not produce the same precipitation because they represent different spatial scales: for example, a substantial subgrid variability may be observed for point measurements compared to grid-averaged precipitation from the NPSs. Meanwhile, the representation of precipitation in NPSs is, for example, sensitive to the horizontal resolution and proximity of complex topography. Point measurements in complex topography can be strongly influenced by orographic wind field distortion (Bromwich et al., 1998;Yang et al., 1999). Both the observational and NPS rainfall determination can be sensitive to air temperature, which is used to discriminate precipitation phase (Huai et al., 2021).

| Verification of rainfall from numerical prediction systems
Agreement between the independent rainfall observations and the NPSs is measured by the daily temporal correlation and the ratio of accumulated precipitation totals (Table 3). The CARRA rainfall simulations agree best with the field data from the relatively wet QAS sites, both in the ratio of totals and according to the temporal correlation. Correlations for the wet QAS sites tend to be high (0.7-0.9), while at the relatively dry KAN sites the correlations vary substantially because of the small data sample and do not produce robust statistics. A spatial difference pattern is evident for all the NPSs: that is, a 10-60% dry difference with the observations at the lowest elevation QAS sites (QAS_L-16 and QAS_L) ( Table 3). The dry difference decreases towards higher elevations (QAS_M) and changes (for MAR and RACMO) to a wet difference at the highest elevation site (QAS_U), while CARRA and NHM-SMAP still have a slight (6-10%) underestimation. There is no average ERA5 difference at QAS_U.
At Narsaq, situated on land $30 km from the ice sheet, using the likely more accurate (more level and wind-shielded) Pluvio2 measurement, we find an inconsistent bias among the NPSs: dry difference for NHM-SMAP and CARRA, no difference for ERA5 and RACMO 1 or 5.5 km and wet difference for MAR (Table 3).

| Discussion of differences
Although verification at just two regions may not represent the simulated accuracy very well, the QAS spatial difference patterns for the various NPSs (Table 3) are consistent with findings for AROME and partly for the ECMWF IFS on the Norwegian coast, that is, too little rainfall and total precipitation at the coast and too much inland (Køltzow et al., 2019).
In addition to the finer grid spacing of the CARRA data compared with the other NPSs, higher agreement with observations for CARRA is likely due to the fact that its simulations are initialized with local observations (not precipitation) within the model domain. ERA5 has some within-domain data assimilation [radiosonde and satellite-derived upper-air humidity and temperature profiles, aircraft-based atmospheric observations (temperature, wind components and specific humidity)], but substantially less than CARRA. The more explicitly modelled interactions in the convection-permitting model physics may also contribute to higher CARRA skill. At the current grid resolution, however, CARRA's skill for T A B L E 3 Comparison of daily rainfall observations with model data. rainfall sensitivity does not appear to depend on the selection of hydrostatic or non-hydrostatic physics. Because the frequency of CARRA, ERA5 and the observed air temperatures above 1.5 C, in which all precipitation should be rain, is equivalent, that is, 52% for the observations, 52% for CARRA and 50% for ERA5, the assessed differences in Table 3 appear to be attributable mainly to simulated precipitation magnitude, not temperature-based phasing.

Model
The difference between nearest neighbour collocation and bilinear interpolation increases (not shown), in general, for higher model grid spacing, suggesting that finer resolution has value.

| CARRA Greenland rainfall climatology
Maximum 31-year average (1991-2021) rainfall in the CARRA data (1207 mm y À1 ) occurs across the southeastern tip of Greenland glaciated area at 60.252 N, 43.271 W at 298 m elevation ( Figure 2). While CAR-RA's annual peak rainfall occurs typically along the southeastern ice sheet, given this study's focus on the south-southwestern ice sheet Qagssimiut lobe, we specify the annual peak rainfall there as well in Figure 2. Along the southeastern ice sheet, CARRA's average annual rainfall above 1000 mm occurs only south of Kap Tobin (Scoreseby Sund, 70.4 N), north of which CARRA ice sheet rainfall is much lower, that is, under (30 mm y À1 ).
The anomalously warm year 2012 has extremely high rainfall amounts in CARRA above 600 mm y À1 along the western ice sheet north of 75 N ( Figure 3). Along western Greenland, low (under 60 mm y À1 ) annual rainfall occurs in CARRA data north of the Thule Air Base (76.5 N). In another anomalously warm year, 2010, rainfall was concentrated instead along the southeast coast, highlighting strongly varying yearto-year spatial patterns.

| Rain fraction of total precipitation
In the high rainfall year 2012 (Figure 3), CARRA rain fraction of total precipitation (f rain ) peaked at 0.59 on the Qagssimiut lobe of the ice sheet where f rain is naturally high given high air temperatures (Figure 4). Along a narrow area (typically under 25 km) for much of the western ice sheet, f rain exceeded 0.3, including the northwesternmost ice sheet (Figure 4). f rain can also be high where total precipitation is extremely low, for example, in sheltered sites in the northwest. The record melt year 2012 had area-averaged f rain that was 1.6 times the 1991-2021 average. In another extremely high temperature year, 2010, f rain was 1.4 times the 1991-2021 average. However, the high-melt year 2019 (Tedesco & Fettweis, 2020) has an average f rain , demonstrating no simple f rain dependence on temperature. See also Section 3.4. On average, from 1991 to 2021, just 3.7% of the Greenland ice area has a rain fraction above 0.1.

| Ice sheet scale rainfall rates
In comparison with past work, the average CARRA rainfall during 1991-2021 is 25.5 Gt y À1 with an inter-annual standard deviation of 7.6 Gt y À1 (or 30%). For comparison, Ettema et al. (2009) found a higher total rainfall in a previous version of RACMO (46 Gt y À1 ) including ice sheet peripheral ice caps and glaciers, although for a different period . Van den Broeke et al. (2016) report that annual average rainfall and rain fraction over the ice sheet from 1991 to 2015 estimated by RACMO2.3p1 are, respectively, 28 Gt y À1 and 3.9%. Over the same period, NHM-SMAP had similar values (28 Gt year À1 and 3.2%). Fettweis et al. (2013) found that annual rainfall over the ice sheet from 1980 to 1999 from MAR v3.5.2 was 28 ± 5 F I G U R E 2 The 1991-2021 average annual CARRA rainfall map with logarithmic colour scale. Ice sheet areas with average rainfall below 1 mm are coloured white. Elevation isolines from 500 to 3000 m appear as dashed lines.
F I G U R E 3 The 2012 total rainfall according to CARRA data. The circle indicates the location of the maximum rainfall occurring on the Qagssimiut lobe, where the QAS PROMICE sites are located. The elevation is indicated as metres above sea level (m a.s.l.).

F I G U R E 4
The 2012 average rain fraction of total precipitation from the CARRA dataset. Areas with total precipitation below 0.005 mm y -1 are coloured white.
Gt y À1 . Over that period, NHM-SMAP was in agreement, with 21 Gt y À1 . Intercomparison at the whole island scale, however, is complicated by the fact that where the land/ ice/ocean masks differ most, around the ice sheet periphery, is where the precipitation fluxes are typically highest.

| Ice sheet and peripheral ice mass rainfall trends
At the ice sheet scale, over the 1991-2021 (31 year) CARRA period, both rainfall and rainfall fraction of total precipitation have increasing trends ( Figure 5, Table 4), which is consistent with earlier findings from RACMO (van den Broeke et al., 2016) and NHM-SMAP data (Niwano et al. 2021). Peak rainfall years are 2010 and 2012, which had high ice sheet-wide June through September air temperatures (and melt). Year 2021 had the highest CARRA rainfall and was examined using field data and ERA5 reanalysis . For the 31 years of CARRA data, snowfall has no trend, consistent with RACMO data (1991-2015) (van den Broeke et al., 2016) and a combination of RACMO and MAR data for 2000 to 2019 . For peripheral ice masses, the 1991-2021 average (± standard deviation) of CARRA rainfall is 1.8 ± 0.5 Gt y À1 , with a significant (1p = 0.874) 31-year increase by 33% (Table 4). A statistical downscaling of RACMO to 1 km resolution by Huai et al. (2022) gave an average total annual rainfall on Greenland ice during 1958-2020 to be 28.6 ± 6.1 Gt y À1 , 2.6 times higher than CARRA (11.4 ± 1.4 Gt y À1 ) over peripheral ice masses. The difference is attributable to the Huai et al. (2022) data using a peripheral glacier mask with 81,400 km 2 ice area (Noël et al., 2017) compared to 42,944 km 2 obtained for CARRA using a different ice masking criterion (see Section 2.3.1).
For snowfall, CARRA simulates an insignificant positive temporal trend (Table 4). For the 1991-2015 period, F I G U R E 5 Snowfall, rainfall and rainfall fraction of total precipitation from CARRA for the ice sheet and peripheral glaciers. Numbers above bars indicate the ice sheet (including peripheral glaciers) mass flux in Gigatonnes. RACMO's total precipitation was 712 ± 64 Gt y À1 (van den Broeke et al., 2016), differing by more than one standard deviation with CARRA (855 ± 69 Gt y À1 ) albeit for a RACMO period that is 6 years shorter. For rainfall, the 1991-2015 RACMO at 28 ± 9 Gt y À1 differs insignificantly with CARRA (25 ± 8 Gt y À1 ). CARRA 1991-2021 rainfall trends above 66% confidence ( Figure 6) occur around the ice sheet, especially for the northwest as shown by Niwano et al. (2021) for NHM-SMAP data. The low-confidence areas account for a net additional 1 Gt y À1 of rain and increase the percent change from +29% rainfall increase (confidence above 66%) to +34%. Only in low-confidence areas do drying trends occur in the CARRA data, signifying how defining rainfall trends is challenged by some areas/years having zero values. The extreme drying trend areas (red areas in Figure 6) have confidence below 0.66.

| Extreme daily rainfall cases
In the 31-year (1991-2021) period, CARRA simulates maximum local ice sheet rainfall up to 448 mm day À1 . Eighteen days had maximum daily local rainfall above 300 mm (Table 5) and occurred not only in summer but also in winter months (see Oltmanns et al., 2019). The average elevation of the extremes in daily rainfall is 419 m with a standard deviation of 224 m. We highlight the challenge to resolve local terrain by coarser model grids, especially ERA5 at $31 km, where terrain-smoothing yields imprecise elevation of local extremes through the smoothing of complex terrain. The high-temperature year 2010 had most (5 of 18) rainfall days above 300 mm, while the other anomalous melt years 2012 Fausto, van As, et al., 2016) and 2019 (Tedesco & Fettweis, 2020) did not register daily extremes over 300 mm. Daily rainfall above 300 mm in CARRA data is purely a sub-Arctic phenomenon given that daily local extremes all occur south of 67 N. Further, the extreme rainfall locations were concentrated only across the southeast and the most southern glaciated areas (Figure 7, Table 5). The 323-mm south-southwestern daily rainfall extreme on 14 September 2017 is examined further in this study given the availability of in situ rainfall data on that date ( Figure 1, Table 1).
A 1991-2021 CARRA maximum daily rainfall integrated over Greenland glaciated area (4.7 Gt day À1 ) occurred on 14 August 2021 (Table 6). On this date, an AR produced rainfall up to the highest areas of the ice sheet. Rainfall was witnessed at the Summit station at 3250 m elevation. Satellite passive microwave recorded an extensive ice sheet wet snow area mainly caused by surface energy budget surplus from turbulent energy fluxes and downward net longwave irradiance .

| Comparison of CARRA and ERA5
during rainfall extremes CARRA has extreme rainfall concentrated more coastward than ERA5 (Figure 8) including other dates (not shown). The difference pattern, with ERA5 wetter inland (Figure 4c), is confirmed by the field data having ERA5 wet difference increasing with elevation sites (Table 3). Nevertheless, CARRA has more inland trace precipitation. On this date (and others, not shown), CARRA's peak localized rainfall exceeds that of ERA5, again suggesting that the finer CARRA resolution is able to resolve higher peaks in addition to sharper horizontal gradients.

| Large-scale atmospheric circulation
For the 14-15 September 2017 episode, ERA5 data are used to illustrate the AR feature in the precipitable water vapour (PWV), where values above 36 mm are found (Figure 9a). At 850 hPa, winds and air temperatures have values exceeding 25 m s À1 and 10 C, respectively (Figure 9b). Indicative of the conversion of moisture to precipitation, the PWV values drop once the AR encountered Greenland (Figure 9a). At this time over the ice F I G U R E 6 Trend in CARRA rainfall. Stippled areas have trend confidence (1-p) below 66%.
sheet, temperatures well above melting (3 to 6 C) are simulated by ERA5 (Figure 9b). In the following, for a more detailed investigation, we turn to CARRA data.

| CARRA meteorology of an atmospheric river reaching Greenland
Mesoscale atmospheric dynamics CARRA data resolve an atmospheric updraft jet 40-140 km offshore in an arc around southern Greenland ( Figure 10). The intense wind pattern driven by a lowpressure system (low 500 hPa heights) to the west over Labrador, Canada, and by the high pressure to the south southeast of Greenland (Figure 9a) contrasts somewhat with tip jets that arise from low pressure east of Greenland (Doyle & Shapiro, 1999;Outten, 2008).
Downstream of the updraft jet, the 925-hPa ($1000 m altitude) winds accelerate from $20 to $30 m s À1 , decelerate upon reaching the southern ice sheet Qagssimiut lobe and continue at high velocity around the southwest coast of Greenland ( Figure 10). This terrain-induced lowlevel blocking induces turning of the winds akin to the barrier wind effect (Harden et al., 2011). The barrier effect becomes less pronounced as the height above ground increases. The strongest updrafts/ downdrafts occur over Greenland mountains and fjords and are examined later.
The CARRA data include linear along-flow updraft regions, which we term 'rapids' (Figure 11). The rapids appear from 750 to 600 hPa in the featured case and at other times on 14 September 2017. Some rapids begin upstream of the offshore updraft jet and also appear $50 km offshore, west of the island, suggesting they are a more generally feature of the AR dynamics. The filaments are 100-200 km in length, 5 to 15 km in width and 3 km deep, flowing 2 km above sea level. At 15:00 UTC and only at the 700 hPa level do the filaments become interrupted along flow (Figure 11).
Atmospheric flow is channelled onto the ice sheet through a $72 km wide (east-west) topographic gap where the coastal mountains opposing the onshore flow T A B L E 5 Ranking of extreme local CARRA daily rainfall cases above 300 mm in from 1991 to 2021, including the ice sheet total rainfall mass flux. have a relatively low maximum height of 330 m. On either side of the terrain gap, the mountain heights are respectively 320 m (470 m) higher to the west (east) (Figure 11). Downstream of the gap, over essentially all of the Qagssimiut ice lobe, vertical motion is upward ( Figure 11) and extends from near the surface to 500 hPa or rouhgly 6 km above ground.

Date
We can now describe mechanisms supporting the updraft jet and rapids. Besides the mass conservation that produces the flow splitting and acceleration around southern Greenland (Figure 10), buoyancy is generated aloft by condensation (diabatic heating). The rapids persist for over 24 hours (not shown).
Starting upstream from the uplift jet, moist (pseudoadiabatic) isentropes incline upward along flow as buoyant (diabatic) uplift is simulated (Figure 12). The updrafting is maintained along the rapid, suggesting continual condensational heating. CARRA data include T A B L E 6 Ranking of Greenland ice sheet daily rainfall totals above 2.5 Gt per day in CARRA data from 1991 to 2021, together with measured amounts of maximum local rainfall and their locations. areas of relative humidity $1% above 100% (supersaturation) (not shown). Confirmation of alongrapid precipitation is confirmed by the presence of rainfall bands underneath the rapids (Figure 13). Vertically extended moist adiabatic isentropes appear in the areas of strongest updraft, that is, from 220 to 260 km on Figure 12. Maximum upward vertical velocities occur from the forced topographic uplift from coastal mountains. The updraft rapid becomes disturbed by downward-propagating gravity waves in the lee of coastal mountains. Over the ice sheet, the moist isentropes incline more steeply than the dry isentropes, indicating condensational buoyancy generation that adds to the terrain-forced uplift. The observed rainfall rates increase with elevation, confirming that the combined forced-topographic and buoyant lifting was generating precipitation. Inclined dry isentropes are coincident with vertically propagating 'gravity waves', which indicate terrain-forcing of vertical motion in a hypothetical state of stable stratification. Again, the stratification gains instability from condensational buoyancy generation. The gravity waves (a.k.a. mountain waves) (Durran, 1990;Menchaca & Durran, 2017;Geldenhuys et al., 2021) propagate vertically more than 3 times the height of the ice sheet and coastal terrain. Once the ice sheet topographic rise stops at $370 km on Figure 12, a downdraft area (and downward gravity wave) appears. Further downstream, the wave oscillates presumably from mass conservation and gravity wave reflection against the stabilitystratified tropopause (above $250 hPa). Yet, further  downstream, where ice sheet topography begins to rise again from 2150 to 2450 m (at $470 km on Figure 12), updrafts appear once again followed by downdraft, further illustrating the gravity wave oscillations. Moist isentropes again align vertically, at 420 km on Figure 12, indicating further buoyancy generation.

Date
F I G U R E 1 2 The 550-km long and 11-km high CARRA vertical wind speed and thermal adiabatic profile denoted in Figure 11. The thermal profile is represented by potential temperature (adiabatic) and irreversible moist-adiabatic (pseudo-adiabatic) isotherms. The CARRA land/ice mask is used to colour ice areas blue and land areas brown.

Mesoscale precipitation patterns
The flow forcing moisture onto the ice sheet produced 3hourly rainfall rates up to 72 mm, peaking at 750 m elevation over the Qagssimiut lobe of the southern ice sheet ( Figure 13). Rainfall intensification is widespread over the Qagssimiut lobe and coincides with the large-scale uplift (see also Figures 11 and 12). Meridional rain bands coincide with the along-flow updraft rapids in the CARRA data. Rainfall is shifted a few kilometres downstream of the updraft jet, presumably by advection. Reduced rainfall and downdraft areas (grey isolines in Figure 13) appear only where the surface slope decreases over the ice sheet, consistent with orographic precipitation theory (Smith & Barstad, 2004).

| Extreme rainfall during the 14 September 2017 atmospheric river
A larger-scale view of the 14 September 2017 extreme rainfall appears in (Figure 14). For this date, CARRA rainfall on the ice sheet totalled 4.4 Gt, the second highest amount in the 1991-2021 period (Table 6). CARRA simulates peak daily rainfall of 323 mm on the Qagssimiut ice lobe, some 30 km west and north of the in situ measurement sites (yellow stars in Figure 11). Rain bands appearing in the 3-hourly data contribute to the formation of meridional (north-south) precipitation patterns visible over the ocean (Figure 14). Local rainfall peaks correspond with terrain opposing air flow. An intensification of rainfall appears upstream of the tip of southern Greenland (see also the updrafts in Figures 11 and 12). The offshore intensification of rainfall drops slightly towards the coast and rises again as the F I G U R E 1 4 Detailed 2.5-km grid CARRA rainfall (colour shaded areas) and 10-m wind streamlines (coloured by temperature) patterns for 14 September 2017 during the second largest Greenland rain event between 1991 and 2021. Surface elevation isolines from 500 to 3000 m appear as dashed lines. air is forced upslope and inland. During this extreme rainfall episode, a north-south transect through the CARRA grid has rainfall at the lowest elevations over Greenland (13 m above sea level) of 71 mm day À1 , increasing to 278 mm day À1 , suggesting that orographic lifting increased the rainfall rate by up to a factor of four.
The on-ice precipitation gauge daily totals were 194 mm (171 mm uncorrected) at the QAS_M site (621 m above sea level) and 156 mm (137 mm uncorrected) at QAS_L (271 m above sea level). Peak hourly rainfall at QAS_M was 21 mm (18 mm uncorrected) from 16 to 17 UTC, and an hour earlier at QAS_L, 17 mm (15 mm uncorrected) was recorded. At 25 m elevation, 35 km east-southeast of the QAS sites in Narsaq, the peak hourly rainfall rate (8.9 mm) was under half the peak hourly rainfall at QAS_L and QAS_M, consistent with the simulated rainfall intensification due to orographic lifting.

| Extreme rainfall in numerical prediction systems and in situ data
Comparing the field measurements with the simulated rainfall, for the $36-h event, we find a 10%-50% dry difference for ERA5 and MAR and a 30%-70% wet difference for NHM-SMAP, RACMO and CARRA (Table 7). Yet, any disagreement with observations for a single event can come from the high spatial variability in the precipitation intensity given, for instance, the narrow 5 to 8 km width of the uplift rapids and the relatively few observation sites available. Consequently, it is somewhat arbitrary how well observations and models agree during such events with complex wind fields.

| Extreme rainfall impact on ice and snow ablation
In the following, we examine the impact of rain on ice ablation rates on the Qagssimiut lobe of the southern ice sheet where the rain gauge observations document the 14-15 September 2017 rainfall event (See Figures 9-13) and another extreme rainfall event on 19 July 2018.
To calculate the heat supply by rainfall at the surface explicitly, rain temperature is taken as the wet-bulb air temperature (T w ), which is assumed to react quickly to the surrounding environment (Anderson et al., 1998). Under high relative humidity conditions that occur during rainfall, T w becomes equivalent to the ambient air temperature, explaining the equality between latent and sensible energy fluxes ( Figure 15). The energy flux calculations for these sites are described in Fausto et al. (2021).
During the featured event, hourly average winds measured by the AWSs exceeded 8 m s À1 at the onset of the extreme rainfall, with air temperatures rising from below 0 C to above 5 C at the highest elevation (lowest temperature) site QAS_U (Figure 15, top row). The combined effect of wind, air temperature and humidity produce high (above 300 W m À2 ) peak net turbulent fluxes at QAS_M (above 150 W m À2 ) and QAS_U (Figure 15, second row). Net longwave irradiance, mostly an energy sink for the ice sheet, was a modest energy source (under 60 W m À2 ). Hourly rainfall follows a pulse shape, peaking above 18 mm h À1 (Figure 15, third row). In Figure 12 (bottom row), the 'A' symbol indicates how, prior to T A B L E 7 Comparison of observations versus models for the 14-15 September 2017 extreme rainfall event. rainfall at QAS_M, we observe a 5-cm convergence of recorded surface height by sonic ranger and independent hydraulic ablation sensor data in a process of the destruction of porous (low density) surface ice 'weathered crust' layer. Surface energy totalling indicates that over the event, rain adds 17% to ice ablation at QAS_M and 14% at QAS_U and improves the agreement with observed ice ablation.
In another rain episode on 19 July 2018, at QAS_U, the calculated rain addition to ablation was +16% of 4.3 cm total ablation over the 24-h period. The magnitude of rain-induced ice ablation is consistent with earlier work Fausto, van As, et al., 2016), and is here confirmed to be realistic by producing a closer match with independently observed ice ablation than the surface energy balance model without rain heat flux ( Figure 15, bottom row). The relatively minor direct impact of rain on ice ablation is consistent with  and earlier studies (Garvelmann et al., 2014;Niwano et al., 2015;Würzer et al., 2016) who find that turbulent sensible and latent heat transfers dominate the surface energy budget. Indirect effects can include snowpack heating through percolation and refreeze and a subsequent melt-albedo feedback initiated by heat and rainfall . We further note that net longwave irradiance is a larger energy source in the cases examined here.

| Summary of this study
In this study, we pursued a comprehensive evaluation of rainfall over the Greenland ice sheet. The initial step involved introducing a new dataset of rainfall (a) (b) F I G U R E 1 5 PROMICE QAS_M and QAS_U station surface energy budget (SEB) calculations after Fausto et al. (2021), air temperature ice ablation during a heavy rainfall event on 14-15 September 2017. Ice ablation (shown as negative values) is recorded either by an acoustic surface height sensor or a hydraulic barometer. With no QAS_U rain measurements available during the event, here QAS_L rainfall is plotted in the right column for convenience. The acoustic surface height (SR-50) on stakes data from QAS_U or QAS_L failed before the event.
measurements obtained from two regions in Greenland.
The observation regions span from terrestrial environments up to an elevation of 893 m above sea level on the ice sheet. Accuracy of the field data was maximized by employing a method for undercatch correction and the exclusion of potentially erroneous data. Subsequently, the rainfall observations were compared with simulations from five numerical weather prediction systems (NPSs) by evaluating temporal correlations and averaged wet or dry differences, while considering the elevation dependence. After evaluating the relative accuracy of the NPSs, we examined the 'added value' of the fine (2.5 km) horizontal resolution of the EU Copernicus Climate Change Service (C3S) Arctic Regional ReAnalysis (CARRA) data compared to the 31-km horizontal resolution ERA5 reanalysis data. We then presented a 31-year (1991-2021) CARRA Greenland ice sheet rainfall climatology and analyses of daily rainfall extremes, with meteorological context provided from both ERA5 and CARRA. The precise locations of the highest daily rainfall cases were presented along with tabulating the dates and amounts for the top 18 daily ice sheet rainfall episodes. We then examined what the fine 2.5-km resolution data could reveal in terms of mesoscale atmospheric dynamics, thermodynamics and precipitation during an extreme rainfall episode that coincided with the field observations. Finally, we quantified the heat impact of two rainfall extremes on snow and ice melt using a surface energy and mass budget approach applied to field data.

| Comparison of numerical weather prediction systems and observations
Of the five evaluated NPSs, the CARRA rainfall simulations showed the closest magnitude and temporal correlation to the available field data. The closer agreement with CARRA is likely due to its unique assimilation within the model domain of satellite and aircraft observations, air temperature and barometric pressure from roughly 40 on-ice automatic weather stations (AWS). All NPSs exhibited a 10%-60% dry difference compared to the observations, which decreased towards higher elevations. For the southern ice sheet, CARRA has a 17% wet bias at 271 m elevation, decreasing to 6% at the uppermost observation site (893 m elevation). Although the verification at two regions does not fully represent the overall NPS accuracy, the pattern in the spatial difference is consistent with earlier findings for the Norwegian coast (Køltzow et al., 2019). The difference between CARRA and ERA5 rainfall in this study follows a similar pattern, which is likely due to the coarser grid resolutions and subsequent ERA5 smoothing of the sharply rising Greenland terrain, resulting in excess inland moisture delivery. Although the improved accuracy of CARRA makes it a credible tool for evaluating and understanding rainfall in all areas of Greenland where in situ data are not available, the agreement among the NPSs for extreme daily totals is inconsistent. Single extreme cases are likely more affected by spatial complexity in precipitation placement exhibited by all NPSs.

| Rainfall climatology
CARRA data from 1991 to 2021 indicated that rainfall accounted for over 50% of total precipitation for narrow areas of the southern and western ice sheet, including the northwest. Local daily rainfall extremes above 300 mm occurred 18 times in 31 years of data. For the 31 years that was examined, the peak daily rainfall rate was 448 mm, occurring along the southeastern Greenland coast. Most rainfall extremes were near the southeast tip of Greenland, and all daily rainfall cases above 300 mm occurred south of the Arctic circle. Extreme rainfall dates occurred throughout the year, with most cases coinciding with high air temperature. CARRA annual rainfall totals for the ice sheet and its peripheral glaciers were in general agreement with previous estimates from NPS studies. CARRA simulated positive trends in rainfall by 39% on average and standard deviation of 25 ± 8 Gt y À1 , respectively, and also a 33% increase in rainfall fraction of precipitation over the 1991 to 2021 period. Similarly, for peripheral ice masses, CARRA rainfall increased significantly (1p = 0.90) by 37% in the 31-year period with an average of 1.8 ± 0.5 Gt y À1 . Snowfall on the Greenland ice sheet and peripheral ice masses in CARRA data showed no trend in this period, consistent with earlier studies.

| Mesoscale dynamics and precipitation for an atmospheric river reaching Greenland
The CARRA data provides new details of atmospheric dynamics associated with orographic precipitation. During an extreme rainfall episode on 14 September 2017, strong onshore flow caused flow-splitting and acceleration around southern Greenland within $2 km of sea level. An offshore updraft jet 40-140 km offshore formed in an arc around southern Greenland, with updraft velocity intensified by buoyancy generation from condensational (diabatic) heating. The same buoyancy updraft mechanism was maintained in long (100-200 km) and narrow (5-15 km) and vertically extensive (3 km deep) along-flow 'rapids' flowing 2 km above sea level at high speed (above 30 m s -1 ) for more than 24 hours. The rapids produce north-south oriented rainfall bands. Rapids occurred upstream and at a distance from the updraft jet, suggesting that they may be a general phenomenon ocurring in ARs.
Surface topographic channelling focuses atmospheric flow entry to the southern Greenland ice sheet, with large-scale forced and buoyant uplift contributing to extreme (over 300 mm) daily rainfall at 750 m elevation on 14 September 2017. Coastal mountains amplify gravity waves. Buoyancy generation from condensation enhanced uplift over the ice sheet and amplified the gravity waves, which oscillated downstream of the mountains in concert with ice sheet terrain variations, and presumably from gravity wave reflection against the stably stratified tropopause. In a detailed view into orographic precipitation, the CARRA data resolved gravity waves and forced uplift onto the ice sheet, which combined with diabatic enhancement of buoyancy, to increase the rainfall rate by up to a factor of four.

| Heat impact of rainfall on snow and ice melt
Rainfall in two extreme cases was estimated to have enhanced ice melt rates directly by 16 ± 4%, a modest amount that is consistent with earlier findings. Indirect effects can include snowpack heating through percolation and refreeze and a subsequent melt-albedo feedback initiated by heat and rainfall .

| Recommendations for future work
To support Greenland ice mass balance assessments, future research should further determine the accuracy of precipitation calculations using a larger number and spatial density of measurement sites. The effort can assist in establishing the level of confidence in NPS simulation of rainfall spatial heterogeneity, especially for extremes.
Further study of CARRA data can illuminate the interplay of forced uplift and buoyancy effects on gravity waves and their relation to orographic precipitation, surface heat transfer, aviation hazards, strong wind events, and so on.
programme (INTAROS, grant no. 727890) and the Danish Ministry of Climate, Energy and Utilities via The Programme for Monitoring of the Greenland Ice Sheet (PROMICE). The C3S Arctic regional reanalysis project has been supported by EU C3S contract 2017/ C3S_322_Lot2_METNO/SC2. B. Noël was funded by the Fonds de la Recherche Scientifique de Belgique (F.R.S.-FNRS). M. Niwano was supported by the Japan Society for the Promotion of Science through Grants-in-Aid for Scientific Research number JP17KK0017 and JP21H03582. M. van den Broeke acknowledges support of the Netherlands Earth System Science Centre (NESSC).