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
28#include <cfloat>
29
30#include <iostream>
31using std::cout;
32using std::endl;
33
34#include <algorithm>
35using std::stable_sort;
36
37#include <utility>
38using std::make_pair;
39using std::pair;
40
41#include "cme.h"
42#include "simulation.h"
43#include "coast.h"
44
45//===============================================================================================================================
47//===============================================================================================================================
48bool bPolygonLengthPairCompare(const pair<int, double> &prLeft, const pair<int, double> &prRight)
49{
50 // Sort in ascending order (i.e. most concave first)
51 return prLeft.second < prRight.second;
52}
53
54//===============================================================================================================================
56//===============================================================================================================================
58{
59 // Do this for each coast
60 for (int nCoast = 0; nCoast < static_cast<int>(m_VCoast.size()); nCoast++)
61 {
62 int nNumPolygons = m_VCoast[nCoast].nGetNumPolygons();
63
64 // Create a vector of pairs: the first value of the pair is the profile index, the second is the seaward length of that profile
65 vector<pair<int, double>> prVPolygonLength;
66
67 for (int nPoly = 0; nPoly < nNumPolygons; nPoly++)
68 {
69 CGeomCoastPolygon const* pPolygon = m_VCoast[nCoast].pGetPolygon(nPoly);
70 double dSeawardLength = pPolygon->dGetLength();
71 prVPolygonLength.push_back(make_pair(nPoly, dSeawardLength));
72 }
73
74 // Sort this pair vector in ascending order, so that the polygons with the shortest length (i.e. the most concave polygons) are first
75 sort(prVPolygonLength.begin(), prVPolygonLength.end(), bPolygonLengthPairCompare);
76
77 // Do this for every coastal polygon in sequence of coastline concavity
78 for (int n = 0; n < nNumPolygons; n++)
79 {
80 int nThisPoly = prVPolygonLength[n].first;
81
82 CGeomCoastPolygon* pPolygon = m_VCoast[nCoast].pGetPolygon(nThisPoly);
83
84 // Calculate the average breaking wave height and angle along this polygon's segment of coastline
85 int nStartNormal = pPolygon->nGetUpCoastProfile();
86 int nEndNormal = pPolygon->nGetDownCoastProfile();
87 int nCoastStartPoint = m_VCoast[nCoast].pGetProfile(nStartNormal)->nGetCoastPoint();
88 int nCoastEndPoint = m_VCoast[nCoast].pGetProfile(nEndNormal)->nGetCoastPoint();
89 int nCoastPoints = 0;
90 int nActiveZonePoints = 0;
91
92 double dAvgBreakingWaveHeight = 0;
93 double dAvgBreakingWaveAngle = 0;
94 double dAvgDeepWaterWavePeriod = 0;
95 double dAvgFluxOrientation = 0;
96 double dAvgBreakingDepth = 0;
97 double dAvgBreakingDist = 0;
98
99 // 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
100 for (int nCoastPoint = nCoastStartPoint; nCoastPoint < nCoastEndPoint - 1; nCoastPoint++)
101 {
102 nCoastPoints++;
103 dAvgFluxOrientation += m_VCoast[nCoast].dGetFluxOrientation(nCoastPoint);
104
105 double dThisBreakingWaveHeight = m_VCoast[nCoast].dGetBreakingWaveHeight(nCoastPoint);
106
107 if (! bFPIsEqual(dThisBreakingWaveHeight, DBL_NODATA, TOLERANCE))
108 {
109 // We are in the active zone
110 nActiveZonePoints++;
111 dAvgBreakingWaveHeight += dThisBreakingWaveHeight;
112
113 double dThisBreakingWaveAngle = m_VCoast[nCoast].dGetBreakingWaveAngle(nCoastPoint);
114 double dThisDeepWaterWavePeriod = m_VCoast[nCoast].dGetCoastDeepWaterWavePeriod(nCoastPoint);
115
116 dAvgBreakingWaveAngle += dThisBreakingWaveAngle;
117 dAvgDeepWaterWavePeriod += dThisDeepWaterWavePeriod;
118
119 dAvgBreakingDepth += m_VCoast[nCoast].dGetDepthOfBreaking(nCoastPoint);
120
121 dAvgBreakingDist += (m_VCoast[nCoast].nGetBreakingDistance(nCoastPoint) * m_dCellSide);
122 }
123 }
124
125 // Safety check
126 if (nCoastPoints == 0)
127 nCoastPoints = 1;
128
129 // Calc the averages
130 dAvgFluxOrientation /= nCoastPoints;
131
132 if (nActiveZonePoints > 0)
133 {
134 // Only calculate sediment movement if the polygon has at least one coast point in its coastline segment which is in the active zone
135 dAvgBreakingWaveHeight /= nActiveZonePoints;
136 dAvgBreakingWaveAngle /= nActiveZonePoints;
137 dAvgDeepWaterWavePeriod /= nActiveZonePoints;
138 dAvgBreakingDepth /= nActiveZonePoints;
139 dAvgBreakingDist /= nActiveZonePoints;
140
141 // Get the coast handedness, and (based on the average tangent) calculate the direction towards which a coastline-normal profile points
142 int nSeaHand = m_VCoast[nCoast].nGetSeaHandedness();
143 double dNormalOrientation;
144
145 if (nSeaHand == RIGHT_HANDED)
146 dNormalOrientation = dKeepWithin360(dAvgFluxOrientation - 90);
147
148 else
149 dNormalOrientation = dKeepWithin360(dAvgFluxOrientation + 90);
150
151 // 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)
152 double dThetaBr = dNormalOrientation - dAvgBreakingWaveAngle;
153
154 if (dThetaBr > 270)
155 dThetaBr = dAvgBreakingWaveAngle + 360.0 - dNormalOrientation;
156
157 else if (dThetaBr < -270)
158 dThetaBr = dNormalOrientation + 360.0 - dAvgBreakingWaveAngle;
159
160 bool bDownCoast = true;
161
162 if (dThetaBr < 0)
163 bDownCoast = false;
164
165 // And save the direction of sediment movement in the polygon object
166 pPolygon->SetDownCoastThisIter(bDownCoast);
167
168 // Now that we have the direction of sediment movement, normalize dThetaBr to be always +ve so that subsequent calculations are clearer
169 dThetaBr = tAbs(dThetaBr);
170
171 // Safety check: not sure why we need this, but get occasional big values
172 dThetaBr = tMin(dThetaBr, 90.0);
173
174 // Calculate the immersed weight of sediment transport
175 double dImmersedWeightTransport = 0;
176
178 {
179 /*
180 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:
181
182 Qls = Kls * Hb^(5/2) * sin(2 * αb)
183
184 where Kls is a transport coefficient which varies between 0.4 to 0.79
185 */
186 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);
187 }
188
190 {
191 /*
192 Use the Kamphuis (1990) equation to estimate the immersive weight transport of sand in kg/s:
193
194 Qls = 2.33 * (Tp^(1.5)) * (tanBeta^(0.75)) * (d50^(-0.25)) * (Hb^2) * (sin(2 * αb)^(0.6))
195
196 where:
197
198 Tp = peak wave period
199 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
200 d50 = median particle size in surf zone (m)
201 */
202
203 if (dAvgBreakingDist > 0)
204 {
205 double dD50 = pPolygon->dGetAvgUnconsD50();
206
207 if (dD50 > 0)
208 {
209 double dBeachSlope = dAvgBreakingDepth / dAvgBreakingDist;
210
211 // Note that we use a calibration constant here (m_dKamphuis)
212 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);
213 }
214 }
215 }
216
217 // Convert from immersed weight rate to bulk volumetric (sand and voids) transport rate (m3/s)
218 double dSedimentVol = m_dInmersedToBulkVolumetric * dImmersedWeightTransport;
219
220 // Convert to a volume during this timestep
221 dSedimentVol *= (m_dTimeStep * 3600);
222
223 // Convert to a depth in m
224 double dSedimentDepth = dSedimentVol / m_dCellArea;
225
226 // If this is just a tiny depth, do nothing
227 if (dSedimentDepth < SEDIMENT_ELEV_TOLERANCE)
228 dSedimentDepth = 0;
229
230 // 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;
231
232 // Store the potential erosion value for this polygon
233 pPolygon->AddPotentialErosion(-dSedimentDepth);
234 // LogStream << "\tPotential erosion on polygon " << nThisPoly << " -dSedimentDepth = " << -dSedimentDepth << endl;
235 }
236
237 // else
238 // LogStream << m_ulIter << ": polygon = " << nThisPoly << " NOT IN ACTIVE ZONE dAvgFluxOrientation = " << dAvgFluxOrientation << endl;
239 }
240 }
241}
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:849
int m_nBeachErosionDepositionEquation
Definition simulation.h:543
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
Definition simulation.h:854
double m_dBreakingWaveHeightDepthRatio
Breaking wave height-to-depth ratio.
Definition simulation.h:789
double m_dSeaWaterDensity
Density of sea water in kg/m**3.
Definition simulation.h:734
double m_dKLS
Transport parameter KLS in the CERC equation.
Definition simulation.h:843
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:707
double m_dCellArea
Area of a cell (in external CRS units)
Definition simulation.h:692
double m_dCellSide
Length of a cell side (in external CRS units)
Definition simulation.h:689
double m_dKamphuis
Transport parameter for the Kamphuis equation.
Definition simulation.h:846
This file contains global definitions for CoastalME.
T tMin(T a, T b)
Definition cme.h:1251
double const TOLERANCE
Definition cme.h:811
int const UNCONS_SEDIMENT_EQUATION_KAMPHUIS
Definition cme.h:785
double const PI
Definition cme.h:793
bool bFPIsEqual(const T d1, const T d2, const T dEpsilon)
Definition cme.h:1284
double const SEDIMENT_ELEV_TOLERANCE
Definition cme.h:813
double const DBL_NODATA
Definition cme.h:822
int const RIGHT_HANDED
Definition cme.h:509
int const UNCONS_SEDIMENT_EQUATION_CERC
Definition cme.h:784
T tAbs(T a)
Definition cme.h:1261
Contains CRWCoast definitions.
bool bPolygonLengthPairCompare(const pair< int, double > &prLeft, const pair< int, double > &prRight)
Function used to sort polygon length values. If the first argument must be ordered before the second,...
Contains CSimulation definitions.