1. Introduction
Evaluation of the contribution of ice masses to sea-level rise requires knowledge of surface mass balance (SMB). Although SMB can be measured directly by site-based methods such as snow pit stratigraphy, snow/firn core stratigraphy, sonic depth data and stake arrays, such measurements are localised and logistically demanding. In contrast, studies assessing overall ice-sheet mass change rely on spatially extensive satellite-based measurements of, for example, gravity and surface elevation change (e.g. Davis and others, Reference Davis, Li, McConnell, Frey and Hanna2005). However, this method depends on the accurate constraint of other contributing processes such as firn densification, ice deformation and isostatic rebound (e.g. Spikes and others, Reference Spikes, Csatho, Hamilton and Whillans2003). The first two of these processes reduce to vertical strain within the snow, firn and ice column. Such strain measurements can also be used to improve model-based ice core chronologies, to reconstruct age–depth relationships in the absence of clear annual layering (e.g. Hawley and others, Reference Hawley, Waddington, Morse, Dunbar and Zielinski2002), and to guide depth-density functions (Kingslake and others, Reference Kingslake, Martín, Arthern, Corr and King2016). Horizontal gradients in vertical strain rate are also responsible for the development of isochrone arches in the stratigraphy of ice rises (Raymond, Reference Raymond1983; Matsuoka and others, Reference Matsuoka2015). Despite this utility, however, vertical strain rate profiles are difficult to measure directly and such measurements are correspondingly scarce. Thus, the overwhelming majority of SMB reconstructions from measurements of annual layer thicknesses have, to date, included model-based approximations of vertical strain (e.g. Thomas and others, Reference Thomas, Marshall and McConnell2008), assumed to be either constant (Nye, Reference Nye1963) or constant near the surface and to decrease linearly at depth (Dansgaard and Johnsen, Reference Dansgaard and Johnsen1969), when they are not guided by a full-Stokes model taking into account the Raymond effect (e.g. Drews and others, Reference Drews2015; Philippe and others, Reference Philippe2016). Thus, as pointed out by Kingslake and others (Reference Kingslake2014), the non-linearity of vertical strain, particularly near ice divides, is not fully captured by current ice-flow models.
Various methods have been used to measure vertical strain directly. The ‘coffee-can’ method (CCM) (Hamilton and others, Reference Hamilton, Whillans and Morgan1998) involves fixing an anchor at depth in the firn or ice column and recording that anchor's vertical displacement over time relative to a reference point fixed at or near the surface. The technique normally involves measuring the length of pre-tensioned cable retrieved over a pulley wheel mounted on a surface reference datum. Cable retrieval is measured either during repeated site visits or automatically at higher temporal resolution by a depth encoder attached to the surface pulley (e.g. Arthern and others, Reference Arthern, Vaughan, Rankin, Mulvaney and Thomas2010). Fibre-optic cable stretch has also been recorded for this purpose (e.g. Zumberge and others, Reference Zumberge2002). Although valuable, the technique is prone to error from factors such as pole bend or tilt, mechanical friction, influence of wind or falling snow on the cable, lateral strain within the installation borehole and anchor slippage. All such ‘surface-to-marker’ measurements also provide only a single net strain rate averaged over the distance separating the marker from the surface reference.
There is therefore a need for precise measurements of distributed strain between closely spaced markers covering the full depth range of interest. While distributed fibre-optic analysis has been used to measure discretised strain along terrestrial surfaces (e.g. Moore and others, Reference Moore, Gischig, Button and Loew2010), and has great potential, the ability of the technique to record distributed vertical strain in ice boreholes has yet to be demonstrated. At least three other techniques have gone some way to this end. First, phase-sensitive radio-echo sounding (pRES) has recently been developed with the capacity to record the depth of radar layers at high spatial resolution (e.g. Kingslake and others, Reference Kingslake2014; Martín and others, Reference Martín, Mulvaney, Gudmundsson and Corr2015; Young and others, Reference Young2019). Repeat or autonomous pRES surveys at a given location can therefore be used to difference the thicknesses of individual layers to produce vertical profiles of spatially discrete vertical strain and velocity across the full depth of radar wave penetration. For example, Kingslake and others (Reference Kingslake2014) measured englacial vertical velocities within ice divides up to 900 m thick and studied spatial variations in vertical strain rates along km-long transects across Berkner Island, Roosevelt Island, Fletcher Promontory and Adelaide Island, Antarctica. However, the application of pRES to millimetre-resolution layer differencing is hampered by several factors, including (i) temporal changes in material properties that influence radar wave velocity and thus inferred depth; (ii) limits to radar wave penetration and (iii) slight variations in antenna geometry, and hence ray-path alignment, between individual surveys. Second, adapting the technique of Rogers and LaChapelle (Reference Rogers and LaChapelle1974), Raymond and others (Reference Raymond1996) and Hawley and others (Reference Hawley, Waddington, Lamorey and Taylor2004) calculated vertical strain by measuring repeatedly the locations of metal markers installed at different depths within boreholes at Dyer Plateau, Antarctica, and Siple Dome, Antarctica, respectively. While these studies did yield vertical strain data, the technique was found to be logistically demanding and potentially subject to errors associated with collar emplacement and slippage. The technique also provides a limited number of observations and cannot easily be applied to irregularly walled boreholes such as those drilled by hot water. Third, Hawley and Waddington (Reference Hawley and Waddington2011) measured vertical strain from repeated measurement of natural layers recorded in optical logs of a borehole drilled to a depth of ~30 m at Summit, Greenland. These authors reported firn compaction using a downward-looking borehole camera on the basis of three logs repeated annually. The resulting vertical velocity profile matched closely that predicted by a compaction model based on measured densities and an assumption of steady surface accumulation over the ~70-year time period concerned. Although representing a notable advance, the use of a downward-looking camera that was not always centred in the borehole resulted in a degree of feature blurring. Stratigraphic logs were therefore low-pass filtered, with a threshold of 0.075 m, and the luminosity traces were differenced by a cross-correlation function applied separately to several 0.5 m-long reference sub-sections. Although somewhat limited by the camera technology available at the time, the repeat optical logging approach of Hawley and Waddington (Reference Hawley and Waddington2011) has the capacity to yield high resolution firn compaction data. Importantly, borehole optical televiewing (OPTV) – whereby the camera images horizontally around the borehole rather than directionally down the borehole (Hubbard and others, Reference Hubbard2008) – overcomes the key technological limitations of directional borehole imaging.
Here, we apply repeat OPTV logging to a ~100 m-deep borehole (IC12) drilled electromechanically into the dome of Derwael Ice Rise (DIR), East Antarctica. These data are used to reconstruct the ice rise's vertical strain profile at millimetre resolution. We combine these data with GNSS-based measurements of surface vertical velocity to calculate DIR's absolute vertical velocity profile. Finally, we compare these data with CCM data recorded over three depth ranges and the output of a full-Stokes ice-flow model.
2. Field site and methods
2.1. Field site, ice coring and surface positioning
The study site is located at the crest of ~550 m-thick DIR (Drews and others, Reference Drews2015) in coastal Dronning Maud land, East Antarctica (Fig. 1). A ~120 m-long ice core (‘IC12’) was recovered from the site by electro-mechanical corer in December 2012. Annual layer counting based on δ18O and δD, major ion concentration and electrical conductivity was used to date the core, the base of which was dated to ~1759 ad (Philippe and others, Reference Philippe2016). Philippe and others (Reference Philippe2016) also measured ice-core density directly, supplemented by OPTV-based estimates following Hubbard and others (Reference Hubbard2013). Combining these density data with annual layer thicknesses allowed DIR's long-term accumulation record to be reconstructed, reported by Philippe and others (Reference Philippe2016) as 0.52 ± 0.01 m ice equivalent a−1 over the full time period concerned. This record also showed an increasing trend, particularly since the mid-20th century, such that mean accumulation at DIR for the period 1992–2011 had increased to 0.70 ± 0.01 m ice equivalent a−1.
A GNSS station was installed within a few tens of metres of IC12 to measure daily surface position from late 2012 to early 2016. The antenna position, initially anchored at a depth of 1.85 m, was determined using Bernese 5.2 software (Dach and others, Reference Dach, Lutz, Walser and Fridez2015). Processing was based on GNSS position referenced, via long-term and overlapping GPS and GLONASS data (Dow and others, Reference Dow, Neilan and Rizos2009), to the global ITRF2014 (Altamimi and others, Reference Altamimi, Rebischung, Métivier and Collilieux2016). Vertical velocity was calculated by stacking the daily solutions using CATREF (Altamimi and others, Reference Altamimi, Sillard and Boucher2007) and modelling the annual signal, yielding a mean vertical surface velocity through 2013 and 2014 of −1.380 ± 0.009 m a−1.
An octagonal strain network was established using eight markers located along a circle of 2 km radius around IC12. These markers were positioned using differential GNSS in 2012 and 2013. Lateral divergence was calculated from orthogonal horizontal strain rates ($\dot{\varepsilon }_{xx}$ and $\dot{\varepsilon }_{yy}$) between two pairs of markers, located on perpendicular axes crossing at the dome. Using the principle of mass conservation ($-\dot{\varepsilon }_{zz} = \dot{\varepsilon }_{xx} + \dot{\varepsilon }_{yy}$) and setting strain as zero at the ice–bed interface yielded a depth-averaged vertical strain rate measured at the surface of −0.0022 ± 0.0002 a−1.
2.2. Borehole optical televiewing and layer thickness change
OPTV logs record an accurately-scaled and -orientated colour image of a complete borehole wall (Hubbard and others, Reference Hubbard, Roberson, Samyn and Merton-Lyn2008). The centralised OPTV probe used for this study illuminates the wall with an outward-facing circular array of 72 white LEDs, the luminosity of which can be controlled by the operator. The reflected image is recorded by a charge-coupled device camera that records the lateral 360°-pixel row of the adjacent wall. The resulting OPTV logs, recorded for this study initially in December 2012 and again in December 2014, have a lateral resolution of ~1.5 mm and a vertical resolution of 1 mm in 2012 and 2 mm in 2014, the former of which was averaged to 2 mm for consistency. The OPTV logs used herein include the uppermost 99.86 m of IC12 because the ~20 m-long section below this was scored by cutters during drilling in 2012, artificially degrading the log's luminosity trace as described by Hubbard and others (Reference Hubbard2013). OPTV logs were analysed using the Borehole and Ice Feature Annotation Tool (BIFAT) (Malone and others, Reference Malone, Hubbard, Merton-Lyn, Worthington and Zwiggelaar2013) and commercial WellCAD software. Once acquired, the 2014 log was shifted vertically to the common datum of the borehole top in 2012. We calculate OPTV luminosity as the mean brightness intensity (0–255, no unit) recorded around each pixel row.
We use two approaches to calculate the change in length (thickness) of multiple individual layers making up the firn/ice column (and hence of spatially distributed vertical strain), both based on matching the luminosity traces of the two OPTV logs. First, we quantify the vertical displacement of 60 manually-identified reference markers between 2012 and 2014. Each marker is defined based on being (i) clearly identified as a single-pixel (i.e., 2 mm) peak or trough in the luminosity trace of both logs, and (ii) separated by at least 1 m, optimising the length over which vertical strain is calculated. These 60 markers, plus the surface, define 60 separate layers of mean thickness 1.67 m and sd 0.95 m. The vertical velocity of each marker is calculated relative to the surface, where vertical velocity was measured by GNSS (Section 2.1). Second, we use an amplitude cross-correlation method (XCM) to calculate the vertical displacement of the firn/ice column at incremental 1 m depth segments. The displacement of each segment between 2012 and 2014 was quantified within a ±2 m search window to the nearest 2 mm by identifying the highest value in the cross-correlation matrix, which indicates the closest luminosity match.
Although the depth encoder linked to the OPTV's cable has a precision of ~0.04 mm, the system is also subject to uncertainty related to system set-up in the field. To evaluate this, we logged IC12 twice on the same day, fully dismantling and reassembling equipment between the logs. Comparing the apparent displacement of each of the study's 60 markers yields a mean difference of 0.0 mm and a depth-independent, normal distribution (p-value = 0.364) of residuals of sd 3.78 mm (Fig. 2). The most likely causes of such mm-scale differences are (i) OPTV probe wobble, (ii) variability in cable stretch between the two logs and (iii) cable slip over the sheave wheel during either or both logs. Hereafter, we take this (1 sd) length of 3.78 mm to represent the uncertainty associated with repeated OPTV-based measurements of marker depth.
2.3. Coffee-can method (CCM)
The CCM involves installing an anchor at a known depth within a borehole and recording its displacement through time relative to a fixed, near-surface datum (Section 1). Anchor slippage, identified as an unfeasibly large strain, is not uncommon in CCM experiments and two of the five installed coffee cans were omitted from the current study for this reason. All anchors were installed in separate boreholes, each drilled by hot water within 100 m of IC12. The three anchors that held firm were installed at depths of 12.00, 21.70 and 42.35 m and their depths were recorded ~12 months apart in December 2012 and December 2013. In order to enable comparison with the OPTV-derived strain rates at the highest possible spatial resolution, the resulting CCM data are expressed as vertical strain averaged over the distance separating successive anchors. Considering the nature of the experimental set-up and borehole conditions at DIR we adopt a CCM depth uncertainty value of 0.2 m.
2.4. Vertical strain rate and submergence velocity
For the layer differencing method, the vertical strain rate $\dot{\varepsilon }_{zz}\lpar z\rpar \;$of each layer (i.e. the space between adjacent markers), l, is calculated as the final layer thickness, l 2, minus the initial layer thickness, l 1, per l 1 and unit time separating the measurements, dt, in this case 720 days (Eqn (1)):
Relative vertical velocity is calculated as the depth difference of each marker between two measurements (2014 cf. 2012) divided by the time separating those measurements (720 days). Since depth is measured relative to the 2012 surface, the absolute velocity of each marker is given by subtracting its local (upwards) velocity from the (downwards) velocity measured at the surface (−1.38 m a−1; Section 2.1). For the XCM, velocity is calculated directly from the relative displacement of each section analysed, and strain is calculated by differencing the velocity profile, following Young and others (Reference Young2019). Adopting a 1 sd uncertainty in OPTV-reconstructed layer depth of 3.78 mm (Section 2.2) yields an equivalent uncertainty in vertical velocity of 0.008 m a−1.
Besides presenting the raw vertical strain and vertical velocity data, we derive smoothed depth profiles by using total-variation regularisation when applying spatial derivatives to the measured data points. This technique applies regularisation to the differentiation operator directly (Chartrand, Reference Chartrand2011), preventing the problematic amplification of noise during layer differencing. Herein, we define the regularisation parameter such that the reconstructed values lie within 1 sd of the measured data.
Since the vertical velocity profile reported here is relative to the 2012 surface, we consider a directional (downwards) 50 mm error to account for additional compaction due to the weight of the wooden board covering the 2012 borehole. This uncertainty also affects the uppermost strain rate datum as it is computed between the first englacial marker and the 2012 surface. Below that first measurement, strain rate uncertainty is calculated using standard error propagation to account for the 3.78 mm uncertainty affecting both the 2012 and 2014 marker depths.
Ice-equivalent vertical strain rates are also calculated for comparison with modelled strain rates (Section 2.5) and are derived from ice-equivalent velocities, using the depth-density profile of Hubbard and others (Reference Hubbard2013) for IC12 and assuming that this profile is in steady-state between 2012 and 2014. The 2012 surface was located ~2.6 m beneath the 2014 surface.
2.5. Ice-flow modelling
We compare our empirically-reconstructed strain rates and (ice-equivalent) vertical velocity data with simulated values for DIR using the full-Stokes ice-flow model Elmer/Ice (Gagliardini and others, Reference Gagliardini2013), applied at this field site by Drews and others (Reference Drews2015), accounting for ice anisotropy (with n = 3 in Glen's flow law) and long-term surface lowering (0.03 m a−1 over 3400 years). The modelled vertical velocity profile was scaled to match the vertical velocity reconstructed from the surface GNSS data (Section 2.1).
3. Results
3.1. OPTV logs
The 2012 and 2014 OPTV logs of IC12 (Figs 3a, b respectively) are remarkably similar overall, showing a general decrease in luminosity with depth and the presence of numerous horizontal layers. The former effect has been interpreted as reduced material reflectivity through snow metamorphism to ice (Hubbard and others, Reference Hubbard2013) and the latter as annual or sub-annual layering (Hubbard and Malone, Reference Hubbard and Malone2013). However, it is likely that the anomalously large negative luminosity excursions, for example at depths of 5, 17 and 35 m, correspond to infiltration ice layers resulting from the refreezing of meltwater produced by substantial surface melt events. It is also apparent from the raw OPTV images shown in Figure 3 that these three infiltration ice layers become displaced vertically from each other during the 2-year period separating the logs, illustrating the compression of the firn column.
The infiltration ice layers noted above represent three of the 60 markers used to calculate strain rates and velocities in the uppermost 100 m of IC12. We illustrate this analysis by focusing on a 6.6 m-long section of the OPTV logs, the luminosity profiles of which are presented in Figure 4. In Figure 4a, the luminosity traces of both logs are overlaid after the 2012 log has been raised to match the 2014 log at the top of the section. With this local alignment, the remarkable overall similarity of the two OPTV logs becomes apparent: individual peaks and troughs match almost precisely at millimetre scale, illustrating the capability of the OPTV-based method for matching layers at high spatial resolution. In Figure 4b, the local depth adjustment is removed and the luminosity traces are separated laterally to illustrate the offset of the five markers (of 60 in total) used in our layer-differencing analysis that are located within this section.
3.2. Vertical strain rates
The dashed lines linking the five markers illustrated in Figure 4b indicate a relative displacement of ~1.5 m over 2 years separating the logs. This equates to a vertical strain rate, averaged over the depth range 0–67 m, of −0.01 a−1. Further, even within this 6.6 m-long section, the offset between the two luminosity traces increases progressively with depth, from zero at the local tie-point at a depth of 59.05 m to ~60 mm at the base of the section at a depth of 67.1 m. This equates to a vertical strain rate, averaged over this depth range, of ~−0.004 a−1. Indeed, analysis of multiple markers and variable-offset XCM data allows higher-resolution spatially-distributed (discretised) strain to be calculated across numerous individual depth ranges (Fig. 5). These raw discretised strain rates, plotted as black dots in Figures 5a and 5b, generally decrease with depth, as anticipated based on the increase in density, from ~−0.07 a−1 near the surface to ~−0.002 a−1 towards the base of the borehole. These discretised strain rates also show substantial local variability (Fig. 5b).
Although showing similar overall patterns, comparing Figure 5a with Figure 5b reveals that the two methods of calculating vertical strain from the OPTV luminosity record differ in detail. This reflects the two methods used, with the peak matching method (Figs 4 and 5a) analysing 60 layers of mean length 1.67 m and the XCM (Fig. 5b) analysing 100 layers of average length 1 m (with the latter therefore recording greater variability). Further, the automated XCM returned three anomalous points resulting from the close similarity of certain adjacent peaks in the luminosity record (such that, in these three cases, an incorrect shift yielded a higher correlation than the correct shift). These three points were removed from the analysis.
The results of the three CCM experiments are superimposed in red on the OPTV-based data in Figure 5. These data highlight the longer distances over which the CCM strain rates are measured, the low number of measurements, and the large error bars associated with the method. Nonetheless, the resulting CCM vertical strain rates are consistent with the OPTV-derived rates.
Finally, the vertical strain rate profile predicted by the full-Stokes 2-D flow model is broadly consistent with the discretised OPTV-based strain rates corrected for density and expressed in ice equivalence (Fig. 5c). However, in detail, modelled strain rates are both less variable (e.g. note the larger measured strain values in the uppermost 10 m and in the depth interval 24–45 m) and, on average, lower than the measured rates, particularly in the uppermost half of the borehole. Thus, the mean OPTV-derived (ice equivalent) strain rate below 50 m depth is ~−0.0028 a−1 (Fig. 5c), which is – accounting for uncertainty – indistinguishable from that calculated from surface GNSS (−0.0022 a−1) and predicted by the model (−0.0019 a−1). In contrast, above 50 m depth, the mean OPTV-derived (ice equivalent) strain rate increases by an order of magnitude to −0.0124 a−1 while the GNSS-derived and model-derived strain rates are unchanged.
3.3. Vertical velocities
Figure 6a shows the vertical velocities of the 60 OPTV markers with the regularised velocity profile plotted as a superimposed green line. Vertical velocity decreases from −1.3 m a−1 at the surface to −0.5 m a−1 – consistent with surface accumulation – at the base of the borehole. This velocity decrease generally follows that of strain (Section 3.2), as anticipated.
The results of the three CCM experiments are superimposed (in red) on the OPTV-based data in Figure 6. Similar to the strain rates from which these data are derived, the CCM velocities are averaged over longer distances and subject to greater uncertainty than the OPTV-derived velocities, yet both techniques yield consistent results.
4. Discussion and conclusions
Two OPTV-based luminosity profiles of the IC12 borehole, recovered initially in late-2012 and subsequently in late-2014, reveal remarkable millimetre-scale similarity in the detail of imaged marker layers. Layer differencing and cross correlation of the luminosity traces of both logs allowed a high-resolution, discretised vertical strain profile to be calculated, indicating rates of up to −0.07 a−1 near the surface, decaying to values of ~−0.002 a−1 below ~50 m. The results are consistent with CCM-derived strain rates and local submergence velocities.
DIR's net vertical strain and velocity are influenced by the Raymond effect (Raymond, Reference Raymond1983), as well as its thickness (~550 m) and surface accumulation (~0.52 m ice equivalent a−1) (Lingle and Troshina, Reference Lingle and Troshina1998; Matsuoka and others, Reference Matsuoka2015). Our OPTV-derived strain rates show notable high-resolution variability at all depths, likely related to temporal changes in SMB and to vertical variations in material compressibility (reflecting properties such as density, layering (e.g. Hörhold and others, Reference Hörhold, Kipfstuhl, Wilhelms, Freitag and Frenzel2011) and impurity content (e.g. Freitag and others, Reference Freitag, Kipfstuhl, Laepple and Wilhelms2013)). These strain rates include the effects of firn densification and lateral ice deformation. Isolating these contributions using a density profile for the borehole (Philippe and others, Reference Philippe2016) allowed the resulting ice equivalent strain rates to be compared with output from a full-Stokes model of ice flow for DIR. Notably, the model, which operates at a vertical spacing of 20 m, generally matches our observations closely. However, the high-resolution variability apparent in the OPTV-derived data is not captured by the model's output because the model assumes time-invariant surface accumulation (using only the long-term average), and does not consider firn densification and its associated strain variability. Similarly, the model's under-estimation of strain, particularly in the upper half of IC12, is consistent with a general increase in accumulation at DIR over recent decades, as indicated by analysis of the IC12 ice core (Philippe and others, Reference Philippe2016).
Assuming that the vertical velocity profile is in steady state, that profile can be used to reconstruct initial annual layer thicknesses and hence SMB for a reasonable range of initial snow densities. Thus, in the absence of independent mass-balance measurements (which are available in this instance from the IC12 core; Philippe and others, Reference Philippe2016), OPTV-based layer-differencing has the capacity to reconstruct SMB without the need for a model or an independent density profile (e.g. Schwerzmann and others, Reference Schwerzmann2006). The age–depth relationship can also be deduced directly from the vertical velocity profile (e.g. Hawley and others, Reference Hawley, Waddington, Morse, Dunbar and Zielinski2002) to improve model-based ice core chronologies and provide useful validation for ice-core dating where annual layering may not be clear. Indeed, deviations from the steady-state condition could potentially be investigated through detailed comparison of such chronologies.
Besides exploring the causes of high-resolution variations in strain reported herein, future work could focus on assessing the potential of this technique for strain rate measurements in boreholes drilled by hot water. If layering is still clearly identifiable, rapid access to multiple boreholes would allow the measurements to be extended laterally and/or vertically by drilling deeper, possibly to the bed (reducing the uncertainty associated with a moving surface reference). Logging two closely-spaced boreholes drilled at different times would introduce additional uncertainty, but it would also avoid the need to locate and excavate a single buried borehole.
Finally, in terms of technique evaluation, repeat OPTV logging provides an accurate way of measuring discretised vertical strain and velocity at high spatial resolution. However, repeat OPTV logging is logistically demanding and should be considered as an option alongside emerging techniques with similar – but currently unrealised – potential such as pRES and distributed fibre-optic strain.
Acknowledgments
This paper forms a contribution to the Belgian Research Programme on the Antarctic (Belgian Federal Science Policy Office Projects SD/SA/06A and BR/165/A2). Borehole drilling and OPTV instrumentation was supported by a Higher Education Funding Council for Wales Capital Equipment grant to Aberystwyth University and UK Natural Environment Research Council Grant (NE/J013544/1) to BH. The authors wish to thank the International Polar Foundation for logistical support in the field. MP was partly funded by a grant from Fonds David et Alice Van Buuren. RD was partially supported by the DFG Emmy Noether Grant DR 822/3-1. The authors thank Sophie Berger for compiling Figure 1. The data presented herein are available from https://figshare.com/articles/Data_sets/11799090/1.