CoastalME (Coastal Modelling Environment)
Simulates the long-term behaviour of complex coastlines
Loading...
Searching...
No Matches
do_beach_potential_erosion.cpp
Go to the documentation of this file.
1
12
13/* ==============================================================================================================================
14
15 This file is part of CoastalME, the Coastal Modelling Environment.
16
17 CoastalME is free software; you can redistribute it and/or modify it under the terms of the GNU General Public License as published by the Free Software Foundation; either version 3 of the License, or (at your option) any later version.
18
19 This program is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more details.
20
21 You should have received a copy of the GNU General Public License along with this program; if not, write to the Free Software Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
22
23==============================================================================================================================*/
24#include <assert.h>
25
26#include <cmath>
27#include <cfloat>
28
29// #include <iostream>
30// using std::cout;
31// using std::endl;
32
33#include <algorithm>
34using std::sort;
35
36#include <utility>
37using std::make_pair;
38using std::pair;
39
40#include "cme.h"
41#include "simulation.h"
42#include "coast.h"
43
44namespace
45{
46 //===============================================================================================================================
48 //===============================================================================================================================
49 bool bPolygonLengthPairCompare(const pair<int, double> &prLeft, const pair<int, double> &prRight)
50 {
51 // Sort in ascending order (i.e. most concave first)
52 return prLeft.second < prRight.second;
53 }
54}
55
56//===============================================================================================================================
58//===============================================================================================================================
60{
61 // Do this for each coast
62 for (int nCoast = 0; nCoast < static_cast<int>(m_VCoast.size()); nCoast++)
63 {
64 int const nNumPolygons = m_VCoast[nCoast].nGetNumPolygons();
65
66 // Create a vector of pairs: the first value of the pair is the profile index, the second is the seaward length of that profile
67 vector<pair<int, double>> prVPolygonLength;
68
69 for (int nPoly = 0; nPoly < nNumPolygons; nPoly++)
70 {
71 CGeomCoastPolygon const* pPolygon = m_VCoast[nCoast].pGetPolygon(nPoly);
72 double const dSeawardLength = pPolygon->dGetLength();
73 prVPolygonLength.push_back(make_pair(nPoly, dSeawardLength));
74 }
75
76 // Sort this pair vector in ascending order, so that the polygons with the shortest length (i.e. the most concave polygons) are first
77 sort(prVPolygonLength.begin(), prVPolygonLength.end(), bPolygonLengthPairCompare);
78
79 // Do this for every coastal polygon in sequence of coastline concavity
80 for (int n = 0; n < nNumPolygons; n++)
81 {
82 int const nThisPoly = prVPolygonLength[n].first;
83
84 CGeomCoastPolygon* pPolygon = m_VCoast[nCoast].pGetPolygon(nThisPoly);
85
86 // Calculate the average breaking wave height and angle along this polygon's segment of coastline
87 int const nStartNormal = pPolygon->nGetUpCoastProfile();
88 int const nEndNormal = pPolygon->nGetDownCoastProfile();
89 int const nCoastStartPoint = m_VCoast[nCoast].pGetProfile(nStartNormal)->nGetCoastPoint();
90 int const nCoastEndPoint = m_VCoast[nCoast].pGetProfile(nEndNormal)->nGetCoastPoint();
91 int nCoastPoints = 0;
92 int nActiveZonePoints = 0;
93
94 double dAvgBreakingWaveHeight = 0;
95 double dAvgBreakingWaveAngle = 0;
96 double dAvgDeepWaterWavePeriod = 0;
97 double dAvgFluxOrientation = 0;
98 double dAvgBreakingDepth = 0;
99 double dAvgBreakingDist = 0;
100
101 // Calculate the average tangent to the polygon's coast segment, the average wave breaking height, the average depth of breaking, and the average distance of breaking, for this coast segment
102 for (int nCoastPoint = nCoastStartPoint; nCoastPoint < nCoastEndPoint - 1; nCoastPoint++)
103 {
104 nCoastPoints++;
105 dAvgFluxOrientation += m_VCoast[nCoast].dGetFluxOrientation(nCoastPoint);
106
107 double const dThisBreakingWaveHeight = m_VCoast[nCoast].dGetBreakingWaveHeight(nCoastPoint);
108
109 if (! bFPIsEqual(dThisBreakingWaveHeight, DBL_NODATA, TOLERANCE))
110 {
111 // We are in the active zone
112 nActiveZonePoints++;
113 dAvgBreakingWaveHeight += dThisBreakingWaveHeight;
114
115 double const dThisBreakingWaveAngle = m_VCoast[nCoast].dGetBreakingWaveAngle(nCoastPoint);
116 double const dThisDeepWaterWavePeriod = m_VCoast[nCoast].dGetCoastDeepWaterWavePeriod(nCoastPoint);
117
118 dAvgBreakingWaveAngle += dThisBreakingWaveAngle;
119 dAvgDeepWaterWavePeriod += dThisDeepWaterWavePeriod;
120
121 dAvgBreakingDepth += m_VCoast[nCoast].dGetDepthOfBreaking(nCoastPoint);
122
123 dAvgBreakingDist += (m_VCoast[nCoast].nGetBreakingDistance(nCoastPoint) * m_dCellSide);
124 }
125 }
126
127 // Safety check
128 if (nCoastPoints == 0)
129 nCoastPoints = 1;
130
131 // Calc the averages
132 dAvgFluxOrientation /= nCoastPoints;
133
134 if (nActiveZonePoints > 0)
135 {
136 // Only calculate sediment movement if the polygon has at least one coast point in its coastline segment which is in the active zone
137 dAvgBreakingWaveHeight /= nActiveZonePoints;
138 dAvgBreakingWaveAngle /= nActiveZonePoints;
139 dAvgDeepWaterWavePeriod /= nActiveZonePoints;
140 dAvgBreakingDepth /= nActiveZonePoints;
141 dAvgBreakingDist /= nActiveZonePoints;
142
143 // Get the coast handedness, and (based on the average tangent) calculate the direction towards which a coastline-normal profile points
144 int const nSeaHand = m_VCoast[nCoast].nGetSeaHandedness();
145 double dNormalOrientation;
146
147 if (nSeaHand == RIGHT_HANDED)
148 dNormalOrientation = dKeepWithin360(dAvgFluxOrientation - 90);
149
150 else
151 dNormalOrientation = dKeepWithin360(dAvgFluxOrientation + 90);
152
153 // Determine dThetaBr, the angle between the coastline-normal orientation and the breaking wave orientation (the direction FROM which the waves move). This tells us whether the sediment movement is up-coast (-ve) or down-coast (+ve)
154 double dThetaBr = dNormalOrientation - dAvgBreakingWaveAngle;
155
156 if (dThetaBr > 270)
157 dThetaBr = dAvgBreakingWaveAngle + 360.0 - dNormalOrientation;
158
159 else if (dThetaBr < -270)
160 dThetaBr = dNormalOrientation + 360.0 - dAvgBreakingWaveAngle;
161
162 bool bDownCoast = true;
163
164 if (dThetaBr < 0)
165 bDownCoast = false;
166
167 // And save the direction of sediment movement in the polygon object
168 pPolygon->SetDownCoastThisIter(bDownCoast);
169
170 // Now that we have the direction of sediment movement, normalize dThetaBr to be always +ve so that subsequent calculations are clearer
171 dThetaBr = tAbs(dThetaBr);
172
173 // Safety check: not sure why we need this, but get occasional big values
174 dThetaBr = tMin(dThetaBr, 90.0);
175
176 // Calculate the immersed weight of sediment transport
177 double dImmersedWeightTransport = 0;
178
180 {
181 /*
182 Use the CERC equation (Komar and Inman, 1970; USACE, 1984), this describes the immersive weight transport of sand (i.e. sand transport in suspension). Depth-integrated alongshore volumetric sediment transport is a function of breaking wave height Hb and angle αb:
183
184 Qls = Kls * Hb^(5/2) * sin(2 * αb)
185
186 where Kls is a transport coefficient which varies between 0.4 to 0.79
187 */
188 dImmersedWeightTransport = m_dKLS / (16 * pow(m_dBreakingWaveHeightDepthRatio, 0.5)) * m_dSeaWaterDensity * pow(m_dG, 1.5) * pow(dAvgBreakingWaveHeight, 2.5) * sin((PI / 180) * 2 * dThetaBr);
189 }
190
192 {
193 /*
194 Use the Kamphuis (1990) equation to estimate the immersive weight transport of sand in kg/s:
195
196 Qls = 2.33 * (Tp^(1.5)) * (tanBeta^(0.75)) * (d50^(-0.25)) * (Hb^2) * (sin(2 * αb)^(0.6))
197
198 where:
199
200 Tp = peak wave period
201 tanBeta = beach slope, defined as the ratio of the water depth at the breaker line and the distance from the still water beach line to the breaker line
202 d50 = median particle size in surf zone (m)
203 */
204
205 if (dAvgBreakingDist > 0)
206 {
207 double const dD50 = pPolygon->dGetAvgUnconsD50();
208
209 if (dD50 > 0)
210 {
211 double const dBeachSlope = dAvgBreakingDepth / dAvgBreakingDist;
212
213 // Note that we use a calibration constant here (m_dKamphuis)
214 dImmersedWeightTransport = m_dKamphuis * 2.33 * pow(dAvgDeepWaterWavePeriod, 1.5) * pow(dBeachSlope, 0.75) * pow(dD50, -0.25) * pow(dAvgBreakingWaveHeight, 2) * pow(sin((PI / 180) * 2 * dThetaBr), 0.6);
215 }
216 }
217 }
218
219 // Convert from immersed weight rate to bulk volumetric (sand and voids) transport rate (m3/s)
220 double dSedimentVol = m_dInmersedToBulkVolumetric * dImmersedWeightTransport;
221
222 // Convert to a volume during this timestep
223 dSedimentVol *= (m_dTimeStep * 3600);
224
225 // Convert to a depth in m
226 double dSedimentDepth = dSedimentVol / m_dCellArea;
227
228 // If this is just a tiny depth, do nothing
229 if (dSedimentDepth < SEDIMENT_ELEV_TOLERANCE)
230 dSedimentDepth = 0;
231
232 // LogStream << m_ulIter << ": polygon = " << nThisPoly << " nActiveZonePoints = " << nActiveZonePoints << " dAvgBreakingWaveHeight = " << dAvgBreakingWaveHeight << " dAvgFluxOrientation = " << dAvgFluxOrientation << " dNormalOrientation = " << dNormalOrientation << " dAvgBreakingWaveAngle = " << dAvgBreakingWaveAngle << " potential sediment transport this timestep = " << dSedimentDepth << " m " << (bDownCoast ? "DOWN" : "UP") << " coast" << endl;
233
234 // Store the potential erosion value for this polygon
235 pPolygon->AddPotentialErosion(-dSedimentDepth);
236 // LogStream << "\tPotential erosion on polygon " << nThisPoly << " -dSedimentDepth = " << -dSedimentDepth << endl;
237 }
238
239 // else
240 // LogStream << m_ulIter << ": polygon = " << nThisPoly << " NOT IN ACTIVE ZONE dAvgFluxOrientation = " << dAvgFluxOrientation << endl;
241 }
242 }
243}
Geometry class used for coast polygon objects.
double dGetLength(void) const
Gets the polygon's length.
void AddPotentialErosion(double const)
Adds in potential beach erosion of unconsolidated sediment (all size classes) on this polygon (m_dPot...
void SetDownCoastThisIter(bool const)
Set a flag to say whether sediment movement on this polygon is down-coast this iteration.
int nGetUpCoastProfile(void) const
Return the number of the up-coast profile.
double dGetAvgUnconsD50(void) const
Get the average d50 for unconsolidated sediment on this polygon.
int nGetDownCoastProfile(void) const
Return the number of the down-coast profile.
double m_dG
Gravitational acceleration (m**2/sec)
Definition simulation.h:779
int m_nBeachErosionDepositionEquation
Which beach erosion-deposition equation is used. Possible values are UNCONS_SEDIMENT_EQUATION_CERC an...
Definition simulation.h:498
vector< CRWCoast > m_VCoast
The coastline objects.
static double dKeepWithin360(double const)
Constrains the supplied angle to be within 0 and 360 degrees.
double m_dInmersedToBulkVolumetric
For beach erosion/deposition, conversion from immersed weight to bulk volumetric (sand and voids) tra...
Definition simulation.h:782
double m_dBreakingWaveHeightDepthRatio
Breaking wave height-to-depth ratio.
Definition simulation.h:719
double m_dSeaWaterDensity
Density of sea water in kg/m**3.
Definition simulation.h:668
double m_dKLS
Transport parameter KLS in the CERC equation.
Definition simulation.h:773
void DoAllPotentialBeachErosion(void)
Uses either the CERC equation or the Kamphuis (1990) equation to calculate potential (unconstrained) ...
double m_dTimeStep
The length of an iteration (a time step) in hours.
Definition simulation.h:644
double m_dCellArea
Area of a cell (in external CRS units)
Definition simulation.h:629
double m_dCellSide
Length of a cell side (in external CRS units)
Definition simulation.h:626
double m_dKamphuis
Transport parameter for the Kamphuis equation.
Definition simulation.h:776
This file contains global definitions for CoastalME.
T tMin(T a, T b)
Definition cme.h:1240
double const TOLERANCE
Definition cme.h:810
int const UNCONS_SEDIMENT_EQUATION_KAMPHUIS
Definition cme.h:784
double const PI
Definition cme.h:792
bool bFPIsEqual(const T d1, const T d2, const T dEpsilon)
Definition cme.h:1273
double const SEDIMENT_ELEV_TOLERANCE
Definition cme.h:812
double const DBL_NODATA
Definition cme.h:821
int const RIGHT_HANDED
Definition cme.h:511
int const UNCONS_SEDIMENT_EQUATION_CERC
Definition cme.h:783
T tAbs(T a)
Definition cme.h:1250
Contains CRWCoast definitions.
Contains CSimulation definitions.