Examination of Analytical Shear Stress Predictions for Coastal Dune Evolution

Cecil, Orie; Cohn, Nicholas; Farthing, Matthew; Dutta, Sourav; Trautz, Andrew

doi:https://doi.org/10.5194/egusphere-2024-855

Preprints

https://doi.org/10.5194/egusphere-2024-855

Preprints

04 Apr 2024

| 04 Apr 2024

Examination of Analytical Shear Stress Predictions for Coastal Dune Evolution

Orie Cecil, Nicholas Cohn, Matthew Farthing, Sourav Dutta, and Andrew Trautz

Abstract. Existing process-based models for simulating coastal foredune evolution largely use the same analytical approach for estimating wind induced surface shear stress distributions over spatially variable topography. Originally developed for smooth, low-sloping hills, these analytical models face significant limitations when the topography of interest exhibits large height-to-length ratios and/or steep, localized features. In this work, we utilize computational fluid dynamics (CFD) to examine the error trends of a commonly used analytical shear stress model for a series of idealized two-dimensional dune profiles. It is observed that the prediction error of the analytical model increases as compared to the CFD simulations for increasing height-to-length ratio and localized slope values. Furthermore, we explore two data-driven methodologies for generating alternative shear stress prediction models, namely, symbolic regression and linear, projection-based, non-intrusive reduced order modeling. These alternative modeling strategies demonstrate reduced overall error, but still suffer in their generalizability to broader sets of dune profiles outside of the training data. Finally, the impact of these improvements to aeolian sediment transport fluxes is examined to demonstrate that even modest improvements to the shear stress prediction can have significant impacts to dune evolution simulations over engineering-relevant timescales.

Received: 21 Mar 2024 – Discussion started: 04 Apr 2024

Publisher's note: Copernicus Publications remains neutral with regard to jurisdictional claims made in the text, published maps, institutional affiliations, or any other geographical representation in this preprint. The responsibility to include appropriate place names lies with the authors.

Download & links

Orie Cecil, Nicholas Cohn, Matthew Farthing, Sourav Dutta, and Andrew Trautz

Status: final response (author comments only)

RC1:
'Comment on egusphere-2024-855', Anonymous Referee #1, 09 May 2024

This is a timely article that rigorously quantifies values of shear stress prediction using a range of methods. This topic is important as shear stress is essential to the accurate prediction of sediment flux and sand dune evolution in a range of aeolian environments. The article provides a good starting point for further exploration of the rapid calculation of shear stress over dunes, particularly in the context of landform evolution numerical modelling tools.
There are however some details that require scrutiny before publication:
Some key details missing in the CFD simulations sections (3.2.1.) that require additional scrutiny before publication.
Section 3.2.1. – what were the discretisation schemes used in the CFD?
Section 3.2.1. – What are the computational boundary dimensions (Length and Height). Please make clear whether in this section that all calculations are 2D.
Section 3.2.1. - Where we the ‘dunes’ placed in the context of the upwind and downwind boundaries?
Section 3.2.1 – What z0 value was used in equations 7 to 11? How does this compare to the smallest cell size? (z0 must be below 50% of the bottom cell height as wall function cannot extend above this. See Figure 3 in https://doi.org/10.1016/j.atmosenv.2006.08.019)
Line 226 – What is meant by a uniform base discretization?
Line 228 – Please quantify what a ‘sufficient drop’ is.
Figures 2, 3, 5, 6, 8, 10 and 11 lack text within the legend. It was therefore difficult to confidently assess the validity of these results as well and the descriptions of the data contained within the results section and discussion. This was particularly pertinent to Figure 10 and the discussion in Section 4.3.
Discussion – The authors responsibly highlight that the results only relate two-dimensional structures. Some brief discussion should also be made that they are also relatively simple idealised structures as well.

Citation: https://doi.org/10.5194/egusphere-2024-855-RC1
- AC1:
  'Reply on RC1', Orie Cecil, 25 May 2024
  Thank you for your comments regarding our manuscript.
  
  We apologize for the lack of text within the figure legends. There was an issue with some fonts not being embedded within the figures. This has been rectified and the pre-print file updated.
  
  Concerning the request for expanded details on the computational fluid dynamics (CFD) modeling inputs and methodology (Section 3.2.1), itemized answers to your comments are listed below. These details are straightforward to include in the manuscript, which we can plan to do at the close of the ESurf open comment period.
  
  What were the discretisation schemes used in the CFD?
  The second order accurate least squares method is used for determining the cell-centered gradient values for all variables. Additionally, a cell based cubic limiter is applied for the velocity and turbulence quantities.
  
  A second order bounded linear upwind scheme was used for all convective terms.
  
  Laplacian terms used a second order accurate linear corrected scheme.
  
  What are the computational boundary dimensions (Length and Height). Please make clear whether in this section that all calculations are 2D.
  The domain height is set to be 200 m over the dune crest; therefore, the height ranges from 202.5 m for H/L = 0.1 to 212.5 m for H/L = 0.5. The total length of the domain is 425 m and the characteristic length (i.e. the half-length at half-height) of the dunes is fixed at 25 m. The domain height was determined based on extensive sensitivity studies associated with this work that were completed to minimize the effect of the upper boundary condition on the surface shear stress. Note that the grid size is variable throughout the domain, as described further below in subsequent responses.
  
  Yes, all CFD simulations were 2D in the cross-wind (cross-shore in the case of coastal dunes) and vertical dimensions. An alongshore dimension is not included in these model setups.
  
  Where were the ‘dunes’ placed in the context of the upwind and downwind boundaries?
  All ‘dunes’ have their crest at x = 30 m while the domain is centered around x = 0 m. With an overall domain length of 425 m, this places the upstream boundary 242.5 m upwind of the crest and the downstream boundary 182.5 m from the crest.
  
  What z0 value was used in equations 7 to 11? How does this compare to the smallest cell size?
  z0 was held at 1e-3 m throughout this study. The smallest size of the first cell was 5.138e-3 m for the quartic profile with H/L = 0.5. However, the wall function implementations in OpenFOAM that were used (i.e., atmNutkWallFunction and atmEpsilonWallFunction) are based directly on the aerodynamic roughness length without having to convert to an equivalent sand grain roughness height as typically required for the standard Fluent and CFX rough wall function implementations. This alleviates the burdensome requirement on the first cell height as discussed in section 7.2 of the reference cited as well as in section 2.2 of Parente et al. 2011 (https://doi.org/10.1016/j.jweia.2010.12.017)
  
  What is meant by a uniform base discretization?
  SnappyHexMesh requires a simple base discretisation which is modified by refinement around object edges (the dune profile in this case), snapping points to object edges, and finally prism inflation. The base discretisation for cases was a uniform grid consisting of 1m x 1m cells over the entire domain. The final mesh spacing is spatially variable ranging from 0.0625 m x 0.0625 m near the surface outside of the inflation layer region to 1m x 1m in the far-field at the top of the domain. The total number of cells of the final computational mesh varied across cases but was generally in the range of 390,000 to 400,000.
  
  Please quantify what a ‘sufficient drop’ is.
  Simulations were considered converged when all relative residuals as reported by OpenFOAM had fallen to a value of 1e-8 and the relative iteration-to-iteration change of the integrated drag was no greater than 1e-8.
  
  On the topic of the simple, idealized nature of the considered profiles:
  
  We agree that some additional discussion on the idealized nature of the considered topographies is warranted. The simple topographies considered only serve as a starting point for additional studies and does not capture more complex features such as compound dune shapes and the complex interactions that may occur with potential separation and reattachment over the dune itself.
  
  Thank you for highlighting the need for this additional information on model inputs and discretisations to be added to the paper!
  
  Citation: https://doi.org/10.5194/egusphere-2024-855-AC1
RC2: 'Comment on egusphere-2024-855', Orencio Duran Vinent, 07 Jun 2024

Using CFD to test the validity of linear approximations of the hydrodynamic equations outside their intended parameter range is certainly a necessary, and sometimes laborious, undertaking. However, there are several important technical and conceptual issues with this work that should be resolved before eventual publication.
1. As shown in Charru et al, Annu. Rev. Fluid Mech (2013) (Fig.3a), it is known that the linear approximation obtained by Jackson and Hunt (1975), and Hunt et al. (1988) (Eqs. 1-3 in the manuscript), overpredicts the value of the constant A (Eq. 5a) for L/z0 < 10^8 compared to the full solution of the turbulent boundary layer in the limit of small amplitudes (H/L << 1). For the value L/z0 = 25000 considered in the manuscript, the KSH model predicts A around 5, whereas the full solution would give A around 4. The authors should add this to their description of the KSH approximation and the discussion of model accuracy.
2. There is a mistake in Eqs.5b and 5c, in the actual approximation by Kroy et al. (2002), the term 2*ln(pi/2) is neglected. This term comes from the term 2*ln(kL) in Eq.2, since for a periodic modulation of wavelength lambda we have k = 2*pi/lambda, and by definition L = lambda/4.
2.1. For clarity, in Eq.2 the term ln(k) should be replaced by ln(kL) since all lengths are rescaled by L in this model.
3. The value L/z0 = 25000 used to evaluate the prediction from KSH seems arbitrary. If L = 50m and z0 = 10^-3 m, as seems to be the case for the CFD, then L/z0 = 50000. Furthermore, given the known overprediction of this approximation, the authors should add another prediction using z0=10^-4 m, which is closer to the actual roughness of sand.
4. There is a major confusion with the dune aspect ratio. In Fig.1, the aspect ratio H/L defines L as the dune's toe-crest length, whereas in all models L is defined as the half-length at half-height (the length between a point at half height and the crest, if I understood correctly). This means that there is no obvious way to compare the different simulated cases (H/L) with the relevant values in Fig.1. Here I have several suggestions:
4.1 Please use different symbols for the different lengths, for example keeping L for the half-length at half-height but using L_base for the toe-crest length. You could also add a diagram illustrating the different terms to avoid confusion.
4.2. Given the toe-crest length is around twice the half-length at half-height, you could add another axis in Fig.1 to convert the standard aspect ratio plotted (using toe-crest length) to the new one used for the simulations and all other figures. This is crucial because the simulated range in H/L = [0.01, 0.5] only represents up to H/L_base = 0.25 in Fig.1, thus missing many relevant cases.
4.3. You should consider extending the simulated ratio H/L (using half-length at half-height) all the way up to the avalanche angle (H/L = 1.3 and H/L_base = 0.65) to include the whole relevant range of foredunes' slopes. Also, please add the avalanche slope in Fig.1 for reference.
5. There is a problem with using a Bump as a dune profile to compare to CFD. Because a Bump, as defined in Eq.6d, is not a smooth profile (there is a discontinuity in the slope and/or curvature at the base) is will experience flow separation at the lee-side and flow stagnation at the toe. Therefore, the predicted shear perturbation from KSH is not defined at the base (tau' < -1) and has to be replaced by tau' = -1, which means a vanishing bed shear stress (tau = 0). Therefore, there is no meaning in using the values at these extreme locations to evaluate model accuracy (Fig.3e,f). Furthermore, the CFD predictions for the Bump (Fig.3c,d) don't really show the stagnation and flow separation expected at steep slopes, which suggests an issue with the simulations. The authors should either avoid the Bump case (and restrict only to smooth profiles) or properly discuss the complex physics involved.
6. The symbolic regression result shown in Eq.19 is physically wrong, since a constant term in the relation between the Fourier transforms of the bed shear stress perturbation and the surface implies that the shear perturbation is proportional to the surface itself which doesn't make sense (the flow field is essentially scale-invariant in the limit of large L/z0).
6.1. There is a similar problem with Eq.20, which is quadratic on the wave number k. There are strong physical reasons behind the form of Eq.5a. If the authors want to improve the KSH model, the simplest way is to fit the parameters A and B to the CFD results on smooth profiles.
7. The way to measure the differences between the models, using MSE (Eq.18) is not particularly useful because it amplifies minor horizontal shifts in the simulations. For example, even a perfect prediction shifted horizontally by 1m would lead to a large MSE. You should also include the difference in the maximum (and minimum) shear stress perturbations between the models.
8. The examples in Fig.12 for the predicted transport rate seem completely arbitrary. The linear models have been successfully used to simulate unvegetated dunes in many different situations mainly because the aspect ratio of dunes (using toe-crest length) is always less than ~0.15. Steeper dunes have to be vegetated, but then there is essentially no sand transport. The authors should either use actual dunes shapes or at least the aspect ratios of actual mobile (unvegetated) dunes to compare sand transport rates.
8.1. The transport equation (Eq.21) is wrong for aeolian sand transport. Please use the one in Martin and Kok, Science Advances (2017), which is consistent with several wind tunnel measurements and transport simulations.
8.2. The data in Fig.12f seems inconsistent with Fig.3b. In Fig.3b, both CDF and KSH are quite similar for H/L = 0.4, but they are quite different in Fig.12f for H/L = 0.5.
9. In lines 100-105. It is not true that the KSH model cannot be applied to scarps. In fact, a better version of this model (given by Weng et al. 1991) is already implemented in the Coastal Dune Model which certainly simulates dune recovery after erosion and scarp formation (Duran and Moore, Nature Clim. Change, 2015). In this case, similarly as for the Bump, the shear stress perturbation reaches < -1, and the only physically meaningful condition is tau = 0 (which follows from flow stagnation).
10. In addition to the technical issues mentioned above, I see two important conceptual problems. The manuscript introduces the main limitation of the linear models (e.g. KSH) using the steep slopes of the foredunes or vegetated dunes at the beach. However, they never compare the models using an actual foredune profile. Furthermore, the relative difference of the maximum shear stress perturbation between the KSH and the CFD for a gaussian hill seems to improve for larger H/L (Fig.3a,b) instead of getting worse, which would undercut the main purpose of the manuscript. There is also a consistency issue here because steep foredunes (where model's predictions would be worst) are covered by vegetation (that is why they can be so steep in the first place) and therefore there is no active sand transport on the crest and lee side. On the other hand, active dunes would have milder slopes and thus the model would be more suitable. These things should be properly discussed to help the reader understand the scope and meaning of the actual problems with flow simulations.
11. The second conceptual problem is related to the general assumption behind the manuscript, that somehow there should be a simple shortcut to improve the calculation of the bed shear stress under arbitrary topographic conditions. These models (in particular the 3D version in Weng et al. 1991) are already very powerful since they allowed for the first time physical simulations of 3D dune formation and interactions, formation of dune fields (Duran et al., ESPL, 2010), the formation of parabolic dunes (Duran and Herrmann, PRL, 2006), coastal dunes, dunes on Mars (Parteli and Herrmann, PRL, 2007), etc. This implementation of the Weng et al. (1991) approximation is very fast to compute and works (with the proper extension for the flow-separation region) for arbitrary topographies. Of course, it is not perfect and it is important to evaluate its accuracy, but it should be clearly stated that it represents a actual physical approximation, not a phenomenological fitting, and cannot be easily replaced by a data-based approach that requires many 3D CFD simulations to provide data. I would suggest to focus on using CFD to improve the parameters in KSH as a way to move forward in the field without having to start from scratch with data-based approaches after +40 years of physical insight and meaningful applications.

Citation: https://doi.org/10.5194/egusphere-2024-855-RC2
EC1: 'Comment on egusphere-2024-855', Andreas Baas, 19 Jun 2024

Dear authors,
Thank you for bearing with us while the manuscript was sent out to referees. The paper has now received two detailed reviews both requesting revisions, one minor and the other major. You will have seen that reviewer #2 has provided a particularly comprehensive set of comments and penetrating questions that may pose some food for thought.
I recommend that you consider options for extending the work to address some of the main concerns and questions raised by reviewer #2 (for example, the suggestion to simulate more realistic or smooth profiles), and/or to provide a detailed rebuttal.
Please do not hesitate to request extensions to the revision timetable if needed, we can easily accommodate.

Andreas Baas, King's College London

Handling Editor

Citation: https://doi.org/10.5194/egusphere-2024-855-EC1

Orie Cecil, Nicholas Cohn, Matthew Farthing, Sourav Dutta, and Andrew Trautz

Viewed

Total article views: 385 (including HTML, PDF, and XML)

HTML	PDF	XML	Total	BibTeX	EndNote
307	57	21	385	10	11

HTML: 307
PDF: 57
XML: 21
Total: 385
BibTeX: 10
EndNote: 11

Views and downloads (calculated since 04 Apr 2024)

Month	HTML	PDF	XML	Total
Apr 2024	133	20	9	162
May 2024	98	22	6	126
Jun 2024	76	15	6	97

Cumulative views and downloads (calculated since 04 Apr 2024)

Month	HTML	PDF	XML	Total
Apr 2024	133	20	9	162
May 2024	98	22	6	126
Jun 2024	76	15	6	97

Viewed (geographical distribution)

Total article views: 371 (including HTML, PDF, and XML) Thereof 371 with geography defined and 0 with unknown origin.

Country	#	Views	%

Latest update: 29 Jun 2024

Short summary

Using computational fluid dynamics, we analyze the error trends of an analytical shear stress distribution model used to drive aeolian transport for coastal dunes which are an important line of defense against storm related flooding hazards. We find that compared to numerical simulations, the analytical model results in a net overprediction of the landward migration rate. Additionally, two data-driven approaches are proposed for reducing the error while maintaining computational efficiency.


Total:	0
HTML:	0
PDF:	0
XML:	0