Arctic soil development on a series of marine terraces on central Spitsbergen, Svalbard: a combined geochronology, fieldwork and modelling approach

Soils in Arctic regions currently enjoy attention because of their sensitivity to climate change. It is therefore important to understand the natural processes and rates of development of these soils. Specifically, there is a need to quantify the rates and interactions between various landscapeand soil-forming processes. Soil chronosequences are ideal natural experiments for this purpose. In this contribution, we combine field observations, luminescence dating and soil–landscape modelling to improve and test our understanding of Arctic soil formation. The field site is a Holocene chronosequence of gravelly raised marine terraces in central Spitsbergen. Field observations show that soil–landscape development is mainly driven by weathering, silt translocation, aeolian deposition and rill erosion. Spatial soil variation is mainly caused by soil age, morphological position within a terrace and depth under the surface. Luminescence dating confirmed existing radiocarbon dating of the terraces, which are between ∼ 1.5 and ∼ 13.3 ka old. The soil-landscape evolution model LORICA was used to test our hypothesis that the field-observed processes indeed dominate soil–landscape development. Model results additionally indicated the importance of aeolian deposition as a source of fine material in the subsoil for both sheltered and vegetated trough positions and barren ridge positions. Simulated overland erosion was negligible. Consequently, an un-simulated process must be responsible for creating the observed erosion rills. Dissolution and physical weathering both play a major role. However, using present-day soil observations, the relative contribution of physical and chemical weathering could not be disentangled. Discrepancies between field and model results indicate that soil formation is non-linear and driven by spatially and temporally varying boundary conditions which were not included in the model. To conclude, Arctic soil and landscape development appears to be more complex and less straightforward than could be reasoned from field observations. Published by Copernicus Publications on behalf of the European Geosciences Union. 222 W. M. van der Meij et al.: Arctic soil development on a series of marine terraces on central Spitsbergen


Introduction
Soils in Arctic and boreal landscapes have recently received intense research interest, because the climate in these regions is expected to experience stronger changes than elsewhere (e.g.Arctic Climate Impact Assessment, 2004;Forland et al., 2011;Zwoliński et al., 2008).The effects of this increase are so far only partially understood (e.g.plant community development;Hodkinson et al., 2003).Another point of interest in the area is the poorly constrained Arctic carbon pool and its potential as a carbon sink (e.g.Ping et al., 2008).To provide context to the short-term changes (∼ 100 years) in Arctic and boreal soils that we are currently observing, knowledge on long-term soil development (∼ 10 000 years) is urgently required as a baseline: we need to better constrain the natural (i.e.paraglacial; Ballantyne, 2002;Slaymaker, 2011) processes, rates and feedbacks in the soil-landscape system.With such understanding, meaningful comparisons can be made between short-term rates of change in soils due to changing climate on the one hand and long-term rates of change in soils on the other.
Chronosequences are a popular means to obtain information about natural rates of soil formation (e.g.Birkeland, 1992;Egli et al., 2006;Phillips, 2015;Sommer and Schlichting, 1997).In chronosequences, the only soil-forming factor that is significantly different for all soils is time.Variation in the other soil-forming factors, i.e. landscape position, climate, lithology and organisms, is assumed equal for all soils in the study area.Consequently, variation in soil properties can mainly be attributed to the age of the soil (Vreeken, 1975).In Arctic regions, two paraglacial landscape settings are particularly suitable for chronosequences.Proglacial areas, where glaciers are currently retreating, are often used to compare soils formed at the onset of the recent retreat (∼ 100 years ago) with those formed in very recently exposed glacial parent material.This can illustrate decadal rates of soil formation (Egli et al., 2014;Kabala and Zapart, 2012).Another chronosequence setting is provided by series of marine terraces, also known as raised beaches, which reflect millennial-scale isostatic rebound after the end of the Last Glacial Maximum.Such terraces are ubiquitous in Arctic landscapes (Scheffers et al., 2012).Terrace chronosequences can provide millennial rates of soil formation, which is particularly helpful because natural soil formation in Arctic regions is relatively slow (Bockheim and Ugolini, 1990;Fischer, 1990) and many differences become apparent only after thousands of years.
Several factors nonetheless distort the temporal signal in chronosequences of marine terraces, as occurs more often in long-term chronosequences (e.g.Birkeland, 1990).First, a typical terrace consists of slightly elevated ridge positions and lower trough positions (Pereverzev and Litvinova, 2010) and thus contains altitude differences resulting in different hydrological conditions that affect soil formation (Makaske and Augustinus, 1998;Scheffers et al., 2012).Second, geo-morphic processes may have a different effect on not only ridges and troughs but also terraces at different positions in the landscape (Pereverzev and Litvinova, 2010) -particularly where a marine terrace complex is part of otherwise mountainous topography.Erosion and deposition can occur with different rates on different terrace levels (Strzelecki, 2012).Third, it is difficult to verify whether the composition and particle size distribution of soil parent material (beach deposits) at the onset of soil formation have been the same within and between terrace levels (Mann et al., 1986).In other words, landscape position and composition of parent material may also have played a role in determining the present soil heterogeneity (Temme and Lange, 2014).These complications to chronosequences can lead to a problem of attribution: are observed differences between soils predominantly the result of a difference in age, or are other factors important as well?
The attribution problem can only be solved by using a combination of various methods.Clearly, geochronology is needed to provide accurate dating of the initiation of soil formation, and field and laboratory observations of soils are needed to determine properties of interest.However, in addition to these methods, model simulations of the various effects of time and other soil-forming factors on soil development in a landscape context are essential to determine which differences in soil-forming factors may have caused differences in observations.A combination of these methods is thus needed to study long-term Arctic soil development that is influenced not only by time but also topographical position.
In this study, we focused on soils in a sequence of marine terraces in central Spitsbergen, Svalbard archipelago, to derive natural processes and rates of soil formation in a landscape context (Elster and Rachlewicz, 2012;Rachlewicz et al., 2013;Zwoliński et al., 2013).We first used optically stimulated luminescence (OSL) dating to complement earlier experimental datings of juvenile marine shells on the same series of terraces (Long et al., 2012).Then, we performed field and laboratory analyses to describe soil properties in a variety of locations on the marine terrace complex.Together with dating results, this allowed us to calculate rates of some soil-forming processes.Third, we used these rates to simulate combined soil-landscape development using a spatially distributed soil-landscape evolution model.Soil-landscape modelling has hitherto rarely been used in soil chronosequence studies (but see Sauer et al., 2012).However, by combining the various interacting geomorphic and pedogenic process, it allowed us to test and increase our understanding of interacting soil and landscape shaping processes in the study site.For simulations, we first hypothesized which soilforming processes played a dominant role.Next, the recently developed soil-landscape evolution model LORICA (Temme and Vanwalleghem, 2015) was adapted to reflect this hypothesis.Model inputs and parameters were derived from field observations.Model outputs were compared to observations and conclusions were drawn with regard to the validity of our hypotheses.

Location and geomorphology
Fieldwork was conducted in the Ebba Valley, one of the glacial valleys that enter Petunia Bay in the north tip of the Billefjorden, central Spitsbergen (Svalbard archipelago, Fig. 1).A sequence of six marine terraces is located at the mouth of the valley, bordered by the Ebba River and floodplain to the north, alluvial material to the east and south and by the fjord to the west.Prominent erosion rills and tundra lakes (Mazurek et al., 2012) were excluded from the study area (Fig. 1).The terrace sediments (i.e.soil parent material) dominantly consist of well-rounded gravel and coarse sand of limestone lithology, but gravel and sand from shale, sandstone and mafic intrusions are also found.
The area has been the subject of research for many years (e.g.Gulińska et al., 2003;Kłysz et al., 1988Kłysz et al., , 1989;;Long et al., 2012;Zwoliński et al., 2013).The marine terraces occupy a range of altitudes in the landscape (∼ 1-50 m), due to isostatic rebound after the Last Glacial Maximum.The typical smooth ridge and trough morphology (Makaske and Augustinus, 1998;Scheffers et al., 2012) of terraces was formed by wave-action and sea-level fluctuations.Six terrace levels can be distinguished, each consisting of a smaller series of ridges and intermediate troughs.The oldest marine terrace in the series (terrace 6, Fig. 1) dates back to the Late Pleistocene (Kłysz et al., 1989), yet it is very small and was not sampled in the present study.Terrace levels 1-4 have been dated using an experimental approach of radiocarbon dating of juvenile marine shells (Long et al., 2012).The ages range from 3156 ± 81 to 9718 ± 91 years, suggesting that younger soils might have been flooded again (Strzelecki, 2012).Age increases continuously with increasing elevation.The approximate locations and individual ages of the datings of Long et al. (2012) are displayed in Figs. 1 and 3.
Due to their slightly more sheltered position and lower altitude relative to the smooth ridges, the troughs have denser vegetation.Ridge positions are in general free from vegetation, but can be partly covered by bacterial soil crusts.In an aerial photograph from summer 2009, the barren ridges and terrace edges are characterized by lighter colours, whereas trough positions are characterized by darker colours (Fig. 1).

Arctic soils
Most soils of Spitsbergen have formed in coastal settings.Soils typically have shallow profiles with poorly differentiated genetic horizons of sandy or loamy texture, pH values varying between 7 and 8 and organic carbon contents from 0 to 10 % (Melke and Chodorowski, 2006;Pereverzev, 2012).In some cases, soils have been affected by geomorphic ac- tivity such as cryogenic processes and erosion (Lindner and Marks, 1990).Thickness of the marine deposits is between 1 and 2 m (Zwoliński et al., 2013).Soils formed in those deposits are well developed compared to proglacial soils but are nonetheless mainly described as incompletely developed soils (Cambisols, Cryosols, Leptosols, Regosols; Kabala and Zapart, 2009).In the Ebba Valley, the thickness of the active layer varies between 0.3 and 2.5 m (Gibas et al., 2005), with thaw depths ranging from 0.45 to 1.2 m on the marine terraces (Rachlewicz and Szczuciński, 2008).
The dominantly mentioned soil-forming processes are weathering through frost action and dissolution (Forman and Miller, 1984;Kabala and Zapart, 2009), calcification (Courty et al., 1994;Ugolini, 1986), silt eluviation (Forman and Miller, 1984) and the formation of organic matter (Melke, 2007).Especially the process of silt eluviation is typical for the coarse-grained Arctic soils.Forman and Miller (1984) identified six stages of silt eluviation, which indicate an increasing presence of silt caps on top of clasts for stage 1-4.In stage 5 and 6 the individual silt caps connect and fill the space between the clasts, eventually leading to a matrix-supported soil.The presence of silt caps is associated with coarsegrained and well drained soils (e.g.Locke, 1986;Ugolini et al., 2006), dense vegetation capturing aeolian silt (Burns, 1980, as stated in Forman andMiller, 1984), illuviation by precipitation (Locke, 1986) and vertical frost sorting (Bockheim and Tarnocai, 1998).

Climate
The study area has an average annual temperature of −5 • C, with average temperatures in summer and winter of +6 and −15 • C respectively (Przybylak et al., 2014).The average annual precipitation is 150-200 mm, mainly as snow (Láska et al., 2012;Rachlewicz and Szczuciński, 2008;Rachlewicz et al., 2013).The climatic conditions are more extreme compared to the western coast of Spitsbergen, with warmer summers, colder winters and less precipitation (Przybylak et al., 2014;Rachlewicz, 2009).The climate is classified as an Arctic desert or tundra (ET; Köppen, 1931).
The prevailing wind directions in the Ebba Valley are south or northeast with the strongest winds (> 6 m s −1 ) blowing from the Ebba Glacier in the northeast (Láska et al., 2012).These strong winds in combination with scarce vegetation and a high availability of sediments on the sandur plains lead to active wind erosion and subsequently to downwind accumulation of aeolian sediments.Deposition occurs when wind speed decreases or when the sediments get mixed with falling snow in autumn and winter, also known as niveoaeolian deposition (Rachlewicz, 2010).In the Ebba Valley plants are a good indicator of hydrological and soil characteristics, yet they reflect the cold climate.Hydrophilic species are found in wet trough positions, whereas vascular species were sporadically found on better drained and better developed soils on ridges (Jónsdóttir et al., 2006;Prach et al., 2012).Vegetation cover in the study area is around 30 % (Buchwal et al., 2013).

Luminescence dating
To complement the experimental rebound chronology from Long et al. (2012), we applied OSL dating to samples taken from the sand fraction of the marine sediments from terrace levels 1, 3 and 5 in the study area (Fig. 1).The samples were collected at a depth of 0.27 or 0.57 m from pits dug in the marine sediments (Table 2) and were shielded from light.The fine sand fraction (180-250 µm) of marine sediments can be assumed to be well bleached, due to reworking by wave action in the swash zone (Reimann et al., 2012).Nonetheless, to account for the possibility of insufficient signal resetting, termed partial bleaching, we applied the single-aliquot regenerative-dose (SAR) measurement protocol of Murray and Wintle (2003) and made use of a small-aliquot approach (e.g.Reimann et al., 2012;Rodnight et al., 2006).
Two quantities are determined for OSL dating.First, measurement of the OSL signal on the purified quartz mineral fraction reveals how much ionizing radiation the sample received since the last bleaching event (i.e.prior to burial).Second, this measurement is combined with a measurement of the background radiation level at the sample position.The luminescence age (ka) is then obtained by dividing the amount of radiation received (palaeodose, Gy) by the rate at which this dose accumulates (dose rate, Gy ka −1 ): OSL age = palaeodose/dose rate. (1) The basic principles of OSL dating are reviewed in Aitken (1998) and Preusser et al. (2008).
For dose rate estimation we used high-resolution gamma ray spectrometry.Activity concentrations of 40 K and several nuclides from the uranium and thorium decay chains were measured.Results were combined with information on geographic location and burial history (Prescott and Hutton, 1994), water and organic content history (Aitken, 1998;Madsen et al., 2005).Furthermore, grain-size-dependent attenuation effects were incorporated (Mejdahl, 1979) to calculate the effective dose rate.The total dose rate is listed in Table 2.
For the OSL measurements the three sediment samples were prepared in the Netherlands Centre for Luminescence dating under subdued orange light conditions.The samples were sieved to obtain the 180-250 µm grain size fractions, which were subsequently cleaned using HCl (10 %) and H 2 O 2 (10 %).Grains of different minerals were separated from each other using a heavy liquid (LST).The quartz-rich fraction (ρ > 2.58 g cm −3 ) was then etched with 40 % HF for 45 min to remove remaining feldspar contamination and the outer rim of the quartz grains.The purified quartz fraction was again sieved with a 180 µm mesh to remove particles that had become too small by etching.
To estimate the palaeodose of the samples, the OSL from quartz was measured by applying the SAR protocol.The most light-sensitive and most suitable OSL signal of the quartz grains was selected using the "early background" approach (Cunningham and Wallinga, 2010).To obtain a good estimate of the palaeodose, measurements were repeated on at least 28 subsamples (aliquots) per sample.Each aliquot consisted of 40-70 grains.To test the SAR procedure, a dose recovery experiment was then carried out on four aliquots of each sample.The average recovered dose agreed with the laboratory given dose.The ratio of measured dose divided by given laboratory dose was 0.96 ± 0.02 (n = 11), confirming the suitability of the selected measurement parameters.
The small-aliquot palaeodose distributions were symmetric and moderately scattered.We calculated over-dispersion values, i.e. the scatter in the palaeodose distributions that cannot be explained by the measurement uncertainties (Galbraith et al., 1999), to be between 12 ± 3 and 33 ± 9 %.These over-dispersion values are typical for well-bleached sediments derived from coarse-grained marine deposits (e.g.Reimann et al., 2012).Furthermore, the over-dispersion increased with age, suggesting that it is unlikely that partial bleaching is the source of the unexplained scatter.Therefore, palaeodoses of our samples were derived from the singlealiquot palaeodose distributions by applying the central age model (CAM;Galbraith et al., 1999).

Soil observations
The study area was divided into three equally sized strata based on altitude, which in turn were divided into vegetated (trough) and unvegetated (ridge) substrata using the aerial photograph.Thirty random locations were divided over the six strata according to stratum size, with at least two locations in each stratum (Fig. 1).Eight pits were located on ridge positions and 22 pits on trough positions.Soil profiles were mostly described according to FAO standards (FAO, 2006;IUSS Working Group WRB, 2015).We deviated from the standard horizon designation in three ways, because that allowed easier comparison between different soils and because it better suits the locally observed soil properties.The deviations are as follows.(1) All aeolian horizons were recorded with prefix 1, whereas all horizons developed in marine parent material were recorded with prefix 2, including the ones without aeolian cover.This was done to easily distinguish between different parent materials.Additionally, where present, we always described the aeolian cover as one horizon.The cover had traces of organic matter throughout the whole horizon and was therefore classified as the 1AC horizon.(2) Marine horizons buried below an aeolian cover were not assigned the typical "b" suffix for buried horizons.This facilitated grouping of comparable horizons independent of morphological position.In the same way we did not include the suffix "k" for horizons with secondary carbonates, although these were present in most subsurface marine horizons and partly in marine A horizons.(3) As silt-enriched horizons occurred in most soil profiles, the Bl horizon classification as proposed by Forman and Miller (1984) was used to distinguish these horizons.The lower-case italic L (l) indicates the pedogenic accumulation of silt.This designation was only applied when silt was not only present on top of the clasts but also filled the matrix.These horizons conform to stage 5 or 6 in the silt morphology classification of Forman and Miller (1984).The three deviations mean that two typical soil profiles in the area (e.g.Fig. 4) were described as 1AC-2A-2Bl-2BC and 2A-2Bl-2BC instead of the formal 1Ah-1ACh-2Ahb-2Bkb-2BCkb and Ah-Bk-BCk.
Soil pits were dug until the unaltered parent material was reached or until further digging was not possible.Each major soil horizon was sampled.Bulk density was measured in the field using a 100 cm 3 bulk density ring.For horizons with predominantly gravel, it was not always possible to completely fill the ring by hammering it into the soil.In these cases, the bulk density ring was manually filled up with soil material, which may have led to an underestimation of bulk density.Field bulk density measurements were corrected for the moisture content, which was determined by drying samples overnight at 105 • C. Samples were subsequently drysieved into three grain size fractions: gravel (> 2 mm), sand (2-0.063mm) and silt and clay (< 0.063 mm).Organic matter content was determined by loss on ignition.Samples were heated to 550 • C for 3 h.
A three-factor ANOVA without interactions was used to test the effect of explanatory variables (terrace level, morphological setting and horizon type) on dependent variables (gravel fraction, sand fraction, organic matter fraction and the logarithm of silt fraction).The use of ANOVA was justified, as the Shapiro-Wilk test indicated a normal distribution of the residuals of a linear model between all explanatory variables and the individual dependent variables.The linear model was used to explain which part of the variation in soil properties could be attributed to each of the explanatory variables.

Soilscape model LORICA
The soilscape model LORICA was used to simulate joint soil and landscape development.This raster-based model simulates lateral geomorphic surface processes together with vertical soil development (Temme and Vanwalleghem, 2015, Fig. 2).Transport and change of sediments and soil material are based on a mass balance of various grain size classes.
The model setup that was used in this study contained 10 soil layers in every raster cell of 10 m × 10 m, each with an initial thickness of 0.15 m each.This created an initial thickness of marine sediments (the soil parent material) of 1.5 m.Only three grain size classes were simulated: gravel (> 2 mm), sand (2-0.063mm) and the combined silt and clay class, from now on called silt (< 0.063 mm).
In model simulations, different processes change the mass of material in each grain size class in each soil layer.Using a bulk density pedotransfer function, this change in mass and composition of soil material is translated to a change in layer thickness and a corresponding change in surface altitude.Geomorphic processes are oriented laterally and only affect the top soil layer.Pedogenic processes are oriented vertically and alter material or transport material from one soil layer to another.Some of LORICA's original soil process formulations were adapted to match our hypothesis of the main processes occurring in marine terraces.Some other processes were assumed less relevant, based on the literature and exploratory fieldwork.Hence, they were deactivated for this study.
Chemical weathering was also not activated.However, it is important to note that chemical weathering in the form of dissolution does occur in the marine terraces (Mazurek et al., 2012) and constitutes a source of sand and silt in Arctic soils elsewhere (reported from the west of Spitsbergen; Forman and Miller, 1984;Ugolini, 1986).However, it is not clear to which extent dissolution contributes to in situ weathering on the marine terraces specifically, where physical weathering also plays a dominant role.Only physical weathering was activated in LORICA.Since dissolution mainly focuses on fine material (Courty et al., 1994), a possible overestimation of the finer fractions, relative to the coarse fractions, would be an indicator of the importance, and could possibly hint at the rate of dissolution.

Model framework
A digital elevation model (DEM) with a cell size of 10 m × 10 m served as input landscape.For trough positions, the thickness of the 1AC horizon following from a trend with age was subtracted from the DEM to simulate altitude before aeolian deposition started.A part of the upslope area was included in the DEM to allow simulation of transport of sediments into the study area.Climatic data required by LORICA are annual precipitation (P ) and evapotranspiration (ET).As we did not have data on the palaeoclimate of the study area, we assumed a constant precipitation and evapotranspiration over the entire model run.The same goes for rates and parameters of the simulated processes (Table 1).Annual precipitation is 150-200 mm (Rachlewicz, 2009;Rachlewicz et al., 2013;Strzelecki, 2012).We assumed that a large fraction is lost to infiltration, evaporation and sublimation, leaving 50 mm for overland flow.The initial composition of the marine parent material was derived from field observations and is 95 % gravel with 5 % sand.
To reflect isostatic rebound, a growing part of the landscape was exposed to process calculations as time progressed.Our results from geochronology of the terraces were used to inform this.Simulations started at the time when terrace level 6 was completely above water and progressed with an annual time step.
The activated processes and modifications to them are described below.Where applicable, the calculation of parameter values is also described.

Geomorphic processes
LORICA generates run-off and infiltration by applying precipitation and snow melt to the grid cells.Run-off flows downhill, potentially eroding and collecting sediment on its way.Deposition starts when the amount of transported sediment surpasses the sediment transport capacity of the water.Undeposited sediment is transported out of the study area to the Ebba River and Petunia Bay.Vegetation protection and surface armouring by coarse grains decrease the mass of material that can be eroded.A more extensive explanation of this landscape process is provided in Temme and Vanwalleghem (2015).Standard parameter values were used for almost all parameters describing this process, except for the bio-protection constant.This dimensionless parameter was set from 1 to 0.5 because of the scarce vegetation in the study site.
For aeolian deposition, a simple linear process description was implemented that added a constant amount of aeolian material to all cells in trough positions for every time step.Ridge positions received no aeolian deposition.The aerial photograph (Fig. 1), aggregated to the raster cell size of the input DEM of 10 m, was used to distinguish between ridge and trough positions.
The annual volume of aeolian deposition per cell surface (m 3 m −2 a −1 , or m a −1 ) was calculated by regressing observed aeolian (1AC) horizon thickness to soil age.The bulk density of aeolian deposits, undisturbed by current vegetation, was measured in the field and used to convert the volume to mass.The initial grain size distribution of aeolian de- posits was calculated by extrapolating trends in sand and silt fractions of 1AC horizons with age to time step 0.

Pedogenic processes
Pedotransfer functions are used to estimate unknown variables from readily available soil data (McBratney et al., 2002).Because LORICA's original pedotransfer function for bulk density (BD l , g cm −3 ) is unsuitable for clast-supported soils, we estimated a new pedotransfer function based on the gravel and sand fractions and depth under the soil surface (m).Soil horizons from both marine and aeolian parent material, where bulk density and particle size distribution were known (n = 62), were used to estimate the parameters of this function (Eq.2).The pedotransfer function was validated using leave-one-out cross-validation on the 62 soil horizons (RMSE = 0.183 g cm −3 ).
Physical weathering in LORICA for the various grain size classes i is described as where the change in mass due to physical weathering M pw in layer l is a function of the mass present in the grain size class M i,l , depth below the surface depth l and the median grain size of the fraction size i (Temme and Vanwalleghem, 2015).With parameter C 5 at its standard value of 5, weathering increases with increasing grain size.
Weathering rate C 3 and depth decay parameter C 4 were parameterized from field data.To calculate these parameters, we assumed that a change in gravel fraction in the subsoil is only due to physical weathering.In contrast, topsoil horizons were assumed to be also affected by geomorphic processes.First, weathering rate C of gravel in 2Bl and 2BC horizons was derived from the decay in gravel fraction using with gravel t,l , gravel 0,l and C gravel,l as gravel fraction at time t (-), initial gravel fraction (-) and weathering rate of gravel in horizon l (a −1 ) respectively.The log function in Eq. ( 4) followed from Eq. ( 3), where weathering results in an exponential decay of the mass of a certain grain size class.Second, depth decay parameter C 4 was derived using the differences in weathering rates and average depths between the Bl and BC horizons.
With the depth decay constant C 4 , the weathering rate at the soil surface (C gravel,0 ) was derived, and weathering rate C 3 was calculated using Eq.(3).
Silt translocation was simulated using LORICA's formulation for clay eluviation (Temme and Vanwalleghem, 2015), but a depth decay factor was introduced to better simulate the belly shape of the silt profiles in the soil.
The values for the maximum silt eluviation in a completely silty sediment and the depth decay factor were determined using manual inverse modelling (i.e. through model calibration), using 40 runs with different parameter values.Simulated silt profiles were compared with observed silt profiles for four representative soil profiles in the field.The objective function for calibration was to minimize the average root mean squared error (RMSE) between the modelled and simulated silt fraction for 5 cm thick layers over the entire depth of the profile.

Model validation
Model results were validated using site-and horizon-specific field observations of the gravel, sand and silt fraction, matched to their location in the simulated soilscape.Because of the small amount of observations, observations used for parameterization and calibration were also used in the validation.
The mean prediction error (ME) was calculated to assess a bias between field measurements and model results.The RMSE was calculated to measure the difference.Normalized ME (ME n ) and RMSE (RMSE n ) were calculated by dividing the ME and RMSE by the average observed value (Janssen and Heuberger, 1995).For the mass fractions this was done for every profile, over the depth of observations available for that profile.For the mass content this was done by considering locations on a certain morphological position together.

Geochronology
The three OSL samples taken in the marine sediments show increasing age with increasing terrace level (Fig. 3, Table 2).Datings of the first main terrace level show an age of 4.4 ± 0.2 ka.The highest part of terrace level 3 has been dated to 7.3 ± 0.4 ka.Terrace level 5 dates to 12.8 ± 1.1 ka.
These results support the radiocarbon datings of Long et al. (2012) that covered terrace levels 1 to 4 (Fig. 3).The combined sets of ages (a) show a clear relation with altitude (m), which we approximated with a quadratic trend (Fig. 3): altitude = 2.69 × 10 −7 • age 2 − 0.64. (5) This trend was used to inform isostatic rebound in LOR-ICA.The offset of the trend suggests that the youngest soils are approximately 1500 years old.Following the trend, the youngest soils on terrace level 6 are 13 300 years old.This age was hence used as the start of our model simulations, when terrace level 6 was completely above water.The age ranges of the different terraces increase with altitude (Table 3).However, there is some overlap in ages of different terraces.

Soil types and properties
Although all observed soils can be classified as Cryosols, there are still distinct differences between ridge (8 soils) and trough (22 soils) positions.Ridge positions are generally well drained and usually contain Skeletic Cryosols (observed 7 times).All ridge soils have accumulated secondary carbonates (calcic).Trough soils are typically vegetated and therefore capture aeolian sediment which forms into an aeolian AC horizon, which is rarely present on ridges (Fig. 4).
The aeolian cover displays darker colours due to moderate content of organic matter in 16 trough soils.Younger trough soils are skeletic, due to their thinner aeolian cover and less weathered 2A horizons (6 times).Older soils display an endoskeletic horizon (12 times).For four trough soils, the skeletic properties are absent, due to very thick combined A horizons.Cryoturbation features like frost heave and patterned ground were visible in some trough positions (Turbic Cryosols).
In general, sand and OM fraction decrease towards deeper lying horizons (Table 3, Fig. 5).OM shows small variation, considering the standard deviations.Silt fraction shows the highest values in 2A horizons, and is lower in lower-lying horizons.The decrease in silt fraction in 1AC and 2A hori-Table 3. Average and standard deviation (between brackets) of soil properties of sampled horizons.Horizons with errors in sampling were left out.Counts indicate the amount of observed horizons, the total number of samples can deviate from these numbers (see Fig. 5).The age ranges of the individual horizons were derived from the minimum and maximum altitude per terrace level and their relation with age (Fig. 3).Carbonate content was estimated in the field according to FAO (2006)  zons and increase in silt fraction in 2Bl horizons with increasing terrace level (Table 3, Fig. 5) indicate transport of silt from surface layers to lower-lying layers.This was also evident from field observations, where silt caps were located on clasts in the subsoil and the 2Bl horizons showed an enrichment of silt throughout the whole horizon (Fig. 4).This enrichment was not homogenous, as occasionally bands of higher silt content could be identified inside the 2Bl horizons.Carbonate content also increases with depth  Part of the variation within soil profiles is explained by soil horizon and terrace level (Fig. 5).The three-factor ANOVA confirms the significant effect (P < 0.05) of the soil horizon and terrace level, as well as morphological setting on the variation in gravel and sand fraction.However, terrace level was not a significant explanatory variable for variation in the log-silt and organic matter fraction.Here only soil horizon and morphological setting were significant in explaining part of the observed variation.A linear model involving all three factors resulted in adjusted R 2 of 0.82, 0.84, 0.55 and 0.52 for gravel, sand, log-silt and organic matter fraction respectively.

Process parameters
Slope of the linear regression between 1AC horizon thickness and age is 2.41 × 10 −5 m a −1 (R 2 = 0.33, p value = 0.007).Multiplying this by bulk density of buried aeolian material gives a deposition rate of 0.040 kg m −2 a −1 .Initial sand and silt fraction of the aeolian deposits are 84 and 16 % respectively.
Following the procedure described in Sect.3.3.3,the calculated physical weathering rate of gravel at the surface (C gravel,0 ) is 3.25 × 10 −5 kg kg −1 a −1 .This corresponds to a weathering rate C 3 of 1.01 × 10 −5 kg kg −1 a −1 , when considering the size-dependent correction factor.The corresponding depth decay constant C 4 is −1.63 m −1 , which means that weathering rate decreases by about 80 % per metre under the soil surface.Calibration of silt eluviation resulted in a maximum eluviation of 0.15 kg and a depth decay factor of 5 m −1 .

Simulated landscape and soils
Model results show that the only significant changes in altitude besides uplift are due to aeolian deposition, with a maximum deposition of 0.48 m, divided over 1AC horizons of max ∼ 0.3 m and silt that eluviated from them into lower horizons, contributing the other ∼ 0.15 m.Simulated altitude change due to erosion and sedimentation is negligible, with amounts of several millimetres.Altitude changes are larger on older terraces.There is a clear distinction between changes for trough and ridge positions, because the latter did not receive aeolian input (Fig. 6).
Variation in simulated profile curves of different particle sizes is mainly caused by morphological position of the soils (Fig. 7).Although the general shapes of these profiles correspond to the mean observed profiles, observed profile curves show a larger spread than simulated profile curves.Observed gravel fractions are lower than simulated fractions.Sand and silt fractions and mass were larger in the field than in the model results (Fig. 8).The silt fraction in the topsoils on both ridge and trough positions is higher in the field than in the model results.
Most accurate predictions for sand and silt fractions and contents are for trough positions (Table 4, Fig. 8).For gravel, ridge positions are predicted most accurately.The relatively high RMSE n values indicate that there is a large spread between modelled and observed mass fractions and contents (cf.Fig. 7).On the other hand, ME n s indicate a low bias in some of the predictions.Examples are sand and silt properties in trough positions, gravel and sand properties in ridge positions, and total mass of soil material in all positions.The positive ME n for total mass of the soil shows that the In some places, the morphological position as derived from the aggregated aerial photograph and used in the model differs from field-observed morphological position, due to small-scale variation between ridge and trough positions.These "mixed" positions (Table 4 and Fig. 8) show the highest differences between observations and simulations and cause the largest errors in the validation statistics.Small differences between RMSE and ME for the mixed positions indicate that the largest part of the error is systematic, and a relatively small part is caused by a random error.

Geochronology and isostatic rebound
Our new OSL dates and the existing calibrated radiocarbon results from Long et al. (2012) show comparable results for the ages of marine terraces in our study area.The altitude above current sea level can be rather well approximated by a quadratic equation with age (Eq.5).The uplift ranges from 7.2 mm a −1 13 300 years ago (altitude = 47 m) down to   0.8 mm a −1 1500 years ago (altitude = 0 m), with an average rate of 4.0 mm a −1 (Fig. 3).The age ranges of the different terraces show some overlap (Table 3).This is due to differences in elevation between ridge and trough positions, lowerlying rills and inaccuracies in the DEM.Forman et al. (2004) reviewed uplift rates of studies all over Svalbard.For every reviewed uplift curve, total uplift (m) at 9000, 7000 and 5000 uncalibrated radiocarbon years was supplied.In order to work with calibrated radiocarbon ages, the marine reservoir effect of 440 years (Forman et al., 2004) was added to these uncalibrated ages.Next, the ages were calibrated with CALIB 7.0.4(Stuiver and Reimer, 1993), using the MARINE13 calibration curve (Reimer et al., 2013).Uncertainty in the uncalibrated ages was assumed to be 1 % (Walker, 2005, p. 23).The deviation from the standard marine reservoir effect ( R) was assumed to be 100 ± 39 years (Long et al., 2012).The corresponding calibrated ages are 10 157, 7808 and 5690 yr BP.Because abovementioned assumptions and uncertainty in the estimation of the marine reservoir effect, the following analysis should be considered with caution.
From the calibrated ages, uplift rates over the last 10 157 years, from 10 157 to 7808 years ago and 7808 to 5690 years ago, and from 5690 years to present could be calculated.Two of the reviewed uplift curves were located near the Ebba Valley: Kapp Ekholm and Blomesletta (Péwé et al., 1982;Salvigsen, 1984).For Kapp Ekholm, the uplift rates over these time intervals were 4.4, 10.6, 4.7 and 1.8 mm a −1 .For Blomesletta they were 2.5, 5.5, 2.4 and 1.2 mm a −1 .In comparison, uplift rates for these time intervals for our uplift curve were 2.7, 4.7, 3.8 and 1.4 mm a −1 .Uplift rates in the Ebba Valley are thus very comparable to the Blomesletta uplift curves.The rates found at Kapp Ekholm are much higher.This discrepancy can be due to a large uncertainty in the uplift rates, which is not incorporated here.Nonetheless, the uplift curve from this study can be considered reliable, as the independent OSL ages support the earlier radiocarbon ages.
The fitted uplift curve suggests that the youngest terrace is around 1500 years old.That could indicate that uplift has stagnated or reversed, leading to flooding of lower-lying terraces, as is also suggested by Long et al. (2012) and Strzelecki (2012).This corresponds to renewed glacier growth in response to a cooler climate starting 3000 years ago, which eventually led to the Little Ice Age (Rachlewicz et al., 2013;Svendsen and Mangerud, 1997).The exact age of the youngest soils remains uncertain.However, the silt and organic matter profiles of the soils on terrace level 1 (Table 3, Fig. 5) nonetheless suggest that they have been developing for a longer time, instead of just having recently emerged from the sea.
A well-known disadvantage of radiocarbon dating is that older ages (> 35 cal ka) can easily be underestimated by contamination with younger carbon (Briant and Bateman, 2009).This may be the case with the radiocarbon dating of the highest terrace in the Ebba Valley (> 37 000 years ago; Kłysz et al., 1989) but is unlikely for the Holocene terrace sequence that was studied in this paper.More importantly, radiocarbon ages derived from marine fauna (e.g.shells) need to be corrected for the marine reservoir effect.However, this correction not only requires extra analysis (e.g.Long et al., 2012) but typically also shows a large regional and biological variation and thus potential bias (Forman and Polyak, 1997).Another common problem of radiocarbon dating, especially in our geomorphologically very dynamic setting, is re-working of the dated material (e.g.Long et al., 2012).In this case the radiocarbon age might potentially overestimate the deposition age of the marine terrace.OSL does not suffer from these potential malign effects, as OSL provides direct depositional ages of sand-sized marine deposits and can be used to independently validate the radiocarbon chronology.In our case both data sets agree and thus support each other.OSL ages have in general a larger uncertainty interval than the radiocarbon ages (Fig. 3), because OSL methods typically provide a lower precision than radiocarbon dating.The typical OSL uncertainty is 5 to 10 % for the 1σ confidence interval (∼ 65 %) (Preusser et al., 2008), which was also achieved for the samples under investigation.Furthermore, OSL is a wellestablished method to date recent coastal dynamics (e.g.Ballarini et al., 2003;Reimann et al., 2010) and aeolian activity (e.g.Sevink et al., 2013).Thus, when considering chronosequences with a longer age span and where recent geomorphic activity plays a role, OSL can validate radiocarbon chronologies and is a powerful alternative dating method.

Landscape evolution
It was our intention to use the LORICA model to test our field-informed hypotheses about the combined evolution of soils and landscapes in the study area.The aspects that were not well simulated therefore suggest improvements to process understanding.
The main geomorphic aspect that was not well simulated is the presence of small rills incising into the smooth terrace ridges (Mazurek et al., 2012).These were observed on all terrace levels but were not simulated (Figs. 1 and 6).
Model tests indicate that unrealistically high erodibility values would have to be adopted to simulate the amounts of erosion that lead to rill formation in the gravelly soil material under the dry climate in the study site.This suggests that the process that has led to rill formation is not included in the model.We suggest two possible processes.First, permafrost that is present not far under the surface can act as an impermeable layer.Combination of seeping groundwater and overland flow at ridge escarpments can then disaggregate coarse material and remove fine material (Higgins and Osterkamp, 1990).This seepage erosion occurs in cliffs and riverbanks (Fox et al., 2007;Higgins and Osterkamp, 1990), but has not, to our knowledge, been described for marine terrace sequences.Second, occasional heavy storms and high tides in the period soon after uplift above sea level may have caused temporary flooding of a beach trough that was already protected by a beach ridge.Drainage of the trough after a storm passes can have formed the rills.Both processes fit with the observed absence of a clear relation between rill size and age: the conditions that initiate rill erosion would be most prevalent after limited uplift over sea level.
Although these erosion rills occur in most of the terrace boundaries, most water is currently drained parallel to the ridges, towards the tundra lakes.These again drain to the Ebba River or the Petunia Bay.Because the flow velocity through these lakes is very low, no erosion occurs.
The rate of aeolian deposition was estimated from observations without model simulations.However, a complicating factor is that this was based on present soil properties: the rate of aeolian deposition was based on current thickness of aeolian cover.Ultimately, part of the silt was translocated from the simulated aeolian deposits to deeper layers, resulting in an underestimation of thickness of 1AC horizons.This effect is visible, for instance, in the overestimation of gravel and underestimation of sand in the top parts of profiles in trough positions (Fig. 8).
The coefficient of determination of the regression between thickness of 1AC horizons and age (Sect. 4.3) shows that age only explains 29 % of the variance in 1AC horizon thickness.This indicates that other factors such as the initial topography, wind shadow, hydrological properties, variation in surface cover by vegetation and bacterial soil crusts, and reworking of aeolian sediments (e.g.Paluszkiewicz, 2003), which were not considered in our model, also play a role.The spatial heterogeneity of aeolian deposition is visible in shorter-term measurements in the study area.Deposition during the summer periods of 2012-2014 ranged from 3 to 1713 g m −2 per summer season, depending on morphological position and vegetation cover (Rymer, 2015), while niveo-aeolian deposition in the years 2000-2005 ranged between 70 and 115 g m −2 a −1 (Rachlewicz, 2010).In comparison, aeolian deposition in Hornsund, southern Spitsbergen, was 300-400 g m −2 a −1 for the winter of 1957/1958 (Czeppe, 1966).The higher numbers in these measured ranges are about an order of magnitude larger than the average deposition rates found by our observation of aeolian horizon thickness (40 g m −2 a −1 ).Part of this difference can be caused by reworking of recent deposits by continued aeolian activity, and another by our underestimation of the thickness of aeolian deposits by observing them after some of the silt has eluviated from them.On top of the spatial heterogeneity, there is also a large temporal variation in niveo-aeolian deposition rates (Christiansen, 1998).

Soil formation
Both physical and chemical weathering in the Arctic are driven by moisture availability (Hall et al., 2002).Consequently, weathering occurs at a faster rate in the wetter troughs of the terraces.This expected variation in weathering between different morphological settings was observed in particle size distributions (Fig. 7) but not quantified in the model, due to data limitations.In addition, the ANOVA also indicates the significant role of morphological position on gravel, sand and silt fraction.However, it should be noted that silt is also significantly influenced by aeolian influx.
The physical weathering rate was calculated based on the gravel fraction of the subsurface horizons -not of the surface horizon.This was done in an attempt to exclude the effect of other processes on grain size changes.Nonetheless, chemical weathering in the form of dissolution may have affected the grain size distribution in the subsurface.Dissolution mainly affects the fine fraction, as it has a larger reactive surface area (Courty et al., 1994;Ford and Williams, 1989, p. 28).Less fine material in the subsurface means that the fraction of coarse material is higher.This subsequently results in an underestimation of the physical weathering rate, as this is based on the coarse fraction.Independent observations would be needed to disentangle the effects of physical and chemical weathering.For example, a scanning electron microscope (SEM) could be used to identify surface morphologies related to certain weathering processes (e.g.Andò et al., 2012;Mavris et al., 2012).The mineralogical content of the groundwater, together with known temperatures, could provide constraints on dissolution rates (Dragon and Marciniak, 2010).
Consequently, the calculated weathering rate must be considered a combined physical and chemical weathering rate.This rate of 1.01 × 10 −5 kg kg −1 a −1 (1.01 × 10 −3 % a −1 ) is orders of magnitudes lower than in field weathering experiments of the dominantly chemical weathering of granite and dolomite in Swedish Lapland (Dixon et al., 2001;Thorn et al., 2002).These resulted in a weight loss of 0.121 and 0.326 % a −1 respectively for an experiment of 5 years.Two reasons for this difference present themselves: time-decreasing weathering rates and moisture availability.Weathering rates decrease with time, amongst other things due to precipitation of secondary minerals which slow the dissolution process (Langman et al., 2015;White et al., 1996).These secondary minerals were widely observed in Ebba Valley, partially coating gravel (cf.Courty et al., 1994).It is also likely that long-term weathering rates in Swedish Lapland exceed those found in the Ebba Valley because of the much larger precipitation (1750 mm a −1 ; Dixon et al., 2001).Other weathering studies also indicate the dominant control of moisture in determining both physical and chemical weathering rates (e.g.Egli et al., 2015;Hall et al., 2002;Langston et al., 2011;Matsuoka, 1990;Wu, 2016;Yokoyama and Matsukura, 2006).
It can nonetheless be argued that dissolution is a significant process in our marine terraces.This is also apparent in the field observations.The 2A horizons have a lower CaCO 3 content than 2B, 2Bl and 2BC horizons.This is consistent with observations elsewhere in Spitsbergen of a more active dissolution regime in near-surface horizons compared to subsurface horizons (Courty et al., 1994;Forman and Miller, 1984).Dissolution mainly occurs in the fine soil mass due to its larger reactive surface.Consequently, the fine fraction in 1AC and 2A horizons consists of less CaCO 3 .
Dissolution was not included in model simulations.This should have caused overestimation of fine material in the model output.However, model results instead show lower sand and silt contents than observed (Figs. 7 and 8, Table 4).For trough positions, this is partly due to underestimation of the thickness of 1AC horizons (see above).By this underestimation, the thinner 1AC horizons give way to underlying marine deposits in the model outputs.Their higher gravel content distorts the comparison between observed and simulated profiles, leading to an underestimation of sand content and fraction.Silt properties are not influenced by this error, as the eluviation rate was calibrated using valley positions.As a consequence, silt predictions there show a low error.
This disturbance by a wrongly estimated thickness of the aeolian cover is not present in ridge positions.However, for those positions, sand and silt contents are also underestimated.Simulated silt content shows the largest deviation from field observations (Fig. 8, Table 4).This indicates that there is another source of sand and silt in ridge positions.
One source of silt is in situ weathering of coarse material into finer material (Fahey and Dagesse, 1984;Forman and Miller, 1984).Another source is an ex situ one, namely aeolian deposition (e.g.Locke, 1986).This is observed in trough positions, which display a higher silt fraction throughout the whole profile compared to ridge positions (Fig. 7).The het-erogeneous deposition of aeolian material over trough positions on the area has resulted in a heterogeneous silt source for the subsurface.Hence, the ANOVA shows no significant relation between terrace level and silt fraction.
Deposition occurs partly through entrapment with falling snow (Rachlewicz, 2010).Meltwater from this snow partly infiltrates in the permeable gravelly soils.This downward flow of water can transport part of the silt it had captured, increasing the silt fraction in the subsurface.This process can also act as a source of silt in ridge positions.Material left behind on the surface can be reworked by strong summer winds, leaving no trace of aeolian deposits on these positions.This theory contradicts the observations of Forman and Miller (1984), who claim that silt in their studied marine ridges on western Spitsbergen mainly originates from in situ frost weathering and dissolution.They attribute this to the absence of a major source area, absence of silt on the surface and absence of vegetation to entrap sediments.In the Ebba Valley, the proglacial sandur plain acts as a major sediment source.In addition to that, our observations of 2A horizons, which partly reach the surface, show high amounts of silt compared to lower-lying horizons (Table 3).Due to these different conditions in the Ebba Valley, it is possible that ridge positions also profit from aeolian silt input.This also can explain the deficit of silt in model simulations (Fig. 7).
The enrichment of silt from aeolian material, silt caps located on top of clasts and heterogeneous distribution of silt in the 2Bl horizons suggest that eluviation is responsible for the silt redistribution, instead of frost sorting, as suggested by Bockheim and Tarnocai (1998).Addition of a depth decay function in the simulation of silt translocation proved to result in simulated profiles comparable to the observed profiles (Fig. 7).This decay in eluviation rate with depth represents the limited water flux in the subsurface.Due to a shallow unfrozen soil in the snow melt season, the water cannot reach deeply into the soil and starts to flow laterally.Later in the season, when the active layer is thawed, the water flux is limited, because precipitation is very limited.More detailed descriptions of the silt content throughout the whole profile can help to better understand the silt dynamics in these Arctic environments.
The soil types and properties observed in this study agree with studies in comparable settings on Svalbard.Cryosols associated with stony and gravelly parent material can be found on the uplifted terraces and Turbic Cryosols can be found in loamy material on the southwest coast of Spitsbergen (Szymański et al., 2013).Marine deposits elsewhere are also covered with a layer of aeolian deposits (e.g.Alexanderson et al., 2013).Dissolution is a main weathering process, which precedes or runs parallel with silicate weathering (Mann et al., 1986).The process of translocation of silt resulting from weathering and the formation of silt-enriched horizons was observed in northwestern Spitsbergen (Forman and Miller, 1984;Mann et al., 1986).

Temporal and spatial soil-landscape interactions
The variation in simulated soil profiles is smaller than the variation in observed soil profiles (Fig. 7).We presume that this is due to a variation in boundary conditions that is not captured in the model.Particularly, in reality, the grain size distribution of the original parent material is not constant but can vary spatially due to different sorting and transport mechanisms driven by the large spread in particle size of gravels (Buscombe and Masselink, 2006).Forman and Miller (1984) suggest a gravel fraction of 85 to 90 % for their marine parent material, which is lower than the 95 % used in this study.However, the coarse fraction of the CB horizons studied by Forman and Miller (1984) ranges from 54 to 99 %, indicating a large spatial variation in texture.Also, aeolian deposition rates vary both in space and time, as is discussed in Sect.5.2.
The glacial retreat since the beginning of the Holocene indicates the general warming climate in Svalbard, which also included several colder episodes (Skirbekk et al., 2010;Svendsen and Mangerud, 1997).Precipitation and evaporation rates were also not constant during the Holocene.Courty et al. (1994) show that characteristics of secondary carbonates on subsurface clasts indicate different climatic and biogenic episodes.Although their research was conducted on western Spitsbergen, which has a much wetter climate than central Spitsbergen, it is likely that variations in climate occurred over the complete Svalbard archipelago.In addition to these long-term trends, short-term extreme rainfall events can also significantly alter the landscape, e.g. through mass movements (André, 1995).These changing climatic conditions thus influence rates of soil and landscape formation, and can therefore have a significant effect on soil properties (e.g.Hall et al., 2002;Jiang et al., 2011;Oliva et al., 2014;Van Vliet-Lanoë, 1998).Nevertheless, these changing climate conditions were not included in our model for the Ebba Valley, because slopes are not steep enough for mass movements and, to our knowledge, there was no information on the Holocene climate for the dry Ebba Valley.
The elevation differences between ridges and troughs in the study area are the main driving force for spatial soil heterogeneity.Water accumulates in the relatively sheltered trough positions.Consequently, weathering occurs at a faster rate and there is a higher vegetation density.These plants capture a part of the aeolian sediment that has been deposited with freshly fallen snow.The resulting 1AC horizon, consisting of finer material, holds more water than the marine sediments, resulting in more plant growth.This feedback has resulted in a local aeolian cover of up to 70 cm.Note that this process is not expected to continue over longer timescales, because the aeolian sediments will ultimately grow out of their sheltered and, more importantly, humid positions, becoming susceptible to reworking by wind erosion.
Ultimately, Arctic soil development is not as straightforward as we hypothesized in the beginning of this paper.The interplay between different processes, known and unknown,

Conclusions
This study combined different methods to study soil development on a series of marine terraces in central Spitsbergen.
The analysis of the results led to the following conclusions: -The changes in soil properties of the gravelly soils on the marine terraces can be attributed to different soilforming processes, such as physical (frost action) and chemical weathering (dissolution) and translocation of silt.Dissolution mainly occurs in A horizons developed in marine material.Translocation of silt occurs everywhere in the landscape, following the water flow.
-Optically stimulated luminescence (OSL) results support earlier radiocarbon dates from the area.Moreover, an uplift curve constructed based on both types of dates concurs with a nearby uplift curve, indicating the potential of OSL for measuring uplift rates in this setting.Combining these datings with field observations enabled the calculation of process rates using field observations.
-However, determining historical rates of weathering and aeolian deposition using current soil properties is difficult when multiple processes have influenced those properties.Especially dissolution, which removes material from the soil, distorts the mass balance of soil constituents that was used to calculate rates.
-Simulation of soil development in a landscape context with the soilscape evolution model LORICA was successful in terms of simulating trends in soil properties.However, there were significant discrepancies between field observations and model results.The larger variation in field observations than in model simulations is likely due to spatially and temporally varying boundary conditions that were not included in simulations.More importantly, bias in model outcomes helped to increase our understanding of Arctic soil development in the marine terraces.
-Soil development is heavily influenced by geomorphic processes, mainly aeolian deposition.Deposition acts as a source of fine material, which mainly accumulates in relatively sheltered beach trough positions.However, our results are consistent with suggestions that aeolian silt has also been added to soils in beach ridge positions.Erosion of overland flow plays a minor role compared to erosion by extruding groundwater or by the effects of storms on young terraces.

Figure 1 .
Figure 1.Aerial photograph (summer 2009) of the study area indicating the 6 terrace levels (in white numerals), 3 OSL-dating locations (in roman numerals), 30 sampling locations (in black numerals) and the approximate location of the radiocarbon datings done by Long et al. (2012).Their corresponding dates can be found in Fig. 3. Ridges are recognizable as light (unvegetated) parts of the terraces, while troughs are darker (vegetated).Disturbed areas such as erosion rills, permanently wet depressions and tundra lakes were excluded from the study area.The inset shows the location of the study area on the island of Spitsbergen.

Figure 2 .
Figure 2. Conceptual framework of LORICA as used in this study.

Figure 3 .
Figure 3. Elevation and age of OSL dates (this study) and radiocarbon dates (Long et al., 2012).Ages are displayed with a confidence interval of 2σ .The black line shows the quadratic regression between altitude and age of each dating.The symbols A-G and I-III refer to the respective locations in Fig. 1.

Figure 4 .
Figure 4. Examples of a typical soil found on a trough (left panel) and ridge location (right panel).The prefixed numbers indicate the parent material: aeolian (1) or marine (2).

Figure 5 .
Figure 5. Box plots of observed soil properties on different terrace levels, ordered by main soil horizon.Counts indicate the occurrence of the displayed properties.Number of values can differ per soil property due to measurement errors.

Figure 6 .
Figure 6.Simulated altitude change in the study area.A clear difference is visible between ridge positions (black grid cells) and trough positions (grey shading), due to the absence of aeolian deposition on ridge positions.

Figure 7 .
Figure7.Total variation and mean of particle fractions in observed and simulated profile curves, divided over the morphological setting.For every centimetre along the soil profile depth, the minimum, maximum and mean mass fraction of the various grain sizes for all profiles in the considered morphological setting are displayed.

Figure 8 .
Figure 8. Scatter plots of simulated versus observed mass in kilograms over the total observed depth, for different particle sizes and morphological positions.The black line indicates the 1 : 1 line, which indicates a perfect match between model and field results.
W. M. van der Meij et al.: Arctic soil development on a series of marine terraces on central Spitsbergen

Table 1 .
Settings and parameters used as input for the LORICA model.

Table 2 .
OSL ages with uncertainty of 1σ for 3 samples taken in the study area (Fig.1).

Table 4 .
Normalized RMSE and ME as validation statistics of the model results, ordered per geomorphic position.Mixed positions are locations where the model follows a different geomorphological setting than was observed in the field.This is due to loss of details with aggregation of the ridge-trough map to the cell size of the input DEM.
www.soil-journal.net/2/221/2016/SOIL, 2, 221-240, 2016 together with variations in initial and boundary conditions in soil and landscape development has resulted in a complex soil-landscape system.Additional research would be required to further unravel soil and landscape development in this fragile environment, especially in the context of a changing climate.