Numerical Model of Geochronological Tracers for Deposition and Reworking Applied to the Mississippi Subaqueous Delta

ABSTRACT Birchler, J.J.; Harris, C.K.; Kniskern, T.A., and Sherwood, C.R., 2018. Numerical model of geochronological tracers for deposition and reworking applied to the Mississippi subaqueous delta. In: Shim, J.-S.; Chun, I., and Lim, H.S. (eds.), Proceedings from the International Coastal Symposium (ICS) 2018 (Busan, Republic of Korea). Journal of Coastal Research, Special Issue No. 85, pp. 456–460. Coconut Creek (Florida), ISSN 0749-0208. Measurements of naturally occurring, short-lived radioisotopes from sediment cores on the Mississippi subaqueous delta have been used to infer event bed characteristics such as depositional thicknesses and accumulation rates. Specifically, the presence of Beryllium-7 (7Be) indicates recent riverine-derived terrestrial sediment deposition; while Thorium-234 (234Th) provides evidence of recent suspension in marine waters. Sediment transport models typically represent coastal flood and storm deposition via estimated grain size patterns and deposit thicknesses, however, and do not directly calculate radioisotope activities and profiles, which leads to a disconnect between the numerical model and field observations. Here, observed radioisotopic profiles from the Mississippi subaqueous delta cores were directly related to a numerical model that represented resuspension and deposition using a new approach to account for the behavior of short-lived radioisotopes. Appropriate selection of parameters such as the biodiffusion coefficient, sediment accumulation rate, and radioisotopic source terms enabled a good match between the modeled and observed cores. Comparisons of modelled profiles with geochronological analytical models that estimate accumulation rate and flood layer thickness revealed potential avenues for refining these tools, and highlight the importance of constraining the biodiffusion coefficient.


INTRODUCTION
Naturally occurring, short-lived radioisotopic tracers have been used to characterize sediment deposition in many coastal environments. For example, 7 Be has been used to infer recent deposition of fluvially-derived sediment (Kniskern et al., 2014;Sommerfield, Nittrouer, and Alexander, 1999), while the presence of 234 Th indicates recent resuspension in marine waters . Assessments of flood sediment budgets and deposition rates derived from radioisotopic data are often based on conceptual or analytical models that rely on estimated rates of mixing and burial (e.g. Nittrouer et al., 1984;Palinkas et al., 2005). Meanwhile, coastal sediment transport models typically restrict their calculations to erosion and deposition, grain size distributions, and changes to sediment deposit thickness (i.e. Warner et al., 2008). Validation of such models has often relied on comparison to observed radioisotopic profiles, however. For example, Xu et al. (2016) compared model estimates of deposit thicknesses from Hurricanes Katrina and Rita with an evaluation based on observed geochronological data in the northern Gulf of Mexico from Goñi et al. (2007). Corbett, McKee, and Duncan (2004) related differences in the sediment bed radioisotopic profiles and inferred depositional rates from two sampling periods to seasonal variations in the northern Gulf of Mexico. At their 'near river' site (red triangle in Figure 1), they obtained radioisotope activity profiles of 7 Be and 234 Th with depth in the sediment bed in April and October, 2000.
The seafloor here predominantly consists of mud and is influenced by wave reworking, and freshwater and sediment discharge from the Mississippi River. Bioturbation did not appear to be intense, based on analysis of X-radiographs and lack of macro-fauna obtained during core collection (Corbett, McKee, and Duncan, 2004). From April to October, 2000, the site experienced high river discharge that then decreased, low wave energy, and significant deposition ( Figure 2). The high inventories of 7 Be in October, 2000 were attributed to spring and summer sediment deposition, while the elevated 234 Th inventories were explained by increased wave energy at the end of the period (Corbett, McKee, and Duncan, 2004).

ABSTRACT
Measurements of naturally occurring, short-lived radioisotopes from sediment cores on the Mississippi subaqueous delta have been used to infer event bed characteristics such as depositional thicknesses and accumulation rates. Specifically, the presence of Beryllium-7 ( 7 Be) indicates recent riverine-derived terrestrial sediment deposition; while Thorium-234 ( 234 Th) provides evidence of recent suspension in marine waters. Sediment transport models typically represent coastal flood and storm deposition via estimated grain size patterns and deposit thicknesses, however, and do not directly calculate radioisotope activities and profiles, which leads to a disconnect between the numerical model and field observations. Here, observed radioisotopic profiles from the Mississippi subaqueous delta cores were directly related to a numerical model that represented resuspension and deposition using a new approach to account for the behavior of short-lived radioisotopes. Appropriate selection of parameters such as the biodiffusion coefficient, sediment accumulation rate, and radioisotopic source terms enabled a good match between the modeled and observed cores. Comparisons of modelled profiles with geochronological analytical models that estimate accumulation rate and flood layer thickness revealed potential avenues for refining these tools, and highlight the importance of constraining the biodiffusion coefficient.

INTRODUCTION
Naturally occurring, short-lived radioisotopic tracers have been used to characterize sediment deposition in many coastal environments. For example, 7 Be has been used to infer recent deposition of fluvially-derived sediment (Kniskern et al., 2014;Sommerfield, Nittrouer, and Alexander, 1999), while the presence of 234 Th indicates recent resuspension in marine waters . Assessments of flood sediment budgets and deposition rates derived from radioisotopic data are often based on conceptual or analytical models that rely on estimated rates of mixing and burial (e.g. Nittrouer et al., 1984;Palinkas et al., 2005). Meanwhile, coastal sediment transport models typically restrict their calculations to erosion and deposition, grain size distributions, and changes to sediment deposit thickness (i.e. Warner et al., 2008). Validation of such models has often relied on comparison to observed radioisotopic profiles, however. For example, Xu et al. (2016) compared model estimates of deposit thicknesses from Hurricanes Katrina and Rita with an evaluation based on observed geochronological data in the northern Gulf of Mexico from Goñi et al. (2007). Corbett, McKee, and Duncan (2004) related differences in the sediment bed radioisotopic profiles and inferred depositional rates from two sampling periods to seasonal variations in the northern Gulf of Mexico. At their 'near river' site (red triangle in Figure 1), they obtained radioisotope activity profiles of 7 Be and 234 Th with depth in the sediment bed in April and October, 2000.
The seafloor here predominantly consists of mud and is influenced by wave reworking, and freshwater and sediment discharge from the Mississippi River. Bioturbation did not appear to be intense, based on analysis of X-radiographs and lack of macro-fauna obtained during core collection (Corbett, McKee, and Duncan, 2004). From April to October, 2000, the site experienced high river discharge that then decreased, low wave energy, and significant deposition ( Figure 2). The high inventories of 7 Be in October, 2000 were attributed to spring and summer sediment deposition, while the elevated 234 Th inventories were explained by increased wave energy at the end of the period (Corbett, McKee, and Duncan, 2004 (Birchler et al., submitted). This paper discusses application of the coupled model to the radioisotopic cores for the Mississippi subaqueous delta measured by Corbett, McKee, and Duncan (2004) and use of the model to evaluate the applicability of analytically based estimates of deposition rates.  Warner et al., 2008), within ROMS (Regional Ocean Modeling System; see Haidvogel et al., 2008;Shchepetkin and McWilliams, 2005) was used. Recently, Moriarty et al. (2017) introduced both particulate and dissolved geochemically reactive tracers in the seabed and the water column. Using a similar framework, Birchler et al.
(submitted) developed a model to directly estimate activities of radioisotopes 7 Be and 234 Th. The one-dimensional model was implemented to represent the 'near river' site offshore of the Mississippi delta sampled by Corbett, McKee, and Duncan (2004) (Figure 1). This section summarizes the model configuration.  provides an archive of the model code, input, and output.
The model used 30 uniformly spaced layers in the vertical to represent the 50 m water column, while the sediment bed model had 40 layers that were each initially 0.5 cm thick. Two grain sizes were used to represent mud, with diameters of 0.015 and 0.063 mm; settling velocities of 0.01 and 0.1 cm s -1 ; and critical shear stresses of 0.03 and 0.08 Pa, respectively. The erosion rate parameter (see Warner et al., 2008) was chosen as 5x10 -5 kg m -2 s -1 so that the model produced erodibility curves consistent with those measured near the Mississippi delta by Xu et al. (2014). As the model proceeded, newly delivered river sediment supplied 7 Be. For suspended material, 234 Th was reset to maintain a constant activity, which supplied 234 Th to the bed during cycles of sediment resuspension and deposition. A constant wind speed of 15 m s -1 was applied to create steady currents of about 10 cm s -1 . For input timeseries, hourly wave height and period were taken from buoys 42040 and 42007 (Figure 1) from the National Oceanic and Atmospheric Administration National Data Buoy Center (NDBC, 2013).
One-dimensional (vertical) models of sediment transport typically neglect horizontal flux convergences and divergences that lead to net erosion or deposition. The study location accumulates ~1-3 cm of new sediment per month, however, and at times also supplies sediment that is eroded and redistributed to downstream locations (Corbett, McKee, and Duncan, 2004). Source and sink terms were therefore added as surface tracer fluxes (m s -1 ) to represent Mississippi River sediment delivered to, or eroded sediment removed from, the study site from April to October, 2000. The timing of the surface tracer source was input at a rate proportional to the observed river sediment discharge at Tarbert Landing, MS, obtained from the U.S. Geological Survey (Figure 2a). The surface tracer flux was scaled so that a total of 4.6 cm of sediment was deposited to produce modeled radioisotope profiles that matched those reported in Corbett, McKee, and Duncan (2004). Episodic erosion coincided with the four wave resuspension events where the bed shear stress exceeded 0.1 Pa (Figure 2c). To account for removal and erosion of sediment, roughly 50% of suspended material was removed as a surface tracer flux during timesteps when the bed shear stress exceeded 0.1 Pa. The sediment deposition, erosion, and input radioisotope activity were chosen to replicate observed penetration depths and activity profiles. The main choices made to match the profiles were the biodiffusion coefficients and deposition rates used, and the activity of new sediment (Birchler, 2014).
The model represented April to October, 2000 ( Figure 2). Because the data from October, 2000 had 234 Th observations only for the very surface sample (Corbett, McKee, and Duncan, 2004), the effort to reproduce the observed profiles focused mainly on 7 Be (Figure 3). A low, but non-zero, biodiffusion rate was needed in the model to obtain reasonable penetration depths. Biodiffusion was implemented as described in Sherwood et al. (in press), where Db,max was 0.95 cm 2 yr -1 (3x10 -12 m 2 s -1 ), Db,min was 0.016 cm 2 yr -1 (5x10 -14 m 2 s -1 ) and Zb,max was 3 cm. When biodiffusion was neglected, the surface activity of both radioisotopes was slightly too high, and the penetration depth was too shallow (Birchler, 2014).
The radioisotope activity profiles were initialized with those found in April, 2000 by Corbett, McKee, and Duncan (2004) (Figure 3a). For 7 Be, a source activity of 20 dpm g -1 was needed to match the final observed 7 Be surface activity from October, 2000. A value of 95 dpm g -1 was chosen for the water column activity of 234 Th in order to match the high surface activity observed in October, 2000. Because nearly all of the modeled deposition occurred in the first half of the time period, 7 Be needed a reasonably high input activity as its signal decayed significantly by the simulation's end.
To investigate the role of biodiffusion on radioisotope profiles, two additional models were run that were identical except for the maximum bioturbation rates; the models had bioturbation rates of 0 and 25 cm 2 yr -1 , respectively.

RESULTS
The profiles of 7 Be and 234 Th changed through time from the initial profile due to sediment deposition, decay, and reworking due to waves and bioturbation (Figure 3a-e). Time series of surface activities, inventories, and penetration depths calculated for both radionuclides illustrates their response to sediment erosion and deposition (Figure 2, 4). Peaks in 7 Be and 234 Th surface activity were associated with fluvial delivery, but did not scale directly with discharge, which was largest early in the model (days 0-15) when seabed activities continued to reflect initial conditions ( Figure 2). Instead, peak surface activities co-occurred with a smaller river pulse around day 100, when bed stresses were smaller and net erosion was nil (Figure 4). For the rest of the model, river discharge decreased which slowed the addition of tracer to the seabed, and surface activities of both tracers fell via radioisotopic decay and biodiffusive dilution. Wave-driven erosion that removed high-activity sediments also decreased bed inventory and surface activity late in the model. Penetration depths of 7 Be and 234 Th were initially about 10 cm and 6 cm, respectively (Figure 4c). Penetration depths of 7 Be and 234 Th changed slightly until day 100, after which a large sediment input increased the penetration depth of 7 Be and 234 Th to 10.5 cm and 8 cm, respectively. After day 100, declining river discharge and two erosion events decreased the penetration depths of 7 Be to about 6 cm by the end of the model. Penetration depths only loosely corresponded to the sediment input timeseries because of confounding processes of sediment erosion and biodiffusive mixing.

DISCUSSION
The modeled radioisotope profiles were interpreted using multiple analytical methods that have been applied to sediment core data to estimate the accumulation rate and deposit thickness. These estimated rates were then compared to the modeled accumulation rates and deposit thicknesses.

Accumulation Rate Estimates
For situations where bioturbation and physical mixing appear to be less influential than deposition, short-lived radioisotopic records from sediment cores are often analyzed to estimate seasonal sediment accumulation rates (Sommerfield, Nittrouer, and Alexander, 1999). Two of these methods were used, that differ in whether they account for biodiffusion, to estimate accumulation rate from the final modeled radioisotopic profiles. These were then compared to the actual accumulation rate: the modeled net deposition divided by the duration of the model run (A3). The first analytical method applied a steady-state solution that assumed accumulation exceeds the effect of bioturbation : where A1 represents accumulation rate in cm month -1 ;  is the decay constant of the radionuclide; z is depth in the seabed; Co is radioactivity at the surface; and Cz is radioactivity at depth. The second method accounted for biodiffusion, fitting an advectiondiffusion equation to the final profile : where Db represented the biodiffusion rate in cm 2 yr -1 . Compared to the actual value of ~0.66 cm month -1 , the steadystate analytical estimates of accumulation rate varied from net erosion (-0.33 cm month -1 ) to double the actual value (Table 1). Rates based on 7 Be were more reliable than those from 234 Th. These analytically-derived values were especially unreliable when biodiffusion was significant (Table 1). Even the steadystate model that accounted for biodiffusion using the correct value of Db,max had this problem (Table 1, A2, bottom row). These analytical solutions were imprecise because deposition and erosion during the modeled period produced radioisotopic profiles that did not represent steady state conditions (e.g. Sadler, 1999). During erosive periods, radioisotope-tagged sediment was removed from the surface, and what remained was mixed more deeply into the seabed. Resuspension increased 234 Th inventory via water column replenishment, which in most cases produced higher apparent accumulation rates based on 234 Th.

_________________________________________________________________________________________________
Journal of Coastal Research, Special Issue No. 85, 2018

Deposit Thickness Estimates
Alternatively, episodic seabed elevation changes can be analyzed using non-steady state solutions. Palinkas et al.'s (2005) approach assesses deposit thickness from modeled non-steadystate radioisotope profiles that account for episodic deposition. In this study, the approach was modified to assess relative rates of diffusive bioturbation and net deposition and applied the appropriate advection-diffusion equation (e.g. Nittrouer et al., 1984) to produce an expected 7 Be profile between each model time step (Palinkas et al., 2005). If the expected activity profile did not match the observed, a new layer of sediment was added to the top of the seabed with activity equal to the observed surface layer until the calculated profile matched the observed (Palinkas et al., 2005). To test the method using the modeled profiles, this routine was applied to each daily time interval and each seabed layer independently for incremental deposit-layer thicknesses ranging from one to the total seabed thickness. The two calculated profiles with the least amount of deviation above and below the actual deposit thickness at each time were identified and adjusted to the actual deposit thickness. This reconstructed the depositional history with precision ( Figure 5), but assumed perfect knowledge of the biodiffusion coefficient and radioisotopic profiles at daily and sub-cm resolution. A more likely scenario would be a case where net deposition would be estimated based on sediment cores taken some weeks or months apart. For example, the difference in 7 Be inventories in the radioisotopic profiles in Figure 3 could be used to infer depositional history. Applying the methods of Canuel, Martens, and Benninger (1990), deposit thicknesses were estimated for the four time intervals illustrated in Figure 3 and compared to the actual modeled deposition (Table 2). For the first three intervals, which were short, the analytical estimate agreed to within 10% of the actual deposition. For the final period, however, which covered days 41-210, the analytical approach underestimated the actual deposition by a factor of two (Table 2).  (Canuel, Martens, and Benninger, 1990) from 7 Be profiles in Figure 3;  This analytical method assumes that all of the new deposition occurs instantaneously, and that the activity of newly deposited sediment equals the activity of surficial sediment at the end of the study period. This method underestimated the actual deposition between days 41-210 because the modeled deposition was somewhat gradual then, but the analytical method assumed that much of the deposition occurred early in the time interval.

CONCLUSIONS
A combined sediment transport and radioisotope model was applied to a realistic shelf setting to calculate sediment bed profiles representing short half-life radioisotopes 7 Be and 234 Th that are used to interpret deposition of river derived sediment and sediment resuspension, respectively. This provided an example in a one-dimensional (vertical) model of suspended and seabed sediment of a realistic scenario that evaluated the response of radioisotope profiles to variations in riverine sediment input, storm intensity, and biodiffusion. Success in reproducing the observed profiles measured by Corbett, McKee, and Duncan (2004) for a 50 m site near the Southwest Pass of the Mississippi River was possible with proper specification of the 7 Be and 234 Th activities of fluvial and resuspended sediment; and selection of biodiffusion coefficients, deposition, and erosion. Radioisotopic values such as surface activity, inventory, and penetration depth reflect processes operating at timescales ranging from individual floods or storms, to seasonal. Episodic increases and decreases in these values occurred in response to depositional and erosional events, but these were also gradually modified by mixing within the sediment bed, background deposition, and radioisotopic decay.
The model's accumulation rates and deposit thicknesses were compared to values obtained via analytical methods commonly applied to sediment cores. Steady-state analytical estimates of accumulation rate based on 7 Be were generally more reliable than those based on 234 Th which was enriched via resuspension. For the case using an intense bioturbation rate, the steady-state analytical expressions were imprecise, even Eq. 2 which accounted for biodiffusion. A non-steady state estimate of new deposition reliably calculated the modeled deposit thickness _________________________________________________________________________________________________ Journal of Coastal Research, Special Issue No. 85, 2018 when it was run at temporal and spatial resolutions that greatly exceeded those typically achieved by field studies. When applied at the temporal scales more typical of field studies, however, the analytical method of Canuel, Martens, and Benninger (1990) underestimated deposit thickness due to uncertainties in the depositional history.