The Cosmolian program for simulating aeolian dynamics and its application to central Australia

The wide spatial coverage of sand dunes in continental interiors makes the understanding of their activity and accumulation history valuable for palaeoenvironmental reconstructions and the interpretation of landscape evolution. Nevertheless, the study of aeolian landscape development at the million‐year timescale is hampered by the complex interaction of factors determining dune migration and the inherently self‐destructive nature of their chronostratigraphy, thus limiting the applicability of traditional dating methods. This study presents a standalone program that simulates aeolian transport based on luminescence‐derived chronologies coupled with numerical modelling of the accumulation of cosmogenic nuclides. This integrative approach to modelling the history of aeolian landforms reveals phases of emergence of aeolian sand into the landscape, and provides a data‐based scheme that facilitates the morphodynamical study of aeolian processes over multiple timescales and up to several millions of years. The application of the program for reanalysing previously reported data from the Australian Simpson Desert reveals multiple pulses of sand dispersion into central Australia at 3.8–3.4, 2.9–2.5 and 1.5–1 Ma, corresponding to pronounced changes in climatic conditions and landscape deformation events. The synchronicity of the results with the established environmental framework that would promote the production and aeolian distribution of sand exemplifies the applicability of process‐based modelling in constructing a timeframe of key landscape evolution events in arid environments by studying aeolian deposits. The dependence of the parameters used to determine environmental settings on sand transportation patterns additionally makes the program a powerful tool to further investigate the triggers and mechanisms of aeolian processes.

traditional dating methods. This study presents a standalone program that simulates aeolian transport based on luminescence-derived chronologies coupled with numerical modelling of the accumulation of cosmogenic nuclides. This integrative approach to modelling the history of aeolian landforms reveals phases of emergence of aeolian sand into the landscape, and provides a data-based scheme that facilitates the morphodynamical study of aeolian processes over multiple timescales and up to several millions of years. The application of the program for reanalysing previously reported data from the Australian Simpson Desert reveals multiple pulses of sand dispersion into central Australia at 3.8-3.4, 2.9-2.5 and 1.5-1 Ma, corresponding to pronounced changes in climatic conditions and landscape deformation events. The synchronicity of the results with the established environmental framework that would promote the production and aeolian distribution of sand exemplifies the applicability of process-based modelling in constructing a timeframe of key landscape evolution events in arid environments by studying aeolian deposits. The dependence of the parameters used to determine environmental settings on sand transportation patterns additionally makes the program a powerful tool to further investigate the triggers and mechanisms of aeolian processes.

K E Y W O R D S
aeolian, Australian landscape evolution, coupling cosmogenic nuclides and luminescence, databased modelling, process-oriented modelling, Simpson Desert, cosmogenic nuclides, OSL, dating sand dunes 1 | INTRODUCTION An understanding of sand dunes dynamics has been long and widely sought out, because of its importance in deciphering terrestrial environmental records in arid environments, which are usually devoid of temporally and spatially continuous geological archives (Bagnold, 1936;Bowler, 1976;Fitzsimmons et al., 2013;Grove, 1969;Haynes, 1982;Lancaster, 1997;Yang et al., 2012). The growing interest in arid aeolian landscapes, which currently cover $5% of Earth's land surface, stems mostly from the fact that non-aeolian climatic and palaeoenvironmental terrestrial records are usually rare (Chase, 2009;Telfer & Hesse, 2013). However, the synergistic interaction between multiple factors that govern dune development (e.g. sediment supply, susceptability to deflation and wind transport capacity) results in non-linear responses, characterized by a complexsystem behaviour (Kocurek & Lancaster, 1999;Werner, 1995). Thus, numerical models have been used to describe the relationships between various components that affect dune mobility under different climatic settings (e.g. Ashkenazy et al., 2012;Bailey & Thomas, 2014;Thomas et al., 2005;Wiggs, 2001;Yizhaq et al., 2007).
Nevertheless, such models do not provide chronological constraints over a million-year timescale, thus limiting their application for interpreting long-term processes (Chase, 2009;Fujioka & Chappell, 2011;Thomas & Bailey, 2017;Tooth, 2012). Therefore, the development of new approaches for quantifying aeolian activity over the million-year timescale remains crucial for determining the response of aeolian landscapes to a changing environment.
The reworking of sand during aeolian transport, and the presumed antiquity of desert dunes worldwide, poses a challenge for dating aeolian landforms with commonly used luminescence dating techniques because they (a) require continuous burial and (b) are limited to several hundreds of thousands of years due to the saturation of the luminescence signal (Faershtein et al., 2020;Fitzsimmons et al., 2013;Fujioka & Chappell, 2010;Stone & Fenn, 2020;Telfer & Hesse, 2013). Therefore, the concentrations of cosmogenic nuclides, which are not being reset upon exposure, may be applicable for dating aeolian landscapes at the million-year timescale (Klein et al., 1986;Vermeesch et al., 2010). However, because the production of cosmogenic nuclides depends on the proximity of the accumulating sandy substance to the surface, a 'tailor-made' method which accounts for repeated burial and exposure cycles is required for interpreting the measured concentrations of cosmogenic nuclides in aeolian landforms (Fujioka & Chappell, 2011;Tooth, 2012).
Given these challenges, there is an urgent need for a tool that enables data-based modelling of aeolian transport chronology over the million-year timescale, which is provided in this paper. To address this need, we present a standalone program named Cosmolian, built for running and interpreting an integrated luminescence-cosmogenic nuclide accumulation model. We demonstrate the applicability and capabilities of the program by modelling published data from the Australian continental interior as a case study. The numerical model encompassed in Cosmolian relies on the method presented by Vainer et al. (2018) that provides a statistically based scheme to reconstruct the chronological development of aeolian dunefields. This approach is set to reproduce both 26 Al and 10 Be concentrations measured in sand by simulating the accumulation of cosmogenic nuclides during aeolian transport, where displacement rates are based on luminescencederived chronologies. We incorporated the model into a graphical user-friendly interface for intuitive use that does not require programming knowledge (Figure 1). It allows the users to simulate and analyse their own data and to conceptualize and test a wide range of scenarios in order to investigate the response of aeolian landforms to various boundary conditions, expressing different environmental settings.
Furthermore, Cosmolian's functional versatility allows it to be adapted and applied to any dunefield around the world.
Simulations were performed using Cosmolian with data from the Simpson Desert dunefield located within the continental interior of Australia ( Figure 2) as input (Table 1). This arid landmass, which lacks substantial topographic barriers to climatic systems, accommodates vast dunefields that were subjected to considerable climatic fluctuations with a poorly constrained chronology at the million-year timescale (Fitzsimmons et al., 2013;Fujioka & Chappell, 2010;Kreig et al., 1990). The Simpson Desert dunefield lies within the greater Lake Eyre Basin (LEB) that covers $40% of the continent, and although the spatial coverage of the considered samples is limited (five dunes that are situated along a $20 km transect), the interpretation of the modelling results is extrapolated to the LEB scale as similar climatic conditions and aeolian activity are assumed to operate collectively throughout the basin (Fitzsimmons et al., 2013). The simulated outcomes were compared with the independently established chronological framework of landscape and climate evolution in the LEB in order to estimate the validity of Cosmolian's results and to point to possible triggers capable of initiating and reviving aeolian processes (see the section entitled 'Phases of aeolian activity and their possible triggers'). Such an examination enables the inspection of some of the features derived from the program, as well as testing their compatibility with the current understating of aeolian structures and mechanisms.

| MODEL FRAMEWORK AND PARAMETRIZATION
The Cosmolian program is based on the modelling approach of Vainer et al. (2018) that includes a pre-aeolian step of cosmogenic nuclides build-up, followed by the main simulation stage of aeolian sand migration. To determine realizations of simulations, the program requires the concentrations of in-situ 26 Al and 10 Be measured in aeolian sand as input. The production rates of these radioactive nuclides by cosmic-ray bombardment entering the quartz lattice in sand grains depend mostly on the elevation of the study site, the depth and density of overburden (if the target material is buried), and the geomagnetic latitude field strength (Lal, 1991). Due to the uncertainty in determining these variables (especially for migrating sediments), Cosmolian allows multiple input values for the different parameters to define a range of boundary conditions for simulating the development of the studied aeolian landscape ( Figure 1; Table 1).
Pre-aeolian concentrations of 26 Al and 10 Be that are produced due to erosion at source areas are calculated by assigning the parameters of a presumed source area that the sand was derived from. For the Simpson Desert dunefield, a source area composed of the proportional contribution of terranes deduced by Pell et al. (2000) was constructed. This was done by determining and averaging the geographical provenance terranes based on a 9-s DEM using GIS software (Hutchinson et al., 2008). Additionally, four erosion rates were assigned to cover several possible palaeoclimatic scenarios, following the estimations of Bierman and Turner (1995), Fujioka et al. (2009), Heimsath et al. (2001 and Vasconcelos et al. (2008). During transport, after each 20 cm displacement, the grains reside for a randomly sampled duration, based on a predefined distribution of vertical displacement rates (VDRs). Empirical VDR datasets were constructed by converting luminescence ages of selected dunes from different environmental settings (globally distributed, Australian dunes and Simpson Desert dunes) published by Lancaster et al. (2016). This is done by dividing the luminescence age of the sample by the sample's depth. An additional set of adjusted VDR datasets was constructed by removing the fastest 10% from each empirical VDR dataset, accounting for the possible bias stemming from oversampling young and shallow dune deposits (Chase, 2009;Fitzsimmons et al., 2013;Hesse, 2016). The OSL dose rates and 26 Al and 10 Be concentrations (Table 1) were measured by Fujioka et al. (2009) in sand collected from boreholes reaching the bases of five dunes intercalated with palaeosols in the Simpson Desert ( Figure 2). In the case that several samples were collected within bordering palaeosol horizons of the same age estimation (Fujioka et al., 2009), we averaged cosmogenic nuclide concentrations of samples that yielded the same 26 Al/ 10 Be burial age (within error). Thus, we reduced the initial 12 samples, which were used by Fujioka et al. (2009) to estimate a minimum age of $1 Ma for the onset of a regional phase of dune build-up, into five representative datapoints that were used to evaluate older episodes of sand distribution. Similarly, we averaged OSL dose rates reported from corresponding depths between the sampled dunes.
Simulations were carried out for 23 Ma, which is the oldest presumed age for the aridification of the Australian continent (Vasconcelos et al., 2008). The model was run 10,000 times for each combination of parameters, resulting in 960 ,000 simulations for each of the five data points. Scenarios that resulted in a convergence of simulated and measured cosmogenic nuclide concentrations were then interpreted as plausible transportation histories that were integrated to decipher the temporal development of the dunefield. These results are discussed below, within the chronological framework of environmental and landscaping events.

| Phases of aeolian activity and their possible triggers
The most striking output of the simulations is the distribution of the simulated timespans required until the simulated 26 Al and 10 Be T A B L E 1 Model parameters used to simulate the aeolian landscape development of the Australian Simpson Desert. Geographical settings during the pre-aeolian and aeolian stages were calculated from a 9-s DEM as the averages for the source areas of the sand and the Simpson Desert, respectively. Names of samples indicate the palaeosol horizons described by Fujioka et al. (2009) and are ordered in decreasing stratigraphical order, defined by the palaeosols. Cosmogenic nuclides and OSL data are the average concentrations and the weighted errors of 12 and 25 samples, respectively, reported by Fujioka et al. (2009). Vertical displacement rates were calculated from luminescence ages of noncoastal dunes published by Lancaster et al. (2016), divided by the depth of the sample (see Vainer et al., 2018  concentrations matched those measured in sand samples from various depths (Figure 4). These are interpreted as pointing to the most probable timing of key aeolian activity events that coalesce into three major pulses of sand production and/or exhumation. The events revealed by the model coincide with shifts in regional environmental conditions that could have led to the production and distribution of sand, and agree with previous observations indicating several phases of aeolian activity since the Pliocene (Kreig et al., 1990). We further note that outliers stem from simulations that were run using the most This period also coincides with global climate changes (Lisiecki & Raymo, 2005) that are expressed as humid conditions in the periphery of LEB (Christensen et al., 2017;Fujioka & Chappell, 2010) and intense continental-scale weathering (Vasconcelos et al., 2008). Such conditions would have triggered the production of loose material available for aeolian distribution. Accordingly, initial stripping of soil mantles from silcrete that developed on the south-central Australian tableland is estimated at $4 Ma (Fujioka et al., 2005), concurrent with the oldest recorded formation of regolith in the northeastern parts of the continent (Bird et al., 1990). The exposure of detritus at plausible source areas (Pell et al., 2000), including the contemporary desiccated extensive lakes within the LEB (Alley, 1998), was thus likely to provide loose material (Fitzsimmons et al., 2013) for the earliest aeolian deposits in the LEB.
Revitalization of aeolian activity was modelled by Cosmolian to 2.9-2.5 Ma, coeval with the continental transitional climatic interval that was followed by the establishment of modern dust transportation pathways (Christensen et al., 2017 (Fujioka et al., 2005), as well as the northeastern margins of the LEB, where a second episode of regolith formation occurred (Bird et al., 1990).
Similar to Australia, increased aeolian activity of early Pleistocene age in southeast Asia and southern Africa (Ding et al., 2002;Matmon et al., 2012;Vainer et al., 2018)  The youngest aeolian phase in central Australia, as revealed by Cosmolian, took place between 1.5 and 1 Ma, during which changes in the periodicity and amplitude of orbital cycles were identified (Lisiecki & Raymo, 2005). This happened simultaneously with aeolian sedimentation in southeastern Australia (McLaren et al., 2011) coupled with a more than twofold increase in dust flux input into the Southern Ocean (Martínez-Garcia et al., 2011). Additionally, the transition from pluvial conditions in marginal LEB lakes into dry playas that occurred between 1.6 and 1 Ma (Chen & Barton, 1991;Chivas et al., 1986) could have possibly contributed sediments for aeolian distribution.
The period between 1.5 and 1 Ma was characterized by augmented aridity in central Australia (Fujioka et al., 2009) that could be related to global-scale environmental changes that commenced by $1.6 Ma (Lisiecki & Raymo, 2005). In the Northern Hemisphere, these changes were expressed by declining atmospheric CO 2 concentrations and massive regolith removal that started at $2 Ma and continued for $1 Myr (Willeit et al., 2019). If similar conditions also prevailed in the Southern Hemisphere, they would have enhanced wind erosion and frequent dispersion events of the newly exposed, and previously stabilized material, leading to the construction of the most pronounced aeolian phase in today's landscape (Fujioka et al., 2009).

| Integrating data and process-oriented modelling to quantitatively interpret long-term aeolian dynamics
Pliocene events of aeolian transport that are revealed by Cosmolian were not previously identified by the modelling approach that relied on 26 Al/ 10 Be ratios of the same dataset of Fujioka et al. (2009), which were interpreted as burial ages. This testifies to the importance of using process-based modelling that simulates aeolian transport for constructing a detailed framework of aeolian processes and landforms (Thomas & Bailey, 2017). Furthermore, Fujioka et al. (2009) suggested that the problematic contradiction between dune stratigraphy and apparent 26 Al/ 10 Be burial ages likely stems from the recycling of older aeolian structures. These structures are probably represented by the identified aeolian phases that are older than 1.5 Ma revealed in this current work by our modelling method, and by constraining the results by the concentrations of the nuclides rather than by their ratios.
While models based solely on 26 Al/ 10 Be ratios can provide a minimum deposition age for the deepest sand within a dune, the low likelihood of dating an original dune structure (Chase, 2009;Fitzsimmons et al., 2013;Fujioka et al., 2009;Lancaster, 2008) implies that a model that simulates aeolian processes-such as Cosmolian-would better represent aeolian landscape-shaping events of a large sand body such as the LEB.
Moreover, there is no dependence between sampling depth within the original dunes and its residence time (Figure 4), thus confirming the effective mixing of sand over geological timescales (Fujioka et al., 2009;Matmon et al., 2015;Vainer et al., 2018;Vermeesch et al., 2010). This indicates that stable periods between episodes of dune accumulation, inferred by horizons that contain centimetre-size hematite nodules, could be short (≤10 4 years), which makes dunes well-mixed from the perspective of cosmogenic nuclides accumulation (Matmon et al., 2015;Vainer et al., 2018;Vermeesch et al., 2010). The fact that a single dune preserves several episodes of sand dispersion, inherited from older fossilized structures or parent F I G U R E 5 Probability density functions of simulated durations until the measured concentrations of 26 Al and 10 Be were simultaneously reached for the first time for sample P3-BR, using the empirical vertical distribution rates datasets. Outliers in the results that are marked with grey rectangles and dashed lines are excluded from the ages shown in Figure 4 [Color figure can be viewed at wileyonlinelibrary. com] bedrock, further implies that an aeolian process-based model provides a practical way for deciphering the long-term development of an entire dunefield by sampling a fairly small number of sites (Vainer et al., 2018).
Furthermore, examining how different environmental factors control aeolian transport is feasible by exploring the way that various combinations of input parameters, representing distinct environmental settings, affect aeolian migration patterns. Cosmolian outputssuch as sand exhumation events and burial frequency (number of OSL cycles in Tables S1, S2, S3, S4 and S5)-change in accordance with other output parameters, for instance the depth at which grains are commonly buried (mean dune depth in Tables S1, S2, S3, S4 and S5).
The output parameters describe the response of aeolian structures to the predefined boundary conditions (input parameters) such as the maximum dune height and the VDR dataset, which together represent varied climatic and geological settings. Therefore, inspecting the relationships between these parameters could be used to mark and quantify triggers and patterns of aeolian transport. Finally, studying the sensitivity and the dynamics of dunes in the past could shed light on the way in which dunefields will respond to forthcoming climate changes that can widely affect vast areas in southern Africa and Australia, where dunes cover over 30 and 40% of the surface, respectively (Thomas & Bailey, 2017).

| Required assumptions and limitations when using Cosmolian
The VDR dataset that is used to construct migration rates/retention periods is a fundamental parameter that controls the dynamics, and hence determines the chronological characteristics of the dunes modelled by Cosmolian during the aeolian phase ( Figure 5).
Its distribution should be therefore carefully constructed. Cosmolian comprises two sections that are used for simulation and the analyses of the results, respectively (Figure 1), and can be operated on any computer without licence requirements or prior programming knowldge. The program is available upon request, and a simplified online graphic interface for its analysis is available online at https://yoav-ben-dor.shinyapps.io/Cosmolian_analyzer/.

| CONCLUSIONS
Understanding and quantifying the rates, triggers and controls of aeolian dynamics is a prerequisite for reconstructing and predicting landscape response to climatic and environmental changes. The standalone Cosmolian program resolves the timing of episodic aeolian activity of sand by simulating the accumulation of cosmogenic nuclide concentrations and OSL signal in sand dunes.
The model used by the program is based on analytically derived datasets and is built from a process-oriented perspective. Thus, it provides a robust statistical scheme for reconstructing a chronological framework of aeolian landscape evolution over the million-year scale.
In this paper, Cosmolian was applied with data reported from the Australian Simpson Desert and revealed three key phases of aeolian emergence of sand into central Australia at 3.8-3.4, 2.9-2.5 and 1.5-1 Ma. These results fit the independently established chronological framework of landscape and climate evolution and coincide with prominent environmental changes that would have encouraged the onset and/or intensification of aeolian processes, thus reinforcing their interpretation as the residence time of the aeolian material. The results suggest that present-day dunes comprise a mixture of sand with different transportation histories that can be resolved at the regional dunefield scale by using the program for interpreting data from selected sites.
Cosmolian was designed for intuitive use by users with no programming skills and is available upon request. We encourage the widespread application of the program in order to establish a comprehensive database of aeolian simulations based on analytical data from varied locations and environments. This will allow a broad framework for inspecting the impact of different environmental settings on the initiation and cessation of aeolian transport, and the analysis of the frequency, amplitude and timing of aeolian dynamics.